Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icealb.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icealb.F90	(revision 11843)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icealb.F90	(revision 13662)
@@ -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
@@ -45,5 +46,5 @@
 CONTAINS
 
-   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os )
+   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE ice_alb  ***
@@ -97,7 +98,8 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky
-      !
+      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   pcloud_fra   !  cloud fraction
+      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
@@ -106,4 +108,5 @@
       REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction
       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
       !!---------------------------------------------------------------------
       !
@@ -116,28 +119,30 @@
       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
             DO ji = 1, jpi
-               !                       !--- Specific snow, ice and pond fractions (for now, we prevent melt ponds and snow at the same time)
-               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN
-                  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
+               !
+               !---------------------------------------------!
+               !--- Specific snow, ice and pond fractions ---!
+               !---------------------------------------------!               
+               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      ! 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
+               !
+               !---------------!
+               !--- Albedos ---!
+               !---------------!               
                !                       !--- Bare ice albedo (for hi > 150cm)
                IF( ld_pnd_alb ) THEN
                   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)
@@ -155,16 +160,15 @@
                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
-               palb_os(ji,jj,jl) = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
-               !
-               palb_cs(ji,jj,jl) = palb_os(ji,jj,jl)  &
-                  &                - ( - 0.1010 * palb_os(ji,jj,jl) * palb_os(ji,jj,jl)  &
-                  &                    + 0.1933 * palb_os(ji,jj,jl) - 0.0148 ) * tmask(ji,jj,1)
-               !
+               zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
+               !
+               zalb_cs = zalb_os - ( - 0.1010 * zalb_os * zalb_os  &
+                  &                  + 0.1933 * zalb_os - 0.0148 ) * tmask(ji,jj,1)
+               !
+               ! albedo depends on cloud fraction because of non-linear spectral effects
+               palb_ice(ji,jj,jl) = ( 1._wp - pcloud_fra(ji,jj) ) * zalb_cs + pcloud_fra(ji,jj) * zalb_os
+
             END DO
          END DO
