Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8592)
@@ -211,14 +211,10 @@
 &namalb     !   albedo parameters
 !------------------------------------------------------------------------------
-   nn_ice_alb       =    1            !  parameterization of ice/snow albedo
-                                      !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
-                                      !     1: "home made" based on Brandt et al. (JClim 2005)
-                                      !                         and Grenfell & Perovich (JGR 2004), giving cloud-sky albedo
-                                      !     2: as 1 with melt ponds
-   rn_alb_sdry      =   0.85          !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
-   rn_alb_smlt      =   0.75          !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
-   rn_alb_idry      =   0.60          !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
-   rn_alb_imlt      =   0.50          !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
-   rn_alb_dpnd      =   0.27          !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' )
+   !                                  !                          !  obs range (cloud-sky)
+   rn_alb_sdry      =   0.85          !  dry snow albedo         :  0.85 -- 0.87
+   rn_alb_smlt      =   0.75          !  melting snow albedo     :  0.72 -- 0.82
+   rn_alb_idry      =   0.60          !  dry ice albedo          :  0.54 -- 0.65
+   rn_alb_imlt      =   0.50          !  bare puddled ice albedo :  0.49 -- 0.58
+   rn_alb_dpnd      =   0.27          !  ponded ice albedo       :  0.10 -- 0.30 
 /
 !------------------------------------------------------------------------------
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8592)
@@ -36,8 +36,4 @@
    !
    ! ** albedo namelist (namalb)
-   INTEGER  ::   nn_ice_alb       ! type of albedo scheme: 0: Shine & Henderson-Sellers (JGR 1985)
-   !                                      !                         1: "home made" based on Brandt et al. (JClim 2005)
-   !                                      !                                             and Grenfell & Perovich (JGR 2004)
-   !                                      !                         2: Same as 1 but with melt ponds
    REAL(wp) ::   rn_alb_sdry      ! dry snow albedo
    REAL(wp) ::   rn_alb_smlt      ! melting snow albedo
@@ -53,5 +49,5 @@
 CONTAINS
 
-   SUBROUTINE ice_alb( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
+   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, pafrac_pnd, ph_pnd, palb_cs, palb_os )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE ice_alb  ***
@@ -59,10 +55,6 @@
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
       !!       
-      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
-      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
-      !!                  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:
+      !! ** Method  :   The scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                      and Grenfell & Perovich (JGR 2004)
       !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
       !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
@@ -76,5 +68,21 @@
       !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
       !!
-      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                     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_imlt = 0.50      ! bare puddled ice albedo
+      !!                        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
+      !!                        rn_alb_smlt = 0.72      ! melting snow
+      !!                        rn_alb_idry = 0.65      ! bare frozen ice
+      !!                     Brandt et al 2005 (East Antarctica)
+      !!                        rn_alb_sdry = 0.87      ! dry snow
+      !!                        rn_alb_smlt = 0.82      ! melting snow
+      !!                        rn_alb_idry = 0.54      ! bare frozen ice
+      !!
+      !! ** Note    :   The old parameterization from Shine & Henderson-Sellers (not here anymore) presented several misconstructions
       !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
       !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
@@ -87,218 +95,76 @@
       !!                Grenfell & Perovich 2004, JGR, vol 109 
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice       !  ice surface temperature (Kelvin)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_su        !  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
+      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
       !
       INTEGER  ::   ji, jj, jl                ! dummy loop indices
-      REAL(wp) ::   zswitch, z1_c1, z1_c2     ! local scalar
-      REAL(wp) ::   z1_href_pnd               !   -      -                 
-      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
-      !
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free)
-!! MV MP
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd        ! ponded sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice        ! bare sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw        ! snow-covered sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw      ! relative snow fraction
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice      ! relative ice fraction
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd      ! relative ice fraction (effective)
+      REAL(wp) ::   z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar
+      REAL(wp) ::   z1_href_pnd               ! inverse of the characteristic length scale (Lecomte et al. 2015)
+      REAL(wp) ::   zalb_pnd, zafrac_pnd      ! ponded sea ice albedo & relative pound fraction
+      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
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('icealb')
       !
