Index: NEMO/trunk/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynspg.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/DYN/dynspg.F90	(revision 14007)
@@ -6,4 +6,6 @@
    !! History :  1.0  ! 2005-12  (C. Talandier, G. Madec, V. Garnier)  Original code
    !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) add Bernoulli Head for
+   !!                            wave coupling
    !!----------------------------------------------------------------------
 
@@ -19,4 +21,5 @@
    USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
    USE sbcapr         ! surface boundary condition: atmospheric pressure
+   USE sbcwave,  ONLY : bhd_wave
    USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
@@ -143,4 +146,11 @@
          ENDIF
          !
+         IF( ln_wave .and. ln_bern_srfc ) THEN          !== Add J terms: depth-independent Bernoulli head
+            DO_2D( 0, 0, 0, 0 )
+               spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
+               spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
+            END_2D
+         ENDIF
+         !
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !== Add all terms to the general trend
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
Index: NEMO/trunk/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynvor.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/DYN/dynvor.F90	(revision 14007)
@@ -21,4 +21,5 @@
    !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
    !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T)
    !!----------------------------------------------------------------------
 
@@ -37,5 +38,5 @@
    USE trddyn         ! trend manager: dynamics
    USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
-   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
+   USE sbc_oce,  ONLY : ln_stcor, ln_vortex_force   ! use Stoke-Coriolis force
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -121,17 +122,12 @@
          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
          !
-         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend
          ztrdv(:,:,:) = pvv(:,:,:,Krhs)
          SELECT CASE( nvor_scheme )
          CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
-            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
-            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
-            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
-            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
-            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          END SELECT
          ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
@@ -161,21 +157,43 @@
          CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
                              CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
                              CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_ENE )                        !* energy conserving scheme
                              CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_ENS )                        !* enstrophy conserving scheme
                              CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
+
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_MIX )                        !* mixed ene-ens scheme
                              CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
                              CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
-            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor )        CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )        ! add the Stokes-Coriolis trend
+            IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add vortex force
          CASE( np_EEN )                        !* energy and enstrophy conserving scheme
                              CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          END SELECT
          !
Index: NEMO/trunk/src/OCE/DYN/dynzad.F90
===================================================================
--- NEMO/trunk/src/OCE/DYN/dynzad.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/DYN/dynzad.F90	(revision 14007)
@@ -16,4 +16,5 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   USE sbcwave, ONLY: wsd   ! Surface Waves (add vertical Stokes-drift)
    !
    USE in_out_manager ! I/O manager
@@ -79,5 +80,9 @@
       DO jk = 2, jpkm1                ! Vertical momentum advection at level w and u- and v- vertical
          DO_2D( 0, 1, 0, 1 )              ! vertical fluxes
+          IF( ln_vortex_force ) THEN
+            zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
+          ELSE
             zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
+          ENDIF
          END_2D
          DO_2D( 0, 0, 0, 0 )              ! vertical momentum advection at w-point
Index: NEMO/trunk/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/cpl_oasis3.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/cpl_oasis3.F90	(revision 14007)
@@ -66,5 +66,5 @@
    INTEGER                    ::   nsnd         ! total number of fields sent 
    INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
-   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=60   ! Maximum number of coupling fields
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=62   ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields
Index: NEMO/trunk/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbc_oce.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/sbc_oce.F90	(revision 14007)
@@ -12,4 +12,5 @@
    !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk)
    !!            4.0  ! 2019-03  (F. Lemarié, G. Samson) add compatibility with ABL mode
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) modified wave parameters in namelist
    !!----------------------------------------------------------------------
 
@@ -36,4 +37,5 @@
    LOGICAL , PUBLIC ::   ln_blk         !: bulk formulation
    LOGICAL , PUBLIC ::   ln_abl         !: Atmospheric boundary layer model
+   LOGICAL , PUBLIC ::   ln_wave        !: wave in the system (forced or coupled)
 #if defined key_oasis3
    LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used
@@ -56,13 +58,4 @@
    !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
    !                                             !:  = 2 annual global mean of e-p-r set to zero
-   LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model
-   LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model
-   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
-   LOGICAL , PUBLIC ::   ln_tauwoc       !: true if normalized stress from wave is used
-   LOGICAL , PUBLIC ::   ln_tauw        !: true if ocean stress components from wave is used
-   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
-   !
-   INTEGER , PUBLIC ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
-   !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
    !
@@ -71,5 +64,19 @@
    !                                   !!* namsbc_cpl namelist *
    INTEGER , PUBLIC ::   nn_cats_cpl    !: Number of sea ice categories over which the coupling is carried out
-
+   !
+   !                                   !!* namsbc_wave namelist *
+   LOGICAL , PUBLIC ::   ln_sdw         !: =T 3d stokes drift from wave model
+   LOGICAL , PUBLIC ::   ln_stcor       !: =T if Stokes-Coriolis and tracer advection terms are used
+   LOGICAL , PUBLIC ::   ln_cdgw        !: =T neutral drag coefficient from wave model
+   LOGICAL , PUBLIC ::   ln_tauoc       !: =T if normalized stress from wave is used
+   LOGICAL , PUBLIC ::   ln_wave_test   !: =T wave test case (constant Stokes drift)
+   LOGICAL , PUBLIC ::   ln_charn       !: =T Chranock coefficient from wave model
+   LOGICAL , PUBLIC ::   ln_taw         !: =T wind stress corrected by wave intake
+   LOGICAL , PUBLIC ::   ln_phioc       !: =T TKE surface BC from wave model 
+   LOGICAL , PUBLIC ::   ln_bern_srfc   !: Bernoulli head, waves' inuced pressure
+   LOGICAL , PUBLIC ::   ln_breivikFV_2016 !: Breivik 2016 profile
+   LOGICAL , PUBLIC ::   ln_vortex_force !: vortex force activation
+   LOGICAL , PUBLIC ::   ln_stshear     !: Stoked Drift shear contribution in zdftke
+   !
    !!----------------------------------------------------------------------
    !!           switch definition (improve readability)
@@ -81,18 +88,5 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 5        !: Pure ocean-atmosphere Coupled formulation
    INTEGER , PUBLIC, PARAMETER ::   jp_none    = 6        !: for OPA when doing coupling via SAS module
-
-   !!----------------------------------------------------------------------
-   !!           Stokes drift parametrization definition
-   !!----------------------------------------------------------------------
-   INTEGER , PUBLIC, PARAMETER ::   jp_breivik_2014 = 0     !: Breivik  2014: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
-   INTEGER , PUBLIC, PARAMETER ::   jp_li_2017      = 1     !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016)
-   !  with depth averaged profile
-   INTEGER , PUBLIC, PARAMETER ::   jp_peakfr       = 2     !: Li et al 2017: using the peak wave number read from wave model instead
-   !  of the inverse depth scale
-   LOGICAL , PUBLIC            ::   ll_st_bv2014  = .FALSE. !  logical indicator, .true. if Breivik 2014 parameterisation is active.
-   LOGICAL , PUBLIC            ::   ll_st_li2017  = .FALSE. !  logical indicator, .true. if Li 2017 parameterisation is active.
-   LOGICAL , PUBLIC            ::   ll_st_bv_li   = .FALSE. !  logical indicator, .true. if either Breivik or Li parameterisation is active.
-   LOGICAL , PUBLIC            ::   ll_st_peakfr  = .FALSE. !  logical indicator, .true. if using Li 2017 with peak wave number
-
+   !
    !!----------------------------------------------------------------------
    !!           component definition
