Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8098)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8099)
@@ -282,4 +282,5 @@
    !                                     !!** melt pond namelist (namicemp)
    LOGICAL , PUBLIC ::   ln_pnd           !: activate ponds or not
+   LOGICAL , PUBLIC ::   ln_pnd_rad       !: ponds radiatively active or not
    INTEGER , PUBLIC ::   nn_pnd_scheme    !: type of melt pond scheme:   =0 prescribed, =1 empirical, =2 topographic
    INTEGER , PUBLIC ::   nn_pnd_cpl       !: type of melt pond coupling: =0 passive, =1 full, =2 radiation only, =3 freshwater only
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 8098)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8099)
@@ -95,4 +95,7 @@
       IF(lwm) WRITE ( numoni, namicemp )
       
+      ! Shortcut for radiatively active melt ponds
+      ln_pnd_rad = ln_pnd .AND. ( ( nn_pnd_cpl == 1 ) .OR. ( nn_pnd_cpl == 2 ) )
+
       IF(lwp) THEN                        ! control print
          WRITE(numout,*)
@@ -104,7 +107,9 @@
          WRITE(numout,*)'    Prescribed pond fraction                                    rn_apnd       = ', rn_apnd
          WRITE(numout,*)'    Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd
+         WRITE(numout,*)'    Melt ponds radiatively active                               ln_pnd_rad    = ', ln_pnd_rad
       ENDIF
 
       IF ( .NOT. ln_pnd ) THEN
+         WRITE(numout,*)
          WRITE(numout,*) ' Melt ponds are not activated '
          WRITE(numout,*) ' nn_pnd_scheme, nn_pnd_cpl, rn_apnd and rn_hpnd are set to zero '
@@ -114,4 +119,5 @@
          rn_hpnd       = 0._wp
       ENDIF
+
 
       !
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 8098)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8099)
@@ -209,5 +209,7 @@
       !------------------------------------------------------------------------!
       CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os )    
-      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
+      ! 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
+
       alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
       CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, 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 8098)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8099)
@@ -39,5 +39,5 @@
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
-   REAL(wp) ::   rn_albice
+   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
 
    !!----------------------------------------------------------------------
@@ -48,5 +48,7 @@
 CONTAINS
 
-   SUBROUTINE albedo_ice( pt_ice, ph_ice, ph_snw, pa_ice_cs, pa_ice_os )
+   ! 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 )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_ice  ***
@@ -83,5 +85,11 @@
       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(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
@@ -89,17 +97,76 @@
       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
-      REAL(wp)            ::   zswitch, z1_c1, z1_c2
+      REAL(wp)                            ::   ralb_im, ralb_sf, ralb_sm, ralb_if, ralb_dp ! MV MP 2016
+      REAL(wp)                            ::   zswitch, z1_c1, z1_c2
+      REAL(wp)                            ::   zhref_pnd                                   ! MV MP 2016
       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
       !!---------------------------------------------------------------------
+      ! TO DO
+      ! Update namelist 
+      ! V Change the reference for rn_albice - should be bare ice with no puddles!!!
+      ! ... add 3rd computation
+
 
       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 )
 
       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
+
+      ENDIF
+
+      ! END MV MP 2016
+
       SELECT CASE ( nn_ice_alb )
 
@@ -109,14 +176,12 @@
       CASE( 0 )
        
-         ralb_sf = 0.80       ! dry snow
-         ralb_sm = 0.65       ! melting snow
-         ralb_if = 0.72       ! bare frozen ice
-         ralb_im = rn_albice  ! bare puddled ice 
+         ! 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)
          
-         !  Computation of ice albedo (free of snow)
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
-         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
-         END  WHERE
-      
          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 )
@@ -148,8 +213,8 @@
                   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
             END DO
@@ -163,9 +228,12 @@
       CASE( 1 ) 
 
-         ralb_im = rn_albice  ! bare puddled ice
+!         ! 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
-         ralb_if = 0.60      ! bare frozen ice
+!        ralb_sf = 0.85      ! dry snow
+!        ralb_sm = 0.75      ! melting snow
+!        ralb_if = 0.60      ! bare frozen ice
+!        ralb_dp = 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
 !         ralb_sf = 0.85       ! dry snow
@@ -177,11 +245,9 @@
 !         ralb_if = 0.54      ! bare frozen ice
 ! 
-         !  Computation of ice albedo (free of snow)
+         ! Computation of ice albedo (free of snow)
          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 = ralb_im
-         ELSE WHERE                                              ;   zalb = ralb_if
-         END  WHERE
-         
+
+         ! 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 *  &
@@ -192,8 +258,10 @@
          z1_c1 = 1. / 0.02
          z1_c2 = 1. / 0.03
-         !  Computation of the snow/ice albedo
+         
+         ! Surface albedo (accounting for snow, ice melt ponds), overcast skies
          DO jl = 1, ijpl
             DO jj = 1, jpj
                DO ji = 1, jpi
+
                   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 );
@@ -205,15 +273,67 @@
                   ! Ice/snow 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)
+                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
 
               END DO
             END DO
          END DO
-         ! Effect of the clouds (2d order polynomial)
+
+         ! Clear sky 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.
+!         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) 
 
       END SELECT
       
       CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
+      CALL wrk_dealloc( jpi,jpj,ijpl, zafrac_pnd, zh_pnd ) ! MV MP 2016
       !
    END SUBROUTINE albedo_ice
@@ -248,5 +368,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry , rn_alb_imlt, rn_alb_dpnd
       !!----------------------------------------------------------------------
       !
@@ -268,5 +388,9 @@
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
          WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
-         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
+         WRITE(numout,*) '      albedo of dry snow                                  rn_alb_sdry = ', rn_alb_sdry
+         WRITE(numout,*) '      albedo of melting snow                              rn_alb_smlt = ', rn_alb_smlt
+         WRITE(numout,*) '      albedo of dry ice                                   rn_alb_idry = ', rn_alb_idry
+         WRITE(numout,*) '      albedo of bare ice                                  rn_alb_imlt = ', rn_alb_imlt
+         WRITE(numout,*) '      albedo of ponded ice                                rn_alb_dpnd = ', rn_alb_dpnd
       ENDIF
       !
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 8098)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8099)
@@ -226,5 +226,7 @@
          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, 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
+
          SELECT CASE( kblk )
 
