Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdf_oce.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdf_oce.F90	(revision 14012)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdf_oce.F90	(revision 14016)
@@ -40,4 +40,5 @@
    LOGICAL , PUBLIC ::   ln_zdfswm   !: surface  wave-induced mixing flag
    LOGICAL , PUBLIC ::   ln_zdfiwm   !: internal wave-induced mixing flag
+   LOGICAL , PUBLIC ::   ln_zdfmfc   !: convection: eddy diffusivity Mass Flux Convection
    !                             ! coefficients 
    REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s)
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfmfc.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfmfc.F90	(revision 14016)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfmfc.F90	(revision 14016)
@@ -0,0 +1,475 @@
+ MODULE zdfmfc
+   !!======================================================================
+   !!                       ***  MODULE  zdfmfc  ***
+   !! Ocean physics: Mass-Flux scheme parameterization of Convection:
+   !!                Non-local transport for the convective ocean boundary
+   !!                layer. Subgrid-scale large eddies are represented by a
+   !!                mass-flux contribution (ln_zdfmfc = .TRUE.)
+   !!======================================================================
+   !! History : NEMO  !
+   !!            3.6  !  2016-06  (H. Giordani, R. Bourdallé-Badie)  Original code
+   !!            4.2  !  2020-12  (H. Giordani, R. Bourdallé-Badie)  adapt to NEM04.2
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   tra_mfc       : Compute the Mass Flux and trends of T/S
+   !!   diag_mfc      : Modify diagonal of trazdf Matrix
+   !!   rhs_mfc       : Modify RHS of trazdf Matrix
+   !!   zdf_mfc_init  : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
+   !
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE domvvl         ! ocean space and time domain : variable volume layer
+   USE domzgr
+   USE zdf_oce        ! ocean vertical physics
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state (eos routine)
+   USE zdfmxl         ! mixed layer
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! MPP manager
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_mfc         ! routine called in step module
+   PUBLIC   diag_mfc        ! routine called in trazdf module
+   PUBLIC   rhs_mfc         ! routine called in trazdf module
+   PUBLIC   zdf_mfc_init    ! routine called in nemo module
+   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  edmfa, edmfb, edmfc   !: diagonal term of the matrix.
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  edmftra               !: y term for matrix inversion
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  edmfm               !: y term for matrix inversion
+   !
+   !! ** Namelist  namzdf_edmf  **
+   REAL(wp) ::   rn_cemf           ! entrain of T/S
+   REAL(wp) ::   rn_cwmf           ! detrain of T/S
+   REAL(wp) ::   rn_cent           ! entrain of the convective mass flux
+   REAL(wp) ::   rn_cdet           ! detrain of the convective mass flux
+   REAL(wp) ::   rn_cap            ! Factor of computation for convective area (negative => area constant)
+   REAL(wp) ::   App_max           ! Maximum of the convective area
+   LOGICAL, PUBLIC, SAVE  ::   ln_edmfuv         !: EDMF flag for velocity  !
+   !
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.2 , NEMO Consortium (2018)
+   !! $Id: zdfmfc.F90 13783 2020-20-02 15:30:22Z rbourdal $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_mfc_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_edmf_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( edmfa(jpi,jpj,jpk) , edmfb(jpi,jpj,jpk) , edmfc(jpi,jpj,jpk)      &
+         &      , edmftra(jpi,jpj,jpk,2), edmfm(jpi,jpj,jpk) ,  STAT= zdf_mfc_alloc )
+         !
+      IF( lk_mpp             )   CALL mpp_sum ( 'zdfmfc', zdf_mfc_alloc )
+      IF( zdf_mfc_alloc /= 0 )   CALL ctl_warn('zdf_mfc_alloc: failed to allocate arrays')
+   END FUNCTION zdf_mfc_alloc
+
+
+   SUBROUTINE tra_mfc( kt, Kmm, pts, Krhs )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_mfc  ***
+      !!
+      !! ** Purpose :      Compute a mass flux, depending on surface flux, over
+      !!            the instable part of the water column.
+      !!
+      !! ** Method  :     Compute surface instability and mix tracers until stable level
+      !!           
+      !!
+      !! ** Action  :      Compute convection plume and (ta,sa)-trends for trazdf (EDMF scheme)
+      !!
+      !! References :      
+      !!                   Giordani, Bourdallé-Badie and Madec JAMES 2020
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in)    :: Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztsp         ! T/S of the plume
+      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztse         ! T/S at W point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp          !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp2         !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zapp          !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zedmf         !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zepsT, zepsW  !
+      !
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar, zustar2   !
+      REAL(wp), DIMENSION(jpi,jpj) :: zuws, zvws, zsws, zfnet          !
+      REAL(wp), DIMENSION(jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl
+      REAL(wp), DIMENSION(jpi,jpj) :: zwpsurf            !
+      REAL(wp), DIMENSION(jpi,jpj) :: zop0 , zsp0 !
+      REAL(wp), DIMENSION(jpi,jpj) :: zrwp_0, zrwp2_0  !
+      REAL(wp), DIMENSION(jpi,jpj) :: zapp0           !
+      REAL(wp), DIMENSION(jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro
+      REAL(wp), DIMENSION(jpi,jpj) :: zhcmo          !
+      !
+      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zn2    ! N^2
+      REAL(wp), DIMENSION(jpi,jpj,2  ) ::   zab, zabm1, zabp ! alpha and beta
+     
+      REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value
+
+      REAL(wp) :: zrho, zrhop
+      REAL(wp) :: zcnh, znum, zden, zcoef1, zcoef2
+      REAL(wp) :: zca, zcb, zcd, zrw, zxl, zcdet, zctre
+      REAL(wp) :: zaw, zbw, zxw
+      REAL(wp) :: alpha
+     !
+      INTEGER, INTENT(in   )    ::   kt   ! ocean time-step index      !
+      !
+      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
+      !
+      !------------------------------------------------------------------
+      ! Initialisation of coefficients
+      !------------------------------------------------------------------
+      zca          = 1._wp
+      zcb          = 1._wp
+      zcd          = 1._wp
+
+      !------------------------------------------------------------------
+      ! Surface boundary condition
+      !------------------------------------------------------------------
+      ! surface Stress
+      !--------------------
+      zuws(:,:) = utau(:,:) * r1_rho0 
+      zvws(:,:) = vtau(:,:) * r1_rho0 
+      zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:))
+      zustar(:,:)  = SQRT(zustar2(:,:))
+
+      ! Heat Flux
+      !--------------------
+      zfnet(:,:) = qns(:,:) + qsr(:,:)
+      zfnet(:,:) = zfnet(:,:) / (rho0 * rcp)
+
+      ! Water Flux
+      !---------------------
+      zsws(:,:) = emp(:,:)
+
+      !-------------------------------------------
+      ! Initialisation of prognostic variables
+      !-------------------------------------------
+      zrwp (:,:,:) =  0._wp ; zrwp2(:,:,:) =  0._wp ; zedmf(:,:,:) =  0._wp
+      zph  (:,:)   =  0._wp ; zphm1(:,:)   =  0._wp ; zphpm1(:,:)  =  0._wp
+      ztsp(:,:,:,:)=  0._wp
+
+      ! Tracers inside plume (ztsp) and environment (ztse)
+      ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
+      ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
+      ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
+      ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
+
+      CALL eos( ztse(:,:,1,:) ,  zrautb(:,:) )
+      CALL eos( ztsp(:,:,1,:) ,  zraupl(:,:) )
+
+      !-------------------------------------------
+      ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain)
+      !-------------------------------------------
+      zhcmo(:,:) = e3t(:,:,1,Kmm)
+      zfbuo(:,:)   = 0._wp
+      WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:)   =   &
+         &      grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:)
+
+      zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:))
+      zedmf(:,:,1) = MAX(0., zedmf(:,:,1))
+
+      zwpsurf(:,:) = 2._wp/3._wp*zustar(:,:) + 2._wp/3._wp*ABS(zfbuo(:,:))**(1._wp/3._wp)
+      zwpsurf(:,:) = MAX(1.e-5_wp,zwpsurf(:,:))
+      zwpsurf(:,:) = MIN(1.,zwpsurf(:,:))
+
+      zapp(:,:,:)  = App_max
+      WHERE(zwpsurf .NE. 0.) zapp(:,:,1)   = MIN(MAX(0.,zedmf(:,:,1)/zwpsurf(:,:)), App_max)
+
+      zedmf(:,:,1) = 0._wp 
+      zrwp (:,:,1) = 0._wp 
+      zrwp2(:,:,1) = 0._wp
+      zepsT(:,:,:) = 0.001_wp
+      zepsW(:,:,:) = 0.001_wp
+
+
+      !--------------------------------------------------------------
+      ! Compute plume properties 
+      ! In the same loop on vert. levels computation of:
+      !    - Vertical velocity: zWp
+      !    - Convective Area: zAp
+      !    - Tracers properties inside the plume (if necessary): ztp
+      !---------------------------------------------------------------
+
+      DO jk= 2, jpk
+
+         ! Compute the buoyancy acceleration on T-points at jk-1
+         zrautbm1(:,:) = zrautb(:,:)
+         CALL eos( pts (:,:,jk  ,:,Kmm) ,  zrautb(:,:)   )
+         CALL eos( ztsp(:,:,jk-1,:    ) ,  zraupl(:,:)   )
+
+         zphm1(:,:)  = zphm1(:,:)  + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm)
+         zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:)   * e3t(:,:,jk-1, Kmm)
+         zph(:,:)    = zphm1(:,:)  + grav * zrautb(:,:)   * e3t(:,:,jk  , Kmm)
+         zph(:,:)    = MAX( zph(:,:), zepsilon)
+
+         WHERE(zrautbm1 .NE. 0.) zfbuo(:,:)  =  grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:)
+
+         DO_2D( 0, 0, 0, 0 )
+
+            ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points
+            zrw              =  (gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm)) &
+               &              / (gdept(ji,jj,jk,Kmm) - gdept(ji,jj,jk-1,Kmm))
+            ztse(ji,jj,jk,:) = (pts(ji,jj,jk,:,Kmm) * zrw + pts(ji,jj,jk-1,:,Kmm)*(1._wp - zrw) )*tmask(ji,jj,jk)
+
+            !---------------------------------------------------------------
+            ! Compute the vertical velocity on W-points
+            !---------------------------------------------------------------
+
+            ! Non-hydrostatic pressure terms in the wp2 equation
+            zcnh = 0.2_wp 
+            znum = 0.5_wp  + zcnh - &
+                   (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) &
+                   *e3t(ji,jj,jk-1,Kmm)*0.5_wp   
+            zden = 0.5_wp + zcnh + &
+                   (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) &
+                   *e3t(ji,jj,jk-1,Kmm)*0.5_wp   
+
+            zcoef1 = zca*e3t(ji,jj,jk-1,Kmm) / zden
+            zcoef2 = znum/zden
+
+            ! compute wp2 
+            zrwp2(ji,jj,jk) = zcoef1*zfbuo(ji,jj) &
+                            + zcoef2*zrwp2(ji,jj,jk-1)
+            zrwp2(ji,jj,jk) = MAX ( zrwp2(ji,jj,jk)*wmask(ji,jj,jk) , 0.)
+            zrwp (ji,jj,jk) = SQRT( zrwp2(ji,jj,jk) )
+
+            !----------------------------------------------------------------------------------
+            ! Compute convective area on W-point
+            ! Compute vertical profil of the convective area with mass conservation hypothesis
+            ! If rn_cap negative => constant value on the water column.
+            !----------------------------------------------------------------------------------
+            IF( rn_cap .GT. 0. ) THEN
+
+               zxw = MAX(zrwp(ji,jj,jk-1), zrwp(ji,jj,jk) )
+               IF( zxw > 0. ) THEN
+
+                  zxl = (zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk))/(e3t(ji,jj,jk-1,Kmm)*zxw)
+                  IF (zxl .LT. 0._wp) THEN
+                     zctre  = -1.*rn_cap*zxl 
+                     zcdet  =  0._wp
+                  ELSE
+                     zctre  =  0._wp
+                     zcdet  =  rn_cap*zxl 
+                  END IF
+                     zapp(ji,jj,jk) = zapp(ji,jj,jk-1)*     &
+                     &                (1._wp + (zxl + zctre - zcdet )*e3t(ji,jj,jk-1,Kmm))
+                  ELSE
+                     zapp(ji,jj,jk) = App_max
+                  END IF
+                  zapp(ji,jj,jk) = MIN( MAX(zapp(ji,jj,jk),0.), App_max)
+               ELSE
+                  zapp(ji,jj,jk) = -1. * rn_cap
+               END IF
+
+            ! Compute Mass Flux on W-point
+            zedmf(ji,jj,jk)   = -zapp(ji,jj,jk) * zrwp(ji,jj,jk)* wmask(ji,jj,jk)
+
+            ! Compute Entrainment coefficient
+            IF(rn_cemf .GT. 0.) THEN
+               zxw = 0.5_wp*(zrwp(ji,jj,jk-1)+ zrwp(ji,jj,jk) )
+               zepsT(ji,jj,jk)  =  0.01_wp
+               IF( zxw > 0.  ) THEN
+                  zepsT(ji,jj,jk)  =  zepsT(ji,jj,jk) +                       &
+                                   &  ABS( zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk) )  &
+                                   &  / ( e3t(ji,jj,jk-1,Kmm) * zxw )
+                  zepsT(ji,jj,jk)  = zepsT(ji,jj,jk) * rn_cemf * wmask(ji,jj,jk)
+               ENDIF
+            ELSE
+               zepsT(ji,jj,jk)  = -rn_cemf
+            ENDIF
+
+            ! Compute the detrend coef for velocity (on W-point and not T-points, bug ???)
+            IF(rn_cwmf .GT. 0.) THEN
+               zepsW(ji,jj,jk)  =  rn_cwmf * zepsT(ji,jj,jk)
+            ELSE
+               zepsW(ji,jj,jk)  = -rn_cwmf
+            ENDIF
+
+            !---------------------------------------------------------------
+            ! Compute the plume properties on T-points
+            !---------------------------------------------------------------
+            IF(zrwp (ji,jj,jk) .LT. 1.e-12_wp .AND. zrwp (ji,jj,jk-1) .LT. 1.e-12_wp) THEN
+               ztsp(ji,jj,jk-1,jp_tem) = pts(ji,jj,jk-1,jp_tem,Kmm)
+               ztsp(ji,jj,jk-1,jp_sal) = pts(ji,jj,jk-1,jp_sal,Kmm)
+            ENDIF
+
+            zcoef1 =  (1._wp-zepsT(ji,jj,jk)*(1._wp-zrw)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk ) ) &
+            &       / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) )
+            !
+            zcoef2 =  zepsT(ji,jj,jk)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk)                       &
+            &       / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk))
+            !
+            ztsp(ji,jj,jk,jp_tem) = (zcoef1 * ztsp(ji,jj,jk-1,jp_tem) +  &
+            &                        zcoef2 * ztse(ji,jj,jk  ,jp_tem) )*tmask(ji,jj,jk)
+            ztsp(ji,jj,jk,jp_sal) = (zcoef1 * ztsp(ji,jj,jk-1,jp_sal) +  &
+            &                        zcoef2 * ztse(ji,jj,jk  ,jp_sal) )*tmask(ji,jj,jk)
+
+         END_2D 
+      END DO ! end of loop on jpk
+
+      ! Compute Mass Flux on T-point
+      DO jk=1,jpk-1
+         edmfm(:,:,jk) = (zedmf(:,:,jk+1)  + zedmf(:,:,jk) )*0.5_wp
+      END DO
+      edmfm(:,:,jpk) = zedmf(:,:,jpk) 
+
+      ! Save variable (on T point)
+      CALL iom_put( "mf_Tp" , ztsp(:,:,:,jp_tem) )  ! Save plume temperature
+      CALL iom_put( "mf_Sp" , ztsp(:,:,:,jp_sal) )  ! Save plume salinity
+      CALL iom_put( "mf_mf" , edmfm(:,:,:)       )  ! Save Mass Flux
+      ! Save variable (on W point)
+      CALL iom_put( "mf_wp" , zrwp (:,:,:)       )  ! Save convective velocity in the plume
+      CALL iom_put( "mf_app", zapp (:,:,:)       )  ! Save convective area
+
+      !=================================================================================
+      !  Computation of a tridiagonal matrix and right hand side terms of the linear system
+      !=================================================================================
+      edmfa(:,:,:)     = 0._wp
+      edmfb(:,:,:)     = 0._wp
+      edmfc(:,:,:)     = 0._wp
+      edmftra(:,:,:,:) = 0._wp
+
+      !---------------------------------------------------------------
+      ! Diagonal terms 
+      !---------------------------------------------------------------
+      DO jk=1,jpk-1
+         edmfa(:,:,jk) =  0._wp
+         edmfb(:,:,jk) = -edmfm(:,:,jk  ) / e3w(:,:,jk+1,Kmm)
+         edmfc(:,:,jk) =  edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm)
+      END DO
+      edmfa(:,:,jpk)   = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm)
+      edmfb(:,:,jpk)   =  edmfm(:,:,jpk  ) / e3w(:,:,jpk,Kmm)
+      edmfc(:,:,jpk)   =  0._wp
+
+      !---------------------------------------------------------------
+      ! right hand side term for Temperature
+      !---------------------------------------------------------------
+      DO jk=1,jpk-1
+        edmftra(:,:,jk,1) = - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_tem) / e3w(:,:,jk+1,Kmm) &
+                          & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm)
+      END DO
+      edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) &
+                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_tem) / e3w(:,:,jpk,Kmm)
+                        
+      !---------------------------------------------------------------
+      ! Right hand side term for Salinity
+      !---------------------------------------------------------------
+      DO jk=1,jpk-1
+         edmftra(:,:,jk,2) =  - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_sal) / e3w(:,:,jk+1,Kmm) &
+                           &  + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm)
+      END DO
+      edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) &
+                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_sal) / e3w(:,:,jpk,Kmm)
+      !
+      !
+      CALL lbc_lnk_multi( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.)
+      !
+   END SUBROUTINE tra_mfc
+
+   
+   SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  zdiagi, zdiagd, zdiags  ! inout: tridaig. terms 
+      REAL(wp)                        , INTENT(in   ) ::   p2dt                   ! tracer time-step
+      INTEGER                         , INTENT(in   ) ::   Kaa                    ! ocean time level indices
+
+      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
+
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
+            zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
+            zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
+         END_3D
+
+   END SUBROUTINE diag_mfc
+
+   SUBROUTINE rhs_mfc( zrhs, jjn )
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   zrhs                   ! inout: rhs trend 
+      INTEGER                         , INTENT(in   ) ::   jjn                    ! tracer indices
+
+      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
+
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         zrhs(ji,jj,jk) = zrhs(ji,jj,jk) + edmftra(ji,jj,jk,jjn)
+      END_3D
+
+   END SUBROUTINE rhs_mfc
+
+
+
+   SUBROUTINE zdf_mfc_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mfc_init  ***
+      !!                    
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
+      !!      mass flux
+      !!
+      !! ** Method  :   Read the namzdf_mfc namelist and check the parameters
+      !!      called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namlist namzdf_mfc
+      !!
+      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
+      !!
+      !!----------------------------------------------------------------------
+      !
+      INTEGER ::   jk    ! dummy loop indices
+      INTEGER ::   ios   ! Local integer output status for namelist read
+      REAL(wp)::   zcr   ! local scalar
+      !!
+      NAMELIST/namzdf_mfc/ ln_edmfuv, rn_cemf, rn_cwmf, rn_cent, rn_cdet, rn_cap, App_max
+      !!----------------------------------------------------------
+      !
+      !
+!      REWIND( numnam_ref )              ! Namelist namzdf_mfc in reference namelist : Vertical eddy diffivity mass flux
+      READ  ( numnam_ref, namzdf_mfc, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in reference namelist' )
+
+!      REWIND( numnam_cfg )              ! Namelist namzdf_mfc in configuration namelist : Vertical eddy diffivity mass flux
+      READ  ( numnam_cfg, namzdf_mfc, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in configuration namelist' )
+      IF(lwm) WRITE ( numond, namzdf_mfc )
+
+      IF(lwp) THEN                     !* Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_mfc_init'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_mfc : set eddy diffusivity Mass Flux Convection'
+         WRITE(numout,*) '   Apply mass flux on velocities (Not yet avail.)     ln_edmfuv = ', ln_edmfuv
+         WRITE(numout,*) '   Coeff for entrain/detrain T/S of plume (Neg => cte) rn_cemf  = ', rn_cemf
+         WRITE(numout,*) '   Coeff for entrain/detrain Wp of plume  (Neg => cte) rn_cwmf  = ', rn_cwmf
+         WRITE(numout,*) '   Coeff for entrain/detrain area of plume             rn_cap   = ', rn_cap
+         WRITE(numout,*) '   Coeff for entrain area of plume                     rn_cent  = ', rn_cent
+         WRITE(numout,*) '   Coeff for detrain area of plume                     rn_cdet  = ', rn_cdet
+         WRITE(numout,*) '   Max convective area                                 App_max  = ', App_max
+       ENDIF
+                                     !* allocate edmf arrays
+      IF( zdf_mfc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_edmf_init : unable to allocate arrays' )
+      edmfa(:,:,:)     = 0._wp
+      edmfb(:,:,:)     = 0._wp
+      edmfc(:,:,:)     = 0._wp
+      edmftra(:,:,:,:) = 0._wp
+      !
+   END SUBROUTINE zdf_mfc_init
+
+   !!======================================================================
+ 
+   !!======================================================================
+END MODULE zdfmfc
+
+
+
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfphy.F90	(revision 14012)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfphy.F90	(revision 14016)
@@ -21,4 +21,5 @@
    USE zdfddm         ! vertical physics: double diffusion mixing      
    USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion  
+   USE zdfmfc         ! vertical physics: Mass Flux Convection 
    USE zdfiwm         ! vertical physics: internal wave-induced mixing  
    USE zdfswm         ! vertical physics: surface  wave-induced mixing
@@ -78,4 +79,5 @@
       NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
          &             ln_zdfosm,                                    &     ! type of closure scheme
+         &             ln_zdfmfc,                                    &     ! convection : mass flux
          &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
          &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
@@ -112,4 +114,5 @@
          WRITE(numout,*) '         OSMOSIS-OBL closure (OSM)               ln_zdfosm = ', ln_zdfosm
          WRITE(numout,*) '      convection: '
+         WRITE(numout,*) '         convection mass flux (mfc)              ln_zdfmfc = ', ln_zdfmfc
          WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
          WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
@@ -172,10 +175,15 @@
       IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
       IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' )
+      IF( ln_zdfmfc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfevd' )
+      IF( ln_zdfmfc .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfnpc' )
+      IF( ln_zdfmfc .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfosm' )
       IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
       IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )
+      IF( lk_top    .AND. ln_zdfmfc )   CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' )
       IF(lwp) THEN
          WRITE(numout,*)
          IF    ( ln_zdfnpc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use non penetrative convective scheme'
          ELSEIF( ln_zdfevd ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use enhanced vertical diffusion scheme'
+         ELSEIF( ln_zdfmfc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use Mass Flux scheme'
          ELSE                      ;   WRITE(numout,*) '   ==>>>   convection: no specific scheme used'
          ENDIF
@@ -205,5 +213,7 @@
       ELSE                   ;   l_zdfsh2 = .TRUE.
       ENDIF
-
+      !                          !== Mass Flux Convectiive algorithm  ==!
+      IF( ln_zdfmfc )   CALL zdf_mfc_init       ! Convection computed with eddy diffusivity mass flux
+      !
       !                          !== gravity wave-driven mixing  ==!
       IF( ln_zdfiwm )   CALL zdf_iwm_init       ! internal wave-driven mixing
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfsh2.F90	(revision 14012)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfsh2.F90	(revision 14016)
@@ -6,4 +6,6 @@
    !! History :   -   !  2014-10  (A. Barthelemy, G. Madec)  original code
    !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm
+   !!   NEMO     4.2  !  2020-12  (G. Madec, E. Clementi) add Stokes Drift Shear
+   !                  !           for wave coupling
    !!----------------------------------------------------------------------
 
@@ -13,4 +15,6 @@
    USE oce
    USE dom_oce        ! domain: ocean
+   USE sbcwave        ! Surface Waves (add Stokes shear)
+   USE sbc_oce , ONLY: ln_stshear  !Stoked Drift shear contribution
    !
    USE in_out_manager ! I/O manager
@@ -21,5 +25,5 @@
 
    PUBLIC   zdf_sh2        ! called by zdftke, zdfglf, and zdfric
-   
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -59,22 +63,37 @@
       !!--------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1
-         DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
-            zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
-               &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
-               &         * (   uu(ji,jj,jk-1,Kbb) -   uu(ji,jj,jk,Kbb) ) & 
-               &         / ( e3uw(ji,jj,jk  ,Kmm) * e3uw(ji,jj,jk,Kbb) ) &
-               &         * wumask(ji,jj,jk)
-            zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
-               &         * (   vv(ji,jj,jk-1,Kmm) -   vv(ji,jj,jk,Kmm) ) &
-               &         * (   vv(ji,jj,jk-1,Kbb) -   vv(ji,jj,jk,Kbb) ) &
-               &         / ( e3vw(ji,jj,jk  ,Kmm) * e3vw(ji,jj,jk,Kbb) ) &
-               &         * wvmask(ji,jj,jk)
-         END_2D
+      DO jk = 2, jpkm1                 !* Shear production at uw- and vw-points (energy conserving form)
+         IF ( cpl_sdrftx .AND. ln_stshear )  THEN       ! Surface Stokes Drift available  ===>>>  shear + stokes drift contibution
+            DO_2D( 1, 0, 1, 0 )
+               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) )        &
+                  &         * ( uu (ji,jj,jk-1,Kmm) -   uu (ji,jj,jk,Kmm)    &
+                  &           + usd(ji,jj,jk-1) -   usd(ji,jj,jk) )  &
+                  &         * ( uu (ji,jj,jk-1,Kbb) -   uu (ji,jj,jk,Kbb) )  &
+                  &         / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk)
+               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) )         &
+                  &         * ( vv (ji,jj,jk-1,Kmm) -   vv (ji,jj,jk,Kmm)     &
+                  &           + vsd(ji,jj,jk-1) -   vsd(ji,jj,jk) )   &
+                  &         * ( vv (ji,jj,jk-1,Kbb) -   vv (ji,jj,jk,Kbb) )   &
+                  &/ ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk)
+            END_2D
+         ELSE
+            DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
+               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
+                  &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
+                  &         * (   uu(ji,jj,jk-1,Kbb) -   uu(ji,jj,jk,Kbb) ) & 
+                  &         / ( e3uw(ji,jj,jk  ,Kmm) * e3uw(ji,jj,jk,Kbb) ) &
+                  &         * wumask(ji,jj,jk)
+               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
+                  &         * (   vv(ji,jj,jk-1,Kmm) -   vv(ji,jj,jk,Kmm) ) &
+                  &         * (   vv(ji,jj,jk-1,Kbb) -   vv(ji,jj,jk,Kbb) ) &
+                  &         / ( e3vw(ji,jj,jk  ,Kmm) * e3vw(ji,jj,jk,Kbb) ) &
+                  &         * wvmask(ji,jj,jk)
+            END_2D
+         ENDIF
          DO_2D( 0, 0, 0, 0 )     !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
             p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
                &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
          END_2D
-      END DO 
+      END DO
       !
    END SUBROUTINE zdf_sh2
Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdftke.F90	(revision 14012)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdftke.F90	(revision 14016)
@@ -29,4 +29,6 @@
    !!            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
+   !!            4.2  !  2020-12  (G. Madec, E. Clementi) add wave coupling
+   !                  !           following Couvelard et al., 2019
    !!----------------------------------------------------------------------
 
@@ -58,4 +60,5 @@
    USE prtctl         ! Print control
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE sbcwave        ! Surface boundary waves
 
    IMPLICIT NONE
@@ -68,4 +71,5 @@
    !                      !!** Namelist  namzdf_tke  **
    LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
+   LOGICAL  ::   ln_mxhsw  ! mixing length scale surface value as a fonction of wave height
    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
@@ -81,4 +85,6 @@
    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)
+   INTEGER  ::   nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
+   INTEGER  ::   nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
    REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
@@ -209,6 +215,7 @@
       REAL(wp) ::   zus   , zwlc  , zind       !   -      -
       REAL(wp) ::   zzd_up, zzd_lw             !   -      -
+      REAL(wp) ::   ztaui, ztauj, z1_norm
       INTEGER , DIMENSION(jpi,jpj)     ::   imlc
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3, zWlc2
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
@@ -219,4 +226,6 @@
       zfact2  = 1.5_wp * rn_Dt * rn_ediss
       zfact3  = 0.5_wp         * rn_ediss
+      !
+      zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used
       !
       ! ice fraction considered for attenuation of langmuir & wave breaking
@@ -232,10 +241,9 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      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)
+      DO_2D( 0, 0, 0, 0 )
          en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+         zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
