Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/ice.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/ice.F90	(revision 11586)
@@ -110,4 +110,5 @@
    !! bv_i        |      -      |    relative brine volume        | ???   | 
    !! at_ip       |      -      |    Total ice pond concentration |       |
+   !! hm_ip       |      -      |    Mean ice pond depth          | m     |
    !! vt_ip       |      -      |    Total ice pond vol. per unit area| m |
    !!=====================================================================
@@ -188,4 +189,5 @@
 
    !                                     !!** ice-ponds namelist (namthd_pnd)
+   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F)
    LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012
    LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth
@@ -196,4 +198,6 @@
    !                                     !!** ice-diagnostics namelist (namdia) **
    LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)
+   REAL(wp), PUBLIC ::   rn_icechk_cel    !: rate of ice spuriously gained/lost (at any gridcell)
+   REAL(wp), PUBLIC ::   rn_icechk_glo    !: rate of ice spuriously gained/lost (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)
@@ -330,4 +334,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_ip      !: mean melt pond depth                     [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area     [m]
 
@@ -351,7 +356,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 +367,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 +400,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(15), ii
+      INTEGER :: ierr(16), ii
       !!-----------------------------------------------------------------
       ierr(:) = 0
@@ -440,5 +451,5 @@
 
       ii = ii + 1
-      ALLOCATE( at_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
+      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
 
       ! * Old values of global variables
@@ -461,4 +472,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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icecor.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icecor.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icecor.F90	(revision 11586)
@@ -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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icectl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icectl.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icectl.F90	(revision 11586)
@@ -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,15 @@
    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 year  spuriously gained/lost
+   REAL(wp), PARAMETER ::   zchk_s   = 1.e-4   ! g/m2/s  <=> 1mm of ice per year  spuriously gained/lost (considering s=10g/kg)
+   REAL(wp), PARAMETER ::   zchk_t   = 3.      ! W/m2    <=> 1mm of ice per year  spuriously gained/lost (considering Lf=3e5J/kg)
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -59,5 +68,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 +77,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_glo * zarea ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice
+            IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice
+            IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * 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_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
          ENDIF
          !
@@ -179,5 +165,4 @@
    END SUBROUTINE ice_cons_hsm
 
-
    SUBROUTINE ice_cons_final( cd_routine )
       !!-------------------------------------------------------------------
@@ -188,40 +173,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_glo * zarea ) &
+            &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice
+         IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) &
+            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice
+         !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * 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) ) > zchk_m * rn_icechk_cel )   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) ) > zchk_s * rn_icechk_cel )   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) ) > zchk_t * rn_icechk_cel )   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 +366,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 +384,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 +396,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 +424,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 +439,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 +447,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 +472,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 +479,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 +489,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 +519,4 @@
    END SUBROUTINE ice_ctl
  
-   
    SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 )
       !!-------------------------------------------------------------------
@@ -443,5 +526,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 +565,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 +586,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 +599,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 +697,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 +757,5 @@
       
    END SUBROUTINE ice_prt3D
-
+      
 #else
    !!----------------------------------------------------------------------
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedia.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedia.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedia.F90	(revision 11586)
@@ -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_cel, rn_icechk_glo, ln_icediahsb, ln_icectl, iiceprt, jiceprt  
       !!----------------------------------------------------------------------
       !
@@ -191,8 +191,10 @@
          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 (gridcell)    rn_icechk_cel = ', rn_icechk_cel
+         WRITE(numout,*) '         threshold for conservation (global)      rn_icechk_glo = ', rn_icechk_glo
+         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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rdgrft.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rdgrft.F90	(revision 11586)
@@ -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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg.F90	(revision 11586)
@@ -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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90	(revision 11586)
@@ -121,5 +121,5 @@
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
-      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
       REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
@@ -130,5 +130,4 @@
       REAL(wp) ::   zshear, zdum1, zdum2
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
@@ -244,12 +243,4 @@
       CALL ice_strength
 
-      ! scale factors
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
-            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
-         END DO
-      END DO
-
       ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
       IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
@@ -280,9 +271,6 @@
 
             ! Ocean currents at U-V points
-            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
-               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
-            
-            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
-               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
+            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
 
             ! Coriolis at T points (m*f)
