Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12720)
@@ -49,4 +49,6 @@
           <field id="icehpnd"      long_name="melt pond depth"                                         standard_name="sea_ice_meltpond_depth"                    unit="m" /> 
           <field id="icevpnd"      long_name="melt pond volume"                                        standard_name="sea_ice_meltpond_volume"                   unit="m" /> 
+          <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" /> 
+          <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" /> 
 	  
 	  <!-- heat -->
@@ -287,5 +289,7 @@
           <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        /> 
           <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" /> 
+          <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" /> 
           <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        /> 
+          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        /> 
           <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        />
           <field id="iceage_cat"   long_name="Ice age per category"                              unit="days"    detect_missing_value="true" />
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12720)
@@ -176,5 +176,9 @@
 !------------------------------------------------------------------------------
    ln_pnd           = .false.         !  activate melt ponds or not
-     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+       ln_pnd_lids  = .true.          !  ponds with frozen lids
+       ln_pnd_flush = .true.          !  ponds flushing trhu the ice  
+       rn_apnd_min  =   0.15          !  minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
+       rn_apnd_max  =   0.85          !  maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
      ln_pnd_CST     = .false.         !  activate constant  melt ponds
        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
@@ -206,4 +210,6 @@
    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North
    rn_hpd_ini_s     =   0.05          !        "            "             South
+   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North
+   rn_hld_ini_s     =   0.0           !        "            "             South
    ! -- for ln_iceini_file = T
    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
@@ -217,4 +223,5 @@
    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hld = 'NOT USED'              , -12 ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    cn_dir='./'
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ref	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ref	(revision 12720)
@@ -643,4 +643,5 @@
    bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
+   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds
    rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice
@@ -649,4 +650,5 @@
    rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            --
    rn_ice_hpnd = 0.05         !       --   pond depth                          --
+   rn_ice_hlid = 0.0          !       --   pond lid depth                      --
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/nambdy_dta
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/nambdy_dta	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/nambdy_dta	(revision 12720)
@@ -29,4 +29,5 @@
    bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
+   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds
    rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice
@@ -35,3 +36,4 @@
    rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            --
    rn_ice_hpnd = 0.05         !       --   pond depth                          --
+   rn_ice_hlid = 0.0          !       --   pond lid depth                      --
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namini
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namini	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namini	(revision 12720)
@@ -23,4 +23,6 @@
    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North
    rn_hpd_ini_s     =   0.05          !        "            "             South
+   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North
+   rn_hld_ini_s     =   0.0           !        "            "             South
    ! -- for ln_iceini_file = T
    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
@@ -34,4 +36,5 @@
    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hld = 'NOT USED'              , -12 ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    cn_dir='./'
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_pnd
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_pnd	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_pnd	(revision 12720)
@@ -3,5 +3,9 @@
 !------------------------------------------------------------------------------
    ln_pnd           = .false.         !  activate melt ponds or not
-     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+       ln_pnd_lids  = .true.          !  ponds with frozen lids
+       ln_pnd_flush = .true.          !  ponds flushing trhu the ice  
+       rn_apnd_min  =   0.15          !  minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
+       rn_apnd_max  =   0.85          !  maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
      ln_pnd_CST     = .false.         !  activate constant  melt ponds
        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12720)
@@ -70,4 +70,5 @@
    !! a_ip        |      -      |    Ice pond concentration       |       |
    !! v_ip        |      -      |    Ice pond volume per unit area| m     |
+   !! v_il        |    v_il_1d  |    Ice pond lid volume per area | m     |
    !!                                                                     |
    !!-------------|-------------|---------------------------------|-------|
@@ -85,4 +86,5 @@
    !! t_su        ! t_su_1d     |    Sea ice surface temperature  ! K     |
    !! h_ip        | h_ip_1d     |    Ice pond thickness           | m     |
+   !! h_il        | h_il_1d     |    Ice pond lid thickness       | m     |
    !!                                                                     |
    !! notes: the ice model only sees a bulk (i.e., vertically averaged)   |
@@ -112,4 +114,6 @@
    !! hm_ip       |      -      |    Mean ice pond depth          | m     |
    !! vt_ip       |      -      |    Total ice pond vol. per unit area| m |
+   !! hm_il       |      -      |    Mean ice pond lid depth      | m     |
+   !! vt_il       |      -      |    Total ice pond lid vol. per area | m |
    !!=====================================================================
 
@@ -190,5 +194,9 @@
    !                                     !!** 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_H12       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010)
+   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Allow ponds to have frozen lids
+   LOGICAL,  PUBLIC ::   ln_pnd_flush     !: Allow ponds to flush thru the ice
+   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds
+   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds
    LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth
    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
@@ -331,9 +339,14 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_eff   !: melt pond effective fraction (not covered up by lid) (a_ip/a_i)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_il       !: melt pond lid volume                     [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_il       !: melt pond lid thickness                  [m]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration
    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 gridcell area [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_il      !: mean melt pond lid depth                     [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_il      !: total melt pond lid volume per gridcell area [m]
 
    !!----------------------------------------------------------------------
@@ -448,8 +461,9 @@
 
       ii = ii + 1
-      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
-
-      ii = ii + 1
-      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
+      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  &
+         &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) )
+
+      ii = ii + 1
+      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , hm_il(jpi,jpj) , vt_il(jpi,jpj) , STAT = ierr(ii) )
 
       ! * Old values of global variables
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice1d.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice1d.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice1d.F90	(revision 12720)
@@ -124,8 +124,11 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   o_i_1d        !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !: ice ponds
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_eff_1d   !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_il_1d       !: Ice pond lid
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_il_1d       !:
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s
@@ -157,4 +160,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_il_2d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d
@@ -208,6 +212,6 @@
          &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &    
          &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  &