Index: NEMO/trunk/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcblk.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/sbcblk.F90	(revision 14007)
@@ -314,23 +314,4 @@
          ENDIF
       END DO
-      !
-      IF( ln_wave ) THEN
-         !Activated wave module but neither drag nor stokes drift activated
-         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
-            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
-            !drag coefficient read from wave model definable only with mfs bulk formulae and core
-         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN
-            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
-         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN
-            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
-         ENDIF
-      ELSE
-         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
-            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
-            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
-            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
-            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &
-            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
-      ENDIF
       !
       IF( ln_abl ) THEN       ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient
Index: NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 14007)
@@ -17,4 +17,5 @@
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!            4.2  !  2020-12  (G. Madec, E. Clementi) Charnock coeff from wave model
    !!----------------------------------------------------------------------
 
@@ -31,5 +32,5 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE sbcwave, ONLY   :  cdn_wave ! wave module
+   USE sbcwave, ONLY   : charn ! wave module
 #if defined key_si3 || defined key_cice
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -233,5 +234,10 @@
       u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
 
-      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
+      IF (ln_charn)  THEN          !  Charnock value if wave coupling
+         z0     = charn*u_star*u_star/grav + 0.11_wp*znu_a/u_star
+      ELSE
+         z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
+      ENDIF
+
       z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
 
@@ -302,5 +308,9 @@
          ztmp2  = u_star*u_star
          ztmp1  = znu_a/u_star
-         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
+         IF (ln_charn) THEN     ! Charnock value if wave coupling
+            z0  = MIN( ABS( alpha_M*ztmp1 + charn*ztmp2/grav ) , 0.001_wp)         
+         ELSE
+            z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
+         ENDIF
          z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
          z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
Index: NEMO/trunk/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbccpl.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/sbccpl.F90	(revision 14007)
@@ -8,4 +8,5 @@
    !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
    !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi)  wave coupling updates
    !!----------------------------------------------------------------------
 
@@ -106,20 +107,26 @@
    INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
    INTEGER, PARAMETER ::   jpr_mslp   = 43   ! mean sea level pressure 
-   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig 
-   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux 
-   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1 
-   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2 
+   !**  surface wave coupling  **
+   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig
+   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux
+   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1
+   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2
    INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period
    INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber
-   INTEGER, PARAMETER ::   jpr_tauwoc = 50   ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_wstrf  = 50   ! Stress fraction adsorbed by waves
    INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient
-   INTEGER, PARAMETER ::   jpr_isf    = 52
-   INTEGER, PARAMETER ::   jpr_icb    = 53
-   INTEGER, PARAMETER ::   jpr_wfreq  = 54   ! Wave peak frequency
-   INTEGER, PARAMETER ::   jpr_tauwx  = 55   ! x component of the ocean stress from waves
-   INTEGER, PARAMETER ::   jpr_tauwy  = 56   ! y component of the ocean stress from waves
-   INTEGER, PARAMETER ::   jpr_ts_ice = 57   ! Sea ice surface temp
-
-   INTEGER, PARAMETER ::   jprcv      = 57   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_charn  = 52   ! Chranock coefficient
+   INTEGER, PARAMETER ::   jpr_twox   = 53   ! wave to ocean momentum flux
+   INTEGER, PARAMETER ::   jpr_twoy   = 54   ! wave to ocean momentum flux
+   INTEGER, PARAMETER ::   jpr_tawx   = 55   ! net wave-supported stress
+   INTEGER, PARAMETER ::   jpr_tawy   = 56   ! net wave-supported stress
+   INTEGER, PARAMETER ::   jpr_bhd    = 57   ! Bernoulli head. waves' induced surface pressure
+   INTEGER, PARAMETER ::   jpr_tusd   = 58   ! zonal stokes transport
+   INTEGER, PARAMETER ::   jpr_tvsd   = 59   ! meridional stokes tranmport
+   INTEGER, PARAMETER ::   jpr_isf    = 60
+   INTEGER, PARAMETER ::   jpr_icb    = 61
+   INTEGER, PARAMETER ::   jpr_ts_ice = 62   ! Sea ice surface temp
+
+   INTEGER, PARAMETER ::   jprcv      = 62   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -184,12 +191,12 @@
       &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
    !                                   ! Received from the atmosphere
-   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr,  &
+   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr,  &
       &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
-   ! Send to waves 
+   !                                   ! Send to waves 
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
-   ! Received from waves 
-   TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, sn_rcv_tauwoc, &
-                    sn_rcv_wdrag, sn_rcv_wfreq
+   !                                   ! Received from waves 
+   TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
+      &             sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
    !                                   ! Other namelist parameters
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -274,9 +281,9 @@
          &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  & 
          &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  & 
-         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  &
-         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  &
-         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
-         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
-         &                  sn_rcv_ts_ice
+         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_wstrf ,  &
+         &                  sn_rcv_charn , sn_rcv_taw   , sn_rcv_bhd  , sn_rcv_tusd  , sn_rcv_tvsd,    &
+         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   &
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice 
+
       !!---------------------------------------------------------------------
       !
@@ -319,4 +326,6 @@
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
+         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
+         WRITE(numout,*)'      surface waves:'
          WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')' 
          WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' 
@@ -325,9 +334,7 @@
          WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')' 
          WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')' 