@@ -459,11 +447,8 @@
                   &                  ) * r1_e1e2v(ji,jj)
                !
-               !                !--- u_ice at V point
-               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+               !                !--- ice currents at U-V point
+               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
+               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
                !
-               !                !--- v_ice at U point
-               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
-                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
             END DO
          END DO
@@ -494,12 +479,13 @@
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-                  !
-                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
-                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
@@ -508,5 +494,5 @@
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
-                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
@@ -544,19 +530,20 @@
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-                  !
-                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
-                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
-                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
-                        &            )   * zmsk00x(ji,jj)
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
-                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
@@ -596,19 +583,20 @@
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-                  !
-                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
-                        &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
-                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
-                        &            )   * zmsk00x(ji,jj)
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
-                        &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
@@ -646,12 +634,13 @@
                   !
                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-                  !
-                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
-                        &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
@@ -660,5 +649,5 @@
                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
-                        &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceistate.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceistate.F90	(revision 11586)
@@ -101,5 +101,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
       !!
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
       !--------------------------------------------------------------------
 
@@ -209,8 +209,11 @@
             ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
             !
-            ! ponds
+            ! pond concentration
             IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
-               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
+               &                              * si(jp_ati)%fnow(:,:,1) 
             zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
+            !
+            ! pond depth
             IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
                &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
@@ -238,5 +241,5 @@
                zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
                ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
-               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:)
+               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
                zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
             ELSEWHERE
@@ -248,9 +251,16 @@
                zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
                ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
-               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:)
+               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
                zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
             END WHERE
             !
          ENDIF
+
+         ! make sure ponds = 0 if no ponds scheme
+         IF ( .NOT.ln_pnd ) THEN
+            zapnd_ini(:,:) = 0._wp
+            zhpnd_ini(:,:) = 0._wp
+         ENDIF
+         
          !-------------!
          ! fill fields !
@@ -275,12 +285,16 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
          CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini )
 
          ! allocate temporary arrays
          ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
-            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl) )
+            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
          
          ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
-         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                   zhi_2d, zhs_2d, zai_2d , &
-            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti),   zti_2d, zts_2d, ztsu_2d, zsi_2d )
+         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
+            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   &
+            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
+            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d )
 
          ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
@@ -289,26 +303,18 @@
             zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
          END DO
-         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d  , h_i  )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d  , h_s  )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d  , a_i  )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d  , zti_3d )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d  , zts_3d )
-         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
-         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d  , s_i  )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d )
+         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
 
          ! deallocate temporary arrays
          DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