-         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) ,                   &
-         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   &
+         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , v_il_1d (jpij) , h_il_1d(jpij) , &
+         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , a_ip_eff_1d(jpij) ,                               &
          &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
       !
@@ -226,5 +230,5 @@
       ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  &
          &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  &
-         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      &
+         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , v_il_2d(jpij,jpl) ,  &
          &      STAT=ierr(ii) )
 
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icealb.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icealb.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icealb.F90	(revision 12720)
@@ -96,4 +96,5 @@
       LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
+                                                                  !  This is the effective fraction not covered up by a pond lid
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90	(revision 12720)
@@ -99,6 +99,8 @@
       WHERE( a_ip(:,:,:) >= epsi20 )
          h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
+         h_il(:,:,:) = v_il(:,:,:) / a_ip(:,:,:)
       ELSEWHERE
          h_ip(:,:,:) = 0._wp
+         h_il(:,:,:) = 0._wp
       END WHERE
       !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv.F90	(revision 12720)
@@ -84,10 +84,10 @@
          !                             !-----------------------!
          CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, &
-            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i )
          !                             !-----------------------!
       CASE( np_advPRA )                ! PRATHER scheme        !
          !                             !-----------------------!
          CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, h_i, h_s, h_ip, &
-            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i )
       END SELECT
 
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90	(revision 12720)
@@ -44,4 +44,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume
 
    !! * Substitutions
@@ -55,5 +56,5 @@
 
    SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
-      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                **  routine ice_dyn_adv_pra  **
@@ -81,8 +82,9 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il     ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
       !
-      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
       INTEGER  ::   icycle                  ! number of sub-timestep for the advection
       REAL(wp) ::   zdt                     !   -      -
@@ -93,5 +95,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
-      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl
       REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
       REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
@@ -159,6 +161,7 @@
             END DO
             IF ( ln_pnd_H12 ) THEN
-               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
-               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction
+               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume
+               z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)      ! Melt pond lid volume
             ENDIF
          END DO
@@ -196,4 +199,6 @@
                CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )    !--- melt pond lid volume
+               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 
             ENDIF
             !                                                               !--------------------------------------------!
@@ -227,5 +232,7 @@
                CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
-            ENDIF
+               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
+               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )    !--- melt pond lid volume 
+           ENDIF
             !
          ENDIF
@@ -247,4 +254,5 @@
                pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
                pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+               pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
             ENDIF
          END DO
@@ -263,5 +271,5 @@
          !     Remove negative values (conservation is ensured)
          !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
-         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
          !
          ! --- Make sure ice thickness is not too big --- !
@@ -756,6 +764,7 @@
          &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
          &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
-         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
-         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
+         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
+         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
+         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   &
          !
          &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
@@ -864,4 +873,10 @@
                CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
                CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
+               !                                                     ! melt pond lid volume
+               CALL iom_get( numrir, jpdom_autoglo, 'sxvl' , sxvl  )
+               CALL iom_get( numrir, jpdom_autoglo, 'syvl' , syvl  )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxxvl', sxxvl )
+               CALL iom_get( numrir, jpdom_autoglo, 'syyvl', syyvl )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxyvl', sxyvl )
             ENDIF
             !
@@ -880,4 +895,5 @@
                sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
                sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
+               sxvl  = 0._wp   ;   syvl  = 0._wp   ;   sxxvl  = 0._wp   ;   syyvl  = 0._wp   ;   sxyvl  = 0._wp   ! melt pond lid volume
             ENDIF
          ENDIF
@@ -954,4 +970,10 @@
             CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
             CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
+            !                                                        ! melt pond lid volume
+            CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
+            CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
+            CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
          ENDIF
          !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_umx.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_umx.F90	(revision 12720)
@@ -60,5 +60,5 @@
 
    SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
-      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ice_dyn_adv_umx  ***
@@ -85,4 +85,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -334,4 +335,9 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, pv_ip, zua_ups, zva_ups )
+            ! lid
+            zamsk = 0._wp
+            zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, pv_il, zua_ups, zva_ups )            
          ENDIF
          !
@@ -350,5 +356,5 @@
          ! Remove negative values (conservation is ensured)
          !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
-         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
          !
          ! --- Make sure ice thickness is not too big --- !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rdgrft.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rdgrft.F90	(revision 12720)
@@ -494,6 +494,6 @@
       REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
       REAL(wp)                  ::   airft1, oirft1, aprft1
-      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
-      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
+      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice
       !
       REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
@@ -569,7 +569,9 @@
                   aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
                   vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
+                  vlrdg (ji) = v_il_2d(ji,jl1) * afrdg
                   aprft1     = a_ip_2d(ji,jl1) * afrft
                   aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
                   vprft (ji) = v_ip_2d(ji,jl1) * afrft
+                  vlrft (ji) = v_il_2d(ji,jl1) * afrft
                ENDIF
 
@@ -601,4 +603,5 @@
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
+                  v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji)
                ENDIF
             ENDIF
@@ -697,4 +700,6 @@
                      a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
                         &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
+                     v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + (   vlrdg (ji) * rn_fpndrdg * fvol   (ji)   &
+                        &                                   + vlrft (ji) * rn_fpndrft * zswitch(ji)   )
                   ENDIF
                   
@@ -727,5 +732,5 @@
       !----------------
       ! In case ridging/rafting lead to very small negative values (sometimes it happens)
-      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
       !
    END SUBROUTINE rdgrft_shift
@@ -839,4 +844,5 @@
          CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -865,4 +871,5 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90	(revision 12720)
@@ -45,9 +45,9 @@
    REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
    REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
-   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n
-   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s
+   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n
+   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s
    !
    !                              ! if ln_iceini_file = T
