Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11586)
@@ -563,8 +563,8 @@
             DO ji = 2, jpim1  
                
-			   zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
-			   zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
+               zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
+               zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
 			   
-			   ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
+               ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
                   &                      +    zrhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
                   &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) )
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11586)
@@ -214,11 +214,12 @@
             ! Check that restoring coefficients are between 0 and 1
             !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   &
-            IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   &
+            !IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   &
+            IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp )   &
                &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' )
             !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   &
-            IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   &
+            IF( zcff  - nn_fsbc > 0.001_wp .OR. zcff  < 0._wp )   &
                &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' )
             IF( zcff  > zcff1                    )   &
-               &                   CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' )
+               &                   CALL ctl_stop( 'abl_init : rn_ldyn_max must be smaller than rn_ldyn_min' )
          END IF
       END IF
@@ -234,11 +235,11 @@
          ! Check that restoring coefficients are between 0 and 1
          !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   &
-         IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   &
+         IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp )   &
             &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' )
          !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   &
-         IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   &
+         IF( zcff  - nn_fsbc > 0.001_wp .OR. zcff  < 0._wp )   &
             &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' )
          IF( zcff  > zcff1                    )   &
-            &                   CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' )
+            &                   CALL ctl_stop( 'abl_init : rn_ltra_max must be smaller than rn_ltra_min' )
       END IF
 
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
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdy_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdy_oce.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdy_oce.F90	(revision 11586)
@@ -57,4 +57,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::  h_i    !: Now ice  thickness climatology
       REAL(wp), POINTER, DIMENSION(:,:) ::  h_s    !: now snow thickness
+      REAL(wp), POINTER, DIMENSION(:,:) ::  t_i    !: now ice  temperature
+      REAL(wp), POINTER, DIMENSION(:,:) ::  t_s    !: now snow temperature
+      REAL(wp), POINTER, DIMENSION(:,:) ::  tsu    !: now surf temperature
+      REAL(wp), POINTER, DIMENSION(:,:) ::  s_i    !: now ice  salinity
+      REAL(wp), POINTER, DIMENSION(:,:) ::  aip    !: now ice  pond concentration
+      REAL(wp), POINTER, DIMENSION(:,:) ::  hip    !: now ice  pond depth
 #if defined key_top
       CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
@@ -68,4 +74,5 @@
    !! Namelist variables
    !!----------------------------------------------------------------------
+   !                                                   !!** nambdy **
    LOGICAL, PUBLIC            ::   ln_bdy                   !: Unstructured Ocean Boundary Condition
 
@@ -101,9 +108,12 @@
    INTEGER , DIMENSION(jp_bdy)          ::   nn_ice_dta     !: = 0 use the initial state as bdy dta ; 
                                                             !: = 1 read it in a NetCDF file
-   REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_tem              !: choice of the temperature of incoming sea ice
-   REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_sal              !: choice of the salinity    of incoming sea ice
-   REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_age              !: choice of the age         of incoming sea ice
+   ! 
+   !                                                   !!** nambdy_dta **
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_tem                !: temperature of incoming sea ice
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_sal                !: salinity    of incoming sea ice
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_age                !: age         of incoming sea ice
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_apnd               !: pond conc.  of incoming sea ice
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_hpnd               !: pond thick. of incoming sea ice
    !
-   
    !!----------------------------------------------------------------------
    !! Global variables
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdydta.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdydta.F90	(revision 11586)
@@ -43,5 +43,5 @@
    PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
 
-   INTEGER , PARAMETER ::   jpbdyfld  = 10    ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jpbdyfld  = 16    ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_bdyssh = 1     ! 
    INTEGER , PARAMETER ::   jp_bdyu2d = 2     ! 
@@ -53,9 +53,15 @@
    INTEGER , PARAMETER ::   jp_bdya_i = 8     ! 
    INTEGER , PARAMETER ::   jp_bdyh_i = 9     ! 
-   INTEGER , PARAMETER ::   jp_bdyh_S = 10    ! 
+   INTEGER , PARAMETER ::   jp_bdyh_s = 10    ! 
+   INTEGER , PARAMETER ::   jp_bdyt_i = 11    ! 
+   INTEGER , PARAMETER ::   jp_bdyt_s = 12    ! 
+   INTEGER , PARAMETER ::   jp_bdytsu = 13    ! 
+   INTEGER , PARAMETER ::   jp_bdys_i = 14    ! 
+   INTEGER , PARAMETER ::   jp_bdyaip = 15    ! 
+   INTEGER , PARAMETER ::   jp_bdyhip = 16    ! 
 #if ! defined key_si3
    INTEGER , PARAMETER ::   jpl = 1
 #endif
