Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 8231)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 8233)
@@ -54,4 +54,7 @@
 
    REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow          [kg/m3]
+   ! MV MP 2016
+   REAL(wp), PUBLIC ::   rhofw    = 1000._wp         !: volumic mass of freshwater in melt ponds [kg/m3]
+   ! END MV MP 2016
    REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice
    REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice               [psu]
@@ -176,4 +179,5 @@
          WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m^3'
          WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m^3'
+         WRITE(numout,*) '          density of freshwater (in melt ponds)     = ', rhofw   , ' kg/m^3'
          WRITE(numout,*) '          emissivity of snow or ice                 = ', emic  
          WRITE(numout,*) '          salinity of ice                           = ', sice    , ' psu'
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8231)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8233)
@@ -7,7 +7,8 @@
    !!   NEMO     1.0  ! 2003-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
    !!             -   ! 2004-11  (C. Talandier)  add albedo_init
-   !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0
+   !!             -   ! 2006-06  (M. Vancoppenolle) LIM 3.0
    !!             -   ! 2006-08  (G. Madec)  cleaning for surface module
    !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo
+   !!            4.0  ! 2017-05  (M. Vancoppenolle, O. Lecomte) Melt ponds
    !!----------------------------------------------------------------------
 
@@ -39,5 +40,5 @@
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
-   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt
+   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
 
    !!----------------------------------------------------------------------
@@ -48,5 +49,5 @@
 CONTAINS
 
-   SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pa_ice_cs, pa_ice_os )
+   SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_ice  ***
@@ -58,4 +59,5 @@
       !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
       !!                                                                           and Grenfell & Perovich (JGR 2004)
+      !!                  2: fractional surface-based formulation of scheme 1 (NEW)
       !!                Description of scheme 1:
       !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
@@ -80,16 +82,30 @@
       !!                Brandt et al. 2005, J. Climate, vol 18
       !!                Grenfell & Perovich 2004, JGR, vol 109 
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice      !  sea-ice thickness
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw      !  snow thickness
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs   !  albedo of ice under clear    sky
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
-      !!
-      INTEGER  ::   ji, jj, jl         ! dummy loop indices
-      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
-      REAL(wp)            ::   zswitch, z1_c1, z1_c2
+      !!
+      !!----------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice              !  ice surface temperature (Kelvin)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice              !  sea-ice thickness
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw              !  snow depth
+      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
+      LOGICAL , INTENT(in   )                   ::   ld_pnd              !  melt ponds radiatively active or not
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs           !  albedo of ice under clear    sky
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os           !  albedo of ice under overcast sky
+      !!
+      INTEGER  ::   ji, jj, jl                                           ! dummy loop indices
+      INTEGER  ::   ijpl                                                 ! number of ice categories (3rd dim of ice input arrays)
+      REAL(wp)                            ::   zswitch, z1_c1, z1_c2
+      REAL(wp)                            ::   zhref_pnd                                 
       REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
+!! MV MP
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_pnd                  ! ponded sea ice albedo
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_ice                  ! bare sea ice albedo
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_snw                  ! snow-covered sea ice albedo
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_snw                ! relative snow fraction
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_ice                ! relative ice fraction
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_pnd                ! relative ice fraction (effective)
+      !!
       !!---------------------------------------------------------------------
 
@@ -97,7 +113,23 @@
       
       CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb_pnd, zalb_ice, zalb_snw )
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb_pnd, zafrac_snw, zafrac_ice, zafrac_pnd )
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
-      
+
+      !-----------------------------------------------------
+      !  Snow-free albedo (no ice thickness dependence yet)
+      !-----------------------------------------------------
+      !
+      ! Part common to nn_ice_alb = 0, 1, 2
+      !
+      IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds
+         ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt
+                                                       ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0
+         ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry
+         END  WHERE
+      ENDIF
+
       SELECT CASE ( nn_ice_alb )
 
@@ -105,11 +137,9 @@
       !  Shine and Henderson-Sellers (1985)
       !------------------------------------------
+      ! NB: This parameterization is based on clear sky values
+
       CASE( 0 )
        
-         !  Computation of ice albedo (free of snow)
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt
-         ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry
-         END  WHERE
-      
+         ! Thickness-dependent bare ice albedo
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
@@ -119,9 +149,18 @@
          ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
          END WHERE
-     
+
+         IF ( ld_pnd ) THEN
+            ! Depth-dependent ponded ice albedo
+            zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd ) 
+
+            ! Snow-free ice albedo is a function of pond fraction
+            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ; zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;   END WHERE
+         ENDIF 
+
          DO jl = 1, ijpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! freezing snow
+                  ! Freezing snow
                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
@@ -130,5 +169,5 @@
                      &        +         zswitch   * rn_alb_sdry  
 
-                  ! melting snow
+                  ! Melting snow
                   ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
                   zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
@@ -136,11 +175,11 @@
                       &     +         zswitch   *   rn_alb_smlt 
                   !
-                  ! snow albedo
+                  ! Snow albedo
                   zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
                   zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
                
-                  ! Ice/snow albedo
-                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
-                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
+                  ! Surface albedo
+                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
                   !
                END DO
@@ -153,10 +192,14 @@
       !  New parameterization (2016)
       !------------------------------------------
+      ! NB: This parameterization is based on overcast skies values
+      
       CASE( 1 ) 
 
-! compilation of values from literature
+! compilation of values from literature (reference overcast sky values)
 !        rn_alb_sdry = 0.85      ! dry snow
 !        rn_alb_smlt = 0.75      ! melting snow
 !        rn_alb_idry = 0.60      ! bare frozen ice