-   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read
+   INTEGER , PARAMETER ::   jpfldi = 10          ! maximum number of files to read
    INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
    INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
@@ -59,4 +59,5 @@
    INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
    INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
+   INTEGER , PARAMETER ::   jp_hld = 10          ! index of pnd lid depth    (m)
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
    !   
@@ -98,8 +99,8 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
       REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini, zhlid_ini            !data from namelist or nc file
       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, zaip_2d, zhip_2d
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d
       !--------------------------------------------------------------------
 
@@ -155,6 +156,9 @@
       a_ip     (:,:,:) = 0._wp
       v_ip     (:,:,:) = 0._wp
+      v_il     (:,:,:) = 0._wp
       a_ip_frac(:,:,:) = 0._wp
+      a_ip_eff (:,:,:) = 0._wp
       h_ip     (:,:,:) = 0._wp
+      h_il     (:,:,:) = 0._wp
       !
       ! ice velocities
@@ -216,4 +220,8 @@
                &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
             !
+            ! pond lid depth
+            IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
+               &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
+            !
             zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
             ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
@@ -222,4 +230,5 @@
             zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
             zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
+            zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1)
             !
             ! change the switch for the following
@@ -246,4 +255,5 @@
                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(:,:)
+               zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
             ELSEWHERE
                zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
@@ -256,4 +266,5 @@
                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(:,:)
+               zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
             END WHERE
             !
@@ -264,4 +275,5 @@
             zapnd_ini(:,:) = 0._wp
             zhpnd_ini(:,:) = 0._wp
+            zhlid_ini(:,:) = 0._wp
          ENDIF
          
@@ -290,14 +302,18 @@
          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 )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_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), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
+         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), &
+            &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_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), a_ip_1d(1:npti), h_ip_1d(1:npti), &
-            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_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), h_il_1d(1:npti), &
+            &              zti_2d          , zts_2d          , ztsu_2d        ,                  &
+            &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d )
 
          ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
@@ -315,8 +331,9 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   )
          CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   )
+         CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   )
 
          ! deallocate temporary arrays
          DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
-            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
+            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )
 
          ! calculate extensive and intensive variables
@@ -365,5 +382,7 @@
             a_ip_frac(:,:,:) = 0._wp
          END WHERE
+         a_ip_eff(:,:,:) = a_ip_frac(:,:,:)
          v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+         v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
           
          ! specific temperatures for coupled runs
@@ -466,5 +485,5 @@
       !
       CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
-      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
+      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld
       TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
       !
@@ -473,6 +492,6 @@
          &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
          &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
-         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
-         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
+         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, &
+         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir
       !!-----------------------------------------------------------------------------
       !
@@ -488,5 +507,5 @@
       slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
       slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
-      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd
+      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd   ;   slf_i(jp_hld) = sn_hld
       !
       IF(lwp) THEN                          ! control print
@@ -508,4 +527,5 @@
             WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
             WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
+            WRITE(numout,*) '      initial pnd lid depth in the north-south         rn_hld_ini     = ', rn_hld_ini_n,rn_hld_ini_s
          ENDIF
       ENDIF
@@ -532,5 +552,6 @@
          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' )
+         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
+         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' )
       ENDIF
       !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceitd.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceitd.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceitd.F90	(revision 12720)
@@ -410,4 +410,5 @@
       CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
       CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
+      CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
       CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       DO jl = 1, jpl
@@ -482,4 +483,8 @@
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
                   v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
+                  !                                              
+                  ztrans          = v_il_2d(ji,jl1) * zworka(ji)     ! Pond lid volume
+                  v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans
+                  v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans
                ENDIF
                !
@@ -526,5 +531,5 @@
       ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
       !       because of truncation error ( i.e. 1. - 1. /= 0 )
-      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
 
       ! at_i must be <= rn_amax
@@ -554,4 +559,5 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
       CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
+      CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       DO jl = 1, jpl
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icerst.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icerst.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icerst.F90	(revision 12720)
@@ -132,4 +132,5 @@
       CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip  )
       CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip  )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_il' , v_il  )
       ! Snow enthalpy
       DO jk = 1, nlay_s 
@@ -171,5 +172,5 @@
       INTEGER           ::   jk
       LOGICAL           ::   llok
-      INTEGER           ::   id0, id1, id2, id3, id4   ! local integer
+      INTEGER           ::   id0, id1, id2, id3, id4, id5   ! local integer
       CHARACTER(len=25) ::   znam
       CHARACTER(len=2)  ::   zchar, zchar1
@@ -250,4 +251,13 @@
             v_ip(:,:,:) = 0._wp
          ENDIF
+         a_ip_eff(:,:,:) = a_ip(:,:,:)
+         ! melt pond lids
+         id5 = iom_varid( numrir, 'v_il' , ldstop = .FALSE. )
+         IF( id5 > 0 ) THEN
+            CALL iom_get( numrir, jpdom_autoglo, 'v_il', v_il)
+         ELSE
+            IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds lids output then set it to zero'
+            v_il(:,:,:) = 0._wp
+         ENDIF
          ! fields needed for Met Office (Jules) coupling
          IF( ln_cpl ) THEN
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90	(revision 12720)
@@ -132,5 +132,5 @@
 
       ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd.F90	(revision 12720)
@@ -355,4 +355,6 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            )
@@ -441,4 +443,5 @@
          sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
          v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
+         v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
          oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
          
@@ -461,4 +464,6 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
@@ -515,4 +520,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
          !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90	(revision 12720)
@@ -89,10 +89,14 @@
          IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
             a_ip_frac_1d(ji) = rn_apnd
+            a_ip_eff_1d(ji)  = rn_apnd
             h_ip_1d(ji)      = rn_hpnd    
             a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
