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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12832)
@@ -178,14 +178,13 @@
 &namthd_pnd     !   Melt ponds
 !------------------------------------------------------------------------------
-   ln_pnd           = .false.         !  activate melt ponds or not
-     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
-       rn_hpnd      =   0.05          !     prescribed pond depth,    at Tsu=0 degC
-     ln_pnd_alb     = .false.         !  melt ponds affect albedo or not
+   ln_pnd            = .false.        !  activate melt ponds or not
+      ln_pnd_LEV     = .false.        !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+         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.        !  constant  melt ponds
+         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC
+         rn_hpnd     =   0.05         !     prescribed pond depth,    at Tsu=0 degC
+      ln_pnd_lids    = .true.         !  frozen lids on top of the ponds (only for ln_pnd_LEV)
+      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo
 /
 !------------------------------------------------------------------------------
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namthd_pnd	(revision 12832)
@@ -2,13 +2,12 @@
 &namthd_pnd     !   Melt ponds
 !------------------------------------------------------------------------------
-   ln_pnd           = .false.         !  activate melt ponds or not
-     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
-       rn_hpnd      =   0.05          !     prescribed pond depth,    at Tsu=0 degC
-     ln_pnd_alb     = .false.         !  melt ponds affect albedo or not
+   ln_pnd            = .false.        !  activate melt ponds or not
+      ln_pnd_LEV     = .false.        !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+         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.        !  constant  melt ponds
+         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC
+         rn_hpnd     =   0.05         !     prescribed pond depth,    at Tsu=0 degC
+      ln_pnd_lids    = .true.         !  frozen lids on top of the ponds (only for ln_pnd_LEV)
+      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo
 /
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12832)
@@ -195,7 +195,5 @@
    !                                     !!** ice-ponds namelist (namthd_pnd)
    LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al (2012), 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
+   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010)
    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
@@ -203,4 +201,5 @@
    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
    REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1)
+   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Allow ponds to have frozen lids
    LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo
 
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icealb.F90	(revision 12832)
@@ -107,4 +107,6 @@
       REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction
       REAL(wp) ::   zalb_cs, zalb_os          ! albedo of ice under clear/overcast sky
+      REAL(wp) ::   zhs_min = 0.1_wp          ! minimum snow thickness to fully shade melt ponds and ice 
+      REAL(wp) ::   z1_hsmin
       !!---------------------------------------------------------------------
       !
@@ -120,17 +122,15 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
+               !
                !---------------------------------------------!
                !--- Specific snow, ice and pond fractions ---!
                !---------------------------------------------!               
-               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN   !--- no snow : we prevent melt ponds and snow at the same time (for now)
-                  zafrac_snw = 0._wp
-                  IF( ld_pnd_alb ) THEN
-                     zafrac_pnd = pafrac_pnd(ji,jj,jl)
-                  ELSE
-                     zafrac_pnd = 0._wp
-                  ENDIF
-                  zafrac_ice = 1._wp - zafrac_pnd
+               IF( ph_snw(ji,jj,jl) < zhs_min ) THEN   !--- thin snow partially shades melt ponds and ice
+                  IF( ld_pnd_alb ) THEN   ;   zafrac_pnd = pafrac_pnd(ji,jj,jl)
+                  ELSE                    ;   zafrac_pnd = 0._wp   ;   ENDIF
+                  zafrac_snw = MAX( 0._wp, ( ph_snw(ji,jj,jl) * z1_hsmin ) * ( 1._wp - zafrac_pnd ) ) ! max for roundoff errors
+                  zafrac_ice = MAX( 0._wp, 1._wp - zafrac_pnd - zafrac_snw )                          ! max for roundoff errors
                ELSE
-                  zafrac_snw = 1._wp                  !--- snow : fully "shades" melt ponds and ice
+                  zafrac_snw = 1._wp                   !--- thick snow fully shades melt ponds and ice
                   zafrac_pnd = 0._wp
                   zafrac_ice = 0._wp
@@ -144,6 +144,6 @@
                   zalb_ice = rn_alb_idry
                ELSE