-                                                             ! =F => baroclinic velocities in 3D boundary conditions
+
 !$AGRIF_DO_NOT_TREAT
    TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::   bf   ! structure of input fields (file informations, fields read)
@@ -181,7 +187,14 @@
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
-                        dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
-                        dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%a_i(ib,jl) =  a_i (ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%h_i(ib,jl) =  h_i (ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%h_s(ib,jl) =  h_s (ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%t_i(ib,jl) =  SUM(t_i (ii,ij,:,jl)) * r1_nlay_i * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%t_s(ib,jl) =  SUM(t_s (ii,ij,:,jl)) * r1_nlay_s * tmask(ii,ij,1)
+                        dta_bdy(jbdy)%tsu(ib,jl) =  t_su(ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%s_i(ib,jl) =  s_i (ii,ij,jl) * tmask(ii,ij,1)
+                        ! melt ponds
+                        dta_bdy(jbdy)%aip(ib,jl) =  a_ip(ii,ij,jl) * tmask(ii,ij,1) 
+                        dta_bdy(jbdy)%hip(ib,jl) =  h_ip(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -280,10 +293,43 @@
 
 #if defined key_si3
-         ! ice: convert N-cat fields (input) into jpl-cat (output)
          IF( dta_alias%lneed_ice ) THEN
-            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)
+            ! fill temperature and salinity arrays
+            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem (jbdy)
+            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_s)%fnow(:,1,:) = rice_tem (jbdy)
+            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' )   bf_alias(jp_bdytsu)%fnow(:,1,:) = rice_tem (jbdy)
+            IF( TRIM(bf_alias(jp_bdys_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdys_i)%fnow(:,1,:) = rice_sal (jbdy)
+            IF( TRIM(bf_alias(jp_bdyaip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyaip)%fnow(:,1,:) = rice_apnd(jbdy) * & ! rice_apnd is the pond fraction
+               &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i )
+            IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy)
+            ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
+            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdytsu)%fnow(:,1,:) + 271.15 )
+            ! if T_su is read and not T_s, set T_s = T_su
+            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdytsu)%fnow(:,1,:)
+            ! if T_s is read and not T_su, set T_su = T_s
+            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_s)%fnow(:,1,:)
+            ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
+            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
+               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdyt_s)%fnow(:,1,:) + 271.15 )
+
+            ! make sure ponds = 0 if no ponds scheme
+            IF ( .NOT.ln_pnd ) THEN
+               bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp
+               bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp
+            ENDIF
+            
+            ! convert N-cat fields (input) into jpl-cat (output)
+            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)            
             IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output)
-               CALL ice_var_itd(bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), &
-                  &              dta_alias%h_i               , dta_alias%h_s               , dta_alias%a_i                 )
+               CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), &
+                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , &
+                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &
+                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &
+                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), &
+                  &              dta_alias%t_i                  , dta_alias%t_s                  , &
+                  &              dta_alias%tsu                  , dta_alias%s_i                  , &
+                  &              dta_alias%aip                  , dta_alias%hip )
             ENDIF
          ENDIF
@@ -332,4 +378,5 @@
       !                                                         ! =F => baroclinic velocities in 3D boundary data
       LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta
+      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
       INTEGER                                ::   ipk,ipl       !
       INTEGER                                ::   idvar         ! variable ID
@@ -341,12 +388,13 @@
       LOGICAL                                ::   llneed        !
       LOGICAL                                ::   llread        !
-      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
-      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
-      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_a_i, bn_h_i, bn_h_s      
+      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
+      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
+      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip       
       TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill
       TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias
       !
       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
-      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s
+      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip
+      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd
       NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp
       !!---------------------------------------------------------------------------
@@ -402,4 +450,19 @@
          ENDIF
 
+#if defined key_si3
+         IF( .NOT.ln_pnd ) THEN
+            rn_ice_apnd = 0. ; rn_ice_hpnd = 0.
+            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' )
+         ENDIF
+#endif
+
+         ! temp, salt, age and ponds of incoming ice
+         rice_tem (jbdy) = rn_ice_tem
+         rice_sal (jbdy) = rn_ice_sal
+         rice_age (jbdy) = rn_ice_age
+         rice_apnd(jbdy) = rn_ice_apnd
+         rice_hpnd(jbdy) = rn_ice_hpnd
+         
+         
          DO jfld = 1, jpbdyfld
 
@@ -497,33 +560,57 @@
             !          ice
             ! =====================
+            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. &
+               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. &
+               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip      ) THEN
+               igrd = 1                                                    ! T point
+               ipk = ipl                                                   ! jpl-cat data
+               llneed = dta_bdy(jbdy)%lneed_ice                            ! ice will be needed
+               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
+               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
+            ENDIF
             IF( jfld == jp_bdya_i ) THEN
                cl3 = 'a_i'