-            &        zti_2d, zts_2d, ztsu_2d, zsi_2d )
-
-         ! Melt ponds: distribute uniformely over the categories
-         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN
-            DO jl = 1, jpl
-               a_ip_frac(:,:,jl) = zapnd_ini(:,:)
-               h_ip     (:,:,jl) = zhpnd_ini(:,:)
-               a_ip     (:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
-               v_ip     (:,:,jl) = h_ip     (:,:,jl) * a_ip(:,:,jl)
-            END DO
-         ENDIF
-          
+            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
+
          ! calculate extensive and intensive variables
          CALL ice_var_salprof ! for sz_i
@@ -350,4 +356,12 @@
          END DO
 
+         ! Melt ponds
+         WHERE( a_i > epsi10 )
+            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
+         ELSEWHERE
+            a_ip_frac(:,:,:) = 0._wp
+         END WHERE
+         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+          
          ! specific temperatures for coupled runs
          tn_ice(:,:,:) = t_su(:,:,:)
@@ -512,4 +526,10 @@
       ENDIF
       !
+      IF( .NOT.ln_pnd ) THEN
+         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
+         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
+         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
+      ENDIF
+      !
    END SUBROUTINE ice_istate_init
 
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceitd.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceitd.F90	(revision 11586)
@@ -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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd.F90	(revision 11586)
@@ -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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_do.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_do.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_do.F90	(revision 11586)
@@ -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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_pnd.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_pnd.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_pnd.F90	(revision 11586)
@@ -205,5 +205,5 @@
       INTEGER  ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
       !!-------------------------------------------------------------------
       !
@@ -221,17 +221,19 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namicethd_pnd:'
-         WRITE(numout,*) '      Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
-         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
-         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
-         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
-         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
+         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd
+         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
+         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
+         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
+         WRITE(numout,*) '            Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
+         WRITE(numout,*) '         Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
       ENDIF
       !
       !                             !== set the choice of ice pond scheme ==!
       ioptio = 0
-                                                            nice_pnd = np_pndNO
-      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
-      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
-      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
+      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
+      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
+      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
+      IF( ioptio /= 1 )   &
+         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
       !
       SELECT CASE( nice_pnd )
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90	(revision 11586)
@@ -47,8 +47,8 @@
    !!   ice_var_zapneg    : remove negative ice fields
    !!   ice_var_roundoff  : remove negative values arising from roundoff erros
-   !!   ice_var_itd       : convert N-cat to M-cat
    !!   ice_var_bv        : brine volume
    !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
    !!   ice_var_sshdyn    : compute equivalent ssh in lead
+   !!   ice_var_itd       : convert N-cat to M-cat
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -115,18 +115,13 @@
       !
       ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
-      !
-      ALLOCATE( z1_at_i(jpi,jpj) )
-      WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
-      ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
-      END WHERE
-      !
-      tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
-      WHERE( at_i(:,:)<=epsi20 ); tm_su(:,:) = rt0; END WHERE      
-      !
+
       ! The following fields are calculated for diagnostics and outputs only
       ! ==> Do not use them for other purposes
       IF( kn > 1 ) THEN
          !
-         ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
+         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
+         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
+         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
+         END WHERE
          WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
          ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
@@ -141,4 +136,5 @@
          !         
          !                          ! mean temperature (K), salinity and age
+         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
@@ -158,13 +154,18 @@
          !                           ! put rt0 where there is no ice
          WHERE( at_i(:,:)<=epsi20 )
+            tm_su(:,:) = rt0
             tm_si(:,:) = rt0
             tm_i (:,:) = rt0
             tm_s (:,:) = rt0
          END WHERE
-
-         DEALLOCATE( z1_vt_i , z1_vt_s )
+         !
+         !                           ! mean melt pond depth
+         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)
+         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp
+         END WHERE         
+         !
+         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
+         !
       ENDIF
-      !
-      DEALLOCATE( z1_at_i )
       !
    END SUBROUTINE ice_var_agg
@@ -664,5 +665,5 @@
       WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
       WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
-      IF ( ln_pnd_H12 ) THEN
+      IF( ln_pnd_H12 ) THEN
          WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
          WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
@@ -785,6 +786,6 @@
    !! ** Purpose :  converting N-cat ice to jpl ice categories
    !!-------------------------------------------------------------------
-   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
+   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting 1-cat ice to 1 ice category
@@ -792,6 +793,6 @@
       REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
-      REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
       !!-------------------------------------------------------------------
       ! == thickness and concentration == !
@@ -800,14 +801,16 @@
       pa_i(:) = pati(:)
       !
-      ! == temperature and salinity == !
-      IF( PRESENT( pt_i  ) )   pt_i (:) = ptmi (:)
-      IF( PRESENT( pt_s  ) )   pt_s (:) = ptms (:)
-      IF( PRESENT( pt_su ) )   pt_su(:) = ptmsu(:)
-      IF( PRESENT( ps_i  ) )   ps_i (:) = psmi (:)
+      ! == temperature and salinity and ponds == !
+      pt_i (:) = ptmi (:)
+      pt_s (:) = ptms (:)
+      pt_su(:) = ptmsu(:)
+      ps_i (:) = psmi (:)
+      pa_ip(:) = patip(:)
+      ph_ip(:) = phtip(:)
       
    END SUBROUTINE ice_var_itd_1c1c
 
-   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
+   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting N-cat ice to 1 ice category
@@ -815,6 +818,6 @@
       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
-      REAL(wp), DIMENSION(:)  , INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
@@ -826,5 +829,5 @@
       !
       ! == thickness and concentration == !
-      ALLOCATE( z1_ai(idim) )
+      ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
       !
       pa_i(:) = SUM( pati(:,:), dim=2 )
@@ -838,30 +841,27 @@
       !
       ! == temperature and salinity == !