-                  IF( ph_snw(ji,jj,jl) == 0._wp .AND. pt_su(ji,jj,jl) >= rt0 ) THEN  ;   zalb_ice = rn_alb_imlt
-                  ELSE                                                               ;   zalb_ice = rn_alb_idry   ;   ENDIF
+                  IF( ph_snw(ji,jj,jl) == 0._wp .AND. pt_su(ji,jj,jl) >= rt0 ) THEN   ;   zalb_ice = rn_alb_imlt
+                  ELSE                                                                ;   zalb_ice = rn_alb_idry   ;   ENDIF
                ENDIF
                !                       !--- Bare ice albedo (for hi < 150cm)
@@ -161,9 +161,6 @@
                ENDIF
                !                       !--- Ponded ice albedo
-               IF( ld_pnd_alb ) THEN
-                  zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd ) 
-               ELSE
-                  zalb_pnd = rn_alb_dpnd
-               ENDIF
+               zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd ) 
+               !
                !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
                zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_pra.F90	(revision 12832)
@@ -160,5 +160,5 @@
                z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
             END DO
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV ) THEN
                z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction
                z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume
@@ -196,5 +196,5 @@
             END DO
             !
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV ) THEN
                CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
                CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
@@ -231,5 +231,5 @@
                   &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
             END DO
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV ) THEN
                CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
                CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
@@ -257,5 +257,5 @@
                pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
             END DO
-            IF ( ln_pnd_H12 ) THEN
+            IF ( ln_pnd_LEV ) THEN
                pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
                pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
@@ -668,5 +668,5 @@
                   !                               ! -- check h_ip -- !
                   ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+                  IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
                      zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
                      IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
@@ -869,5 +869,5 @@
             END DO
             !
-            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction
+            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction
                IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
                   CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  )
@@ -911,5 +911,5 @@
             sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
             sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
-            IF( ln_pnd_H12 ) THEN
+            IF( ln_pnd_LEV ) THEN
                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
@@ -979,5 +979,5 @@
          END DO
          !
-         IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction
+         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction
             CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
             CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_adv_umx.F90	(revision 12832)
@@ -325,5 +325,5 @@
          !
          !== melt ponds ==!
-         IF ( ln_pnd_H12 ) THEN
+         IF ( ln_pnd_LEV ) THEN
             ! concentration
             zamsk = 1._wp
@@ -1558,5 +1558,5 @@
                   !                               ! -- check h_ip -- !
                   ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
-                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
+                  IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
                      zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
                      IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rdgrft.F90	(revision 12832)
@@ -565,5 +565,5 @@
                oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft 
 
-               IF ( ln_pnd_H12 ) THEN
+               IF ( ln_pnd_LEV ) THEN
                   aprdg1     = a_ip_2d(ji,jl1) * afrdg
                   aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
@@ -602,5 +602,5 @@
                sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
                oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
-               IF ( ln_pnd_H12 ) THEN
+               IF ( ln_pnd_LEV ) THEN
                   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)
@@ -699,5 +699,5 @@
                   v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
                      &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) )
-                  IF ( ln_pnd_H12 ) THEN
+                  IF ( ln_pnd_LEV ) THEN
                      v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
                         &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
@@ -705,6 +705,6 @@
                         &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
                      IF ( ln_pnd_lids ) THEN
-                        v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + (   vlrdg (ji) * rn_fpndrdg * fvol   (ji)   &
-                           &                                   + vlrft (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
                   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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceitd.F90	(revision 12832)
@@ -304,5 +304,5 @@
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
                a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
-               IF( ln_pnd_H12 )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
+               IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
                h_i_1d(ji) = rn_himin
             ENDIF
@@ -475,5 +475,5 @@
                zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
                !  
-               IF ( ln_pnd_H12 ) THEN
+               IF ( ln_pnd_LEV ) THEN
                   ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
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 12831)
+++ NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90	(revision 12832)
@@ -680,5 +680,5 @@
       WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
       WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
-      IF( ln_pnd_H12 ) THEN
+      IF( ln_pnd_LEV ) THEN
          WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
          WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