-         WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'
-         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauwoc%cldes), ' (', TRIM(sn_rcv_tauwoc%clcat ), ')' 
-         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
-         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 
+         WRITE(numout,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
@@ -351,5 +358,24 @@
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
       ENDIF
-
+      IF( lwp .AND. ln_wave) THEN                        ! control print
+      WRITE(numout,*)'      surface waves:'
+         WRITE(numout,*)'      Significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'
+         WRITE(numout,*)'      Wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'
+         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'
+         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'
+         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'
+         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
+         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
+         WRITE(numout,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift grid u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift grid v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
+         WRITE(numout,*)'      Bernouilli pressure head        = ', TRIM(sn_rcv_bhd%cldes   ), ' (', TRIM(sn_rcv_bhd%clcat  ), ')'
+         WRITE(numout,*)'Wave to ocean momentum flux and Net wave-supported stress = ', TRIM(sn_rcv_taw%cldes ), ' (', TRIM(sn_rcv_taw%clcat ), ')'
+         WRITE(numout,*)'      Surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'
+         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref
+         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor
+         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd
+      ENDIF
       !                                   ! allocate sbccpl arrays
       IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
@@ -629,9 +655,4 @@
          cpl_wper = .TRUE.
       ENDIF
-      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency 
-      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN
-         srcv(jpr_wfreq)%laction = .TRUE.
-         cpl_wfreq = .TRUE.
-      ENDIF
       srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number
       IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN
@@ -639,15 +660,8 @@
          cpl_wnum = .TRUE.
       ENDIF
-      srcv(jpr_tauwoc)%clname = 'O_TauOce'   ! stress fraction adsorbed by the wave
-      IF( TRIM(sn_rcv_tauwoc%cldes ) == 'coupled' )  THEN
-         srcv(jpr_tauwoc)%laction = .TRUE.
-         cpl_tauwoc = .TRUE.
-      ENDIF
-      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction
-      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction
-      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN
-         srcv(jpr_tauwx)%laction = .TRUE.
-         srcv(jpr_tauwy)%laction = .TRUE.
-         cpl_tauw = .TRUE.
+      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave
+      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wstrf)%laction = .TRUE.
+         cpl_wstrf = .TRUE.
       ENDIF
       srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
@@ -656,7 +670,36 @@
          cpl_wdrag = .TRUE.
       ENDIF
-      IF( srcv(jpr_tauwoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) &
-            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
-                                     '(sn_rcv_tauwoc=coupled and sn_rcv_tauw=coupled)' )
+      srcv(jpr_charn)%clname = 'O_Charn'     ! Chranock coefficient
+      IF( TRIM(sn_rcv_charn%cldes ) == 'coupled' )  THEN
+         srcv(jpr_charn)%laction = .TRUE.
+         cpl_charn = .TRUE.
+      ENDIF
+      srcv(jpr_bhd)%clname = 'O_Bhd'     ! Bernoulli head. waves' induced surface pressure
+      IF( TRIM(sn_rcv_bhd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_bhd)%laction = .TRUE.
+         cpl_bhd = .TRUE.
+      ENDIF
+      srcv(jpr_tusd)%clname = 'O_Tusd'     ! zonal stokes transport
+      IF( TRIM(sn_rcv_tusd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tusd)%laction = .TRUE.
+         cpl_tusd = .TRUE.
+      ENDIF
+      srcv(jpr_tvsd)%clname = 'O_Tvsd'     ! meridional stokes tranmport
+      IF( TRIM(sn_rcv_tvsd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tvsd)%laction = .TRUE.
+         cpl_tvsd = .TRUE.
+      ENDIF
+
+      srcv(jpr_twox)%clname = 'O_Twox'     ! wave to ocean momentum flux in the u direction
+      srcv(jpr_twoy)%clname = 'O_Twoy'     ! wave to ocean momentum flux in the v direction
+      srcv(jpr_tawx)%clname = 'O_Tawx'     ! Net wave-supported stress in the u direction
+      srcv(jpr_tawy)%clname = 'O_Tawy'     ! Net wave-supported stress in the v direction
+      IF( TRIM(sn_rcv_taw%cldes ) == 'coupled' )  THEN
+         srcv(jpr_twox)%laction = .TRUE.
+         srcv(jpr_twoy)%laction = .TRUE.
+         srcv(jpr_tawx)%laction = .TRUE.
+         srcv(jpr_tawy)%laction = .TRUE.
+         cpl_taw = .TRUE.
+      ENDIF
       !
       !                                                      ! ------------------------------- !
@@ -1058,5 +1101,4 @@
       !   initialisation of the coupler  !
       ! ================================ !
-
       CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
       
@@ -1071,4 +1113,5 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
+      !
       !
    END SUBROUTINE sbc_cpl_init
@@ -1146,4 +1189,8 @@
          IF( ncpl_qsr_freq /= 0) ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top
          
+         IF ( ln_wave .AND. nn_components == 0 ) THEN
+            ncpl_qsr_freq = 1;
+            WRITE(numout,*) 'ncpl_qsr_freq is set to 1 when coupling NEMO with wave (without SAS) '
+         ENDIF
       ENDIF
       !
@@ -1320,9 +1367,4 @@
          IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
       ! 
-      !                                                      ! ========================= !  
-      !                                                      !    Wave peak frequency    ! 
-      !                                                      ! ========================= !  
-         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)
-      !
       !                                                      ! ========================= ! 
       !                                                      !    Vertical mixing Qiao   !
@@ -1331,6 +1373,6 @@
 
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
-         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
-                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) THEN
+         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. &
+             srcv(jpr_wper)%laction .OR. srcv(jpr_hsig)%laction )   THEN
             CALL sbc_stokes( Kmm )
          ENDIF
@@ -1339,19 +1381,48 @@
       !                                                      ! Stress adsorbed by waves  !
       !                                                      ! ========================= ! 
-      IF( srcv(jpr_tauwoc)%laction .AND. ln_tauwoc ) tauoc_wave(:,:) = frcv(jpr_tauwoc)%z3(:,:,1)
-
-      !                                                      ! ========================= !  
-      !                                                      ! Stress component by waves ! 
-      !                                                      ! ========================= !  
-      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN
-         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1)
-         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1)
-      ENDIF
-
+      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc )  tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)
+      !
       !                                                      ! ========================= ! 
       !                                                      !   Wave drag coefficient   !
       !                                                      ! ========================= ! 
       IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw )   cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
-
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Chranock coefficient    !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_charn)%laction .AND. ln_charn )  charn(:,:) = frcv(jpr_charn)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      ! net wave-supported stress !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tawx)%laction .AND. ln_taw )     tawx(:,:) = frcv(jpr_tawx)%z3(:,:,1)
+      IF( srcv(jpr_tawy)%laction .AND. ln_taw )     tawy(:,:) = frcv(jpr_tawy)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !wave to ocean momentum flux!
+      !                                                      ! ========================= !
+      IF( srcv(jpr_twox)%laction .AND. ln_taw )     twox(:,:) = frcv(jpr_twox)%z3(:,:,1)
+      IF( srcv(jpr_twoy)%laction .AND. ln_taw )     twoy(:,:) = frcv(jpr_twoy)%z3(:,:,1)
+      !                                                      
+      !                                                      ! ========================= !
+      !                                                      !    wave TKE flux at sfc   !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_phioc)%laction .AND. ln_phioc )     phioc(:,:) = frcv(jpr_phioc)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !      Bernoulli head       !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_bhd)%laction .AND. ln_bern_srfc )   bhd_wave(:,:) = frcv(jpr_bhd)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Stokes transport u dir  !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tusd)%laction .AND. ln_breivikFV_2016 )    tusd(:,:) = frcv(jpr_tusd)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Stokes transport v dir  !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tvsd)%laction .AND. ln_breivikFV_2016 )     tvsd(:,:) = frcv(jpr_tvsd)%z3(:,:,1)
+      !
       !  Fields received by SAS when OASIS coupling
       !  (arrays no more filled at sbcssm stage)
Index: NEMO/trunk/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcmod.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/sbcmod.F90	(revision 14007)
@@ -16,4 +16,5 @@
    !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation
    !!            4.0  ! 2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) modified wave forcing and coupling  
    !!----------------------------------------------------------------------
 