-      IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN
-         !
-         ALLOCATE( z1_vi(idim), z1_vs(idim) )
-         !
-         WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
-         ELSEWHERE                                 ;   z1_vi(:) = 0._wp
-         END WHERE
-         WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
-         ELSEWHERE                                 ;   z1_vs(:) = 0._wp
-         END WHERE
-         !
-         IF( PRESENT( pt_i  ) )   pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
-         IF( PRESENT( pt_s  ) )   pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
-         IF( PRESENT( pt_su ) )   pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
-         IF( PRESENT( ps_i  ) )   ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
-         !
-         DEALLOCATE( z1_vi, z1_vs )
-         !
-      ENDIF
-      !
-      DEALLOCATE( z1_ai )
+      WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
+      ELSEWHERE                                 ;   z1_vi(:) = 0._wp
+      END WHERE
+      WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
+      ELSEWHERE                                 ;   z1_vs(:) = 0._wp
+      END WHERE
+      pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+      pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
+      pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
+      ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+
+      ! == ponds == !
+      pa_ip(:) = SUM( patip(:,:), dim=2 )
+      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+      ELSEWHERE                    ;   ph_ip(:) = 0._wp
+      END WHERE
+      !
+      DEALLOCATE( z1_ai, z1_vi, z1_vs )
       !
    END SUBROUTINE ice_var_itd_Nc1c
    
-   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
+   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
       !!-------------------------------------------------------------------
       !!
@@ -894,8 +894,9 @@
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:)  , INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
-      REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
       !
       INTEGER , DIMENSION(4) ::   itest
+      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra
       INTEGER  ::   ji, jk, jl
       INTEGER  ::   idim, i_fill, jl0  
@@ -1010,29 +1011,36 @@
       !
       ! == temperature and salinity == !
-      IF( PRESENT( pt_i  ) ) THEN
-         DO jl = 1, jpl
-            pt_i(:,jl) = ptmi (:)
-         END DO
-      ENDIF
-      IF( PRESENT( pt_s  ) ) THEN
-         DO jl = 1, jpl
-            pt_s (:,jl) = ptms (:)
-         END DO
-      ENDIF
-      IF( PRESENT( pt_su ) ) THEN
-         DO jl = 1, jpl
-            pt_su(:,jl) = ptmsu(:)
-         END DO
-      ENDIF
-      IF( PRESENT( ps_i  ) ) THEN
-         DO jl = 1, jpl
-            ps_i (:,jl) = psmi (:)
-         END DO
-      ENDIF
+      DO jl = 1, jpl
+         pt_i (:,jl) = ptmi (:)
+         pt_s (:,jl) = ptms (:)
+         pt_su(:,jl) = ptmsu(:)
+         ps_i (:,jl) = psmi (:)
+         ps_i (:,jl) = psmi (:)         
+      END DO
+      !
+      ! == ponds == !
+      ALLOCATE( zfra(idim) )
+      ! keep the same pond fraction atip/ati for each category
+      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:)
+      ELSEWHERE                   ;   zfra(:) = 0._wp
+      END WHERE
+      DO jl = 1, jpl
+         pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
+      END DO
+      ! keep the same v_ip/v_i ratio for each category
+      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
+      ELSEWHERE                                 ;   zfra(:) = 0._wp
+      END WHERE
+      DO jl = 1, jpl
+         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
+         END WHERE
+      END DO
+      DEALLOCATE( zfra )
       !
    END SUBROUTINE ice_var_itd_1cMc
 
-   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
+   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
       !!-------------------------------------------------------------------
       !!
@@ -1065,10 +1073,10 @@
       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
-      REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
       !
       INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
       INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin
-      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp
+      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra
       !
       REAL(wp), PARAMETER ::   ztrans = 0.25_wp
@@ -1088,21 +1096,25 @@
          pa_i(:,:) = pati(:,:)
          !
