Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref	(revision 12890)
@@ -56,9 +56,9 @@
    rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
    ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016
-      rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
+      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
                                       !          recommended range: [0.1 ; 0.25]
-      rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
-      rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1]
-      rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??]
+      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
+      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1]
+      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??]
 /
 !------------------------------------------------------------------------------
@@ -110,5 +110,9 @@
 !------------------------------------------------------------------------------
    rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-)
-   rn_blow_s        =   0.66          !  mesure of snow blowing into the leads
+   nn_snwfra        =   0             !  calculate the fraction of ice covered by snow (for zdf and albedo)
+                                      !     = 0  fraction = 1 (if snow) or 0 (if no snow)
+                                      !     = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation]
+                                      !     = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation]
+   rn_snwblow       =   0.66          !  mesure of snow blowing into the leads
                                       !     = 1 => no snow blowing, < 1 => some snow blowing
    nn_flxdist       =  -1             !  Redistribute heat flux over ice categories
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn	(revision 12890)
@@ -10,8 +10,8 @@
    rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
    ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016
-      rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
+      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
                                       !          recommended range: [0.1 ; 0.25]
-      rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
-      rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1]
-      rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??]
+      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
+      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1]
+      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??]
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90	(revision 12890)
@@ -141,8 +141,8 @@
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
    LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016 
-   REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice
-   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home) 
-   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction
-   REAL(wp), PUBLIC ::   rn_tensile       !:    isotropic tensile strength
+   REAL(wp), PUBLIC ::   rn_lf_depfra     !:    fraction of ocean depth that ice must reach to initiate landfast ice
+   REAL(wp), PUBLIC ::   rn_lf_bfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home) 
+   REAL(wp), PUBLIC ::   rn_lf_relax      !:    relaxation time scale (s-1) to reach static friction
+   REAL(wp), PUBLIC ::   rn_lf_tensile    !:    isotropic tensile strength
    !
    !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) **
@@ -163,5 +163,10 @@
    !                                     !!** ice-surface boundary conditions namelist (namsbc) **
                                           ! -- icethd_dh -- !
-   REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice
+   REAL(wp), PUBLIC ::   rn_snwblow       !: coef. for partitioning of snowfall between leads and sea ice
+                                          ! -- icethd_zdf and icealb -- !
+   INTEGER , PUBLIC ::   nn_snwfra        !: calculate the fraction of ice covered by snow
+   !                                      !   = 0  fraction = 1 (if snow) or 0 (if no snow)
+   !                                      !   = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation]
+   !                                      !   = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation]
                                           ! -- icethd -- !
    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90	(revision 12890)
@@ -223,5 +223,5 @@
       NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
          &             rn_ishlat ,                                                           &
-         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
+         &             ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile
       !!-------------------------------------------------------------------
       !
@@ -246,8 +246,8 @@
          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
          WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
-         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
-         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
-         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
+         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_lf_depfra    = ', rn_lf_depfra
+         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_lf_bfr       = ', rn_lf_bfr
+         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lf_relax     = ', rn_lf_relax
+         WRITE(numout,*) '         isotropic tensile strength                          rn_lf_tensile   = ', rn_lf_tensile
          WRITE(numout,*)
       ENDIF
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90	(revision 12890)
@@ -258,5 +258,5 @@
 
       ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
-      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
       ELSE                         ;   zkt = 0._wp
       ENDIF
@@ -329,12 +329,12 @@
                zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
                ! ice-bottom stress at U points
-               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
-               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj)
+               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
                ! ice-bottom stress at V points
-               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
-               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj)
+               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
                ! ice_bottom stress at T points
-               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
-               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
             END DO
          END DO
@@ -506,5 +506,5 @@
                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0
                         &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &           )   * zmsk00y(ji,jj)
@@ -513,5 +513,5 @@
                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
-                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                         &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &            )   * zmsk00y(ji,jj)
@@ -557,5 +557,5 @@
                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0
                         &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                         &           )   * zmsk00x(ji,jj)
@@ -564,5 +564,5 @@
                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
-                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                         &            )   * zmsk00x(ji,jj)
@@ -610,5 +610,5 @@
                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0
                         &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
                         &           )   * zmsk00x(ji,jj)
@@ -617,5 +617,5 @@
                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
-                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &            )   * zmsk00x(ji,jj)
@@ -661,5 +661,5 @@
                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
-                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0
                         &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &           )   * zmsk00y(ji,jj)
@@ -668,5 +668,5 @@
                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