-      !-----------------------------------------------------
-      !  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 )   ;   zalb(:,:,:) = rn_alb_imlt
-         ELSEWHERE                                           ;   zalb(:,:,:) = rn_alb_idry
-         END WHERE
-      ENDIF
-
-      SELECT CASE ( nn_ice_alb )       ! select a parameterization
-      !
-      !              !------------------------------------------
-      CASE( 0 )      !  Shine and Henderson-Sellers (1985)     !   (based on clear sky values)
-         !           !------------------------------------------
-         !
-         !                       ! Thickness-dependent bare ice albedo
-         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_it = zalb
-         ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
-         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )   ;   zalb_it = 0.2467 + 0.7049 * ph_ice              &
-            &                                                                  - 0.8608 * ph_ice * ph_ice     &
-            &                                                                  + 0.3812 * ph_ice * ph_ice * ph_ice
-         ELSEWHERE                                        ;   zalb_it = 0.1    + 3.6    * ph_ice
-         END WHERE
-         !
-         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo
-            z1_href_pnd = 1. / 0.05    ! inverse of the characteristic length scale (Lecomte et al. 2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
-            !                          ! Snow-free ice albedo is a function of pond fraction
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   
-               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
-            END WHERE
-         ENDIF 
-         !
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !                    ! Freezing snow
-                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > ppc2
-                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - ppc1 ) ) )
-                  zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
-                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * pp1_c1  )   &
-                     &    +           zswitch   * rn_alb_sdry  
-                     !
-                  !                    ! Melting snow
-                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > ppc2
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - ppc2 ) )
-                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * pp1_c2 )   &
-                     &      +         zswitch   *   rn_alb_smlt 
-                     !
-                  !                    ! 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
-                  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
+      z1_href_pnd = 0.05
+      z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+      z1_c2 = 1. / 0.05
+      z1_c3 = 1. / 0.02
+      z1_c4 = 1. / 0.03
+      !
+      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
+                  zafrac_pnd = pafrac_pnd(ji,jj,jl)
+                  zafrac_ice = 1._wp - zafrac_pnd
+               ELSE
+                  zafrac_snw = 1._wp      ! Snow fully "shades" melt ponds and ice
+                  zafrac_pnd = 0._wp
+                  zafrac_ice = 0._wp
+               ENDIF
+               !
+               !                       !--- Bare ice albedo (for hi > 150cm)
+               IF( zafrac_pnd > 0._wp ) 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
+               ENDIF
+               !                       !--- Bare ice albedo (for hi < 150cm)
+               IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN      ! 5cm < hi < 150cm
+                  zalb_ice = zalb_ice    + ( 0.18 - zalb_ice   ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice(ji,jj,jl)) )
+               ELSEIF( ph_ice(ji,jj,jl) <= 0.05 ) THEN                               ! 0cm < hi < 5cm
+                  zalb_ice = rn_alb_oce  + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl)
+               ENDIF
+               !
+               !                       !--- Snow-covered ice albedo (freezing, melting cases)
+               IF( pt_su(ji,jj,jl) < rt0_snow ) THEN
+                  zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c3 )
+               ELSE
+                  zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c4 )
+               ENDIF
+               !                       !--- Ponded ice albedo
+               IF( zafrac_pnd > 0._wp ) 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
+               !                       !--- 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
+               !
             END DO
          END DO
-         !
-         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + ppcloud       ! Oberhuber correction for overcast sky
-         !
-         !
-         !           !------------------------------------------
-      CASE( 1 )      !  New parameterization (2016)            !    (based on overcast sky values)
-         !           !------------------------------------------
-         !
-         !                    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
-         !                          rn_alb_smlt = 0.72      ! melting snow
-         !                          rn_alb_idry = 0.65      ! bare frozen ice
-         !                    Brandt et al 2005 (East Antarctica)
-         !                          rn_alb_sdry = 0.87      ! dry snow
-         !                          rn_alb_smlt = 0.82      ! melting snow
-         !                          rn_alb_idry = 0.54      ! bare frozen ice
-         ! 
-         !              !--- 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
-         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 = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
-         END WHERE
-         !
-         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo
-            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
-            !
-            !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )
-               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd
-            END WHERE
-         ENDIF 
-         !
-         !              !--- Overcast sky surface albedo (accounting for snow, ice melt ponds)
-         z1_c1 = 1. / 0.02
-         z1_c2 = 1. / 0.03
-         DO jl = 1, jpl
-            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 (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
-
-                  ! 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)
-
-              END DO
-            END DO
-         END DO
-         !
-         !              !--- 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 ); 
-         !
-         !
-         !           !------------------------------------------
-      CASE( 2 )      !  Fractional surface-based parameterization (2017)
-         !           !------------------------------------------
-         !              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   ! Snow fully "shades" melt ponds
-         ENDIF         
-         !
-         !                       !--- Specific snow fraction (for now, prescribed)
-         WHERE( ph_snw > 0._wp )   ;   zafrac_snw = 1.
-         ELSEWHERE                 ;   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
-         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_ice = zalb       + ( 0.18 - zalb       ) * z1_c1 *  &
-            &                                                                       ( LOG(1.5) - LOG(ph_ice) )
-         ELSEWHERE                                        ;   zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
-         END WHERE
-         !
-         z1_c1 = 1. / 0.02       !--- Snow-covered ice albedo (freezing, melting cases)
-         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 )
-         ELSEWHERE                    ;   zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 )
-         END WHERE
-         !
-         IF ( ld_pnd ) THEN      !--- Depth-dependent ponded ice albedo
-            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd * z1_href_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
+      END DO
+      !
+      palb_cs(:,:,:) = palb_os(:,:,:) - ( - 0.1010 * palb_os(:,:,:) * palb_os(:,:,:) + 0.1933 * palb_os(:,:,:) - 0.0148 )
       !
       IF( nn_timing == 1 )   CALL timing_stop('icealb')