+            h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever
          ELSE
             a_ip_frac_1d(ji) = 0._wp
+            a_ip_eff_1d(ji)  = 0._wp
             h_ip_1d(ji)      = 0._wp    
             a_ip_1d(ji)      = 0._wp
+            h_il_1d(ji)      = 0._wp
          ENDIF
          !
@@ -106,85 +110,220 @@
       !!                ***  ROUTINE pnd_H12  ***
       !!
-      !! ** Purpose    : Compute melt pond evolution
-      !!
-      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds 
-      !!                 and sent to ocean when surface is freezing
-      !!
-      !!                 pond growth:      Vp = Vp + dVmelt
-      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
-      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
-      !!                    with Tp = -2degC
-      !!  
-      !! ** Tunable parameters : (no real expertise yet, ideas?)
+      !! ** Purpose : Compute melt pond evolution
+      !!
+      !! ** Method  : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
+      !!              We  work with volumes and then redistribute changes into thickness and concentration
+      !!              assuming linear relationship between the two. 
+      !!
+      !! ** Action  : - pond growth:      Vp = Vp + dVmelt                                          --- from Holland et al 2012 ---
+      !!                                     dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
+      !!                                        dh_i  = meltwater from ice surface melt
+      !!                                        dh_s  = meltwater from snow melt
+      !!                                        (1-r) = fraction of melt water that is not flushed
+      !!
+      !!              - limtations:       a_ip must not exceed (1-r)*a_i
+      !!                                  h_ip must not exceed 0.5*h_i
+      !!
+      !!              - pond shrinking:
+      !!                       if lids:   Vp = Vp -dH * a_ip
+      !!                                     dH = lid thickness change. Retrieved from this eq.:    --- from Flocco et al 2010 ---
+      !!
+      !!                                                                   rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 
+      !!                                                                      H = lid thickness
+      !!                                                                      Lf = latent heat of fusion
+      !!                                                                      Tp = -2C
+      !!
+      !!                                                                And solved implicitely as:
+      !!                                                                   H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
+      !!
+      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 ---
+      !!
+      !!              - Overflow:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 ---
+      !!                                     perm = permability of sea-ice
+      !!                                     visc = water viscosity
+      !!                                     Hp   = height of top of the pond above sea-level
+      !!                                     Hi   = ice thickness thru which there is flushing
+      !!
+      !!
+      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness
+      !!
+      !!              - effective pond area: to be used for albedo 
+      !!
+      !!              - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
+      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect
+      !!
+      !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
       !! 
-      !! ** Note       : Stolen from CICE for quick test of the melt pond
-      !!                 radiation and freshwater interfaces
-      !!                 Coupling can be radiative AND freshwater
-      !!                 Advection, ridging, rafting are called
-      !!
-      !! ** References : Holland, M. M. et al (J Clim 2012)
-      !!-------------------------------------------------------------------
-      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
-      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            -
-      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
-      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
-      !
-      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
-      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
-      REAL(wp) ::   z1_Tp            ! inverse reference temperature
-      REAL(wp) ::   z1_rhow          ! inverse freshwater density
-      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
-      REAL(wp) ::   zfac, zdum
-      !
-      INTEGER  ::   ji   ! loop indices
-      !!-------------------------------------------------------------------
-      z1_rhow        = 1._wp / rhow 
-      z1_zpnd_aspect = 1._wp / zpnd_aspect
-      z1_Tp          = 1._wp / zTp 
+      !! ** Note       :   mostly stolen from CICE
+      !!
+      !! ** References :   Flocco and Feltham (JGR, 2007)
+      !!                   Flocco et al       (JGR, 2010)
+      !!                   Holland et al      (J. Clim, 2012)
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(nlay_i) ::   zperm          ! Permeability of sea ice
+      !!
+      REAL(wp), PARAMETER ::   zaspect =  0.8_wp      ! pond aspect ratio
+      REAL(wp), PARAMETER ::   zTp     = -2._wp       ! reference temperature
+      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity
+      REAL(wp), PARAMETER ::   zhl_max =  0.015_wp    ! pond lid thickness above which the ponds disappear from the albedo calculation
+      REAL(wp), PARAMETER ::   zhl_min =  0.005_wp    ! pond lid thickness below which the full pond area is used in the albedo calculation
+      !!
+      REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding
+      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes
+      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh
+      REAL(wp) ::   v_ip_max                          ! max pond volume allowed
+      REAL(wp) ::   zdT                               ! zTp-t_su
+      REAL(wp) ::   zsbr                              ! Brine salinity
+      REAL(wp) ::   zfac, zdum                        ! temporary arrays
+      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
+      !!
+      INTEGER  ::   ji, jk                            ! loop indices
+      !!-------------------------------------------------------------------
+      z1_rhow   = 1._wp / rhow 
+      z1_aspect = 1._wp / zaspect
+      z1_Tp     = 1._wp / zTp 
 
       DO ji = 1, npti
-         !                                                        !--------------------------------!
-         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
-            !                                                     !--------------------------------!
-            !--- Remove ponds on thin ice
+         !                                                            !----------------------------------------------------!
+         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
+            !                                                         !----------------------------------------------------!
+            !--- Remove ponds on thin ice or tiny ice fractions
             a_ip_1d(ji)      = 0._wp
             a_ip_frac_1d(ji) = 0._wp
             h_ip_1d(ji)      = 0._wp
