Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg	(revision 13249)
@@ -353,5 +353,4 @@
 &namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
 !-----------------------------------------------------------------------
-      rn_eice     =   0       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4   
 /
 !!======================================================================
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 13249)
@@ -374,5 +374,4 @@
                                !        = 2 add a tke source just at the base of the ML
                                !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T)
-      rn_eice     =   0       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4   
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ref	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ref	(revision 13249)
@@ -681,4 +681,5 @@
    !
    ln_drgimp   = .true.    !  implicit top/bottom friction flag
+      ln_drgice_imp = .false. ! implicit ice-ocean drag
 /
 !-----------------------------------------------------------------------
@@ -1059,6 +1060,6 @@
                               !    = 0 no scaling under sea-ice
                               !    = 1 scaling with constant sea-ice thickness
-                              !    = 2  scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
-                              !    = 3  scaling with maximum sea-ice thickness
+                              !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
+                              !    = 3 scaling with maximum sea-ice thickness
       rn_mxlice   = 10.       ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1)
    rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
@@ -1074,5 +1075,9 @@
                               !        = 0  constant 10 m length scale
                               !        = 1  0.5m at the equator to 30m poleward of 40 degrees
-      rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4   
+   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking
+   !                       !           = 1 weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 weighted by 1-fr_i
+   !                       !           = 3 weighted by 1-MIN(1,4*fr_i)   
 /
 !-----------------------------------------------------------------------
@@ -1087,7 +1092,13 @@
    rn_charn      = 70000.  !  Charnock constant for wb induced roughness length
    rn_hsro       =  0.02   !  Minimum surface roughness
+   rn_hsri       =  0.03   !  Ice-ocean roughness
    rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met>1)
    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3)
-   !                             ! =3 requires ln_wave=T
+   !                       !     = 3 requires ln_wave=T
+   nn_z0_ice     =   1     !  attenutaion of surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover
+   !                       !           = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i
+   !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i)
    nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
    nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SPITZ12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SPITZ12/EXPREF/namelist_cfg	(revision 13249)
@@ -217,4 +217,5 @@
    ln_loglayer = .true.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U|
    ln_drgimp   = .true.   !  implicit top/bottom friction flag
+      ln_drgice_imp = .true. ! implicit ice-ocean drag
 /
 !-----------------------------------------------------------------------
@@ -340,4 +341,19 @@
    nn_havtb    =    1         !  horizontal shape for avtb (=1) or not (=0)
 /
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
+!-----------------------------------------------------------------------
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+      nn_mxlice    = 2        ! type of scaling under sea-ice
+                              !    = 0 no scaling under sea-ice
+                              !    = 1 scaling with constant sea-ice thickness
+                              !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
+                              !    = 3 scaling with maximum sea-ice thickness
+   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking
+   !                       !           = 1 weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 weighted by 1-fr_i
+   !                       !           = 3 weighted by 1-MIN(1,4*fr_i)
+/
 !!======================================================================
 !!                  ***  Diagnostics namelists  ***                   !!
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdrg
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdrg	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdrg	(revision 13249)
@@ -8,3 +8,4 @@
    !
    ln_drgimp   = .true.    !  implicit top/bottom friction flag
+      ln_drgice_imp = .false. ! implicit ice-ocean drag
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namzdf_gls
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namzdf_gls	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namzdf_gls	(revision 13249)
@@ -13,4 +13,9 @@
    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3)
    !                             ! =3 requires ln_wave=T
+   nn_z0_ice     =   1     !  attenutaion of surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover
+   !                       !           = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i
+   !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i)
    nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
    nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namzdf_tke
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namzdf_tke	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namzdf_tke	(revision 13249)
@@ -26,4 +26,8 @@
                               !        = 0  constant 10 m length scale
                               !        = 1  0.5m at the equator to 30m poleward of 40 degrees
-      rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4   
+   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking
+   !                       !           = 1 weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 weighted by 1-fr_i
+   !                       !           = 3 weighted by 1-MIN(1,4*fr_i)
 /
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceupdate.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceupdate.F90	(revision 13249)
@@ -26,4 +26,5 @@
    USE icectl         ! sea-ice: control prints
    USE bdy_oce , ONLY : ln_bdy
+   USE zdfdrg  , ONLY : ln_drgice_imp
    !
    USE in_out_manager ! I/O manager
@@ -321,4 +322,5 @@
       REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
       REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
+      REAL(wp) ::   zflagi                          !   -      -
       !!---------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('ice_update_tau')
@@ -353,4 +355,11 @@
       !
       !                                      !==  every ocean time-step  ==!
+      IF ( ln_drgice_imp ) THEN
+         ! Save drag with right sign to update top drag in the ocean implicit friction 
+         rCdU_ice(:,:) = -r1_rau0 * tmod_io(:,:) * at_i(:,:) * tmask(:,:,1) 
+         zflagi = 0._wp
+      ELSE
+         zflagi = 1._wp
+      ENDIF
       !
       DO jj = 2, jpjm1                                !* update the stress WITHOUT an ice-ocean rotation angle
@@ -362,6 +371,6 @@
                &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji  ,jj+1,1) )
             !                                                   ! linearized quadratic drag formulation
-            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
-            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - zflagi * pu_oce(ji,jj) )
+            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - zflagi * pv_oce(ji,jj) )
             !                                                   ! stresses at the ocean surface
             utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynnxt.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynnxt.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynnxt.F90	(revision 13249)
@@ -34,5 +34,5 @@
    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
    USE domvvl         ! variable volume