@@ -317,5 +183,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namalb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
+      NAMELIST/namalb/ rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
       !!----------------------------------------------------------------------
       !
@@ -334,5 +200,4 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namalb:'
-         WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
          WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry
          WRITE(numout,*) '      albedo of melting snow               rn_alb_smlt = ', rn_alb_smlt
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8592)
@@ -130,5 +130,5 @@
 
       ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os )
+      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8592)
@@ -147,17 +147,13 @@
       END DO
 
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'a_ip'//'_htc'//zchar
-            z2d(:,:) = a_ip(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip
-            znam = 'v_ip'//'_htc'//zchar
-            z2d(:,:) = v_ip(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip
-         END DO
-      ENDIF
-      ! END MV MP 2016
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I2.2)') jl
+         znam = 'a_ip'//'_htc'//zchar
+         z2d(:,:) = a_ip(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip
+         znam = 'v_ip'//'_htc'//zchar
+         z2d(:,:) = v_ip(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip
+      END DO
 
       DO jl = 1, jpl 
@@ -248,17 +244,13 @@
       END DO
 
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'a_ip'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            a_ip(:,:,jl) = z2d(:,:)
-            znam = 'v_ip'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            v_ip(:,:,jl) = z2d(:,:)
-         END DO
-      ENDIF
-      ! END MV MP 2016
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I2.2)') jl
+         znam = 'a_ip'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         a_ip(:,:,jl) = z2d(:,:)
+         znam = 'v_ip'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         v_ip(:,:,jl) = z2d(:,:)
+      END DO
 
       DO jl = 1, jpl 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8592)
@@ -199,5 +199,5 @@
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
       !
       alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8592)
@@ -167,7 +167,7 @@
 !!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
 
-      !-------------------------------------------------------
-      ! Ice thickness, snow thickness, ice salinity, ice age
-      !-------------------------------------------------------
+      !---------------------------------------------------------------
+      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
+      !---------------------------------------------------------------
       !                                            !--- inverse of the ice area
       WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
@@ -178,26 +178,31 @@
       ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
       END WHERE
-      !
-      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
+      !                                           !--- ice thickness
+      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
 
       zhmax    =          hi_max(jpl)
       z1_zhmax =  1._wp / hi_max(jpl)               
-      WHERE( h_i(:,:,jpl) > zhmax )               !--- bound h_i by hi_max (i.e. 99 m) with associated update of ice area
+      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
          h_i  (:,:,jpl) = zhmax
          a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
-         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)          ! NB: v_i always /=0 as h_i > hi_max
+         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
       END WHERE
-
-      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)    !--- snow thickness
-      
-      o_i(:,:,:)  = oa_i(:,:,:) * z1_a_i(:,:,:)    !--- ice age
-
-      IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere)
+      !                                           !--- snow thickness
+      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
+      !                                           !--- ice age      
+      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
+      !                                           !--- pond fraction and thickness      
+      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
+      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
+      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
+      END WHERE
+      !
+      !                                           !---  salinity (with a minimum value imposed everywhere)     
+      IF( nn_icesal == 2 ) THEN
          WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
          ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
          END WHERE
       ENDIF
-
-      CALL ice_var_salprof      ! salinity profile
+      CALL ice_var_salprof   ! salinity profile
 
       !-------------------
