Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 8105)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 8106)
@@ -383,12 +383,18 @@
 
             DO jl = 1, jpl
-               a_ip(:,:,jl) = 0.1 * zswitch(:,:)
-               v_ip(:,:,jl) = 0.1 * zswitch(:,:) * a_ip(:,:,jl)
+
+               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
+               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
+               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
+               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
+
             END DO
 
          ELSE
 
-            a_ip = 0._wp
-            v_ip = 0._wp
+            a_ip(:,:,:)      = 0._wp
+            v_ip(:,:,:)      = 0._wp
+            a_ip_frac(:,:,:) = 0._wp
+            h_ip     (:,:,:) = 0._wp
 
          ENDIF
@@ -418,8 +424,8 @@
          END DO
 
-         ! MV MP 2016
          a_ip(:,:,:)      = 0._wp
          v_ip(:,:,:)      = 0._wp
-         ! END MV MP 2016
+         a_ip_frac(:,:,:) = 0._wp
+         h_ip     (:,:,:) = 0._wp
 
       ENDIF ! ln_limini
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8105)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8106)
@@ -120,5 +120,4 @@
       ENDIF
 
-
       !
    END SUBROUTINE lim_mp_init
@@ -200,8 +199,8 @@
 
        WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 
-          a_ip      = rn_apnd * z1_jpl 
+          a_ip_frac = rn_apnd
           h_ip      = rn_hpnd    
-          v_ip      = a_ip * h_ip 
-          a_ip_frac = a_ip / a_i
+          v_ip      = a_ip_frac * a_i * h_ip 
+          a_ip      = a_ip_frac * a_i
        ELSE WHERE
           a_ip      = 0._wp
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8105)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8106)
@@ -209,6 +209,5 @@
       !------------------------------------------------------------------------!
       CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os )    
-      ! MV PND 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, a_i, ht_s, a_ip, v_ip, zalb_cs, zalb_os, ln_pnd_rad ) ! 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
 
       alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8105)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8106)
@@ -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
+   !!            x.x  ! 2017-05  (M. Vancoppenolle, O. Lecomte) Melt ponds
    !!----------------------------------------------------------------------
 
@@ -48,7 +49,5 @@
 CONTAINS
 
-   ! MV MP 2016
-   ! SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pa_ice_cs, pa_ice_os )
-   SUBROUTINE albedo_ice( pt_ice, ph_ice, pa_ice, ph_snw, pa_pnd, pv_pnd, pa_ice_cs, pa_ice_os, ld_pnd )
+   SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_ice  ***
@@ -60,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)
@@ -82,90 +82,64 @@
       !!                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 depth
-      ! MV MP 2016
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_ice      !  ice fraction
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_pnd      !  melt pond fraction
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pv_pnd      !  melt pond volume
-      LOGICAL , INTENT(in   )                   ::   ld_pnd      !  melt ponds radiatively active 
-      ! END MV MP 2016
-      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)                            ::   ralb_im, ralb_sf, ralb_sm, ralb_if, ralb_dp ! MV MP 2016
+      !!
+      !!----------------------------------------------------------------------
+      !!
+      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)                            ::   ralb_im, ralb_sf, ralb_sm, ralb_if, ralb_dp
       REAL(wp)                            ::   zswitch, z1_c1, z1_c2
-      REAL(wp)                            ::   zhref_pnd                                   ! MV MP 2016
+      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 2016
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_pnd, zh_pnd        ! melt pond fraction and thickness
       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
-      ! END MV MP 2016
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_snw                ! relative snow fraction
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zafrac_ice                ! relative ice fraction
+      !!
       !!---------------------------------------------------------------------
-      ! TO DO
-      ! Update namelist 
-      ! V Change the reference for rn_albice - should be bare ice with no puddles!!!
-      ! ... add 3rd computation
-
+
+      WRITE(numout,*) ' albedo_ice '
+      WRITE(numout,*) ' ld_pnd     : ', ld_pnd
+      WRITE(numout,*) ' nn_ice_alb : ', nn_ice_alb
 
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
       
       CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
-      CALL wrk_alloc( jpi,jpj,ijpl, zafrac_pnd, zh_pnd, zalb_pnd )
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb_pnd, zalb_ice, zalb_snw )
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb_pnd, zafrac_snw, zafrac_ice )
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
 
