Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfdrg.F90	(revision 13268)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfdrg.F90	(revision 13284)
@@ -32,4 +32,5 @@
    USE lib_mpp        ! distributed memory computing
    USE prtctl         ! Print control
+   USE sbc_oce , ONLY : nn_ice 
 
    IMPLICIT NONE
@@ -46,5 +47,5 @@
    LOGICAL          ::   ln_loglayer  ! logarithmic drag: Cd = vkarmn/log(z/z0)
    LOGICAL , PUBLIC ::   ln_drgimp    ! implicit top/bottom friction flag
-
+   LOGICAL , PUBLIC ::   ln_drgice_imp ! implicit ice-ocean drag 
    !                                 !!* Namelist namdrg_top & _bot: TOP or BOTTOM coefficient namelist *
    REAL(wp)         ::   rn_Cd0       !: drag coefficient                                           [ - ]
@@ -231,5 +232,5 @@
       INTEGER   ::   ios, ioptio   ! local integers
       !!
-      NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp
+      NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp, ln_drgice_imp
       !!----------------------------------------------------------------------
       !
@@ -254,4 +255,5 @@
          WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer
          WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp
+         WRITE(numout,*) '      implicit ice-ocean drag                   ln_drgice_imp  =', ln_drgice_imp
       ENDIF
       !
@@ -264,4 +266,9 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' )
       !
+      IF ( ln_drgice_imp.AND.(.NOT.ln_drgimp) ) & 
+         &                CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires ln_drgimp=T' )
+      !
+      IF ( ln_drgice_imp.AND.( nn_ice /=2 ) ) &
+         &  CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires si3' )
       !
       !                     !==  BOTTOM drag setting  ==!   (applied at seafloor)
@@ -274,6 +281,10 @@
       !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities)
       !
-      IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting
-         ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) )
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN              ! Ocean cavities: top friction setting
+         ALLOCATE( rCdU_top(jpi,jpj) )
+      ENDIF
+      !
+      IF( ln_isfcav ) THEN
+         ALLOCATE( rCd0_top(jpi,jpj))
          CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in
             &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfgls.F90	(revision 13268)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfgls.F90	(revision 13284)
@@ -54,4 +54,5 @@
    INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
    INTEGER  ::   nn_z0_met         ! Method for surface roughness computation
+   INTEGER  ::   nn_z0_ice         ! Roughness accounting for sea ice
    INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
    INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
@@ -62,4 +63,5 @@
    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
+   REAL(wp) ::   rn_hsri           ! Ice ocean roughness
    REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
 
@@ -151,4 +153,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves 
       REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves)
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before
@@ -166,4 +169,11 @@
       ustar2_bot (:,:) = 0._wp
 
+      SELECT CASE ( nn_z0_ice )
+      CASE( 0 )   ;   zice_fra(:,:) = 0._wp
+      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp )
+      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:)
+      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp )
+      END SELECT
+      
       ! Compute surface, top and bottom friction at T-points
       DO jj = 2, jpjm1              !==  surface ocean friction
@@ -211,4 +221,7 @@
       END SELECT
       !
+      ! adapt roughness where there is sea ice
+      zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
+      !
       DO jk = 2, jpkm1              !==  Compute dissipation rate  ==!
          DO jj = 1, jpjm1
@@ -305,5 +318,5 @@
       CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2) 
       ! First level
-      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  )
+      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  )
       zd_lw(:,:,1) = en(:,:,1)
       zd_up(:,:,1) = 0._wp
@@ -311,6 +324,6 @@
       ! 
       ! One level below
-      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   &
-         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   )
+      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) &
+         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   )
       zd_lw(:,:,2) = 0._wp 
       zd_up(:,:,2) = 0._wp
@@ -321,5 +334,5 @@
       !
       ! Dirichlet conditions at k=1
-      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  )
+      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  )
       zd_lw(:,:,1) = en(:,:,1)
       zd_up(:,:,1) = 0._wp
@@ -331,5 +344,5 @@
       zd_lw(:,:,2) = 0._wp
       zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