@@ -242,8 +247,4 @@
       vt_s (:,:) = SUM( v_s, dim=3 )
       at_i (:,:) = SUM( a_i, dim=3 )
-
-! MV MP 2016
-! probably should resum for melt ponds ???
-
       !
    END SUBROUTINE ice_var_glo2eqv
Index: anches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/albedooce.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/albedooce.F90	(revision 8591)
+++ 	(revision )
@@ -1,49 +1,0 @@
-MODULE albedooce
-   !!======================================================================
-   !!                       ***  MODULE  albedooce  ***
-   !! Ocean forcing:  bulk thermohaline forcing of the ocean
-   !!=====================================================================
-   !! History :
-   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ocean and ice albedos
-   !!----------------------------------------------------------------------
-   !!   albedo_oce    : albedo for ocean (clear and overcast skies)
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constants
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   albedo_oce   ! routine called by sbccpl
-  
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
-   !! $Id: albedooce.F90 8268 2017-07-03 15:01:04Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE albedo_oce( pa_oce_os , pa_oce_cs )
-      !!----------------------------------------------------------------------
-      !!               ***  ROUTINE albedo_oce  ***
-      !! 
-      !! ** Purpose :   Computation of the albedo of the ocean
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_os   !  albedo of ocean under overcast sky
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
-      !!
-      REAL(wp) ::   zcoef 
-      REAL(wp) ::   rmue = 0.40    !  cosine of local solar altitude
-      !!----------------------------------------------------------------------
-      !
-      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
-      pa_oce_cs(:,:) = zcoef 
-      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
-      !
-   END SUBROUTINE albedo_oce
-
-   !!======================================================================
-END MODULE albedooce
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90	(revision 8592)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90	(revision 8592)
@@ -0,0 +1,49 @@
+MODULE ocealb
+   !!======================================================================
+   !!                       ***  MODULE  ocealb  ***
+   !! Ocean forcing:  bulk thermohaline forcing of the ocean
+   !!=====================================================================
+   !! History :
+   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ocean and ice albedos
+   !!----------------------------------------------------------------------
+   !!   oce_alb    : albedo for ocean (clear and overcast skies)
+   !!----------------------------------------------------------------------
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   oce_alb   ! routine called by sbccpl
+  
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
+   !! $Id: ocealb.F90 8268 2017-07-03 15:01:04Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE oce_alb( palb_os , palb_cs )
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE oce_alb  ***
+      !! 
+      !! ** Purpose :   Computation of the albedo of the ocean
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   palb_os   !  albedo of ocean under overcast sky
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   palb_cs   !  albedo of ocean under clear sky
+      !!
+      REAL(wp) ::   zcoef 
+      REAL(wp) ::   rmue = 0.40    !  cosine of local solar altitude
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
+      palb_cs(:,:) = zcoef 
+      palb_os(:,:) = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      !
+   END SUBROUTINE oce_alb
+
+   !!======================================================================
+END MODULE ocealb
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8591)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8592)
@@ -32,5 +32,5 @@
    USE geo2ocean      ! 
    USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
-   USE albedooce      ! 
+   USE ocealb         ! 
    USE eosbn2         ! 
    USE sbcrnf, ONLY : l_rnfcpl
@@ -177,5 +177,5 @@
    TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
@@ -201,5 +201,5 @@
       ierr(:) = 0
       !
-      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
+      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
       
 #if ! defined key_lim3 && ! defined key_cice
@@ -735,7 +735,7 @@
       !     2. receiving mixed oce-ice solar radiation 
       IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
-         CALL albedo_oce( zaos, zacs )
+         CALL oce_alb( zaos, zacs )
          ! Due to lack of information on nebulosity : mean clear/overcast sky
-         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
+         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
       ENDIF
 
@@ -1884,5 +1884,5 @@
 !       ( see OASIS3 user guide, 5th edition, p39 )
          zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
-            &            / (  1.- ( albedo_oce_mix(:,:  ) * ziceld(:,:)       &
+            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
             &                     + palbi         (:,:,1) * picefr(:,:) ) )
       END SELECT
@@ -2051,5 +2051,5 @@
                    ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
                 ELSEWHERE
-                   ztmp1(:,:) = albedo_oce_mix(:,:)
+                   ztmp1(:,:) = alb_oce_mix(:,:)
                 END WHERE
              CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
@@ -2079,5 +2079,5 @@
 
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
-         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
+         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
          DO jl=1,jpl
             ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
