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 13006)
+++ 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,*)