@@ -54,4 +55,5 @@
    USE usrdef_sbc     ! user defined: surface boundary condition
    USE closea         ! closed sea
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    !
    USE prtctl         ! Print control                    (prt_ctl routine)
@@ -70,5 +72,6 @@
 
    INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
-
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -99,7 +102,6 @@
          &             nn_ice   , ln_ice_embd,                                       &
          &             ln_traqsr, ln_dm2dc ,                                         &
-         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_apr_dyn,                  &
-         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc , ln_stcor  ,      &
-         &             ln_tauw  , nn_lsm, nn_sdrift
+         &             ln_rnf   , nn_fwb     , ln_ssr   , ln_apr_dyn,                &
+         &             ln_wave  , nn_lsm
       !!----------------------------------------------------------------------
       !
@@ -133,4 +135,5 @@
          WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk
          WRITE(numout,*) '         ABL          formulation                   ln_abl        = ', ln_abl
+         WRITE(numout,*) '         Surface wave (forced or coupled)           ln_wave       = ', ln_wave
          WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : '
          WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
@@ -150,32 +153,5 @@
          WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
          WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
-         WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
-         WRITE(numout,*) '               Stokes drift corr. to vert. velocity ln_sdw        = ', ln_sdw
-         WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift
-         WRITE(numout,*) '               wave modified ocean stress           ln_tauwoc     = ', ln_tauwoc
-         WRITE(numout,*) '               wave modified ocean stress component ln_tauw       = ', ln_tauw
-         WRITE(numout,*) '               Stokes coriolis term                 ln_stcor      = ', ln_stcor
-         WRITE(numout,*) '               neutral drag coefficient (CORE,NCAR) ln_cdgw       = ', ln_cdgw
-      ENDIF
-      !
-      IF( .NOT.ln_wave ) THEN
-         ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false.
-      ENDIF 
-      IF( ln_sdw ) THEN
-         IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) &
-            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
-      ENDIF
-      ll_st_bv2014  = ( nn_sdrift==jp_breivik_2014 )
-      ll_st_li2017  = ( nn_sdrift==jp_li_2017 )
-      ll_st_bv_li   = ( ll_st_bv2014 .OR. ll_st_li2017 )
-      ll_st_peakfr  = ( nn_sdrift==jp_peakfr )
-      IF( ln_tauwoc .AND. ln_tauw ) &
-         CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
-                                  '(ln_tauwoc=.true. and ln_tauw=.true.)' )
-      IF( ln_tauwoc ) &
-         CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauwoc=.true.)' )
-      IF( ln_tauw ) &
-         CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', &
-                              'This will override any other specification of the ocean stress' )
+      ENDIF
       !
       IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case)
@@ -357,5 +333,22 @@
       IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc, Kbb, Kmm )   ! CICE initialization
       !
-      IF( ln_wave     )   CALL sbc_wave_init                     ! surface wave initialisation
+      IF( ln_wave     ) THEN
+                          CALL sbc_wave_init                     ! surface wave initialisation
+      ELSE
+                          IF(lwp) WRITE(numout,*)
+                          IF(lwp) WRITE(numout,*) '   No surface waves : all wave related logical set to false'
+                          ln_sdw       = .false.
+                          ln_stcor     = .false.
+                          ln_cdgw      = .false.
+                          ln_tauoc     = .false.
+                          ln_wave_test = .false.
+                          ln_charn     = .false.
+                          ln_taw       = .false.
+                          ln_phioc     = .false.
+                          ln_bern_srfc = .false.
+                          ln_breivikFV_2016 = .false.
+                          ln_vortex_force = .false.
+                          ln_stshear  = .false.
+      ENDIF
       !
    END SUBROUTINE sbc_init
@@ -380,4 +373,5 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      INTEGER  ::   jj, ji          ! dummy loop argument
       !
       LOGICAL ::   ll_sas, ll_opa   ! local logical
@@ -412,9 +406,8 @@
       !
       IF( .NOT.ll_sas )   CALL sbc_ssm ( kt, Kbb, Kmm )  ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
-      IF( ln_wave     )   CALL sbc_wave( kt, Kmm )       ! surface waves
-
       !
       !                                            !==  sbc formulation  ==!
       !                                                   
+      !
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
       !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
@@ -423,4 +416,9 @@
       CASE( jp_blk     )
          IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OPA-SAS coupling: SAS receiving fields from OPA
+!!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+         IF( ln_wave )   THEN
+             IF ( lk_oasis )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OPA-wave coupling
+             CALL sbc_wave ( kt, Kmm )
+         ENDIF
                                CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean
                                !
@@ -436,5 +434,33 @@
       IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing
       !
-      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves 
+      IF( ln_wave .AND. ln_tauoc )  THEN            ! Wave stress reduction
+         DO_2D( 0, 0, 0, 0)
+            utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji-1,jj) ) * 0.5_wp
+            vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj-1) ) * 0.5_wp
+         END_2D
+         !
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         !
+         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+         !
+         IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   &
+            &                                'If not requested select ln_tauoc=.false.' )
+         !
+      ELSEIF( ln_wave .AND. ln_taw ) THEN                  ! Wave stress reduction
+         utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:)
+         vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:)
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         !
+         DO_2D( 0, 0, 0, 0)
+             taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
+         END_2D
+         !
+         IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   &
+            &                                'If not requested select ln_taw=.false.' )
+         !
+      ENDIF
+      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. )
       !
       !                                            !==  Misc. Options  ==!
Index: NEMO/trunk/src/OCE/SBC/sbcwave.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcwave.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/SBC/sbcwave.F90	(revision 14007)
@@ -9,15 +9,17 @@
    !!             -   !  2016-12  (G. Madec, E. Clementi) update Stoke drift computation
    !!                                                    + add sbc_wave_ini routine
+   !!            4.2  !  2020-12  (G. Madec, E. Clementi) updates, new Stoke drift computation 
+   !!                                                    according to Couvelard et al.,2019
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   sbc_stokes    : calculate 3D Stokes-drift velocities
-   !!   sbc_wave      : wave data from wave model in netcdf files 
+   !!   sbc_wave      : wave data from wave model: forced (netcdf files) or coupled mode
    !!   sbc_wave_init : initialisation fo surface waves 
    !!----------------------------------------------------------------------
-   USE phycst         ! physical constants 
+   USE phycst         ! physical constants
    USE oce            ! ocean variables
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
-   USE zdf_oce,  ONLY : ln_zdfswm
+   USE dom_oce        ! ocean domain variables
+   USE sbc_oce        ! Surface boundary condition: ocean fields
    USE bdy_oce        ! open boundary condition variables
    USE domvvl         ! domain: variable volume layers
@@ -26,5 +28,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing library
-   USE fldread	       ! read input fields
+   USE fldread        ! read input fields
 
    IMPLICIT NONE
@@ -32,19 +34,21 @@
 
    PUBLIC   sbc_stokes      ! routine called in sbccpl