-            !                                                     !--------------------------------!
-         ELSE                                                     ! Case ice thickness >= rn_himin !
-            !                                                     !--------------------------------!
-            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
-            !
-            ! available meltwater for melt ponding [m, >0] and fraction
-            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
-            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
-            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
-            !
-            !--- Pond gowth ---!
-            ! v_ip should never be negative, otherwise code crashes
-            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
-            !
-            ! melt pond mass flux (<0)
+            h_il_1d(ji)      = 0._wp
+            !
+            ! clem: problem with conservation or not ?
+            !                                                         !--------------------------------!
+         ELSE                                                         ! Case ice thickness >= rn_himin !
+            !                                                         !--------------------------------!
+            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
+            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
+            !
+            !------------------!
+            ! case ice melting !
+            !------------------!
+            !
+            !--- available meltwater for melt ponding ---!
+            zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
+            zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) in H12 = fraction of melt water that is not flushed
+            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 
+            !
+            !--- overflow ---!
+            ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
+            !    a_ip_max = zfr_mlt * a_i
+            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
+            v_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
+            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) )
+
+            ! If pond depth exceeds half the ice thickness then reduce the pond volume
+            !    h_ip_max = 0.5 * h_i
+            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
+            v_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
+            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) )
+            
+            !--- Pond growing ---!
+            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
+            !
+            !--- Lid melting ---!
+            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
+            !
+            !--- mass flux ---!
             IF( zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
+               zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1]
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
                !
-               ! adjust ice/snow melting flux to balance melt pond flux (>0)
-               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
+               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux
                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
             ENDIF
+
+            !-------------------!
+            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
+            !-------------------!
+            !
+            zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
             !
             !--- Pond contraction (due to refreezing) ---!
-            v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
-            !
-            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
-            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
-            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
+            IF( ln_pnd_lids ) THEN
+               !
+               !--- Lid growing and subsequent pond shrinking ---! 
+               zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
+                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
+               
+               ! Lid growing
+               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
+               
+               ! Pond shrinking
+               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
+
+            ELSE
+               ! Pond shrinking
+               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
+            ENDIF
+            !
+            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
+            ! v_ip     = h_ip * a_ip
+            ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip)
+            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) )
             a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
-            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
+            h_ip_1d(ji)      = zaspect * a_ip_frac_1d(ji)
+
+            !---------------!            
+            ! Pond flushing !
+            !---------------!
+            IF( ln_pnd_flush ) THEN
+               ! height of top of the pond above sea-level
+               zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_frac_1d(ji) ) ) * r1_rau0
+
+               ! Calculate the permeability of the ice (Assur 1958)
+               DO jk = 1, nlay_i
+                  zsbr = - 1.2_wp                                  &
+                     &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
+                     &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
+                     &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 ! clem: error here the factor was 0.01878 instead of 0.0178 (cf Flocco 2010)
+                  zperm(jk) = MAX( 0._wp, 3.e-08_wp * (sz_i_1d(ji,jk) / zsbr)**3 )
+               END DO
+
+               ! Do the drainage using Darcy's law
+               zdv_flush   = -MINVAL(zperm(:)) * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
+               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
+               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
+               
+               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
+               a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) )
+               a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
+               h_ip_1d(ji)      = zaspect * a_ip_frac_1d(ji)
+
+            ENDIF
+
+            !--- Corrections and lid thickness ---!
+            IF( ln_pnd_lids ) THEN
+               !--- remove ponds if lids are much larger than ponds ---!
+               IF ( v_il_1d(ji) > v_ip_1d(ji) * 10._wp ) THEN
+                  a_ip_1d(ji)      = 0._wp
+                  a_ip_frac_1d(ji) = 0._wp
+                  h_ip_1d(ji)      = 0._wp
+                  v_il_1d(ji)      = 0._wp
+               ENDIF
+               !--- retrieve lid thickness from volume ---!
+               IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
+               ELSE                              ;   h_il_1d(ji) = 0._wp ;
+               ENDIF
+            ENDIF
             !
          ENDIF
+         
       END DO
+
+      !-------------------------------------------------!            
+      ! How much melt pond is exposed to the atmosphere !
+      !-------------------------------------------------!            
+      ! Calculate the melt pond effective area (used for albedo)
+      WHERE    ( h_il_1d(1:npti) <= zhl_min )   ;   a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti)       ! lid is very thin.  Expose all the pond
+      ELSEWHERE( h_il_1d(1:npti) >= zhl_max )   ;   a_ip_eff_1d(1:npti) = 0._wp                      ! lid is very thick. Cover all the pond up with ice and snow
+      ELSEWHERE                                 ;   a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti) * &   ! lid is in between. Expose part of the pond
+         &                                                                ( h_il_1d(1:npti) - zhl_min ) / ( zhl_max - zhl_min )
+      END WHERE
       !
    END SUBROUTINE pnd_H12
@@ -205,5 +344,7 @@
       INTEGER  ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_lids, ln_pnd_flush, rn_apnd_min, rn_apnd_max, &
+         &                          ln_pnd_CST, rn_apnd, rn_hpnd, &
+         &                          ln_pnd_alb
       !!-------------------------------------------------------------------
       !
@@ -221,10 +362,14 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namicethd_pnd:'
-         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
+         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
+         WRITE(numout,*) '         Evolutive  melt pond fraction and depth                  ln_pnd_H12   = ', ln_pnd_H12
+         WRITE(numout,*) '            Melt ponds can have frozen lids                       ln_pnd_lids  = ', ln_pnd_lids
+         WRITE(numout,*) '            Allow ponds to flush thru the ice                     ln_pnd_flush = ', ln_pnd_flush
+         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
+         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
+         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
       !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceupdate.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceupdate.F90	(revision 12720)
@@ -185,5 +185,5 @@
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
       !
       alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90	(revision 12720)
@@ -113,4 +113,5 @@
       at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
       vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
+      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 )
       !
       ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
@@ -161,6 +162,6 @@
          !
          !                           ! mean melt pond depth