-                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
                         &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
                         &            )   * zmsk00y(ji,jj)
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90	(revision 12890)
@@ -199,17 +199,17 @@
                si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
                si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
-            ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
-               si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
-            ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
-               si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
-            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s
-               si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
-            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i
-               si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
-            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su
-               si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
-            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i
-               si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
             ENDIF
+            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
+               &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
+            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
+               &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
+            IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s
+               &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
+            IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i
+               &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
+            IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su
+               &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
+            IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i
+               &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
             !
             ! pond concentration
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90	(revision 12890)
@@ -276,5 +276,5 @@
       INTEGER ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namsbc/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate
+      NAMELIST/namsbc/ rn_cio, rn_snwblow, nn_snwfra, nn_flxdist, ln_cndflx, ln_cndemulate
       !!-------------------------------------------------------------------
       !
@@ -293,5 +293,6 @@
          WRITE(numout,*) '   Namelist namsbc:'
          WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio        = ', rn_cio
-         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s     = ', rn_blow_s
+         WRITE(numout,*) '      fraction of ice covered by snow (options 0,1,2)  nn_snwfra     = ', nn_snwfra
+         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_snwblow    = ', rn_snwblow
          WRITE(numout,*) '      Multicategory heat flux formulation              nn_flxdist    = ', nn_flxdist
          WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_dh.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_dh.F90	(revision 12890)
@@ -13,6 +13,5 @@
    !!----------------------------------------------------------------------
    !!   ice_thd_dh        : vertical sea-ice growth and melt
-   !!   ice_thd_snwblow   : distribute snow fall between ice and ocean
-  !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -20,4 +19,5 @@
    USE ice1D          ! sea-ice: thermodynamics variables
    USE icethd_sal     ! sea-ice: salinity profiles
+   USE icevar         ! for CALL ice_var_snwblow
    !
    USE in_out_manager ! I/O manager
@@ -29,9 +29,4 @@
 
    PUBLIC   ice_thd_dh        ! called by ice_thd
-   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbccpl and here
-
-   INTERFACE ice_thd_snwblow
-      MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
-   END INTERFACE
 
    !!----------------------------------------------------------------------
@@ -186,5 +181,5 @@
       ! Snow precipitation
       !-------------------
-      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
+      CALL ice_var_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing
 
       zdeltah(1:npti,:) = 0._wp
@@ -636,33 +631,4 @@
    END SUBROUTINE ice_thd_dh
 
-
-   !!--------------------------------------------------------------------------
-   !! INTERFACE ice_thd_snwblow
-   !!
-   !! ** Purpose :   Compute distribution of precip over the ice
-   !!
-   !!                Snow accumulation in one thermodynamic time step
-   !!                snowfall is partitionned between leads and ice.
-   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
-   !!                but because of the winds, more snow falls on leads than on sea ice
-   !!                and a greater fraction (1-at_i)^beta of the total mass of snow 
-   !!                (beta < 1) falls in leads.
-   !!                In reality, beta depends on wind speed, 
-   !!                and should decrease with increasing wind speed but here, it is 
-   !!                considered as a constant. an average value is 0.66
-   !!--------------------------------------------------------------------------
-!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
-   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
-      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
-      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
-      pout = ( 1._wp - ( pin )**rn_blow_s )
-   END SUBROUTINE ice_thd_snwblow_2d
-
-   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
-      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
-      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
-      pout = ( 1._wp - ( pin )**rn_blow_s )
-   END SUBROUTINE ice_thd_snwblow_1d
-
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90	(revision 12890)
@@ -94,12 +94,11 @@
       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
       REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
-      REAL(wp) ::   zhs_min   =  0.1_wp       ! minimum snow thickness for conductivity calculation 
+      REAL(wp) ::   zhs_ssl   =  0.03_wp      ! surface scattering layer in the snow 
+      REAL(wp) ::   zhi_ssl   =  0.10_wp      ! surface scattering layer in the ice
       REAL(wp) ::   ztmelts                   ! ice melting temperature
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
       REAL(wp) ::   zcpi                      ! Ice specific heat
       REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
-      REAL(wp) ::   zfac                      ! dummy factor
-      !
-      REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only)
+      !
       REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
@@ -150,28 +149,14 @@
       !------------------
       DO ji = 1, npti