-   PUBLIC   sbc_wstress     ! routine called in sbcmod 
    PUBLIC   sbc_wave        ! routine called in sbcmod
    PUBLIC   sbc_wave_init   ! routine called in sbcmod
    
    ! Variables checking if the wave parameters are coupled (if not, they are read from file)
-   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_tauwoc = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_hsig          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_phioc         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_sdrftx        = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_sdrfty        = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wper          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wnum          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wstrf         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wdrag         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_charn         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_taw           = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_bhd           = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tusd          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tvsd          = .FALSE.
 
    INTEGER ::   jpfld    ! number of files to read for stokes drift
@@ -53,22 +57,32 @@
    INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point
    INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
-   INTEGER ::   jp_wfr   ! index of wave peak frequency         (1/s)    at T-point
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_cd      ! structure of input fields (file informations, fields read) Drag Coefficient
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sd      ! structure of input fields (file informations, fields read) Stokes Drift
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_wn      ! structure of input fields (file informations, fields read) wave number for Qiao
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauwoc  ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauw    ! structure of input fields (file informations, fields read) ocean stress components from wave model
-
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:  
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x, tauw_y      !:  
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp.
-
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauoc   ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
+
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave        !: Neutral drag coefficient at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw             !: Significant Wave Height at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wmp             !: Wave Mean Period at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wnum            !: Wave Number at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave      !: stress reduction factor  at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d           !: Surface Stokes Drift module at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd          !: barotropic stokes drift divergence
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd    !: surface Stokes drift velocities at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd, vsd, wsd   !: Stokes drift velocities at u-, v- & w-points, resp.u
+!
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   charn           !: charnock coefficient at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tawx            !: Net wave-supported stress, u
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tawy            !: Net wave-supported stress, v
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   twox            !: wave-ocean momentum flux, u
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   twoy            !: wave-ocean momentum flux, v
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wavex     !: stress reduction factor  at, u component
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wavey     !: stress reduction factor  at, v component
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   phioc           !: tke flux from wave model
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   KZN2            !: Kz*N2
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   bhd_wave        !: Bernoulli head. wave induce pression
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tusd, tvsd      !: Stokes drift transport
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   ZMX             !: Kz*N2
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -88,28 +102,79 @@
       !!                2014 (DOI: 10.1175/JPO-D-14-0020.1)
       !!
-      !! ** Method  : - Calculate Stokes transport speed 
-      !!              - Calculate horizontal divergence 
-      !!              - Integrate the horizontal divergenze from the bottom 
-      !! ** action  
+      !! ** Method  : - Calculate the horizontal Stokes drift velocity (Breivik et al. 2014)
+      !!              - Calculate its horizontal divergence
+      !!              - Calculate the vertical Stokes drift velocity
+      !!              - Calculate the barotropic Stokes drift divergence
+      !!
+      !! ** action  : - tsd2d         : module of the surface Stokes drift velocity
+      !!              - usd, vsd, wsd : 3 components of the Stokes drift velocity
+      !!              - div_sd        : barotropic Stokes drift divergence
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: Kmm ! ocean time level index
       INTEGER  ::   jj, ji, jk   ! dummy loop argument
       INTEGER  ::   ik           ! local integer 
-      REAL(wp) ::  ztransp, zfac, zsp0
-      REAL(wp) ::  zdepth, zsqrt_depth,  zexp_depth, z_two_thirds, zsqrtpi !sqrt of pi
-      REAL(wp) ::  zbot_u, zbot_v, zkb_u, zkb_v, zke3_u, zke3_v, zda_u, zda_v
-      REAL(wp) ::  zstokes_psi_u_bot, zstokes_psi_v_bot
-      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd     ! 2D workspace
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zstokes_psi_u_top, zstokes_psi_v_top ! 2D workspace
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh                              ! 3D workspace
-      !!---------------------------------------------------------------------
-      !
-      ALLOCATE( ze3divh(jpi,jpj,jpkm1) )   ! jpkm1 -> avoid lbc_lnk on jpk that is not defined
+      REAL(wp) ::  ztransp, zfac, ztemp, zsp0, zsqrt, zbreiv16_w
+      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v, sdtrp
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd ! 2D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh, zInt_w                  ! 3D workspace
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( ze3divh(jpi,jpj,jpkm1) ) ! jpkm1 -> avoid lbc_lnk on jpk that is not defined
+      ALLOCATE( zInt_w(jpi,jpj,jpk) )
       ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) )
+      zk_t    (:,:) = 0._wp
+      zk_u    (:,:) = 0._wp
+      zk_v    (:,:) = 0._wp
+      zu0_sd  (:,:) = 0._wp
+      zv0_sd  (:,:) = 0._wp
+      ze3divh (:,:,:) = 0._wp
+
       !
       ! select parameterization for the calculation of vertical Stokes drift
       ! exp. wave number at t-point
-      IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) )
+      IF( ln_breivikFV_2016 ) THEN
+      ! Assumptions :  ut0sd and vt0sd are surface Stokes drift at T-points
+      !                sdtrp is the norm of Stokes transport
+      !
+         zfac = 0.166666666667_wp
+         DO_2D( 1, 1, 1, 1 ) ! In the deep-water limit we have ke = ||ust0||/( 6 * ||transport|| )
+            zsp0          = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) !<-- norm of Surface Stokes drift
+            tsd2d(ji,jj)  = zsp0
+            IF( cpl_tusd .AND. cpl_tvsd ) THEN  !stokes transport is provided in coupled mode
+               sdtrp      = SQRT( tusd(ji,jj)*tusd(ji,jj) + tvsd(ji,jj)*tvsd(ji,jj) )  !<-- norm of Surface Stokes drift transport
+            ELSE 
+               ! Stokes drift transport estimated from Hs and Tmean 
+               sdtrp      = 2.0_wp * rpi / 16.0_wp *                             &
+                   &        hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
+            ENDIF
+            zk_t (ji,jj)  = zfac * zsp0 / MAX ( sdtrp, 0.0000001_wp ) !<-- ke = ||ust0||/( 6 * ||transport|| )
+         END_2D
+      !# define zInt_w ze3divh
+         DO_3D( 1, 1, 1, 1, 1, jpk ) ! Compute the primitive of Breivik 2016 function at W-points
+            zfac             = - 2._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm)  !<-- zfac should be negative definite
+            ztemp            = EXP ( zfac )
+            zsqrt            = SQRT( -zfac )
+            zbreiv16_w       = ztemp - SQRT(rpi)*zsqrt*ERFC(zsqrt) !Eq. 16 Breivik 2016
+            zInt_w(ji,jj,jk) = ztemp - 4._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) * zbreiv16_w
+         END_3D
+!
+         DO jk = 1, jpkm1
+            zfac = 0.166666666667_wp
+            DO_2D( 1, 1, 1, 1 ) !++ Compute the FV Breivik 2016 function at T-points
+               zsp0          = zfac / MAX(zk_t (ji,jj),0.0000001_wp)
+               ztemp         = zInt_w(ji,jj,jk) - zInt_w(ji,jj,jk+1)
+               zu0_sd(ji,jj) = ut0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk)
+               zv0_sd(ji,jj) = vt0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk)
+            END_2D
+            DO_2D( 1, 0, 1, 0 ) ! ++ Interpolate at U/V points
+               zfac          =  1.0_wp / e3u(ji  ,jj,jk,Kmm)
+               usd(ji,jj,jk) =  0.5_wp * zfac * ( zu0_sd(ji,jj)+zu0_sd(ji+1,jj) ) * umask(ji,jj,jk)
+               zfac          =  1.0_wp / e3v(ji  ,jj,jk,Kmm)
+               vsd(ji,jj,jk) =  0.5_wp * zfac * ( zv0_sd(ji,jj)+zv0_sd(ji,jj+1) ) * vmask(ji,jj,jk)
+            END_2D
+         ENDDO
+      !# undef zInt_w
+      !
+      ELSE
          zfac = 2.0_wp * rpi / 16.0_wp
          DO_2D( 1, 1, 1, 1 )