-      ralb_sf = rn_alb_sdry ! dry snow
-      ralb_sm = rn_alb_smlt ! melting snow
-      ralb_if = rn_alb_idry ! bare frozen ice
-      ralb_im = rn_alb_imlt ! bare puddled ice 
-      ralb_dp = rn_alb_dpnd ! deep pond albedo   ! MV MP 2016
-      
-      ! MV MP 2016
-      !------------------------------------------------
-      !  Pond fraction and thickness, snow-free albedo
-      !------------------------------------------------
-      ! If melt ponds are radiatively active ( ld_pnd = .TRUE. )
-      ! Then - melt pond fraction and thickness are assumed non-zero in the albedo computation 
-      !      - the snow-free albedo is the weighted mean of bare ice (ralb_if) and ponded-ice albedos (zalb_pnd)
-
-      ! Otherwise, we just specify the bare ice albedo
-      ! for cold ice and melting ice, the latter implicitly assuming a specified fraction of melt ponds
-
-      IF ( ld_pnd ) THEN
-
-         ! ponded ice fraction and thickness
-         WHERE ( pa_ice > 1.e-10_wp ) ;   zafrac_pnd = pa_pnd / pa_ice ;   zh_pnd = pv_pnd / pa_ice
-         ELSE WHERE                   ;   zafrac_pnd = 0._wp           ;   zh_pnd = 0._wp 
-         END WHERE
-          
-         ! ponded ice albedo, including thickness dependence:
-         zhref_pnd = 0.05         ! characteristic length scale for albedo dependence of ponded ice, Lecomte et al (2015)
-         zalb_pnd  = ralb_dp - ( ralb_dp - ralb_if ) * EXP( - zh_pnd / zhref_pnd )
-
-         ! weighted bare / ponded ice albedo (no ice thickness correction)
-         WHERE     ( ph_snw == 0._wp )  ;   zalb = ralb_if * ( 1. - zafrac_pnd ) + zalb_pnd * zafrac_pnd
-         END  WHERE
-
-      ELSE ! No melt ponds or radiatively inactive
-      
-         ! ponded ice fraction and thickness
-         zafrac_pnd = 0._wp
-         zh_pnd     = 0._wp
-
-         ! snow-free ice albedo
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im    ! prescribed bare ice albedo with assumed contribution of melt ponds
-         ELSE WHERE                                              ;   zalb = ralb_if
-         END  WHERE
+      ralb_sf = rn_alb_sdry    ! dry snow
+      ralb_sm = rn_alb_smlt    ! melting snow
+      ralb_if = rn_alb_idry    ! bare frozen ice
+      ralb_im = rn_alb_imlt    ! bare puddled ice (assumes implicit contribution of melt ponds)
+      ralb_dp = rn_alb_dpnd    ! deep pond albedo
+
+      !-----------------------------------------------------
+      !  Snow-free albedo (no ice thickness dependence yet)
+      !-----------------------------------------------------
+      !
+      ! Part common to nn_ice_alb = 0, 1, 2
+      !
+
+      zalb = ralb_if         
+      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 = ralb_im  ;  END WHERE  
+                                                       ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0
 
       ENDIF
-
-      ! END MV MP 2016
 
       SELECT CASE ( nn_ice_alb )
@@ -174,14 +148,9 @@
       !  Shine and Henderson-Sellers (1985)
       !------------------------------------------
+      ! NB: This parameterization is based on clear sky values
+
       CASE( 0 )
        
-         ! Reference (clear sky) values
-         !ralb_sf = 0.80       ! dry snow
-         !ralb_sm = 0.65       ! melting snow
-         !ralb_if = 0.72       ! bare frozen ice
-         !ralb_im = rn_albice  ! bare ice 
-         !ralb_dp = 0.30       ! deep ponded ice albedo, clear sky value, Lecomte et al (2015) MV MP 2016
-         !                     ! early melt ponds 0.24, late melt ponds 0.13 Grenfell & Perovich (2004)
-         
+         ! 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 )
@@ -191,9 +160,23 @@
          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  = ralb_dp - ( ralb_dp - 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 
+
+         ! Note that this way of coding things means that the thickness
+         ! dependence of bare ice albedo is also altering the ponded ice albedo
+         ! I'm not sure that this is valid, but at least, this maintains
+         ! consistency with previous formulation (see below)
      
          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
                   !                                        !  freezing snow        
@@ -203,15 +186,15 @@
                      &        +         zswitch   * ralb_sf  
 
-                  ! 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 ) )
                   zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
                       &     +         zswitch   *   ralb_sm 
