Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdftke.F90	(revision 13277)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdftke.F90	(revision 13899)
@@ -28,5 +28,5 @@
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
    !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only 
-   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition (ln_drg)
+   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition
    !!----------------------------------------------------------------------
 
@@ -68,8 +68,8 @@
    !                      !!** 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]
-   INTEGER  ::      nn_mxlice ! type of scaling under sea-ice
-   REAL(wp) ::      rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1)
    INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
    REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
@@ -79,11 +79,10 @@
    REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
    REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
-   LOGICAL  ::   ln_drg    ! top/bottom friction forcing flag 
    INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
    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)
@@ -200,29 +199,42 @@
       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 / rho0       ! Local constant initialisation
-      zfact1 = -.5_wp * rn_Dt 
-      zfact2 = 1.5_wp * rn_Dt * rn_ediss
-      zfact3 = 0.5_wp       * rn_ediss
+      zbbrau  = rn_ebb / rho0       ! Local constant initialisation
+      zbbirau = 3.75_wp / rho0
+      zfact1  = -.5_wp * rn_Dt 
+      zfact2  = 1.5_wp * rn_Dt * 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_2D_00_00
+      !
+      DO_2D( 0, 0, 0, 0 )         ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+!! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly
+!!       one way around would be to increase zbbirau 
+!!          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)
          en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
       END_2D
@@ -236,7 +248,7 @@
       ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
       !
-      IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE
-         !
-         DO_2D_00_00
+      IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE
+         !
+         DO_2D( 0, 0, 0, 0 )        ! bottom friction
             zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
@@ -246,6 +258,6 @@
             en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
          END_2D
-         IF( ln_isfcav ) THEN       ! top friction
-            DO_2D_00_00
+         IF( ln_isfcav ) THEN
+            DO_2D( 0, 0, 0, 0 )     ! top friction
                zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
@@ -274,26 +286,24 @@
          zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
-         DO_3DS_11_11( jpkm1, 2, -1 )
-            zus  = zcof * taum(ji,jj)
+         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! Last w-level at which zpelc>=0.5*us*us 
+            zus = zcof * taum(ji,jj)          !      with us=0.016*wind(starting from jpk-1)
             IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
          END_3D
          !                               ! finite LC depth
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
          END_2D
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             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) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
          END_2D
-         DO_3D_00_00( 2, jpkm1 )
-            IF ( zfr_i(ji,jj) /= 0. ) THEN               
-               ! vertical velocity due to LC   
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
+            IF ( zus3(ji,jj) /= 0._wp ) THEN               
                IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
                   !                                           ! vertical velocity due to LC
-                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i
+                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
                ENDIF
             ENDIF
@@ -309,6 +319,6 @@
       !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
-         DO_3D_00_00( 2, jpkm1 )
+      IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             !                             ! local Richardson number
             IF (rn2b(ji,jj,jk) <= 0.0_wp) then
@@ -322,14 +332,12 @@
       ENDIF
       !         
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en
          zcof   = zfact1 * tmask(ji,jj,jk)
          !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
          !                                   ! eddy coefficient (ensure numerical stability)
          zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
-            &          /    (  e3t(ji,jj,jk  ,Kmm)   &
-            &                * e3w(ji,jj,jk  ,Kmm)  )
+            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
          zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
-            &          /    (  e3t(ji,jj,jk-1,Kmm)   &
-            &                * e3w(ji,jj,jk  ,Kmm)  )
+            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
          !
          zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
@@ -344,20 +352,20 @@
       END_3D
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                          ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
       END_2D
-      DO_3D_00_00( 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
       END_3D
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
       END_2D
-      DO_3DS_00_00( jpk-2, 2, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke
          en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
       END_3D
@@ -368,19 +376,20 @@
 !!gm BUG : in the exp  remove the depth of ssh !!!
 !!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
-      
-      
+      !
+      ! penetration is partly switched off below sea-ice if nn_eice/=0
+      !
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
             en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
-               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
          END_3D
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             jk = nmln(ji,jj)
             en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
-               &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
          END_2D
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -389,5 +398,5 @@
             zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
             en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
-               &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
          END_3D
       ENDIF
@@ -451,10 +460,10 @@
       zmxlm(:,:,:)  = rmxl_min    
       zmxld(:,:,:)  = rmxl_min
-      !
+      ! 
      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
          !
          zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
 #if ! defined key_si3 && ! defined key_cice
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
             zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
          END_2D
@@ -463,27 +472,31 @@
          !
          CASE( 0 )                      ! No scaling under sea-ice
-            DO_2D_00_00
+            DO_2D( 0, 0, 0, 0 )
                zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
             END_2D
             !
-         CASE( 1 )                           ! scaling with constant sea-ice thickness
-            DO_2D_00_00
-               zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1)
+         CASE( 1 )                      ! scaling with constant sea-ice thickness
+            DO_2D( 0, 0, 0, 0 )
+               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_2D
             !
-         CASE( 2 )                                 ! scaling with mean sea-ice thickness
-            DO_2D_00_00
+         CASE( 2 )                      ! scaling with mean sea-ice thickness
+            DO_2D( 0, 0, 0, 0 )
 #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_2D
             !
-         CASE( 3 )                                 ! scaling with max sea-ice thickness
-            DO_2D_00_00
+         CASE( 3 )                      ! scaling with max sea-ice thickness
+            DO_2D( 0, 0, 0, 0 )
                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_2D
             !
@@ -491,5 +504,5 @@
 #endif
          !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
          END_2D
@@ -500,5 +513,5 @@
 
       !
-      DO_3D_00_00( 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zrn2 = MAX( rn2(ji,jj,jk), rsmall )
          zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
@@ -515,5 +528,5 @@
       ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
       CASE ( 0 )           ! bounded by the distance to surface and bottom
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
             &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
@@ -526,5 +539,5 @@
          !
       CASE ( 1 )           ! bounded by the vertical scale factor
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -533,9 +546,9 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom :
             zmxlm(ji,jj,jk) =   &
                &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3DS_00_00( jpkm1, 2, -1 )
+         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface :
             zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -544,13 +557,13 @@
          !
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup
             zmxld(ji,jj,jk) =    &
                &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3DS_00_00( jpkm1, 2, -1 )
+         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
             zmxlm(ji,jj,jk) =   &
                &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3D_00_00( 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
             zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
@@ -564,5 +577,5 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
          zsqen = SQRT( en(ji,jj,jk) )
          zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
@@ -573,6 +586,6 @@
       !
       !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
-         DO_3D_00_00( 2, jpkm1 )
+      IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
          END_3D
@@ -610,6 +623,6 @@
          &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
          &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
-         &                 nn_pdl  , ln_drg   , ln_lc    , rn_lc,      &
-         &                 nn_etau , nn_htau  , rn_efr   , rn_eice  
+         &                 nn_pdl  , ln_lc    , rn_lc    ,             &
+         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  
       !!----------------------------------------------------------------------
       !
@@ -637,11 +650,18 @@
          WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
          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
-         ENDIF         
-         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
-         WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg
+            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
@@ -649,12 +669,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   '
-         IF( ln_drg ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
-            WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top
-            WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot
-         ENDIF
+         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      
          WRITE(numout,*)
          WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
@@ -737,8 +758,8 @@
             !
             IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
-               CALL iom_get( numror, jpdom_autoglo, 'en'   , en   , ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'en'   , en   , ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'dissl', dissl, ldxios = lrxios )
             ELSE                                          ! start TKE from rest
                IF(lwp) WRITE(numout,*)