-         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)
-         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp
+         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)   ;   hm_il(:,:) = vt_il(:,:) / at_ip(:,:)
+         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp                     ;   hm_il(:,:) = 0._wp
          END WHERE         
          !
@@ -221,4 +222,8 @@
       WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
       ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
+      END WHERE
+      !                                           !--- pond lid thickness      
+      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_il(:,:,:) = v_il(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
+      ELSEWHERE                            ;   h_il(:,:,:) = 0._wp
       END WHERE
       !
@@ -289,4 +294,5 @@
       sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
       v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+      v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
       !
    END SUBROUTINE ice_var_eqv2glo
@@ -533,4 +539,5 @@
                a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
                v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+               v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj)
                !
             END DO
@@ -555,5 +562,5 @@
 
 
-   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_zapneg ***
@@ -570,4 +577,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -636,10 +644,10 @@
       WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
-      !                                                        but it does not change conservation, so keep it this way is ok
+      WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok
       !
    END SUBROUTINE ice_var_zapneg
 
 
-   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_roundoff ***
@@ -654,4 +662,5 @@
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -668,4 +677,5 @@
          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
+         WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 )    pv_il(1:npti,:)   = 0._wp   ! v_il must be >= 0
       ENDIF
       !
@@ -786,6 +796,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, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting 1-cat ice to 1 ice category
@@ -793,6 +803,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)    ::   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), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !!-------------------------------------------------------------------
       ! == thickness and concentration == !
@@ -808,9 +818,10 @@
       pa_ip(:) = patip(:)
       ph_ip(:) = phtip(:)
+      ph_il(:) = phtil(:)
       
    END SUBROUTINE ice_var_itd_1c1c
 
-   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 )
+   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting N-cat ice to 1 ice category
@@ -818,6 +829,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)    ::   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), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
@@ -854,6 +865,10 @@
       ! == 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
+      WHERE( pa_ip(:) /= 0._wp )
+         ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+         ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+      ELSEWHERE
+         ph_ip(:) = 0._wp
+         ph_il(:) = 0._wp
       END WHERE
       !
@@ -862,6 +877,6 @@
    END SUBROUTINE ice_var_itd_Nc1c
    
-   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 )
+   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !!
@@ -885,6 +900,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)    ::   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), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti
@@ -997,10 +1012,19 @@
          END WHERE
       END DO
+      ! keep the same v_il/v_i ratio for each category
+      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
+      ELSEWHERE                                 ;   zfra(:) = 0._wp
+      END WHERE
+      DO jl = 1, jpl
+         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+         ELSEWHERE                       ;   ph_il(:,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, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !!
@@ -1033,6 +1057,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)    ::   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), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
@@ -1063,4 +1087,5 @@
          pa_ip(:,:) = patip(:,:)
          ph_ip(:,:) = phtip(:,:)
+         ph_il(:,:) = phtil(:,:)
          !                              ! ---------------------- !
       ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
@@ -1068,6 +1093,6 @@
          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(:,:)  )
+            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
+            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  )
          !                              ! ---------------------- !
       ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
@@ -1075,6 +1100,6 @@
          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)  )
+            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
+            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  )
          !                              ! ----------------------- !
       ELSE                              ! input cat /= output cat !
@@ -1218,4 +1243,15 @@
             END WHERE
          END DO
+         ! keep the same v_il/v_i ratio for each category
+         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
+            zfra(:) = SUM( phtil(:,:) * 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_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
+            END WHERE
+         END DO
          DEALLOCATE( zfra )
          !
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icewri.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icewri.F90	(revision 12720)
@@ -116,4 +116,6 @@
       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
+      IF( iom_use('icehlid' ) )   CALL iom_put( 'icehlid', hm_il  * zmsk00      )                                           ! melt pond lid depth
+      IF( iom_use('icevlid' ) )   CALL iom_put( 'icevlid', vt_il  * zmsk00      )                                           ! melt pond lid total volume per unit area
       ! salt
       IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity
@@ -162,6 +164,8 @@
       IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume
       IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories
-      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond frac for categories
+      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories
+      IF( iom_use('icehlid_cat' ) )   CALL iom_put( 'icehlid_cat' ,   h_il         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories
       IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories
+      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories
       IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories
 
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdy_oce.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdy_oce.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdy_oce.F90	(revision 12720)
@@ -63,4 +63,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::  aip    !: now ice  pond concentration
       REAL(wp), POINTER, DIMENSION(:,:) ::  hip    !: now ice  pond depth
+      REAL(wp), POINTER, DIMENSION(:,:) ::  hil    !: now ice  pond lid depth
 #if defined key_top
       CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
@@ -115,4 +116,5 @@
    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
+   REAL(wp), DIMENSION(jp_bdy) ::   rice_hlid               !: pond lid thick. of incoming sea ice
    !
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdydta.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdydta.F90	(revision 12720)
@@ -43,5 +43,5 @@
    PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
 
-   INTEGER , PARAMETER ::   jpbdyfld  = 16    ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jpbdyfld  = 17    ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_bdyssh = 1     ! 
    INTEGER , PARAMETER ::   jp_bdyu2d = 2     ! 
@@ -60,4 +60,5 @@
    INTEGER , PARAMETER ::   jp_bdyaip = 15    ! 
    INTEGER , PARAMETER ::   jp_bdyhip = 16    ! 
+   INTEGER , PARAMETER ::   jp_bdyhil = 17    ! 
 #if ! defined key_si3
    INTEGER , PARAMETER ::   jpl = 1
@@ -189,4 +190,5 @@
                         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) 
+                        dta_bdy(jbdy)%hil(ib,jl) =  h_il(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -294,4 +296,5 @@
                &                                                                         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( TRIM(bf_alias(jp_bdyhil)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhil)%fnow(:,1,:) = rice_hlid(jbdy)
 
             ! if T_i is read and not T_su, set T_su = T_i
@@ -318,4 +321,5 @@
                bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp
                bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp
+               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp
             ENDIF
             
@@ -323,12 +327,12 @@
             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                  , &
-                  &              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 )
+               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,:), & ! in
+                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & ! out
+                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &                                  ! in (optional)
+                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &                                  ! in     -
+                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), bf_alias(jp_bdyhil)%fnow(:,1,:), & ! in     -
+                  &              dta_alias%t_i                  , dta_alias%t_s                  , &                                  ! out    -
+                  &              dta_alias%tsu                  , dta_alias%s_i                  , &                                  ! out    -
+                  &              dta_alias%aip                  , dta_alias%hip                  , dta_alias%hil )                    ! out    -
             ENDIF
          ENDIF