-
-         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead
-         !!IF( h_s_1d(ji) < zhs_min ) THEN
-         !!  isnow(ji) = h_s_1d(ji) / zhs_min
-         !!  zh_s(ji) = zhs_min * r1_nlay_s
-         !!ELSE
-         !!  isnow(ji) = 1.0_wp
-         !!  zh_s(ji) = h_s_1d(ji) * r1_nlay_s
-         !!END IF
-         isnow(ji) = za_s_fra(ji) !!clem: 2 variables for the same thing
-         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
-         
-         
-         ! layer thickness
-         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
-
+         zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction
+         zh_i(ji) = MAX( zhi_ssl , h_i_1d(ji) ) * r1_nlay_i ! set a minimum ice thickness for conduction
       END DO
       !
-      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
-      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
+      WHERE( h_i_1d(1:npti) >= zhi_ssl )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
+      ELSEWHERE                            ;   z1_h_i(1:npti) = 0._wp ! put 0 if ice thick < scattering layer
       END WHERE
       !
-      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
-      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
+      WHERE( h_s_1d(1:npti) >= zhs_ssl )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
+      ELSEWHERE                            ;   z1_h_s(1:npti) = 0._wp ! put 0 if snw thick < scattering layer
       END WHERE
       !
@@ -198,5 +183,5 @@
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) )
             !                             ! radiation absorbed by the layer-th snow layer
             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
@@ -204,9 +189,9 @@
       END DO
       !
-      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
+      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * (1._wp - za_s_fra(1:npti))
       DO jk = 1, nlay_i 
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
             !                             ! radiation absorbed by the layer-th ice layer
             zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
@@ -216,5 +201,5 @@
       qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
       !
-      iconv    = 0          ! number of iterations
+      iconv = 0          ! number of iterations
       !
       l_T_converged(:) = .FALSE.
@@ -243,5 +228,5 @@
                DO ji = 1, npti
                   ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
-                     &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+                     &                    MIN( -epsi10, 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 )
                END DO
             END DO
@@ -251,13 +236,13 @@
             DO ji = 1, npti
                ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
-                  &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
+                  &                            - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
                ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
-                  &                           - 0.011_wp * ( t_bo_1d(ji) - rt0 )
+                  &                            - 0.011_wp * ( t_bo_1d(ji) - rt0 )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1, npti
-                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
-                     &                        MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
-                     &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
+                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /       &
+                     &                         MIN( -epsi10, 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 ) &
+                     &                        - 0.011_wp * ( 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 )
                END DO
             END DO
@@ -304,7 +289,7 @@
          DO ji = 1, npti   ! Snow-ice interface
             IF ( .NOT. l_T_converged(ji) ) THEN
-               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
-               IF( zfac > epsi10 ) THEN
-                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
+               IF( h_s_1d(ji) >= zhs_ssl ) THEN
+                  zkappa_s(ji,nlay_s) =   zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / &
+                     &                  ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) )
                ELSE
                   zkappa_s(ji,nlay_s) = 0._wp
@@ -324,5 +309,5 @@
          DO ji = 1, npti   ! Snow-ice interface
             IF ( .NOT. l_T_converged(ji) ) &
-               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
+               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji))
          END DO
          !
@@ -558,5 +543,5 @@
                   IF( t_su_1d(ji) < rt0 ) THEN
                      t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
-                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
+                        &          ( za_s_fra(ji) * t_s_1d(ji,1) + (1._wp - za_s_fra(ji)) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
                   ENDIF
                ENDIF
@@ -790,6 +775,6 @@
          !
          DO ji = 1, npti
-            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
-               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
+            qcn_ice_top_1d(ji) =  -           za_s_fra(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
+               &                  - ( 1._wp - za_s_fra(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
          END DO
          !
@@ -806,7 +791,7 @@
          DO ji = 1, npti
             t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature
-               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &
-               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) &
-               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 )
+               &                +          za_s_fra(ji)  * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &
+               &                + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) &
+               &          ) / MAX( epsi10, za_s_fra(ji)  * zkappa_s(ji,0) * zg1s + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 )
             t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su
          END DO
@@ -866,11 +851,12 @@
       !--------------------------------------------------------------------
       ! effective conductivity and 1st layer temperature (needed by Met Office)
+      ! this is a conductivity at mid-layer, hence the factor 2
       DO ji = 1, npti
-         IF( h_i_1d(ji) > 0.1_wp ) THEN
+         IF( h_i_1d(ji) >= zhi_ssl ) THEN   ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0)
             cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
          ELSE
-            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp ! cnd_ice is capped by: cond_i/0.1m
+            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl
          ENDIF
-         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
+         t1_ice_1d(ji) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1)
       END DO
       !
@@ -886,8 +872,7 @@
       DO ji = 1, npti         
          !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-         IF( h_s_1d(ji) >= zhs_min ) THEN !!clem change that