@@ -128,23 +193,11 @@
             zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
          END_2D
-      ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
-         DO_2D( 1, 1, 1, 1 )
-            zk_t(ji,jj) = ( 2.0_wp * rpi * wfreq(ji,jj) ) * ( 2.0_wp * rpi * wfreq(ji,jj) ) / grav
-         END_2D
-         DO_2D( 1, 0, 1, 0 )
-            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
-            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
-            !
-            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
-            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
-         END_2D
-      ENDIF
-      !
+
       !                       !==  horizontal Stokes Drift 3D velocity  ==!
-      IF( ll_st_bv2014 ) THEN
+
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zdep_u = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji+1,jj,jk,Kmm) )
             zdep_v = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji,jj+1,jk,Kmm) )
-            !                          
+            !
             zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
             zkh_v = zk_v(ji,jj) * zdep_v
@@ -156,45 +209,4 @@
             vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
          END_3D
-      ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN
-         ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) )
-         DO_2D( 1, 0, 1, 0 )
-            zstokes_psi_u_top(ji,jj) = 0._wp
-            zstokes_psi_v_top(ji,jj) = 0._wp
-         END_2D
-         zsqrtpi = SQRT(rpi)
-         z_two_thirds = 2.0_wp / 3.0_wp
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! exp. wave number & Stokes drift velocity at u- & v-points
-            zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth
-            zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth
-            zkb_u  = zk_u(ji,jj) * zbot_u                             ! 2 * k * bottom depth
-            zkb_v  = zk_v(ji,jj) * zbot_v                             ! 2 * k * bottom depth
-            !
-            zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u(ji,jj,jk,Kmm))     ! 2k * thickness
-            zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v(ji,jj,jk,Kmm))     ! 2k * thickness
-
-            ! Depth attenuation .... do u component first..
-            zdepth      = zkb_u
-            zsqrt_depth = SQRT(zdepth)
-            zexp_depth  = EXP(-zdepth)
-            zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
-                 &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
-                 &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
-            zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u
-            zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot
-
-            !         ... and then v component
-            zdepth      =zkb_v
-            zsqrt_depth = SQRT(zdepth)
-            zexp_depth  = EXP(-zdepth)
-            zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
-                 &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
-                 &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
-            zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v
-            zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot
-            !
-            usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
-            vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
-         END_3D
-         DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top )
       ENDIF
 
@@ -235,63 +247,28 @@
       CALL iom_put( "vstokes",  vsd  )
       CALL iom_put( "wstokes",  wsd  )
-      !
-      DEALLOCATE( ze3divh )
+!      !
+      DEALLOCATE( ze3divh, zInt_w )
       DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
       !
    END SUBROUTINE sbc_stokes
-
-
-   SUBROUTINE sbc_wstress( )
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE sbc_wstress  ***
-      !!
-      !! ** Purpose :   Updates the ocean momentum modified by waves
-      !!
-      !! ** Method  : - Calculate u,v components of stress depending on stress
-      !!                model 
-      !!              - Calculate the stress module
-      !!              - The wind module is not modified by waves 
-      !! ** action  
-      !!---------------------------------------------------------------------
-      INTEGER  ::   jj, ji   ! dummy loop argument
-      !
-      IF( ln_tauwoc ) THEN
-         utau(:,:) = utau(:,:)*tauoc_wave(:,:)
-         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
-         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
-      ENDIF
-      !
-      IF( ln_tauw ) THEN
-         DO_2D( 1, 0, 1, 0 )
-            ! Stress components at u- & v-points
-            utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
-            vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
-            !
-            ! Stress module at t points
-            taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
-         END_2D
-         CALL lbc_lnk_multi( 'sbcwave', utau(:,:), 'U', -1.0_wp , vtau(:,:), 'V', -1.0_wp , taum(:,:) , 'T', -1.0_wp )
-      ENDIF
-      !
-   END SUBROUTINE sbc_wstress
-
-
+!
+!
    SUBROUTINE sbc_wave( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_wave  ***
       !!
-      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
-      !!
-      !! ** Method  : - Read namelist namsbc_wave
-      !!              - Read Cd_n10 fields in netcdf files 
-      !!              - Read stokes drift 2d in netcdf files 
-      !!              - Read wave number in netcdf files 
-      !!              - Compute 3d stokes drift using Breivik et al.,2014
-      !!                formulation
-      !! ** action  
+      !! ** Purpose :   read wave parameters from wave model in netcdf files
+      !!                or from a coupled wave mdoel
+      !!
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
       INTEGER, INTENT(in   ) ::   Kmm  ! ocean time index
       !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave : update the read waves fields'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
       !
       IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==!
@@ -300,40 +277,44 @@
       ENDIF
 
-      IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN  !==  Wave induced stress  ==!
-         CALL fld_read( kt, nn_fsbc, sf_tauwoc )         ! read wave norm stress from external forcing
-         tauoc_wave(:,:) = sf_tauwoc(1)%fnow(:,:,1) * tmask(:,:,1)
-      ENDIF
-
-      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==!
-         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid)
-         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1) * tmask(:,:,1)
-         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) * tmask(:,:,1)
-      ENDIF
-
-      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==! 
+      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read stress reduction factor due to wave from external forcing
+         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) * tmask(:,:,1)
+      ELSEIF ( ln_taw .AND. cpl_taw ) THEN
+         IF (kt < 1) THEN ! The first fields gave by OASIS have very high erroneous values ....
+            twox(:,:)=0._wp
+            twoy(:,:)=0._wp
+            tawx(:,:)=0._wp
+            tawy(:,:)=0._wp
+            tauoc_wavex(:,:) = 1._wp
+            tauoc_wavey(:,:) = 1._wp
+         ELSE
+            tauoc_wavex(:,:) = abs(twox(:,:)/tawx(:,:))
+            tauoc_wavey(:,:) = abs(twoy(:,:)/tawy(:,:))
+         ENDIF
+      ENDIF
+
+      IF ( ln_phioc .and. cpl_phioc .and.  kt == nit000 ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave : PHIOC from wave model'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+
+      IF( ln_sdw .AND. .NOT. cpl_sdrftx)  THEN       !==  Computation of the 3d Stokes Drift  ==! 
          !
          IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
             CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
+            !                                            ! NB: test case mode, not read as jpfld=0
             IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) * tmask(:,:,1)  ! significant wave height
             IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) * tmask(:,:,1)  ! wave mean period