-      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
+      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
           &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
 !!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)
@@ -582,5 +595,6 @@
          zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
          zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
-         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
+         zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &
+            &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
          zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &
             &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.)
@@ -855,8 +869,8 @@
       REAL(wp)::   zcr   ! local scalar
       !!
-      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
-         &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
-         &            rn_crban, rn_charn, rn_frac_hs,        &
-         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
+      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       &
+         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   &
+         &            rn_crban, rn_charn, rn_frac_hs,              &
+         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, &
          &            nn_stab_func, nn_clos
       !!----------------------------------------------------------
@@ -886,8 +900,18 @@
          WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
          WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
+         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice
+         SELECT CASE( nn_z0_ice )
+         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking'
+         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )'
+         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)'
+         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
+         CASE DEFAULT
+            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3')
+         END SELECT
          WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
          WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
          WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
          WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
+         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri
          WRITE(numout,*)
          WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfphy.F90	(revision 13268)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfphy.F90	(revision 13284)
@@ -28,4 +28,5 @@
    USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test)
    USE sbcrnf         ! surface boundary condition: runoff variables
+   USE sbc_ice        ! sea ice drag
 #if defined key_agrif
    USE agrif_oce_interp   ! interpavm
@@ -252,4 +253,14 @@
       ENDIF
       !
+#if defined key_si3
+      IF ( ln_drgice_imp) THEN
+         IF ( ln_isfcav ) THEN
+            rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:)
+         ELSE
+            rCdU_top(:,:) = rCdU_ice(:,:)
+         ENDIF
+      ENDIF
+#endif
+      ! 
       !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
       !
Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdftke.F90	(revision 13268)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdftke.F90	(revision 13284)
@@ -46,4 +46,10 @@
    USE zdfmxl         ! vertical physics: mixed layer
    !
+#if defined key_si3
+   USE ice, ONLY: hm_i, h_i
+#endif
+#if defined key_cice
+   USE sbc_ice, ONLY: h_i
+#endif
    USE in_out_manager ! I/O manager
    USE iom            ! I/O manager library
@@ -62,4 +68,6 @@
    !                      !!** Namelist  namzdf_tke  **
    LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
+   INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3)
+   REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice
    INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
    REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
@@ -74,7 +82,7 @@
    INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
    REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
-   REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4   
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
    REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells
+   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)   
 
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
@@ -190,31 +198,41 @@
       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
       !
-      INTEGER ::   ji, jj, jk              ! dummy loop arguments
+      INTEGER ::   ji, jj, jk                  ! dummy loop arguments
       REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars
       REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient
-      REAL(wp) ::   zbbrau, zri                ! local scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3     !   -         -
-      REAL(wp) ::   ztx2  , zty2  , zcof       !   -         -
-      REAL(wp) ::   ztau  , zdif               !   -         -
-      REAL(wp) ::   zus   , zwlc  , zind       !   -         -
-      REAL(wp) ::   zzd_up, zzd_lw             !   -         -
+      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars
+      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      -
+      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      -
+      REAL(wp) ::   ztau  , zdif               !   -      -
+      REAL(wp) ::   zus   , zwlc  , zind       !   -      -
+      REAL(wp) ::   zzd_up, zzd_lw             !   -      -
       INTEGER , DIMENSION(jpi,jpj)     ::   imlc
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
       !
-      zbbrau = rn_ebb / rau0       ! Local constant initialisation
-      zfact1 = -.5_wp * rdt 
-      zfact2 = 1.5_wp * rdt * rn_ediss
-      zfact3 = 0.5_wp       * rn_ediss
+      zbbrau  =  rn_ebb / rau0       ! Local constant initialisation
+      zbbirau =  3.75_wp / rau0
+      zfact1  = -0.5_wp * rdt 
+      zfact2  =  1.5_wp * rdt * rn_ediss
+      zfact3  =  0.5_wp       * rn_ediss
+      !
+      ! ice fraction considered for attenuation of langmuir & wave breaking
+      SELECT CASE ( nn_eice )
+      CASE( 0 )   ;   zice_fra(:,:) = 0._wp
+      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp )
+      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:)
+      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp )
+      END SELECT
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     !  Surface/top/bottom boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      
+      !
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + &
+               &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1)
          END DO
       END DO