+         zd_lw(ji,jj,1) = 1._wp  
+         zd_up(ji,jj,1) = 0._wp
       END_2D
       !
@@ -274,8 +282,45 @@
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke (Axell JGR 2002)
          !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
          !
-         !                        !* total energy produce by LC : cumulative sum over jk
+         !                       !* Langmuir velocity scale
+         !
+         IF ( cpl_sdrftx )  THEN       ! Surface Stokes Drift available
+            !                                ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2    with u* = (taum/rho0)^1/2
+            !                                ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s
+            !                                ! more precisely, it is the dot product that must be used :
+            !                                !     1/2  (W_lc)^2 = MAX( u* u_s + v* v_s , 0 )   only the positive part
+!!gm  ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
+!!gm  ! so we will overestimate the LC velocity....   !!gm I will do the work if !LC have an effect !
+            DO_2D( 0, 0, 0, 0 )
+!!XC                  zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )  )
+                  zWlc2(ji,jj) = 0.5_wp *  ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )
+            END_2D
+!
+!  Projection of Stokes drift in the wind stress direction
+!
+            DO_2D( 0, 0, 0, 0 )
+                  ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
+                  ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
+                  z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1)
+                  zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
+            END_2D
+         CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. )
+!
+         ELSE                          ! Surface Stokes drift deduced from surface stress
+            !                                ! Wlc = u_s   with u_s = 0.016*U_10m, the surface stokes drift  (Axell 2002, Eq.44)
+            !                                ! using |tau| = rho_air Cd |U_10m|^2 , it comes:
+            !                                ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2   and thus:
+            !                                ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag )
+            zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag
+            DO_2D( 1, 1, 1, 1 )
+               zWlc2(ji,jj) = zcof * taum(ji,jj)
+            END_2D
+            !
+         ENDIF
+         !
+         !                       !* Depth of the LC circulation  (Axell 2002, Eq.47)
+         !                             !- LHS of Eq.47
          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm)
          DO jk = 2, jpk