-            IF( jp_wfr > 0 )   wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1) * tmask(:,:,1)  ! Peak wave frequency
             IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) * tmask(:,:,1)  ! 2D zonal Stokes Drift at T point
             IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) * tmask(:,:,1)  ! 2D meridional Stokes Drift at T point
          ENDIF
          !
-         ! Read also wave number if needed, so that it is available in coupling routines
-         IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN
-            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing
-            wnum(:,:) = sf_wn(1)%fnow(:,:,1) * tmask(:,:,1)
-         ENDIF
-           
-         ! Calculate only if required fields have been read
-         ! In coupled wave model-NEMO case the call is done after coupling
+         IF( jpfld == 4 .OR. ln_wave_test )   &
+            &      CALL sbc_stokes( Kmm )                 ! Calculate only if all required fields are read
+            !                                            ! or in wave test case
+         !  !                                            ! In coupled case the call is done after (in sbc_cpl)
+      ENDIF
          !
-         IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &
-           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes( Kmm )
-         !
-      ENDIF
-      !
    END SUBROUTINE sbc_wave
 
@@ -343,12 +324,9 @@
       !!                     ***  ROUTINE sbc_wave_init  ***
       !!
-      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
+      !! ** Purpose :   Initialisation fo surface waves
       !!
       !! ** Method  : - Read namelist namsbc_wave
-      !!              - Read Cd_n10 fields in netcdf files 
-      !!              - Read stokes drift 2d in netcdf files 
-      !!              - Read wave number in netcdf files 
-      !!              - Compute 3d stokes drift using Breivik et al.,2014
-      !!                formulation
+      !!              - create the structure used to read required wave fields
+      !!                (its size depends on namelist options)
       !! ** action  
       !!---------------------------------------------------------------------
@@ -357,129 +335,185 @@
       !!
       CHARACTER(len=100)     ::  cn_dir                            ! Root directory for location of drag coefficient files
-      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i            ! array of namelist informations on the fields to read
       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
-                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, &
-                             &   sn_tauwoc, sn_tauwx, sn_tauwy     ! informations about the fields to be read
-      !
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, &
-                             sn_wnum, sn_tauwoc, sn_tauwx, sn_tauwy
-      !!---------------------------------------------------------------------
+                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc    ! informations about the fields to be read
+      !
+      NAMELIST/namsbc_wave/ cn_dir, sn_cdg, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc,   &
+         &                  ln_cdgw, ln_sdw, ln_tauoc, ln_stcor, ln_charn, ln_taw, ln_phioc,     &
+         &                  ln_wave_test, ln_bern_srfc, ln_breivikFV_2016, ln_vortex_force, ln_stshear
+      !!---------------------------------------------------------------------
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave_init : surface waves in the system'
+         WRITE(numout,*) '~~~~~~~~~~~~~ '
+      ENDIF
       !
       READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_wave in reference namelist' )
-         
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_wave in reference namelist')
+
       READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist' )
       IF(lwm) WRITE ( numond, namsbc_wave )
       !
-      IF( ln_cdgw ) THEN
-         IF( .NOT. cpl_wdrag ) THEN
-            ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist namsbc_wave'
+         WRITE(numout,*) '      Stokes drift                                  ln_sdw = ', ln_sdw
+         WRITE(numout,*) '      Breivik 2016                       ln_breivikFV_2016 = ', ln_breivikFV_2016
+         WRITE(numout,*) '      Stokes Coriolis & tracer advection terms    ln_stcor = ', ln_stcor
+         WRITE(numout,*) '      Vortex Force                         ln_vortex_force = ', ln_vortex_force
+         WRITE(numout,*) '      Bernouilli Head Pressure                ln_bern_srfc = ', ln_bern_srfc
+         WRITE(numout,*) '      wave modified ocean stress                  ln_tauoc = ', ln_tauoc
+         WRITE(numout,*) '      neutral drag coefficient (CORE bulk only)    ln_cdgw = ', ln_cdgw
+         WRITE(numout,*) '      charnock coefficient                        ln_charn = ', ln_charn
+         WRITE(numout,*) '      Stress modificated by wave                    ln_taw = ', ln_taw
+         WRITE(numout,*) '      TKE flux from wave                          ln_phioc = ', ln_phioc
+         WRITE(numout,*) '      Surface shear with Stokes drift           ln_stshear = ', ln_stshear
+         WRITE(numout,*) '      Test with constant wave fields          ln_wave_test = ', ln_wave_test
+      ENDIF
+
+      !                                ! option check
+      IF( .NOT.( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_charn) )   &
+         &     CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')
+      IF( ln_cdgw .AND. ln_blk )   &
+         &     CALL ctl_stop( 'drag coefficient read from wave model NOT available yet with aerobulk package')
+      IF( ln_stcor .AND. .NOT.ln_sdw )   &
+         &     CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+
+      !                             !==  Allocate wave arrays  ==!
+      ALLOCATE( ut0sd (jpi,jpj)    , vt0sd (jpi,jpj) )
+      ALLOCATE( hsw   (jpi,jpj)    , wmp   (jpi,jpj) )
+      ALLOCATE( wnum  (jpi,jpj) )
+      ALLOCATE( tsd2d (jpi,jpj)    , div_sd(jpi,jpj)    , bhd_wave(jpi,jpj)     )
+      ALLOCATE( usd   (jpi,jpj,jpk), vsd   (jpi,jpj,jpk), wsd     (jpi,jpj,jpk) )
+      ALLOCATE( tusd  (jpi,jpj)    , tvsd  (jpi,jpj)    , ZMX     (jpi,jpj,jpk) )
+      usd   (:,:,:) = 0._wp
+      vsd   (:,:,:) = 0._wp
+      wsd   (:,:,:) = 0._wp
+      hsw     (:,:) = 0._wp
+      wmp     (:,:) = 0._wp
+      ut0sd   (:,:) = 0._wp
+      vt0sd   (:,:) = 0._wp
+      tusd    (:,:) = 0._wp
+      tvsd    (:,:) = 0._wp
+      bhd_wave(:,:) = 0._wp
+      ZMX   (:,:,:) = 0._wp
+!
+      IF( ln_wave_test ) THEN       !==  Wave TEST case  ==!   set uniform waves fields
+         jpfld    = 0                   ! No field read
+         ln_cdgw  = .FALSE.             ! No neutral wave drag input
+         ln_tauoc = .FALSE.             ! No wave induced drag reduction factor
+         ut0sd(:,:) = 0.13_wp * tmask(:,:,1)   ! m/s
+         vt0sd(:,:) = 0.00_wp                  ! m/s
+         hsw  (:,:) = 2.80_wp                  ! meters
+         wmp  (:,:) = 8.00_wp                  ! seconds
+         !
+      ELSE                          !==  create the structure associated with fields to be read  ==!
+         IF( ln_cdgw ) THEN                       ! wave drag
+            IF( .NOT. cpl_wdrag ) THEN
+               ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+               !
+                                      ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
+               IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
+            ENDIF
+            ALLOCATE( cdn_wave(jpi,jpj) )
+            cdn_wave(:,:) = 0._wp
+         ENDIF
+         IF( ln_charn ) THEN                     ! wave drag
+            IF( .NOT. cpl_charn ) THEN
+               CALL ctl_stop( 'STOP', 'Charnock based wind stress can be used in coupled mode only' )
+            ENDIF
+            ALLOCATE( charn(jpi,jpj) )
+            charn(:,:) = 0._wp
+         ENDIF
+         IF( ln_taw ) THEN                     ! wind stress
+            IF( .NOT. cpl_taw ) THEN
+               CALL ctl_stop( 'STOP', 'wind stress from wave model can be used in coupled mode only, use ln_cdgw instead' )
+            ENDIF
+            ALLOCATE( tawx(jpi,jpj) )
+            ALLOCATE( tawy(jpi,jpj) )
+            ALLOCATE( twox(jpi,jpj) )
+            ALLOCATE( twoy(jpi,jpj) )
+            ALLOCATE( tauoc_wavex(jpi,jpj) )
+            ALLOCATE( tauoc_wavey(jpi,jpj) )
+            tawx(:,:) = 0._wp
+            tawy(:,:) = 0._wp
+            twox(:,:) = 0._wp
+            twoy(:,:) = 0._wp
+            tauoc_wavex(:,:) = 1._wp
+            tauoc_wavey(:,:) = 1._wp
+         ENDIF
+
+         IF( ln_phioc ) THEN                     ! TKE flux
+            IF( .NOT. cpl_phioc ) THEN
+                CALL ctl_stop( 'STOP', 'phioc can be used in coupled mode only' )
+            ENDIF
+            ALLOCATE( phioc(jpi,jpj) )
+            phioc(:,:) = 0._wp
+         ENDIF
+
+         IF( ln_tauoc ) THEN                    ! normalized wave stress into the ocean
+            IF( .NOT. cpl_wstrf ) THEN
+               ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' )
+               !
+                                       ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
+               IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
+            ENDIF
+            ALLOCATE( tauoc_wave(jpi,jpj) )
+            tauoc_wave(:,:) = 0._wp
+         ENDIF
+
+         IF( ln_sdw ) THEN                      ! Stokes drift
+            ! 1. Find out how many fields have to be read from file if not coupled
+            jpfld=0
+            jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0
+            IF( .NOT. cpl_sdrftx ) THEN
+               jpfld  = jpfld + 1
+               jp_usd = jpfld
+            ENDIF
+            IF( .NOT. cpl_sdrfty ) THEN
+               jpfld  = jpfld + 1
+               jp_vsd = jpfld
+            ENDIF
+            IF( .NOT. cpl_hsig ) THEN
+               jpfld  = jpfld + 1
+               jp_hsw = jpfld
+            ENDIF
+            IF( .NOT. cpl_wper ) THEN
+               jpfld  = jpfld + 1
+               jp_wmp = jpfld
+            ENDIF
+            ! 2. Read from file only the non-coupled fields 
+            IF( jpfld > 0 ) THEN
+               ALLOCATE( slf_i(jpfld) )
+               IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
+               IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
+               IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
+               IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
+               ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+               !
+               DO ifpr= 1, jpfld
+                  ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
+                  IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
+               END DO
+               !
+               CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
+            ENDIF
             !
