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 12853)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icealb.F90	(revision 12854)
@@ -14,7 +14,8 @@
    !!   ice_alb_init   : initialisation of albedo computation
    !!----------------------------------------------------------------------
-   USE ice, ONLY: jpl ! sea-ice: number of categories
    USE phycst         ! physical constants
    USE dom_oce        ! domain: ocean
+   USE ice, ONLY: jpl ! sea-ice: number of categories
+   USE icevar         ! sea-ice: operations
    !
    USE in_out_manager ! I/O manager
@@ -100,4 +101,5 @@
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_ice     !  albedo of ice
       !
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: za_s_fra   ! ice fraction covered by snow
       INTEGER  ::   ji, jj, jl                ! dummy loop indices
       REAL(wp) ::   z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar
@@ -107,6 +109,4 @@
       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
       !!---------------------------------------------------------------------
       !
@@ -119,4 +119,6 @@
       z1_c4 = 1. / 0.03
       !
+      CALL ice_var_snwfra( ph_snw, za_s_fra )   ! calculate ice fraction covered by snow
+      !
       DO jl = 1, jpl
          DO jj = 1, jpj
@@ -126,14 +128,11 @@
                !--- Specific snow, ice and pond fractions ---!
                !---------------------------------------------!               
-               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
+               zafrac_snw = za_s_fra(ji,jj,jl)
+               IF( ld_pnd_alb ) THEN
+                  zafrac_pnd = MIN( pafrac_pnd(ji,jj,jl), 1._wp - zafrac_snw ) ! make sure (a_ip_eff + a_s_fra) <= 1
                ELSE
-                  zafrac_snw = 1._wp                   !--- thick snow fully shades melt ponds and ice
                   zafrac_pnd = 0._wp
-                  zafrac_ice = 0._wp
-               ENDIF
+               ENDIF
+               zafrac_ice = MAX( 0._wp, 1._wp - zafrac_pnd - zafrac_snw ) ! max for roundoff errors
                !
                !---------------!
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12853)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12854)
@@ -130,4 +130,5 @@
       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
+      REAL(wp), DIMENSION(jpij)            ::   za_s_fra    ! ice fraction covered by snow 
       !
       ! Mono-category
@@ -143,4 +144,6 @@
       END DO
 
+      CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) )   ! calculate ice fraction covered by snow
+      
       !------------------
       ! 1) Initialization
@@ -149,12 +152,15 @@
 
          ! If the snow thickness drops below zhs_min then reduce the snow fraction instead
-         IF( h_s_1d(ji) < zhs_min ) THEN
-           isnow(ji) = h_s_1d(ji) / zhs_min
-           zh_s(ji) = zhs_min * r1_nlay_s
-         ELSE
-           isnow(ji) = 1.0_wp
-           zh_s(ji) = h_s_1d(ji) * r1_nlay_s
-         END IF
-
+         !!IF( h_s_1d(ji) < zhs_min ) THEN
+         !!  isnow(ji) = h_s_1d(ji) / zhs_min
+         !!  zh_s(ji) = zhs_min * r1_nlay_s
+         !!ELSE
+         !!  isnow(ji) = 1.0_wp
+         !!  zh_s(ji) = h_s_1d(ji) * r1_nlay_s
+         !!END IF
+         isnow(ji) = za_s_fra(ji) !!clem: 2 variables for the same thing
+         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
+         
+         
          ! layer thickness
          zh_i(ji) = h_i_1d(ji) * r1_nlay_i
@@ -166,5 +172,5 @@
       END WHERE
       !
-      WHERE( zh_s(1:npti) > 0._wp   )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
+      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
       ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
       END WHERE
@@ -881,5 +887,5 @@
          !--- Snow-ice interfacial temperature (diagnostic SIMIP)
          zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-         IF( h_s_1d(ji) >= zhs_min ) THEN
+         IF( h_s_1d(ji) >= zhs_min ) THEN !!clem change that
             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
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 12853)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90	(revision 12854)
@@ -51,4 +51,5 @@
    !!   ice_var_sshdyn    : compute equivalent ssh in lead
    !!   ice_var_itd       : convert N-cat to M-cat
+   !!   ice_var_snwfra    : fraction of ice covered by snow
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -77,7 +78,12 @@
    PUBLIC   ice_var_sshdyn
    PUBLIC   ice_var_itd
+   PUBLIC   ice_var_snwfra
 
    INTERFACE ice_var_itd
       MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
+   END INTERFACE
+
+   INTERFACE ice_var_snwfra
+      MODULE PROCEDURE ice_var_snwfra_1d, ice_var_snwfra_2d, ice_var_snwfra_3d
    END INTERFACE
 
@@ -187,5 +193,5 @@
       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), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i, z1_a_ip
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i, z1_a_ip, za_s_fra
       !!-------------------------------------------------------------------
 
@@ -234,4 +240,8 @@
          &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min )
       END WHERE
+      !
+      CALL ice_var_snwfra( h_s, za_s_fra )           ! calculate ice fraction covered by snow
+      a_ip_eff = MIN( a_ip_eff, 1._wp - za_s_fra )   ! make sure (a_ip_eff + a_s_fra) <= 1
+      !
       !                                           !---  salinity (with a minimum value imposed everywhere)     
       IF( nn_icesal == 2 ) THEN
@@ -1268,4 +1278,43 @@
    END SUBROUTINE ice_var_itd_NcMc
 
+   !!-------------------------------------------------------------------
+   !! INTERFACE ice_var_snwfra
+   !!
+   !! ** Purpose :  fraction of ice covered by snow
+   !!
+   !! ** Method  :  In absence of proper snow model on top of sea ice,
+   !!               we argue that snow does not cover the whole ice because
+   !!               of wind blowing...
+   !!                
+   !! ** Arguments : ph_s: snow thickness
+   !!                
+   !! ** Output    : pa_s_fra: fraction of ice covered by snow
+   !!
+   !!-------------------------------------------------------------------
+   SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra )
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ph_s        ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
+      !!-------------------------------------------------------------------
+      pa_s_fra(:,:,:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:,:,:) )
+   END SUBROUTINE ice_var_snwfra_3d
+
+   SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra )
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
+      !!-------------------------------------------------------------------
+      pa_s_fra(:,:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:,:) )
+   END SUBROUTINE ice_var_snwfra_2d
+
+   SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra )
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness
+      REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
+      !!-------------------------------------------------------------------
+      pa_s_fra(:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:) )
+   END SUBROUTINE ice_var_snwfra_1d
+   
+
 #else
    !!----------------------------------------------------------------------