@@ -283,10 +328,9 @@
                &        MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
          END DO
-         !                        !* finite Langmuir Circulation depth
-         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
+         !
+         !                             !- compare LHS to RHS of Eq.47
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
-         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
+         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
+            IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk
          END_3D
          !                               ! finite LC depth
@@ -294,7 +338,8 @@
             zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
          END_2D
+         !
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
          DO_2D( 0, 0, 0, 0 )
-            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
+            zus = SQRT( 2. * zWlc2(ji,jj) )             ! Stokes drift
             zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
          END_2D
@@ -351,12 +396,43 @@
             &                                ) * wmask(ji,jj,jk)
       END_3D
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface boundary condition on tke if
+      !                     !  coupling with waves
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF ( cpl_phioc .and. ln_phioc )  THEN
+         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves 
+
+         CASE ( 0 ) ! Dirichlet BC
+            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
+               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
+            END_2D
+
+         CASE ( 1 ) ! Neumann BC
+            DO_2D( 0, 0, 0, 0 )
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
+               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
+               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
+               zdiag(ji,jj,1) = 1._wp
+               zd_lw(ji,jj,2) = 0._wp
+            END_2D
+
+         END SELECT
+
+      ENDIF
+      !
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO_3D( 0, 0, 0, 0, 3, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+      DO_3D( 0, 0, 0, 0, 2, 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( 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( 0, 0, 0, 0, 3, jpkm1 )
+!XC : commented to allow for neumann boundary condition
+!      DO_2D( 0, 0, 0, 0 )
+!         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( 0, 0, 0, 0, 2, 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
@@ -460,56 +536,63 @@
       zmxlm(:,:,:)  = rmxl_min    
       zmxld(:,:,:)  = rmxl_min
+      !
+      IF(ln_sdw .AND. ln_mxhsw) THEN
+         zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
+         ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6
+         zcoef       = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff
+         zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
+      ELSE
       ! 
-     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
-         !
-         zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
+         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( 0, 0, 0, 0 )                  ! No sea-ice
-            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
-         END_2D
+            DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
+               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
+            END_2D
 #else
-         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
-         !
-         CASE( 0 )                      ! No scaling under sea-ice
+            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
+            !
+            CASE( 0 )                      ! No scaling under sea-ice
+               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( 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( 0, 0, 0, 0 )
+#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_2D
+               !
+            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._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
+               END_2D
+               !
+            END SELECT
+#endif
+            !
             DO_2D( 0, 0, 0, 0 )
-               zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
+               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
             END_2D
             !
-         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( 0, 0, 0, 0 )
-#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_2D
-            !
-         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._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
-                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
-            END_2D
-            !
-         END SELECT
-#endif
-         !
-         DO_2D( 0, 0, 0, 0 )
-            zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
-         END_2D
-         !
-      ELSE
-         zmxlm(:,:,1) = rn_mxl0
-      ENDIF
-
+         ELSE
+            zmxlm(:,:,1) = rn_mxl0
+         ENDIF
+      ENDIF
       !
       DO_3D( 0, 0, 0, 0, 2, jpkm1 )
@@ -624,5 +707,6 @@
          &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
          &                 nn_pdl  , ln_lc    , rn_lc    ,             &
-         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  
+         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  ,  &   
+         &                 nn_bc_surf, nn_bc_bot, ln_mxhsw
       !!----------------------------------------------------------------------
       !
@@ -666,4 +750,10 @@
          WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
          WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
+         IF ( cpl_phioc .and. ln_phioc )  THEN
+            SELECT CASE( nn_bc_surf)             ! Type of scaling under sea-ice
+            CASE( 0 )   ;   WRITE(numout,*) '  nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
+            CASE( 1 )   ;   WRITE(numout,*) '  nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
+            END SELECT
+         ENDIF
          WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
          WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
