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 8426)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8486)
@@ -4,15 +4,14 @@
    !! Ocean forcing:  bulk thermohaline forcing of the ice
    !!=====================================================================
-   !! History :
-   !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos
-   !!----------------------------------------------------------------------
-   !!   ice_alb    : albedo for ice (clear and overcast skies)
-   !!   alb_init   : initialisation of albedo computation
+   !! History :  4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos
    !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
-   !!----------------------------------------------------------------------
-   USE ice, ONLY : jpl
+   !!   'key_lim3'                                    LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_alb       : albedo for ice (clear and overcast skies)
+   !!   alb_init      : initialisation of albedo computation
+   !!----------------------------------------------------------------------
+   USE ice     , ONLY : jpl   ! number of ice category
    USE phycst         ! physical constants
    !
@@ -24,18 +23,22 @@
    PRIVATE
 
-   PUBLIC   ice_alb   ! routine called icestp.F90
-
-   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   PUBLIC   ice_alb   ! routine called in iceforcing.F90 and iceupdate.F90
+
+   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+
+!!gm real variable name starting with a "c" NOT DOCTOR !!!!
    INTEGER  ::   albd_init = 0       ! control flag for initialization  
-   REAL(wp) ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0)
-   REAL(wp) ::   c2        = 0.10    !  "        "
-   REAL(wp) ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   c2        = 0.10    !  "        "
+   REAL(wp) , PARAMETER ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   r1_c1 = 1. / c1
+   REAL(wp) , PARAMETER ::   r1_c2 = 1. / c2
  
-   !                             !!* namelist namsbc_alb
+   !                             !!* namelist namsbc_alb *
    INTEGER  ::   nn_ice_alb
    REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2010)
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
    !! $Id: icealb.F90 8268 2017-07-03 15:01:04Z clem $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -76,30 +79,27 @@
       !!                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
-      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
-      REAL(wp) ::   zswitch, z1_c1, z1_c2
-      REAL(wp) ::   zhref_pnd                                 
+      !!----------------------------------------------------------------------
+      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
+      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)
+      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), 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)
       !!---------------------------------------------------------------------
 
@@ -120,14 +120,11 @@
       ENDIF
 
-      SELECT CASE ( nn_ice_alb )
-
-      !------------------------------------------
-      !  Shine and Henderson-Sellers (1985)
-      !------------------------------------------
-      ! NB: This parameterization is based on clear sky values
-
-      CASE( 0 )
-       
-         ! Thickness-dependent bare ice albedo
+      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
          ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
@@ -137,35 +134,37 @@
          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
+         !
+         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_ice )   
+               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
+            END WHERE
          ENDIF 
-
+         !
+!!gm: optimization ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 can be computed one for all 
+!!gm  before the DO-loop below
          DO jl = 1, jpl
             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 ) ) )
-                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
-                     &                           + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) / c1  )   &
-                     &        +         zswitch   * rn_alb_sdry  
-
-                  ! Melting snow
+                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
+                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1  )   &
+                     &    +           zswitch   * rn_alb_sdry  
+                     !
+                  !                    ! 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 ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) / c2 )   &
-                      &     +         zswitch   *   rn_alb_smlt 
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_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
                   !
-                  ! 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)
@@ -174,53 +173,50 @@
             END DO
          END DO
-
+         !
          pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
-
-      !------------------------------------------
-      !  New parameterization (2016)
-      !------------------------------------------
-      ! NB: This parameterization is based on overcast skies values
-      
-      CASE( 1 ) 
-
-! 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 
+         !
+         !
+         !           !------------------------------------------
+      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
-            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
+         !              !--- 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_ice )
+               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
-         
-         ! Overcast sky surface albedo (accounting for snow, ice melt ponds)
          DO jl = 1, jpl
             DO jj = 1, jpj
@@ -241,69 +237,62 @@
             END DO
          END DO
-
-         ! Clear sky surface 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 ); 
-
-      !---------------------------------------------------
-      !  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 *  &
+         !
+         !
+         !           !------------------------------------------
+      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 ;  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 = rn_alb_oce + ( 0.18 - rn_alb_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 )
+         ELSE WHERE                                       ;  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 );
+         ELSE WHERE                 ; 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
@@ -311,4 +300,5 @@
    END SUBROUTINE ice_alb
 
+
    SUBROUTINE alb_init
       !!----------------------------------------------------------------------
@@ -320,4 +310,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
+      !!
       NAMELIST/namicealb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
       !!----------------------------------------------------------------------
@@ -353,9 +344,4 @@
    !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   !
-   SUBROUTINE ice_alb     ! Dummy routine
-      WRITE(*,*) 'ice_alb: You should not have seen this print! error?'
-   END SUBROUTINE ice_alb
 #endif
 