-               igrd = 1                                                    ! T point
-               ipk = ipl                                                   ! 
-               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%a_i will be needed
-               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
-               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
-               bn_alias => bn_a_i                                          ! alias for ssh structure of nambdy_dta 
-               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
-           ENDIF
+               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for a_i structure of bdy number jbdy
+               bn_alias => bn_a_i                                          ! alias for a_i structure of nambdy_dta 
+            ENDIF
             IF( jfld == jp_bdyh_i ) THEN
                cl3 = 'h_i'
-               igrd = 1                                                    ! T point
-               ipk = ipl                                                   ! 
-               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_i will be needed
-               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
-               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
-               bn_alias => bn_h_i                                          ! alias for ssh structure of nambdy_dta 
-               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
+               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for h_i structure of bdy number jbdy
+               bn_alias => bn_h_i                                          ! alias for h_i structure of nambdy_dta 
             ENDIF
             IF( jfld == jp_bdyh_s ) THEN
                cl3 = 'h_s'
-               igrd = 1                                                    ! T point
-               ipk = ipl                                                   ! 
-               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_s will be needed
-               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
-               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
-               bn_alias => bn_h_s                                          ! alias for ssh structure of nambdy_dta 
-               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
+               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for h_s structure of bdy number jbdy
+               bn_alias => bn_h_s                                          ! alias for h_s structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyt_i ) THEN
+               cl3 = 't_i'
+               bf_alias => bf(jp_bdyt_i,jbdy:jbdy)                         ! alias for t_i structure of bdy number jbdy
+               bn_alias => bn_t_i                                          ! alias for t_i structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyt_s ) THEN
+               cl3 = 't_s'
+               bf_alias => bf(jp_bdyt_s,jbdy:jbdy)                         ! alias for t_s structure of bdy number jbdy
+               bn_alias => bn_t_s                                          ! alias for t_s structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdytsu ) THEN
+               cl3 = 'tsu'
+               bf_alias => bf(jp_bdytsu,jbdy:jbdy)                         ! alias for tsu structure of bdy number jbdy
+               bn_alias => bn_tsu                                          ! alias for tsu structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdys_i ) THEN
+               cl3 = 's_i'
+               bf_alias => bf(jp_bdys_i,jbdy:jbdy)                         ! alias for s_i structure of bdy number jbdy
+               bn_alias => bn_s_i                                          ! alias for s_i structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyaip ) THEN
+               cl3 = 'aip'
+               bf_alias => bf(jp_bdyaip,jbdy:jbdy)                         ! alias for aip structure of bdy number jbdy
+               bn_alias => bn_aip                                          ! alias for aip structure of nambdy_dta 
+            ENDIF
+            IF( jfld == jp_bdyhip ) THEN
+               cl3 = 'hip'
+               bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy
+               bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta 
             ENDIF
 
@@ -542,11 +629,11 @@
 
                ! associate the pointer and get rid of the dimensions with a size equal to 1
-               IF( jfld == jp_bdyssh           ) dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1)
-               IF( jfld == jp_bdyu2d           ) dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1)
-               IF( jfld == jp_bdyv2d           ) dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1)
-               IF( jfld == jp_bdyu3d           ) dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:)
-               IF( jfld == jp_bdyv3d           ) dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:)
-               IF( jfld == jp_bdytem           ) dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:)
-               IF( jfld == jp_bdysal           ) dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdyssh )        dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1)
+               IF( jfld == jp_bdyu2d )        dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1)
+               IF( jfld == jp_bdyv2d )        dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1)
+               IF( jfld == jp_bdyu3d )        dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdyv3d )        dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdytem )        dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:)
+               IF( jfld == jp_bdysal )        dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:)
                IF( jfld == jp_bdya_i ) THEN
                   IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:)