@@ -377,5 +381,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 
+      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid
       INTEGER                                ::   ipk,ipl       !
       INTEGER                                ::   idvar         ! variable ID
@@ -389,11 +393,11 @@
       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(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, bn_hil       
       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, 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/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil
+      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid
       NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp
       !!---------------------------------------------------------------------------
@@ -452,6 +456,6 @@
 #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' )
+            rn_ice_apnd = 0. ; rn_ice_hpnd = 0. ; rn_ice_hlid = 0.
+            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 & rn_ice_hlid = 0 when no ponds' )
          ENDIF
 #endif
@@ -463,5 +467,6 @@
          rice_apnd(jbdy) = rn_ice_apnd
          rice_hpnd(jbdy) = rn_ice_hpnd
-         
+         rice_hlid(jbdy) = rn_ice_hlid
+
          
          DO jfld = 1, jpbdyfld
@@ -562,5 +567,5 @@
             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
+               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip .OR. jfld == jp_bdyhil ) THEN
                igrd = 1                                                    ! T point
                ipk = ipl                                                   ! jpl-cat data
@@ -613,4 +618,9 @@
                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
+            IF( jfld == jp_bdyhil ) THEN
+               cl3 = 'hil'
+               bf_alias => bf(jp_bdyhil,jbdy:jbdy)                         ! alias for hil structure of bdy number jbdy
+               bn_alias => bn_hil                                          ! alias for hil structure of nambdy_dta 
             ENDIF
 
@@ -682,4 +692,9 @@
                   ENDIF
                ENDIF
+               IF( jfld == jp_bdyhil ) THEN
+                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hil => bf_alias(1)%fnow(:,1,:)
+                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hil(iszdim,jpl) )
+                  ENDIF
+               ENDIF
             ENDIF
 
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdyice.F90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/BDY/bdyice.F90	(revision 12720)
@@ -94,8 +94,8 @@
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
             ! exchange 3d arrays
-            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. &
-                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
-                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                &
-                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      )
+            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1.                 &
+               &                        , s_i , 'T', 1., t_su, 'T', 1., v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.  &
+               &                        , a_ip, 'T', 1., v_ip, 'T', 1., v_il, 'T', 1.                                &
+               &                        , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
             ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk
             CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
@@ -163,4 +163,5 @@
             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
+            h_il(ji,jj,  jl) = ( h_il(ji,jj,  jl) * zwgt1 + dta%hil(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid depth
             !
             sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
@@ -170,4 +171,5 @@
                a_ip(ji,jj,jl) = 0._wp
                h_ip(ji,jj,jl) = 0._wp
+               h_il(ji,jj,jl) = 0._wp
             ENDIF
             !
@@ -231,4 +233,5 @@
                a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl)
                h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl)
+               h_il(ji,jj,  jl) = h_il(ib,jb,  jl)
                !
                sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
@@ -274,4 +277,5 @@
                ENDIF
                v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl)
+               v_il(ji,jj,jl) = h_il(ji,jj,jl) * a_ip(ji,jj,jl)
                !
             ELSE   ! no ice at the boundary
@@ -281,6 +285,4 @@
                h_s (ji,jj,  jl) = 0._wp
                oa_i(ji,jj,  jl) = 0._wp
-               a_ip(ji,jj,  jl) = 0._wp
-               v_ip(ji,jj,  jl) = 0._wp
                t_su(ji,jj,  jl) = rt0
                t_s (ji,jj,:,jl) = rt0
@@ -288,7 +290,7 @@
 
                a_ip_frac(ji,jj,jl) = 0._wp