-                  !
-                  ! 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
                
-                  ! surface albedo
+                  ! 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)
@@ -226,9 +209,9 @@
       !  New parameterization (2016)
       !------------------------------------------
+      ! NB: This parameterization is based on overcast skies values
+      
       CASE( 1 ) 
 
 !         ! Reference (overcast sky) values
-!       ralb_im = rn_albice  ! bare ice
-! compilation of values from literature
 !        ralb_sf = 0.85      ! dry snow
 !        ralb_sm = 0.75      ! melting snow
@@ -245,31 +228,41 @@
 !         ralb_if = 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
 
-         ! Accounting for the ice-thickness dependency (this has not been carefully checked with melt ponds)
-         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  = ralb_dp - ( ralb_dp - 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
          
-         ! Surface albedo (accounting for snow, ice melt ponds), overcast skies
+         ! 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 = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
                   zalb_sm = ralb_sm - ( ralb_sm - 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) 
@@ -279,61 +272,65 @@
          END DO
 
-         ! Clear sky albedo
+         ! 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 ); 
 
-!     !------------------------------------------
-!     !  Surface-based parameterization (2017)
-!     !------------------------------------------
-!      CASE( 2 ) 
-!
-!         !-----------------------------------------
-!         ! Snow, bare ice and ponded ice fractions 
-!         !-----------------------------------------
-!         ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
-!
-!         !--- Snow (for now, prescribed)
-!         WHERE     ( pa_ice > epsi10 .OR. ph_snw <= epsi06 )  ;  zafrac_snw = 0.
-!         ELSE WHERE                                           ;  zafrac_snw = 1.
-!         END WHERE
-!
-!         !--- Melt ponds 
-!         ! for now, done at the beginning
-!
-!         !--- Ice
-!         WHERE     ( pa_ice > epsi10 )   ;  zafrac_ice = 1. - zafrac_snw - zafrac_pnd
-!         ELSE WHERE                      ;  zafrac_ice = 0.
+      !---------------------------------------------------
+      !  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
+ 
+          !-----------------------------------------
+          ! Snow, bare ice and ponded ice fractions 
+          !-----------------------------------------
+          ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
+ 
+          !--- 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 - pafrac_pnd
+!          WHERE     ( ph_ice > 1.0e-6_wp ) ;  zafrac_ice = 1. - zafrac_snw - pafrac_pnd
+!         ELSE WHERE                       ;  zafrac_ice = 0.
 !         ENDWHERE
-!
-!         !--------------------------------------------------
-!         ! Snow-covered, pond-covered, and bare ice albedos
-!         !--------------------------------------------------
-!         ! zalb_snw, zalb_pnd, zalb_ice
-!         k
-!
-!
-!         
-!
-!          
-!
-!
-!
-!      ENDch grid cell has specific fractions of snow, ponds and ice
-!      ãs = 1 if hs > 0; 0 otherwise
-!      ãp = a_pnd / a_i
-!      ãi = ( 1 - ãs - ãp ) 
-!      
-!      2) Compute albedo for each of the three surface types
-!      αi = f(hi,Tsu)
-!      αs = f(hi, hs, Tsu)
-!      αp = f(hp, hi?, Tsu?)
-!      
-!      3) Surface albedo = weighted mean
-!      
-!      α = αs * ãs + αp * ãp + αi * ( 1 - ãs - ãp) 
+ 
+          !--------------------------------------------------
+          ! 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 = ralb_sf - ( ralb_sf - zalb_ice ) * EXP( - ph_snw * z1_c1 );
+          ELSE WHERE                 ; zalb_snw = ralb_sm - ( ralb_sm - zalb_ice ) * EXP( - ph_snw * z1_c2 );
+          END WHERE
+
+          ! 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  = ralb_dp - ( ralb_dp - zalb_ice ) * EXP( - ph_pnd / zhref_pnd ) 
+
+          ! Surface albedo is weighted mean of snow, ponds and bare ice contributions
+          pa_ice_os = zafrac_snw * zalb_snw  +  pafrac_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, zafrac_pnd, zh_pnd ) ! MV MP 2016
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb_pnd, zalb_ice, zalb_snw )
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb_pnd, zafrac_snw, zafrac_ice )
       !
    END SUBROUTINE albedo_ice
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8105)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8106)
@@ -226,6 +226,5 @@
          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, a_i, ht_s, a_ip, v_ip, zalb_cs, zalb_os, ln_pnd_rad ) ! 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( kblk )