-            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
+         IF( h_s_1d(ji) >= zhs_ssl ) THEN
+            t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / &
+               &          ( rn_cnd_s * zh_i(ji)                + ztcond_i(ji,1) * zh_s(ji)                )
          ELSE
             t_si_1d(ji) = t_su_1d(ji)
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90	(revision 12890)
@@ -52,4 +52,5 @@
    !!   ice_var_itd       : convert N-cat to M-cat
    !!   ice_var_snwfra    : fraction of ice covered by snow
+   !!   ice_var_snwblow   : distribute snow fall between ice and ocean
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -79,4 +80,5 @@
    PUBLIC   ice_var_itd
    PUBLIC   ice_var_snwfra
+   PUBLIC   ice_var_snwblow
 
    INTERFACE ice_var_itd
@@ -86,4 +88,8 @@
    INTERFACE ice_var_snwfra
       MODULE PROCEDURE ice_var_snwfra_1d, ice_var_snwfra_2d, ice_var_snwfra_3d
+   END INTERFACE
+
+   INTERFACE ice_var_snwblow
+      MODULE PROCEDURE ice_var_snwblow_1d, ice_var_snwblow_2d
    END INTERFACE
 
@@ -1009,5 +1015,4 @@
          pt_su(:,jl) = ptmsu(:)
          ps_i (:,jl) = psmi (:)
-         ps_i (:,jl) = psmi (:)         
       END DO
       !
@@ -1059,5 +1064,5 @@
       !!
       !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
-       !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
+      !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
       !!              
       !!               3) Expand the filling to the empty cat between jlmin and jlmax 
@@ -1293,27 +1298,72 @@
    !!-------------------------------------------------------------------
    SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra )
-      !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ph_s        ! snow thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
-      !!-------------------------------------------------------------------
-      pa_s_fra(:,:,:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:,:,:) )
+      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
+         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
+         ELSEWHERE             ; pa_s_fra = 0._wp
+         END WHERE
+      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
+         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
+      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
+         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
+      ENDIF
    END SUBROUTINE ice_var_snwfra_3d
 
    SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra )
-      !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
-      !!-------------------------------------------------------------------
-      pa_s_fra(:,:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:,:) )
+      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
+         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
+         ELSEWHERE             ; pa_s_fra = 0._wp
+         END WHERE
+      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
+         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
+      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
+         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
+      ENDIF
    END SUBROUTINE ice_var_snwfra_2d
 
    SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra )
-      !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness
       REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
-      !!-------------------------------------------------------------------
-      pa_s_fra(:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:) )
+      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
+         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
+         ELSEWHERE             ; pa_s_fra = 0._wp
+         END WHERE
+      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
+         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
+      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
+         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
+      ENDIF
    END SUBROUTINE ice_var_snwfra_1d
    
+   !!--------------------------------------------------------------------------
+   !! INTERFACE ice_var_snwblow
+   !!
+   !! ** Purpose :   Compute distribution of precip over the ice
+   !!
+   !!                Snow accumulation in one thermodynamic time step
+   !!                snowfall is partitionned between leads and ice.
+   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
+   !!                but because of the winds, more snow falls on leads than on sea ice
+   !!                and a greater fraction (1-at_i)^beta of the total mass of snow 
+   !!                (beta < 1) falls in leads.
+   !!                In reality, beta depends on wind speed, 
+   !!                and should decrease with increasing wind speed but here, it is 
+   !!                considered as a constant. an average value is 0.66
+   !!--------------------------------------------------------------------------
+!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
+   SUBROUTINE ice_var_snwblow_2d( pin, pout )
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
+      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
+      pout = ( 1._wp - ( pin )**rn_snwblow )
+   END SUBROUTINE ice_var_snwblow_2d
+
+   SUBROUTINE ice_var_snwblow_1d( pin, pout )
+      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
+      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
+      pout = ( 1._wp - ( pin )**rn_snwblow )
+   END SUBROUTINE ice_var_snwblow_1d
 
 #else
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90	(revision 12890)
@@ -47,5 +47,5 @@
 #if defined key_si3
    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
-   USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icevar         ! for CALL ice_var_snwblow
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
@@ -890,5 +890,5 @@
       ! --- evaporation minus precipitation --- !
       zsnw(:,:) = 0._wp
-      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
+      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
       emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90	(revision 12889)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90	(revision 12890)
@@ -41,5 +41,5 @@
 #endif
 #if defined key_si3
-   USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icevar         ! for CALL ice_var_snwblow
 #endif
    !
@@ -1730,5 +1730,5 @@
 
       ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
-      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw )
+      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw )
       
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