-                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
-            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( cdn_wave(jpi,jpj) )
-      ENDIF
-
-      IF( ln_tauwoc ) THEN
-         IF( .NOT. cpl_tauwoc ) THEN
-            ALLOCATE( sf_tauwoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwoc
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+            ! 3. Wave number (only needed for Qiao parametrisation, ln_zdfqiao=T)
+            IF( .NOT. cpl_wnum ) THEN
+               ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' )
+                                      ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
+               IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            ENDIF
             !
-                                     ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   )
-            IF( sn_tauwoc%ln_tint )  ALLOCATE( sf_tauwoc(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_tauwoc, (/ sn_tauwoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( tauoc_wave(jpi,jpj) )
-      ENDIF
-
-      IF( ln_tauw ) THEN
-         IF( .NOT. cpl_tauw ) THEN
-            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' )
-            !
-            ALLOCATE( slf_j(2) )
-            slf_j(1) = sn_tauwx
-            slf_j(2) = sn_tauwy
-                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   )
-                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   )
-            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) )
-            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( tauw_x(jpi,jpj) )
-         ALLOCATE( tauw_y(jpi,jpj) )
-      ENDIF
-
-      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled
-         jpfld=0
-         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
-         IF( .NOT. cpl_sdrftx ) THEN
-            jpfld  = jpfld + 1
-            jp_usd = jpfld
-         ENDIF
-         IF( .NOT. cpl_sdrfty ) THEN
-            jpfld  = jpfld + 1
-            jp_vsd = jpfld
-         ENDIF
-         IF( .NOT. cpl_hsig  .AND. ll_st_bv_li  ) THEN
-            jpfld  = jpfld + 1
-            jp_hsw = jpfld
-         ENDIF
-         IF( .NOT. cpl_wper  .AND. ll_st_bv_li  ) THEN
-            jpfld  = jpfld + 1
-            jp_wmp = jpfld
-         ENDIF
-         IF( .NOT. cpl_wfreq .AND. ll_st_peakfr ) THEN
-            jpfld  = jpfld + 1
-            jp_wfr = jpfld
-         ENDIF
-
-         ! Read from file only the non-coupled fields 
-         IF( jpfld > 0 ) THEN
-            ALLOCATE( slf_i(jpfld) )
-            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
-            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
-            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
-            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
-            IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
-
-            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
-            !
-            DO ifpr= 1, jpfld
-               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
-               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
-            END DO
-            !
-            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
-         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     )
-         ALLOCATE( wfreq(jpi,jpj) )
-         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
-         ALLOCATE( div_sd(jpi,jpj) )
-         ALLOCATE( tsd2d (jpi,jpj) )
-
-         ut0sd(:,:) = 0._wp
-         vt0sd(:,:) = 0._wp
-         hsw(:,:) = 0._wp
-         wmp(:,:) = 0._wp
-
-         usd(:,:,:) = 0._wp
-         vsd(:,:,:) = 0._wp
-         wsd(:,:,:) = 0._wp
-         ! Wave number needed only if ln_zdfswm=T
-         IF( .NOT. cpl_wnum ) THEN
-            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' )
-                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
-            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( wnum(jpi,jpj) )
+         ENDIF
+         !
       ENDIF
       !
Index: NEMO/trunk/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- NEMO/trunk/src/OCE/ZDF/zdfsh2.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/ZDF/zdfsh2.F90	(revision 14007)
@@ -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/trunk/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/trunk/src/OCE/ZDF/zdftke.F90	(revision 14006)
+++ NEMO/trunk/src/OCE/ZDF/zdftke.F90	(revision 14007)
@@ -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