@@ -562,4 +649,34 @@
                   IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:)
                   ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdyt_i ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_i => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_i(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdyt_s ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_s => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_s(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdytsu ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%tsu => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%tsu(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdys_i ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%s_i => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%s_i(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdyaip ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%aip => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%aip(iszdim,jpl) )
+                  ENDIF
+               ENDIF
+               IF( jfld == jp_bdyhip ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hip => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hip(iszdim,jpl) )
                   ENDIF
                ENDIF
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyice.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyice.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyice.F90	(revision 11586)
@@ -63,4 +63,5 @@
       IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       !
       CALL ice_var_glo2eqv
@@ -109,6 +110,7 @@
       !
       ! controls
+      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
+      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing
       !
@@ -152,8 +154,22 @@
             zwgt  = idx%nbw(i_bdy,jgrd)
             zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd)
-            a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction 
-            h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth 
-            h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
-
+            a_i (ji,jj,  jl) = ( a_i (ji,jj,  jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  concentration 
+            h_i (ji,jj,  jl) = ( h_i (ji,jj,  jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  depth 
+            h_s (ji,jj,  jl) = ( h_s (ji,jj,  jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
+            t_i (ji,jj,:,jl) = ( t_i (ji,jj,:,jl) * zwgt1 + dta%t_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  temperature
+            t_s (ji,jj,:,jl) = ( t_s (ji,jj,:,jl) * zwgt1 + dta%t_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow temperature
+            t_su(ji,jj,  jl) = ( t_su(ji,jj,  jl) * zwgt1 + dta%tsu(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Surf temperature
+            s_i (ji,jj,  jl) = ( s_i (ji,jj,  jl) * zwgt1 + dta%s_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  salinity
+            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration
+            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth
+            !
+            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
+            !
+            ! make sure ponds = 0 if no ponds scheme
+            IF( .NOT.ln_pnd ) THEN
+               a_ip(ji,jj,jl) = 0._wp
+               h_ip(ji,jj,jl) = 0._wp
+            ENDIF
+            !
             ! -----------------
             ! Pathological case
@@ -170,5 +186,5 @@
             h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
             h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos ) 
-
+            !
          ENDDO
       ENDDO
@@ -206,33 +222,22 @@
             IF( a_i(ib,jb,jl) > 0._wp ) THEN   ! there is ice at the boundary
                !
-               a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration
-               h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice
-               h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw
-               !
-               SELECT CASE( jpbound )
-                  !
-               CASE( 0 )   ! velocity is inward
-                  !
-                  oa_i(ji,jj,  jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age
-                  a_ip(ji,jj,  jl) = 0._wp                            ! pond concentration
-                  v_ip(ji,jj,  jl) = 0._wp                            ! pond volume
-                  t_su(ji,jj,  jl) = rn_ice_tem(jbdy)                 ! temperature surface
-                  t_s (ji,jj,:,jl) = rn_ice_tem(jbdy)                 ! temperature snw
-                  t_i (ji,jj,:,jl) = rn_ice_tem(jbdy)                 ! temperature ice
-                  s_i (ji,jj,  jl) = rn_ice_sal(jbdy)                 ! salinity
-                  sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy)                 ! salinity profile
-                  !
-               CASE( 1 )   ! velocity is outward
-                  !
-                  oa_i(ji,jj,  jl) = oa_i(ib,jb,  jl) ! age
-                  a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) ! pond concentration
-                  v_ip(ji,jj,  jl) = v_ip(ib,jb,  jl) ! pond volume
-                  t_su(ji,jj,  jl) = t_su(ib,jb,  jl) ! temperature surface
-                  t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw
-                  t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice
-                  s_i (ji,jj,  jl) = s_i (ib,jb,  jl) ! salinity
-                  sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile
-                  !
-               END SELECT
+               a_i (ji,jj,  jl) = a_i (ib,jb,  jl)
+               h_i (ji,jj,  jl) = h_i (ib,jb,  jl)
+               h_s (ji,jj,  jl) = h_s (ib,jb,  jl)
+               t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl)
+               t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl)
+               t_su(ji,jj,  jl) = t_su(ib,jb,  jl)
+               s_i (ji,jj,  jl) = s_i (ib,jb,  jl)
+               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl)
+               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl)
+               !
+               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
+               !
+               ! ice age
+               IF    ( jpbound == 0 ) THEN  ! velocity is inward
+                  oa_i(ji,jj,jl) = rice_age(jbdy) * a_i(ji,jj,jl)
+               ELSEIF( jpbound == 1 ) THEN  ! velocity is outward
+                  oa_i(ji,jj,jl) = oa_i(ib,jb,jl)
+               ENDIF
                !
                IF( nn_icesal == 1 ) THEN     ! if constant salinity
@@ -259,4 +264,12 @@
                END DO
                !
+               ! melt ponds
+               IF( a_i(ji,jj,jl) > epsi10 ) THEN
+                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl)
+               ELSE
+                  a_ip_frac(ji,jj,jl) = 0._wp
+               ENDIF
+               v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl)
+               !
             ELSE   ! no ice at the boundary
                !
@@ -270,4 +283,9 @@
                t_s (ji,jj,:,jl) = rt0
                t_i (ji,jj,:,jl) = rt0 
+
+               a_ip_frac(ji,jj,jl) = 0._wp
+               h_ip     (ji,jj,jl) = 0._wp
+               a_ip     (ji,jj,jl) = 0._wp
+               v_ip     (ji,jj,jl) = 0._wp
                
                IF( nn_icesal == 1 ) THEN     ! if constant salinity
@@ -372,5 +390,5 @@
                      jj    = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
                      zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd)
-                     !    ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨     !  ¨¨¨¨ïce¨¨¨(jj+1)¨¨     ! ¨¨¨¨¨¨ö¨¨¨¨(jj+1)       
+                     !                         !      ice   (jj+1)       !       o    (jj+1)
                      !       ^    (jj  )       !       ^    (jj  )       !       ^    (jj  )       
                      !      ice   (jj  )       !       o    (jj  )       !       o    (jj  )       
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyini.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyini.F90	(revision 11586)
@@ -67,5 +67,4 @@
          &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
          &             cn_ice, nn_ice_dta,                                     &
-         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
          &             ln_vol, nn_volctl, nn_rimwidth
          !
@@ -94,7 +93,4 @@
       cn_ice         (2:jp_bdy) = cn_ice         (1)
       nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
-      rn_ice_tem     (2:jp_bdy) = rn_ice_tem     (1)
-      rn_ice_sal     (2:jp_bdy) = rn_ice_sal     (1)
-      rn_ice_age     (2:jp_bdy) = rn_ice_age     (1)
       REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
       READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
@@ -328,8 +324,4 @@
             CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
             END SELECT
-            WRITE(numout,*)
-            WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)         
-            WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)         
-            WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)         
          ENDIF
 #else
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diacfl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diacfl.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diacfl.F90	(revision 11586)
@@ -29,9 +29,4 @@
    REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
 
-!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
-!!gm          I don't understand why.
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
-!!gm end
-
    PUBLIC   dia_cfl       ! routine called by step.F90
    PUBLIC   dia_cfl_init  ! routine called by nemogcm
@@ -55,8 +50,8 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER                ::   ji, jj, jk                            ! dummy loop indices
-      REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
-      INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace
-!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
+      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
+      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars
+      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
       !!----------------------------------------------------------------------
       !
@@ -71,5 +66,5 @@
       DO jk = 1, jpk       ! calculate Courant numbers
          DO jj = 1, jpj
-            DO ji = 1, fs_jpim1   ! vector opt.
+            DO ji = 1, jpi
                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
@@ -111,5 +106,5 @@
       !                    ! write out to file
       IF( lwp ) THEN
-         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
@@ -172,8 +167,4 @@
       rCw_max = 0._wp
       !
-!!gm required to work
-      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
-!!gm end
-      !      
    END SUBROUTINE dia_cfl_init
 
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11586)
@@ -71,10 +71,9 @@
 #endif
    INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