+               a_ip     (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
+               h_il     (ji,jj,jl) = 0._wp
                
                IF( nn_icesal == 1 ) THEN     ! if constant salinity
@@ -306,4 +308,6 @@
                e_s (ji,jj,:,jl) = 0._wp
                e_i (ji,jj,:,jl) = 0._wp
+               v_ip(ji,jj,  jl) = 0._wp
+               v_il(ji,jj,  jl) = 0._wp
 
             ENDIF
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/LBC/lbc_lnk_multi_generic.h90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/LBC/lbc_lnk_multi_generic.h90	(revision 12719)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/LBC/lbc_lnk_multi_generic.h90	(revision 12720)
@@ -15,26 +15,30 @@
 #endif
 
-   SUBROUTINE ROUTINE_MULTI( cdname                                                                             &
-      &                    , pt1, cdna1, psgn1, pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4, cdna4, psgn4   &
-      &                    , pt5, cdna5, psgn5, pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8, cdna8, psgn8   &
-      &                    , pt9, cdna9, psgn9, pt10, cdna10, psgn10, pt11, cdna11, psgn11                      &
+   SUBROUTINE ROUTINE_MULTI( cdname                                                                               &
+      &                    , pt1 , cdna1 , psgn1 , pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4 , cdna4 , psgn4   &
+      &                    , pt5 , cdna5 , psgn5 , pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8 , cdna8 , psgn8   &
+      &                    , pt9 , cdna9 , psgn9 , pt10, cdna10, psgn10, pt11, cdna11, psgn11, pt12, cdna12, psgn12  &
+      &                    , pt13, cdna13, psgn13, pt14, cdna14, psgn14, pt15, cdna15, psgn15, pt16, cdna16, psgn16  &
       &                    , kfillmode, pfillval, lsend, lrecv, ihlcom )
       !!---------------------------------------------------------------------
-      CHARACTER(len=*)   ,                   INTENT(in   ) :: cdname  ! name of the calling subroutine
-      ARRAY_TYPE(:,:,:,:)          , TARGET, INTENT(inout) :: pt1     ! arrays on which the lbc is applied
-      ARRAY_TYPE(:,:,:,:), OPTIONAL, TARGET, INTENT(inout) :: pt2  , pt3  , pt4  , pt5  , pt6  , pt7  , pt8  , pt9  , pt10  , pt11
-      CHARACTER(len=1)                     , INTENT(in   ) :: cdna1   ! nature of pt2D. array grid-points
-      CHARACTER(len=1)   , OPTIONAL        , INTENT(in   ) :: cdna2, cdna3, cdna4, cdna5, cdna6, cdna7, cdna8, cdna9, cdna10, cdna11
-      REAL(wp)                             , INTENT(in   ) :: psgn1   ! sign used across the north fold
-      REAL(wp)           , OPTIONAL        , INTENT(in   ) :: psgn2, psgn3, psgn4, psgn5, psgn6, psgn7, psgn8, psgn9, psgn10, psgn11
-      INTEGER            , OPTIONAL        , INTENT(in   ) :: kfillmode   ! filling method for halo over land (default = constant)
-      REAL(wp)           , OPTIONAL        , INTENT(in   ) :: pfillval    ! background value (used at closed boundaries)
-      LOGICAL, DIMENSION(4), OPTIONAL      , INTENT(in   ) :: lsend, lrecv   ! indicate how communications are to be carried out
-      INTEGER            , OPTIONAL        , INTENT(in   ) :: ihlcom         ! number of ranks and rows to be communicated
+      CHARACTER(len=*)     ,                   INTENT(in   ) ::   cdname  ! name of the calling subroutine
+      ARRAY_TYPE(:,:,:,:)            , TARGET, INTENT(inout) ::   pt1     ! arrays on which the lbc is applied
+      ARRAY_TYPE(:,:,:,:)  , OPTIONAL, TARGET, INTENT(inout) ::   pt2   , pt3   , pt4   , pt5   , pt6   , pt7   , pt8   , pt9  , &
+         &                                                        pt10  , pt11  , pt12  , pt13  , pt14  , pt15  , pt16
+      CHARACTER(len=1)                       , INTENT(in   ) ::   cdna1   ! nature of pt2D. array grid-points
+      CHARACTER(len=1)     , OPTIONAL        , INTENT(in   ) ::   cdna2 , cdna3 , cdna4 , cdna5 , cdna6 , cdna7 , cdna8 , cdna9, &
+         &                                                        cdna10, cdna11, cdna12, cdna13, cdna14, cdna15, cdna16
+      REAL(wp)                               , INTENT(in   ) ::   psgn1   ! sign used across the north fold
+      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   psgn2 , psgn3 , psgn4 , psgn5 , psgn6 , psgn7 , psgn8 , psgn9, &
+         &                                                        psgn10, psgn11, psgn12, psgn13, psgn14, psgn15, psgn16
+      INTEGER              , OPTIONAL        , INTENT(in   ) ::   kfillmode   ! filling method for halo over land (default = constant)
+      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries)
+      LOGICAL, DIMENSION(4), OPTIONAL        , INTENT(in   ) ::   lsend, lrecv   ! indicate how communications are to be carried out
+      INTEGER              , OPTIONAL        , INTENT(in   ) ::   ihlcom         ! number of ranks and rows to be communicated
       !!
       INTEGER                          ::   kfld        ! number of elements that will be attributed
-      PTR_TYPE         , DIMENSION(11) ::   ptab_ptr    ! pointer array
-      CHARACTER(len=1) , DIMENSION(11) ::   cdna_ptr    ! nature of ptab_ptr grid-points
-      REAL(wp)         , DIMENSION(11) ::   psgn_ptr    ! sign used across the north fold boundary
+      PTR_TYPE         , DIMENSION(16) ::   ptab_ptr    ! pointer array
+      CHARACTER(len=1) , DIMENSION(16) ::   cdna_ptr    ! nature of ptab_ptr grid-points
+      REAL(wp)         , DIMENSION(16) ::   psgn_ptr    ! sign used across the north fold boundary
       !!---------------------------------------------------------------------
       !
@@ -55,6 +59,11 @@
       IF( PRESENT(psgn10) )   CALL ROUTINE_LOAD( pt10, cdna10, psgn10, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
       IF( PRESENT(psgn11) )   CALL ROUTINE_LOAD( pt11, cdna11, psgn11, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn12) )   CALL ROUTINE_LOAD( pt12, cdna12, psgn12, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn13) )   CALL ROUTINE_LOAD( pt13, cdna13, psgn13, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn14) )   CALL ROUTINE_LOAD( pt14, cdna14, psgn14, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn15) )   CALL ROUTINE_LOAD( pt15, cdna15, psgn15, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
+      IF( PRESENT(psgn16) )   CALL ROUTINE_LOAD( pt16, cdna16, psgn16, ptab_ptr, cdna_ptr, psgn_ptr, kfld )
       !
-      CALL lbc_lnk_ptr    ( cdname,               ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ihlcom )
+      CALL lbc_lnk_ptr( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ihlcom )
       !
    END SUBROUTINE ROUTINE_MULTI