-   USE bdy_oce   , ONLY: ln_bdy
+   USE bdy_oce , ONLY : ln_bdy
    USE bdydta         ! ocean open boundary conditions
    USE bdydyn         ! ocean open boundary conditions
@@ -48,4 +48,5 @@
    USE prtctl         ! Print control
    USE timing         ! Timing
+   USE zdfdrg ,  ONLY : ln_drgice_imp, rCdU_top
 #if defined key_agrif
    USE agrif_oce_interp
@@ -99,4 +100,5 @@
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
@@ -354,4 +356,40 @@
       ENDIF
       !
+      IF ( iom_use("utau") ) THEN
+         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
+            ALLOCATE(zutau(jpi,jpj)) 
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  jk = miku(ji,jj) 
+                  zutau(ji,jj) = utau(ji,jj) & 
+                  &  + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk) 
+               END DO
+            END DO
+            CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
+            CALL iom_put(  "utau", zutau(:,:) )
+            DEALLOCATE(zutau)
+         ELSE
+            CALL iom_put(  "utau", utau(:,:) )
+         ENDIF
+      ENDIF
+      !
+      IF ( iom_use("vtau") ) THEN
+         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
+            ALLOCATE(zvtau(jpi,jpj))
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  jk = mikv(ji,jj)
+                  zvtau(ji,jj) = vtau(ji,jj) &
+                  &  + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
+               END DO
+            END DO
+            CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
+            CALL iom_put(  "vtau", zvtau(:,:) )
+            DEALLOCATE(zvtau)
+         ELSE
+            CALL iom_put(  "vtau", vtau(:,:) )
+         ENDIF
+      ENDIF
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynspg_ts.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynspg_ts.F90	(revision 13249)
@@ -1465,5 +1465,5 @@
       !                    !==  Set the barotropic drag coef.  ==!
       !
-      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
          
          DO jj = 2, jpjm1
@@ -1528,5 +1528,5 @@
       !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
       !
-      IF( ln_isfcav ) THEN
+      IF( ln_isfcav.OR.ln_drgice_imp ) THEN
          !
          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynzdf.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/DYN/dynzdf.F90	(revision 13249)
@@ -141,5 +141,5 @@
             END DO
          END DO
-         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
+         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
             DO jj = 2, jpjm1        
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -258,5 +258,5 @@
             END DO
          END DO
-         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
+         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
@@ -423,5 +423,5 @@
             END DO
          END DO
-         IF ( ln_isfcav ) THEN
+         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbc_ice.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbc_ice.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbc_ice.F90	(revision 13249)
@@ -70,4 +70,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: sea surface freezing temperature            [degC]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rCdU_ice       !: ice-ocean drag at T-point (<0)               [m/s]
 #endif
 
@@ -132,5 +133,5 @@
          &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce   (jpi,jpj)     ,   &
          &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce    (jpi,jpj)     ,   &
-         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , STAT= ierr(2) )
+         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , rCdU_ice   (jpi,jpj)     , STAT= ierr(2) )
 #endif
 
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcmod.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcmod.F90	(revision 13249)
@@ -564,7 +564,4 @@
       ENDIF
       !
-      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice)
-      CALL iom_put( "vtau", vtau )   ! j-wind stress
-      !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
          CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask )
Index: /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfdrg.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfdrg.F90	(revision 13249)
@@ -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_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp
+      NAMELIST/namdrg/ ln_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/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfgls.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfgls.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfgls.F90	(revision 13249)
@@ -53,4 +53,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
@@ -61,4 +62,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) 
 
@@ -150,4 +152,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
@@ -165,4 +168,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          
@@ -206,4 +216,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
@@ -300,5 +313,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
@@ -306,6 +319,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
@@ -316,5 +329,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
@@ -326,5 +339,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)
@@ -577,5 +590,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.)
@@ -850,8 +864,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
       !!----------------------------------------------------------
@@ -881,8 +895,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/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfphy.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfphy.F90	(revision 13249)
@@ -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/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdftke.F90
===================================================================
--- /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdftke.F90	(revision 13248)
+++ /NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdftke.F90	(revision 13249)
@@ -83,7 +83,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)
@@ -199,31 +199,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
@@ -257,5 +267,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
@@ -295,6 +305,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         
@@ -302,11 +311,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
@@ -408,9 +417,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
@@ -421,5 +430,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
@@ -434,5 +443,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
@@ -499,10 +508,10 @@
       zmxlm(:,:,:)  = rmxl_min    
       zmxld(:,:,:)  = rmxl_min
-      !
-     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
+      ! 
+      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
          !
          zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
 #if ! defined key_si3 && ! defined key_cice
-         DO jj = 2, jpjm1
+         DO jj = 2, jpjm1                     ! No sea-ice
             DO ji = fs_2, fs_jpim1
                zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
@@ -510,4 +519,5 @@
          END DO
 #else
+
          SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
          !
@@ -519,28 +529,32 @@
             END DO
             !
-         CASE( 1 )                           ! scaling with constant sea-ice thickness
+         CASE( 1 )                      ! scaling with constant sea-ice thickness
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1)
+                  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
+         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. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1)
+                  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. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
+                  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
+         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. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
+                  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
@@ -704,5 +718,5 @@
          &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
          &                 nn_pdl  , ln_drg   , ln_lc    , rn_lc,      &
-         &                 nn_etau , nn_htau  , rn_efr   , rn_eice  
+         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  
       !!----------------------------------------------------------------------
       !
@@ -737,4 +751,12 @@
             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,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg
@@ -744,6 +766,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( ln_drg ) THEN
             WRITE(numout,*)