-   INTEGER ::          nb_T, ndim_bT                 ! grid_T file
+   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
    INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
    INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
    INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
-   INTEGER ::   ndim_A, ndim_hA                      ! ABL file   
-   INTEGER ::   nid_A, nz_A, nh_A                    ! grid_ABL file   
+   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file   
    INTEGER ::   ndex(1)                              ! ???
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
@@ -216,4 +215,6 @@
       ENDIF
 
+      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
+      !
       CALL iom_put( "woce", wn )                   ! vertical velocity
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
@@ -226,4 +227,6 @@
          IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
       ENDIF
+      !
+      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
 
       CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
@@ -415,4 +418,5 @@
          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
+         !
 	  dia_wri_alloc = MAXVAL(ierr)
       CALL mpp_sum( 'diawri', dia_wri_alloc )
@@ -460,5 +464,4 @@
          ninist = 0
       ENDIF
-	  
       !
       IF( nn_write == -1 )   RETURN   ! we will never do any output
@@ -488,5 +491,5 @@
       ijmi = 1      ;      ijma = jpj
       ipk = jpk
-	  IF(ln_abl) ipka = jpkam1
+      IF(ln_abl) ipka = jpkam1
 
       ! define time axis
@@ -724,6 +727,6 @@
        
          clmx ="l_max(only(x))"    ! max index on a period
-         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
+!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
+!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
 #if defined key_diahth
          CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
@@ -891,6 +894,6 @@
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
-      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
-      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
+!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
+!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
 
 #if defined key_diahth
@@ -907,5 +910,9 @@
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
-      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
+      IF( ln_zad_Aimp ) THEN
+         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current
+      ELSE
+         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current
+      ENDIF
       CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
       CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
@@ -969,5 +976,9 @@
       CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
       CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