@@ -248,5 +266,5 @@
                   zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
                      &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
-                  en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &
+                  en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &     
                      &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
                END DO
@@ -286,6 +304,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
-               zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
-               IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0.
+               zus3(ji,jj) = ( 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
             END DO
          END DO         
@@ -293,11 +310,11 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF ( zfr_i(ji,jj) /= 0. ) THEN               
+                  IF ( zus3(ji,jj) /= 0. ) THEN               
                      ! vertical velocity due to LC   
                      IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
                         !                                           ! vertical velocity due to LC
-                        zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i
+                        zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )
                         !                                           ! TKE Langmuir circulation source term
-                        en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
+                        en(ji,jj,jk) = en(ji,jj,jk) + rdt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
                      ENDIF
                   ENDIF
@@ -399,9 +416,9 @@
       
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
-         DO jk = 2, jpkm1                       ! rn_eice =0 ON below sea-ice, =4 OFF when ice fraction > 0.25
+         DO jk = 2, jpkm1                       ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -412,5 +429,5 @@
                jk = nmln(ji,jj)
                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                  &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
          END DO
@@ -425,5 +442,5 @@
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -477,5 +494,5 @@
       REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
       REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
-      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      -
+      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
       !!--------------------------------------------------------------------
@@ -490,13 +507,66 @@
       zmxlm(:,:,:)  = rmxl_min    
       zmxld(:,:,:)  = rmxl_min
-      !
+      ! 
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
+         !
          zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
-         DO jj = 2, jpjm1
+#if ! defined key_si3 && ! defined key_cice
+         DO jj = 2, jpjm1                     ! No sea-ice
             DO ji = fs_2, fs_jpim1
-               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
-            END DO
-         END DO
-      ELSE 
+               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
+            END DO
+         END DO
+#else
+
+         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
+         !
+         CASE( 0 )                      ! No scaling under sea-ice
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
+               END DO
+            END DO
+            !
+         CASE( 1 )                      ! scaling with constant sea-ice thickness
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
+               END DO
+            END DO
+            !
+         CASE( 2 )                      ! scaling with mean sea-ice thickness
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+#if defined key_si3
+                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
+#elif defined key_cice
+                  zmaxice = MAXVAL( h_i(ji,jj,:) )
+                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
+#endif
+               END DO
+            END DO
+            !
+         CASE( 3 )                      ! scaling with max sea-ice thickness
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zmaxice = MAXVAL( h_i(ji,jj,:) )
+                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
+               END DO
+            END DO
+            !
+         END SELECT
+#endif
+         !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
+            END DO
+         END DO
+         !
+      ELSE
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
@@ -643,8 +713,9 @@
       INTEGER ::   ios
       !!
-      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,  &
-         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,  &
-         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc,      &
-         &                 nn_etau , nn_htau  , rn_efr , rn_eice  
+      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  &
+         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
+         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
+         &                 nn_pdl  , ln_lc    , rn_lc,                 &
+         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  
       !!----------------------------------------------------------------------
       !
@@ -675,4 +746,17 @@
          WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
          WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
+         IF( ln_mxl0 ) THEN
+            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice
+            IF( nn_mxlice == 1 ) &
+            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice
+            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
+            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice'
+            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness'
+            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness'
+            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness'
+            CASE DEFAULT
+               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4')
+            END SELECT
+         ENDIF
          WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
          WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
@@ -680,6 +764,13 @@
          WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
          WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
-         WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice
-         WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   '
+         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice
+         SELECT CASE( nn_eice ) 
+         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking'
+         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )'
+         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)'
+         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
+         CASE DEFAULT
+            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3')
+         END SELECT      
          IF( .NOT.ln_drg_OFF ) THEN
             WRITE(numout,*)