+!        rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
+!                                ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
 ! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
 !        rn_alb_sdry = 0.85      ! dry snow
@@ -168,31 +211,40 @@
 !        rn_alb_idry = 0.54      ! bare frozen ice
 ! 
-         !  Computation of ice albedo (free of snow)
+         ! Computation of snow-free ice albedo 
          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
          z1_c2 = 1. / 0.05
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = rn_alb_imlt
-         ELSE WHERE                                              ;   zalb = rn_alb_idry
-         END  WHERE
-         
-         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
-         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+
+         ! Accounting for the ice-thickness dependency
+         WHERE     ( 1.5  < ph_ice                     )        ;  zalb_it = zalb
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )        ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
             &                                                                     ( LOG(1.5) - LOG(ph_ice) )
-         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
+         ELSE WHERE                                             ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
          END WHERE
+
+         IF ( ld_pnd ) THEN
+            ! Depth-dependent ponded ice albedo
+            zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd ) 
+
+            ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
+            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;  zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;  END WHERE
+         ENDIF 
 
          z1_c1 = 1. / 0.02
          z1_c2 = 1. / 0.03
-         !  Computation of the snow/ice albedo
+         
+         ! Overcast sky surface albedo (accounting for snow, ice melt ponds)
          DO jl = 1, ijpl
             DO jj = 1, jpj
                DO ji = 1, jpi
+                  ! Snow depth dependence of snow albedo
                   zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
                   zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
 
-                  ! snow albedo
+                  ! Snow albedo (MV I guess we could use rt0 instead of rt0_snow)
                   zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
                   zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
 
-                  ! Ice/snow albedo   
+                  ! Surface albedo   
                   zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
                   pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
@@ -201,10 +253,75 @@
             END DO
          END DO
-         ! Effect of the clouds (2d order polynomial)
+
+         ! Clear sky surface albedo
          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+
+      !---------------------------------------------------
+      !  Fractional surface-based parameterization (2017)
+      !---------------------------------------------------
+      CASE( 2 ) 
+ 
+      ! MV: I propose this formulation that is more elegant, and more easy to expand towards
+      !     varying pond and snow fraction.
+      !     Formulation 1 has issues to handle ponds and snow together that
+      !     can't easily be fixed. This one handles it better, I believe.
+
+          !-----------------------------------------
+          ! Snow, bare ice and ponded ice fractions 
+          !-----------------------------------------
+          ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
+
+          !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time)
+          zafrac_pnd = 0._wp
+          IF ( ld_pnd ) THEN  
+             WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds
+          ENDIF         
+
+          !--- Specific snow fraction (for now, prescribed)
+          WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1.
+          ELSE WHERE                       ;  zafrac_snw = 0.
+          END WHERE
+ 
+          !--- Specific ice fraction 
+          zafrac_ice = 1. - zafrac_snw - zafrac_pnd
+ 
+          !--------------------------------------------------
+          ! Snow-covered, pond-covered, and bare ice albedos
+          !--------------------------------------------------
+          ! Bare ice albedo
+          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+          z1_c2 = 1. / 0.05
+          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb
+          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+            &                                                                       ( LOG(1.5) - LOG(ph_ice) )
+          ELSE WHERE                                       ;  zalb_ice = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
+          END WHERE
+
+          ! Snow-covered ice albedo (freezing, melting cases)
+          z1_c1 = 1. / 0.02
+          z1_c2 = 1. / 0.03
+         
+          WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 );
+          ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 );
+          END WHERE
+
+          ! Depth-dependent ponded ice albedo
+          IF ( ld_pnd ) THEN
+             zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd / zhref_pnd ) 
+          ELSE
+             zalb_pnd  = rn_alb_dpnd
+          ENDIF
+
+          ! Surface albedo is weighted mean of snow, ponds and bare ice contributions
+          pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice
+          
+          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 )
 
       END SELECT
       
       CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb_pnd, zalb_ice, zalb_snw )
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb_pnd, zafrac_snw, zafrac_ice, zafrac_pnd )
       !
    END SUBROUTINE albedo_ice
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8231)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8233)
@@ -48,4 +48,8 @@
    USE limvar          ! Ice variables switch
    USE limctl          !
+   ! MV MP 2016
+   USE limmp
+   ! END MV MP 2016
+
    USE limistate       ! LIM initial state
    USE limthd_sal      ! LIM ice thermodynamics: salinity
@@ -208,5 +212,6 @@
          CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
          
-                                      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+                                      CALL albedo_ice( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos MV MP 2016
+
          SELECT CASE( ksbc )
             CASE( jp_usr )   ;        CALL usrdef_sbc_ice_flx( kt ) ! user defined formulation
@@ -231,4 +236,9 @@
          ! --- zap this if no ice thermo --- !
          IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics      
+
+         ! MV MP 2016
+         IF ( ln_pnd )                CALL lim_mp( kt )         ! -- Melt ponds
+         ! END MV MP 2016
+
          IF( ln_limthd )              CALL lim_update2( kt )    ! -- Corrections
          ! ---
@@ -311,5 +321,9 @@
       !
       CALL lim_thd_sal_init            ! set ice salinity parameters
-      !
+       
+      ! MV MP 2016
+      CALL lim_mp_init                 ! set melt ponds parameters
+      ! END MV MP 2016
+
       IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation
       !                                ! Initial sea-ice state
@@ -620,4 +634,8 @@
       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
       wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
+
+      ! MV MP 2016
+      wfx_pnd(:,:) = 0._wp   ;   wfx_snw_sum(:,:) = 0._wp
+      ! END MV MP 2016
       
       hfx_thd(:,:) = 0._wp   ;