-      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
+      IF( ln_zad_Aimp ) THEN
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity
+      ELSE
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity
+      ENDIF
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/dommsk.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/dommsk.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/dommsk.F90	(revision 11586)
@@ -100,5 +100,4 @@
          &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
          &             cn_ice, nn_ice_dta,                                     &
-         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
          &             ln_vol, nn_volctl, nn_rimwidth
       !!---------------------------------------------------------------------
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domvvl.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domvvl.F90	(revision 11586)
@@ -327,6 +327,6 @@
       END DO
       !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
-         !                                                            ! ------baroclinic part------ !
+      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
+         !                                                               ! ------baroclinic part------ !
          ! I - initialization
          ! ==================
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domwri.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domwri.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domwri.F90	(revision 11586)
@@ -162,8 +162,14 @@
       CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
       !															             ! vertical mesh
-      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
-      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
-      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
-      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3t_1d', e3t_1d, ktype = jp_r8  )    !    ! scale factors
+      CALL iom_rstput( 0, 0, inum, 'e3w_1d', e3w_1d, ktype = jp_r8  )
+      
+      CALL iom_rstput( 0, 0, inum, 'e3t_0' , e3t_0 , ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3uw_0', e3uw_0, ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3vw_0', e3vw_0, ktype = jp_r8  )
       !
       CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/dynhpg.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/dynhpg.F90	(revision 11586)
@@ -37,5 +37,5 @@
    USE trd_oce         ! trends: ocean variables
    USE trddyn          ! trend manager: dynamics
-!jc   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
+   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
    !
    USE in_out_manager  ! I/O manager
@@ -338,4 +338,5 @@
       REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
+      REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv
       !!----------------------------------------------------------------------
       !
@@ -346,6 +347,6 @@
       ENDIF
 
-      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level
-!jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv )
+      ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level
+      CALL zps_hde( kt, jpts, tsn, zgtsu, zgtsv, rhd, zgru , zgrv )
 
       ! Local constant initialization
@@ -385,5 +386,5 @@
       END DO
 
-      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
+      ! partial steps correction at the last level  (use zgru & zgrv computed in zpshde.F90)
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
@@ -395,5 +396,5 @@
                ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
                zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
-                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
+                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj)
                ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
             ENDIF
@@ -401,5 +402,5 @@
                va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
                zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
-                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
+                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj)
                va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
             ENDIF
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/sshwzv.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/sshwzv.F90	(revision 11586)
@@ -9,4 +9,5 @@
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
    !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!            4.0  !  2018-12  (A. Coward) add mixed implicit/explicit advection
    !!----------------------------------------------------------------------
 
@@ -279,5 +280,7 @@
       !!            :   wi      : now vertical velocity (for implicit treatment)
       !!
-      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
+      !!              implicit scheme for vertical advection in oceanic modeling. 
+      !!              Ocean Modelling, 91, 38-69.
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
@@ -286,5 +289,5 @@
       REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars
       REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
-      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
+      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
       REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
       REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
@@ -300,20 +303,40 @@
       ENDIF
       !
-      !
-      DO jk = 1, jpkm1            ! calculate Courant numbers
-         DO jj = 2, jpjm1
-            DO ji = 2, fs_jpim1   ! vector opt.
-               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
-               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   &  ! 2*rdt and not r2dt (for restartability)
-                  &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
-                  &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
-                  &                               * r1_e1e2t(ji,jj)                                                 &
-                  &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
-                  &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
-                  &                               * r1_e1e2t(ji,jj)                                                 &
-                  &                             ) * z1_e3t
+      ! Calculate Courant numbers
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, fs_jpim1   ! vector opt.
+                  z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+                  ! 2*rdt and not r2dt (for restartability)
+                  Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )                       &  
+                     &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk) + un_td(ji  ,jj,jk), 0._wp ) -   &
+                     &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                                     &
+                     &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk) + vn_td(ji,jj  ,jk), 0._wp ) -   &
+                     &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                                     &
+                     &                             ) * z1_e3t
+               END DO
             END DO
          END DO
-      END DO
+      ELSE
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, fs_jpim1   ! vector opt.
+                  z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+                  ! 2*rdt and not r2dt (for restartability)
+                  Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   & 
+                     &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
+                     &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                 &
+                     &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
+                     &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                     &                               * r1_e1e2t(ji,jj)                                                 &
+                     &                             ) * z1_e3t
+               END DO
+            END DO
+         END DO
+      ENDIF
       CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
       !
@@ -346,5 +369,5 @@
                   wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
                   !
-                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
+                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
                END DO
             END DO
@@ -353,5 +376,5 @@
       ELSE
          ! Fully explicit everywhere