-         ! == temperature and salinity == !
-         IF( PRESENT( pt_i  ) )   pt_i (:,:) = ptmi (:,:)
-         IF( PRESENT( pt_s  ) )   pt_s (:,:) = ptms (:,:)
-         IF( PRESENT( pt_su ) )   pt_su(:,:) = ptmsu(:,:)
-         IF( PRESENT( ps_i  ) )   ps_i (:,:) = psmi (:,:)
+         ! == temperature and salinity and ponds == !
+         pt_i (:,:) = ptmi (:,:)
+         pt_s (:,:) = ptms (:,:)
+         pt_su(:,:) = ptmsu(:,:)
+         ps_i (:,:) = psmi (:,:)
+         pa_ip(:,:) = patip(:,:)
+         ph_ip(:,:) = phtip(:,:)
          !                              ! ---------------------- !
       ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
          !                              ! ---------------------- !
-         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:) )
-!!         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), &
-!!            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) )
+         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
+            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), &
+            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &
+            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  )
          !                              ! ---------------------- !
       ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
          !                              ! ---------------------- !
-         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) )
-!!         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), &
-!!            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )         
+         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
+            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), &
+            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &
+            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  )
          !                              ! ----------------------- !
       ELSE                              ! input cat /= output cat !
@@ -1194,47 +1206,57 @@
          ! == temperature and salinity == !
          !
-         IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN
-            !
-            ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
-            !
-            WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
-            ELSEWHERE                                               ;   z1_ai(:) = 0._wp
+         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
+         !
+         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
+         ELSEWHERE                                               ;   z1_ai(:) = 0._wp
+         END WHERE
+         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
+         ELSEWHERE                                               ;   z1_vi(:) = 0._wp
+         END WHERE
+         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
+         ELSEWHERE                                               ;   z1_vs(:) = 0._wp
+         END WHERE
+         !
+         ! fill all the categories with the same value
+         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+         DO jl = 1, jpl
+            pt_i (:,jl) = ztmp(:)
+         END DO
+         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
+         DO jl = 1, jpl
+            pt_s (:,jl) = ztmp(:)
+         END DO
+         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
+         DO jl = 1, jpl
+            pt_su(:,jl) = ztmp(:)
+         END DO
+         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+         DO jl = 1, jpl
+            ps_i (:,jl) = ztmp(:)
+         END DO
+         !
+         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
+         !
+         ! == ponds == !
+         ALLOCATE( zfra(idim) )
+         ! keep the same pond fraction atip/ati for each category
+         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
+         ELSEWHERE                                   ;   zfra(:) = 0._wp
+         END WHERE
+         DO jl = 1, jpl
+            pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
+         END DO
+         ! keep the same v_ip/v_i ratio for each category
+         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
+            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
+         ELSEWHERE
+            zfra(:) = 0._wp
+         END WHERE
+         DO jl = 1, jpl
+            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
             END WHERE
-            WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
-            ELSEWHERE                                               ;   z1_vi(:) = 0._wp
-            END WHERE
-            WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
-            ELSEWHERE                                               ;   z1_vs(:) = 0._wp
-            END WHERE
-            !
-            ! fill all the categories with the same value
-            IF( PRESENT( pt_i  ) ) THEN
-               ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
-               DO jl = 1, jpl
-                  pt_i (:,jl) = ztmp(:)
-               END DO
-            ENDIF
-            IF( PRESENT( pt_s  ) ) THEN
-               ztmp(:) =  SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
-               DO jl = 1, jpl
-                  pt_s (:,jl) = ztmp(:)
-               END DO
-            ENDIF
-            IF( PRESENT( pt_su ) ) THEN
-               ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
-               DO jl = 1, jpl
-                  pt_su(:,jl) = ztmp(:)
-               END DO
-            ENDIF
-            IF( PRESENT( ps_i  ) ) THEN
-               ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
-               DO jl = 1, jpl
-                  ps_i (:,jl) = ztmp(:)
-               END DO
-            ENDIF
-            !
-            DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
-            !
-         ENDIF
+         END DO
+         DEALLOCATE( zfra )
          !
       ENDIF
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icewri.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icewri.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icewri.F90	(revision 11586)
@@ -114,4 +114,5 @@
       ! melt ponds
       IF( iom_use('iceapnd' ) )   CALL iom_put( 'iceapnd', at_ip  * zmsk00      )                                           ! melt pond total fraction
+      IF( iom_use('icehpnd' ) )   CALL iom_put( 'icehpnd', hm_ip  * zmsk00      )                                           ! melt pond depth
       IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area
       ! salt