-         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient
+         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
          wi    (:,:,:) = 0._wp
       ENDIF
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbrst.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbrst.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbrst.F90	(revision 11586)
@@ -131,6 +131,4 @@
       CALL iom_get( ncid, jpdom_unknown, 'kount' , zdata(:) )
       num_bergs(:) = INT(zdata(:))
-      ! Close file
-      CALL iom_close( ncid )
       !
 
@@ -146,4 +144,6 @@
       IF( lwp )   WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, icb_rst_read: there were',ibergs_in_file,   &
          &                                    ' bergs in the restart file and', jn,' bergs have been read'
+      ! Close file
+      CALL iom_close( ncid )
       !
       ! Confirm that all areas have a suitable base for assigning new iceberg
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbstp.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbstp.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbstp.F90	(revision 11586)
@@ -86,5 +86,5 @@
       !                                   !* write out time
       ll_verbose = .FALSE.
-      IF( nn_verbose_write > 0 .AND. MOD( kt-1 , nn_verbose_write ) == 0 )   ll_verbose = ( nn_verbose_level >= 0 )
+      IF( nn_verbose_write > 0 .AND. MOD( kt-1 , nn_verbose_write ) == 0 )   ll_verbose = ( nn_verbose_level > 0 )
       !
       IF( ll_verbose )   WRITE(numicb,9100) nktberg, ndastp, nsec_day
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11586)
@@ -30,6 +30,6 @@
 #if defined key_iomput
    USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1
-   USE trc_oce  , ONLY :   nn_dttrc                  ! frequency of step on passive tracers
-   USE icb_oce  , ONLY :   nclasses, class_num       ! iceberg classes
+   USE trc_oce  , ONLY :   nn_dttrc        !  !: frequency of step on passive tracers
+   USE icb_oce  , ONLY :   nclasses, class_num       !  !: iceberg classes
 #if defined key_si3
    USE ice      , ONLY :   jpl
@@ -230,4 +230,5 @@
           za_bnds(2,:) = ght_abl(2:jpka  ) + e3w_abl(2:jpka)
           CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds )
+
           CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) )
 # if defined key_si3
@@ -1960,8 +1961,10 @@
       INTEGER :: icnr, jcnr                             ! Offset such that the vertex coordinate (i+icnr,j+jcnr)
       !                                                 ! represents the bottom-left corner of cell (i,j)
-      REAL(wp), DIMENSION(4,jpi,jpj,2) ::   z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
-      REAL(wp), DIMENSION(  jpi,jpj  ) ::   z_fld       ! Working array to determine where to rotate cells
-      REAL(wp), DIMENSION(4        ,2) ::   z_rot       ! Lat/lon working array for rotation of cells
-      !!----------------------------------------------------------------------
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_fld       ! Working array to determine where to rotate cells
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_rot       ! Lat/lon working array for rotation of cells
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( z_bnds(4,jpi,jpj,2), z_fld(jpi,jpj), z_rot(4,2)  )
       !
       ! Offset of coordinate representing bottom-left corner
@@ -2043,4 +2046,6 @@
       CALL iom_set_domain_attr("grid_"//cdgrd, bounds_lat = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,1),(/ 4,ni*nj /)),           &
           &                                    bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 )
+      !
+      DEALLOCATE( z_bnds, z_fld, z_rot ) 
       !
    END SUBROUTINE set_grid_bounds
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom_nf90.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom_nf90.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom_nf90.F90	(revision 11586)
@@ -196,6 +196,6 @@
       CHARACTER(len=*)     , INTENT(in   )           ::   cdvar    ! name of the variable
       INTEGER              , INTENT(in   )           ::   kiv   ! 
-      INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of the dimensions
-      INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions
+      INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of each dimension
+      INTEGER              , INTENT(  out), OPTIONAL ::   kndims   ! number of dimensions
       LOGICAL              , INTENT(  out), OPTIONAL ::   lduld    ! true if the last dimension is unlimited (time)
       !
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/LBC/mppini.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/LBC/mppini.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/LBC/mppini.F90	(revision 11586)
@@ -167,5 +167,4 @@
            &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
            &             cn_ice, nn_ice_dta,                                     &
-           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
            &             ln_vol, nn_volctl, nn_rimwidth
       NAMELIST/nammpp/ jpni, jpnj, ln_nnogather, ln_listonly
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11586)
@@ -121,5 +121,6 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
 
-   INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read
+   !INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read
+   INTEGER , PUBLIC            ::   jpfld         !: maximum number of files to read
    INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   !: index of 10m wind velocity (i-component) (m/s)    at T-point
    INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   !: index of 10m wind velocity (j-component) (m/s)    at T-point
@@ -171,5 +172,6 @@
       !!
       CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
+      !TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i                 ! array of namelist informations on the fields to read
       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
@@ -216,4 +218,8 @@
       !                                   !* set the bulk structure
       !                                      !- store namelist information in an array
+      IF( ln_blk ) jpfld = 9
+      IF( ln_abl ) jpfld = 11
+      ALLOCATE( slf_i(jpfld) )
+      !
       slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
       slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
@@ -221,5 +227,7 @@
       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
       slf_i(jp_slp ) = sn_slp
-      slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj  
+      IF( ln_abl ) THEN
+        slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj  
+      END IF
       !
       !                                      !- allocate the bulk structure
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traadv_fct.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traadv_fct.F90	(revision 11586)
@@ -21,4 +21,5 @@
    USE diaar5         ! AR5 diagnostics
    USE phycst  , ONLY : rau0_rcp
+   USE zdf_oce , ONLY : ln_zad_Aimp
    !
    USE in_out_manager ! I/O manager
@@ -86,4 +87,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
+      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
       !!----------------------------------------------------------------------
       !
@@ -97,4 +100,5 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
+      ll_zAimp = .FALSE.
       IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
       IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
@@ -116,4 +120,22 @@
       !
       zwi(:,:,:) = 0._wp        
+      !
+      ! If adaptive vertical advection, check if it is needed on this PE at this time
+      IF( ln_zad_Aimp ) THEN
+         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
+      END IF
+      ! If active adaptive vertical advection, build tridiagonal matrix
+      IF( ll_zAimp ) THEN
+         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
+                  zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk)
+                  zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t_a(ji,jj,jk)
+                  zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      END IF
       !
       DO jn = 1, kjpt            !==  loop over the tracers  ==!
@@ -169,4 +191,29 @@
             END DO
          END DO
+         
+         IF ( ll_zAimp ) THEN
+            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
+            !
+            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
+                        &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         END IF
          !                
          IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
@@ -277,4 +324,30 @@
             zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
          ENDIF
+         !         
+         IF ( ll_zAimp ) THEN
+            DO jk = 1, jpkm1     !* trend and after field with monotonic scheme
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     !                             ! total intermediate advective trends
+                     ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                        &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                        &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+                     ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
+            !
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
          !
          CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. )
@@ -289,11 +362,35 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
-                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk)
+                  zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF ( ll_zAimp ) THEN
+            !
+            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
+                        &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF         
          !
          IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
@@ -318,4 +415,7 @@
       END DO                     ! end of tracer loop
       !
+      IF ( ll_zAimp ) THEN
+         DEALLOCATE( zwdia, zwinf, zwsup )
+      ENDIF
       IF( l_trd .OR. l_hst ) THEN 
          DEALLOCATE( ztrdx, ztrdy, ztrdz )
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traldf_iso.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traldf_iso.F90	(revision 11586)
@@ -289,5 +289,5 @@
          !!----------------------------------------------------------------------
          !
-         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp
+         ztfw(fs_2:1,:,:) = 0._wp     ;     ztfw(jpi:fs_jpim1,:,:) = 0._wp   ! avoid to potentially manipulate NaN values
          !
          ! Vertical fluxes
@@ -323,5 +323,5 @@
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
             DO jk = 2, jpkm1       
-               DO jj = 1, jpjm1
+               DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
@@ -336,5 +336,5 @@
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
                DO jk = 2, jpkm1 
-                  DO jj = 1, jpjm1
+                  DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
@@ -346,5 +346,5 @@
             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
                DO jk = 2, jpkm1 
-                  DO jj = 1, jpjm1
+                  DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      &
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/step.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/step.F90	(revision 11586)
@@ -165,15 +165,4 @@
                             CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
                             
-!!jc: fs simplification
-!!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636) 
-!!                                         but ensures reproductible results
-!!                                         with previous versions using split-explicit free surface          
-            IF( ln_zps .AND. .NOT. ln_isfcav )                               &
-               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
-               &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level
-            IF( ln_zps .AND.       ln_isfcav )                                          &
-               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
-               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
-!!jc: fs simplification
                             
                          ua(:,:,:) = 0._wp            ! set dynamics trends to zero
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/stpctl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/stpctl.F90	(revision 11413)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/stpctl.F90	(revision 11586)
@@ -96,5 +96,5 @@
             IF( ln_zad_Aimp ) THEN
                istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  )
-               istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  )
+               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  )
             ENDIF
             istatus = NF90_ENDDEF(idrun)
@@ -123,5 +123,5 @@
       IF( ln_zad_Aimp ) THEN
          zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max
-         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max
+         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max
       ENDIF
       !
