Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5983)
@@ -253,4 +253,6 @@
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
    ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
+   ln_tauoc= .false.       !  Activate ocean stress modified by external wave induced stress (T => fill namsbc_wave)
+   ln_stcor= .false.       !  Activate Stokes Coriolis term (T => fill namsbc_wave)
    nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
                            !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
@@ -915,4 +917,5 @@
    ln_zdfexp   = .false.   !  time-stepping: split-explicit (T) or implicit (F) time stepping
    nn_zdfexp   =    3            !  number of sub-timestep for ln_zdfexp=T
+   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010)
 /
 !-----------------------------------------------------------------------
@@ -1240,5 +1243,8 @@
    sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
-   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_swh      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_tauoc    -  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
 !
    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5983)
@@ -41,4 +41,5 @@
    USE sbcapr          ! surface boundary condition: atmospheric pressure
    USE dynadv, ONLY: ln_dynadv_vec
+   USE sbcwave, ONLY: usd2d, vsd2d
 #if defined key_agrif
    USE agrif_opa_interp ! agrif
@@ -422,4 +423,20 @@
       ENDIF  
       !
+      ! Add Stokes Coriolis if defined
+      IF ( ln_stcor ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+
+               zy1 = ff(ji  ,jj-1) * ( vsd2d(ji  ,jj-1) + vsd2d(ji+1,jj-1) )
+               zy2 = ff(ji  ,jj  ) * ( vsd2d(ji  ,jj  ) + vsd2d(ji+1,jj  ) )
+               zx1 = ff(ji-1,jj  ) * ( usd2d(ji-1,jj  ) + usd2d(ji-1,jj+1) )
+               zx2 = ff(ji  ,jj  ) * ( usd2d(ji  ,jj  ) + usd2d(ji  ,jj+1) )
+
+               zu_frc(ji,jj) = zu_frc(ji,jj) + 0.25 * ( zy1 + zy2 ) * hur(ji,jj)
+               zv_frc(ji,jj) = zv_frc(ji,jj) - 0.25 * ( zx1 + zx2 ) * hvr(ji,jj)
+           END DO
+         END DO
+      ENDIF
+      !
       IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
          IF (ln_bt_fw) THEN
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 5983)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 5983)
@@ -0,0 +1,83 @@
+MODULE dynstcor
+   !!======================================================================
+   !!                       ***  MODULE  dynstcor  ***
+   !! Ocean dynamics: Stokes-Coriolis effect.
+   !!
+   !!======================================================================
+   !! History :  0.1  !  2012-10  (Oyvind Breivik)
+   !!----------------------------------------------------------------------
+   !! dyn_stcor        : Add the Stokes-Coriolis forcing to the momentum equation
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! distributed memory computing
+   USE prtctl           ! Print control
+   USE phycst
+   USE lbclnk
+   USE wrk_nemo         ! Memory Allocation
+   USE sbcmod           ! Access to ln_stcor (sbc_oce) and wave parameters (sbc_wave)
+   USE sbcwave          ! Wave module
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dyn_stcor      ! routine called by step.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2010)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+ CONTAINS
+
+   SUBROUTINE dyn_stcor( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_stcor ***
+      !!
+      !! ** Purpose:  Add Stokes-Coriolis forcing to horizontal momentum equation.
+      !!
+      !! ** History:  3.4  !   2012-10  (O. Breivik)   Initial version
+      !!              3.6  !   2014-10  (R. Benshila) 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt  ! ocean time-step index
+      !!
+      INTEGER  ::  ji, jj, jk        ! dummy loop indices
+      REAL(wp) :: zx1, zx2, zy1, zy2
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dyn_stcor')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_stcor : time stepping'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+
+      !
+      ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities 
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zy1 = ff(ji  ,jj-1) * ( vsd3d(ji  ,jj-1,jk) + vsd3d(ji+1,jj-1,jk) )
+               zy2 = ff(ji  ,jj  ) * ( vsd3d(ji  ,jj  ,jk) + vsd3d(ji+1,jj  ,jk) )
+               zx1 = ff(ji-1,jj  ) * ( usd3d(ji-1,jj  ,jk) + usd3d(ji-1,jj+1,jk) )
+               zx2 = ff(ji  ,jj  ) * ( usd3d(ji  ,jj  ,jk) + usd3d(ji  ,jj+1,jk) )
+
+               ua(ji,jj,jk) = ua(ji,jj,jk) + 0.25 * ( zy1 + zy2 ) *umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( zx1 + zx2 ) *vmask(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dynst_cor')
+      !
+   END SUBROUTINE dyn_stcor
+
+END MODULE dynstcor
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5983)
@@ -9,4 +9,5 @@
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
    !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) add wave contribution to surface vertical velocity 
    !!----------------------------------------------------------------------
 
@@ -38,4 +39,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE sbcwave,  ONLY: usd2dt, vsd2dt,wsd3d
 
    IMPLICIT NONE
@@ -160,4 +162,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::  sshnu, sshnv, dsshnu, dsshnv
       !!----------------------------------------------------------------------
       !
@@ -206,5 +209,40 @@
          END DO
       ENDIF
-
+!
+!     In case ln_wave and ln_sdw the surface vertical velocity is modified
+!     accounting for Sokes Drift velocity
+      IF ( ln_wave .AND. ln_sdw )  THEN
+       ALLOCATE(sshnu(jpi,jpj),sshnv(jpi,jpj),dsshnu(jpi,jpj),dsshnv(jpi,jpj))
+       sshnu (:,:) = 0._wp
+       sshnv (:,:) = 0._wp
+       dsshnu(:,:) = 0._wp
+       dsshnv(:,:) = 0._wp
+       ! sshn interpolated on U-V grid
+       !---------------------------------
+       DO jj = 1 , jpjm1
+         DO ji = 1 , jpim1
+           sshnu(ji,jj) =  0.5 * ( 2. - umask(ji,jj,1) ) *               &
+                        &        ( sshn(ji  ,jj) * tmask(ji  ,jj,1)      &
+                        &        + sshn(ji+1,jj) * tmask(ji+1,jj,1) )
+           sshnv(ji,jj) =  0.5 * ( 2. - vmask(ji,jj,1) ) *               &
+                        &        ( sshn(ji,jj  ) * tmask(ji,jj  ,1)      &
+                        &        + sshn(ji,jj+1) * tmask(ji,jj+1,1) )
+         ENDDO
+       ENDDO
+       ! Compute d(ssh)/dx  and d(ssh)/dy  
+       !---------------------------------
+       DO jj = 1 , jpjm1
+         DO ji = 1 , jpim1
+           dsshnu(ji,jj) = ( sshnu(ji+1,jj) - sshnu(ji,jj) ) / e1u(ji,jj)
+           dsshnv(ji,jj) = ( sshnv(ji,jj+1) - sshnv(ji,jj) ) / e2v(ji,jj)
+         ENDDO
+       ENDDO
+       ! Compute the surface vertical velocity accounting for the Stokes Drift 
+       !---------------------------------------------------------------------
+       wn(:,:,1) = wn(:,:,1) + usd2dt(:,:) * dsshnu(:,:)     &
+                 &           + vsd2dt(:,:) * dsshnv(:,:)      &
+                 &           - ( wsd3d (:,:,1) ) * tmask(:,:,1)
+      ENDIF
+      !
 #if defined key_bdy
       IF( lk_bdy ) THEN
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5983)
@@ -65,4 +65,6 @@
    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_tauoc       !: true if normalized stress from wave is used
+   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
    !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5983)
@@ -88,5 +88,5 @@
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
-         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
+         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx  , nn_components, ln_cpl 
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -360,5 +360,16 @@
                              CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
       END SELECT
-
+      IF (ln_wave .AND. ln_tauoc) THEN
+            utau(:,:) = utau(:,:)*tauoc_wave(:,:)
+            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
+            taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+      !
+            SELECT CASE( nsbc )
+            CASE(  0,1,2,3,5,-1 )  ;
+                IF(lwp) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. &
+                        & If not requested select ln_tauoc=.false'
+            END SELECT
+      !
+      END IF
       IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
 
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 5983)
@@ -6,4 +6,5 @@
    !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
    !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
+   !! History :  3.6  !2014-09  (Clementi E, Oddo P)New Stokes Drift Computation
    !!----------------------------------------------------------------------
 
@@ -12,5 +13,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! 
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
+   USE sbc_oce	      ! Surface boundary condition: ocean fields
    USE bdy_oce        !
    USE domvvl         !
@@ -19,6 +20,7 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing library
-   USE fldread	       ! read input fields
+   USE fldread	      ! read input fields
    USE wrk_nemo       !
+   USE phycst         ! physical constants 
 
    IMPLICIT NONE
@@ -27,15 +29,20 @@
    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
    
-   INTEGER , PARAMETER ::   jpfld  = 3   ! maximum number of files to read for srokes drift
-   INTEGER , PARAMETER ::   jp_usd = 1   ! index of stokes drift  (i-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_vsd = 2   ! index of stokes drift  (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_wn  = 3   ! index of wave number                 (1/m)    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
-
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:)   :: cdn_wave 
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:,:) :: usd3d, vsd3d, wsd3d 
-   REAL(wp),         ALLOCATABLE, DIMENSION (:,:)   :: usd2d, vsd2d, uwavenum, vwavenum 
+   INTEGER , PARAMETER ::   jpfld  = 4           ! number of files to read for stokes drift
+   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
+   INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
+   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave hight      (m)      at T-point
+   INTEGER , PARAMETER ::   jp_wmp = 4           ! index of mean wave period            (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_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: usd2d,vsd2d
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: swh,wmp,wnum
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: usd2dt,vsd2dt,tsd2d
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: tauoc_wave
 
    !! * Substitutions
@@ -58,19 +65,27 @@
       !!              - 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 monochromatic
+      !!              - Read wave number in netcdf files 
+      !!              - Compute 3d stokes drift using Breivik et al.,2014
+      !!                formulation
       !! ** action  :   
       !!---------------------------------------------------------------------
+      USE zdf_oce,  ONLY : ln_zdfqiao
+
       INTEGER, INTENT( in  ) ::   kt       ! ocean time step
       !
       INTEGER                ::   ierror   ! return error code
       INTEGER                ::   ifpr, jj,ji,jk 
-      INTEGER                ::   ios     ! Local integer output status for namelist read
+      INTEGER                ::   ios      ! Local integer output status for namelist read
+
+      REAL(wp)                       ::  ztransp,zsp0, zk, zus,zvs
+      REAL(wp), DIMENSION(jpi,jpj)   ::  zfac 
+
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
       CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
-      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read
+      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
+                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
       REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace
       !!
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc
       !!---------------------------------------------------------------------
       !
@@ -97,7 +112,20 @@
             cdn_wave(:,:) = 0.0
          ENDIF
+!
+         IF ( ln_tauoc ) THEN
+            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
+            !
+                                   ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
+            IF( sn_cdg%ln_tint )   ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            ALLOCATE( tauoc_wave(jpi,jpj) )
+            tauoc_wave(:,:) = 0.0
+        ENDIF
+!
          IF ( ln_sdw ) THEN
-            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
-            ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd;
+            slf_i(jp_swh) = sn_swh ; 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: unable to allocate sf_wave structure' )
             !
@@ -106,10 +134,22 @@
                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', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj) , vsd2d(jpi,jpj) , uwavenum(jpi,jpj) , vwavenum(jpi,jpj) )
+            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),usd2dt(jpi,jpj),vsd2dt(jpi,jpj))
             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
-            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;    uwavenum(:,:) = 0._wp
-            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ;    vwavenum(:,:) = 0._wp
-            wsd3d(:,:,:) = 0._wp
+            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
+            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;    usd2dt(:,:) = 0._wp ;
+            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ;    vsd2dt(:,:) = 0._wp ;
+            wsd3d(:,:,:) = 0._wp   ;
+            swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ;
+            IF ( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==!
+               ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: 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' )
+               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) )
+               wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp
+            ENDIF
          ENDIF
       ENDIF
@@ -119,25 +159,58 @@
          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
       ENDIF
-      !
-      IF ( ln_sdw )  THEN              !==  Computation of the 3d Stokes Drift  ==!
-         !
-         CALL fld_read( kt, nn_fsbc, sf_sd )    !* read drag coefficient from external forcing
-         !
-         !
-         CALL wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
-         !                                      !* distribute it on the vertical
-         DO jk = 1, jpkm1
-            zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * fsdept_n(:,:,jk) )
-            zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * fsdept_n(:,:,jk) )
-         END DO
-         !                                      !* interpolate the stokes drift from t-point to u- and v-points
-         DO jk = 1, jpkm1
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                   usd3d(ji,jj,jk) = 0.5_wp * ( zusd_t(ji  ,jj,jk) + zusd_t(ji+1,jj,jk) ) * umask(ji,jj,jk)
-                   vsd3d(ji,jj,jk) = 0.5_wp * ( zvsd_t(ji  ,jj,jk) + zvsd_t(ji,jj+1,jk) ) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+!
+      IF ( ln_tauoc ) THEN             !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing
+         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1)
+      ENDIF
+!
+      IF ( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==! 
+         !
+         CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
+         swh(:,:)    = sf_sd(jp_swh)%fnow(:,:,1)  ! significant wave height
+         wmp(:,:)    = sf_sd(jp_wmp)%fnow(:,:,1)  ! wave mean period
+         usd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)  ! 2D zonal Stokes Drift
+         vsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)  ! 2D meridional Stokes Drift
+         !
+         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014
+         !(DOI: 10.1175/JPO-D-14-0020.1)==! 
+         !
+         CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv )
+          DO jk = 1, jpk
+             DO jj = 1, jpj
+                DO ji = 1, jpi
+               ! On T grid
+               ! Stokes transport speed estimated from Hs and Tmean
+               ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp))
+               ! Stokes surface speed
+               zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 +  sf_sd(jp_vsd)%fnow(ji,jj,1)**2)
+               ! Wavenumber scale
+               zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
+               ! Depth attenuation
+               zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk))
+                END DO
+             END DO
+!
+             DO jj = 1, jpjm1
+                DO ji = 1, jpim1
+                 ! Into the U and V Grid 
+                 zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
+                 &                                + zfac(ji+1,jj) * tmask(ji+1,jj,1) )
+!
+                 zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
+                 &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) )
+!
+                 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( usd2dt(ji,jj) * tmask(ji,jj,1) &
+                 &                                             +  usd2dt(ji+1,jj) * tmask(ji+1,jj,1) )
+!
+                 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vsd2dt(ji,jj) * tmask(ji,jj,1) &
+                 &                                              + vsd2dt(ji,jj+1) * tmask(ji,jj+1,1) )
+!
+                 usd3d(ji,jj,jk) = usd2d(ji,jj)*zus
+                 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*zvs
+                END DO
+             END DO
+          END DO
+         !
          CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
          CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
@@ -171,6 +244,19 @@
          ENDIF
 #endif
-         CALL wrk_dealloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
-         ! 
+!         CALL wrk_dealloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
+         CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv )
+!
+        IF ( ln_zdfqiao )  THEN
+          wnum(:,:) = sf_wn(1)%fnow(:,:,1)
+         ! Calculate the module of the stokes drift on T grid
+         !-------------------------------------------------
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+                tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0  +     &
+                &               (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5
+            END DO
+         END DO
+        ENDIF
+      !
       ENDIF
       !
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5983)
@@ -9,4 +9,5 @@
    !!            3.7  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes 
    !!             -   !  2014-12  (G. Madec) suppression of cross land advection option
+   !!            3.6  !  2015-06  (E. Clementi) Addition of Stokes drift in case of wave coupling
    !!----------------------------------------------------------------------
 
@@ -34,6 +35,10 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-
-   USE diaptr          ! Poleward heat transport 
+   USE sbcwave        ! wave module
+   USE sbc_oce        ! surface boundary condition: ocean
+
+   USE diaptr         ! Poleward heat transport 
+   USE sbcwave        ! wave module
+   USE sbc_oce        ! surface boundary condition: ocean
 
    IMPLICIT NONE
@@ -95,4 +100,8 @@
       !
       !                                          ! set time step
+      zun(:,:,:) = 0.0
+      zvn(:,:,:) = 0.0
+      zwn(:,:,:) = 0.0
+      !    
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
          r2dtra(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
@@ -102,4 +111,14 @@
       !
       !                                         !==  effective transport  ==!
+      IF (ln_wave .AND. ln_sdw) THEN
+         DO jk = 1, jpkm1
+            zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) *      &
+                        &  ( un(:,:,jk) + usd3d(:,:,jk) )                     !eulerian transport + Stokes Drift
+            zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) *      &
+                        &  ( vn(:,:,jk) + vsd3d(:,:,jk) )
+            zwn(:,:,jk) = e1e2t(:,:) *                    &
+                        &  ( wn(:,:,jk) + wsd3d(:,:,jk) )
+         END DO
+      ELSE
       DO jk = 1, jpkm1
          zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
@@ -107,4 +126,5 @@
          zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
       END DO
+      ENDIF
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 5983)
@@ -35,4 +35,5 @@
    INTEGER , PUBLIC ::   nn_npc      !: non penetrative convective scheme call  frequency
    INTEGER , PUBLIC ::   nn_npcp     !: non penetrative convective scheme print frequency
+   LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag
 
 
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 5983)
@@ -51,6 +51,7 @@
       INTEGER ::   ioptio, ios       ! local integers
       !!
-      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,   &
-         &              ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp
+      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  &
+         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp,       &
+         &        ln_zdfqiao
       !!----------------------------------------------------------------------
 
@@ -81,4 +82,5 @@
          WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc
          WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp
+         WRITE(numout,*) '      Qiao formulation flag               ln_zdfqiao=', ln_zdfqiao
       ENDIF
 
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 5983)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 5983)
@@ -0,0 +1,130 @@
+MODULE zdfqiao
+   !!======================================================================
+   !!                       ***  MODULE  zdfqiao  ***
+   !! Qiao module      : vertical mixing enhancement due to surface waves 
+   !!======================================================================
+   !! History :  3.6  !  2014-10  (E. Clementi)  Original code
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   qiao_init      
+   !!   zdf_qiao        : compute Qiao parameters
+   !!----------------------------------------------------------------------
+
+   USE iom             ! I/O manager library
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE zdf_oce
+   USE sbcwave         ! wave module
+   USE dom_oce
+   
+   !!----------------------------------------------------------------------
+   !!   qiao_init       : compute QBv: Qiao terms to be added to vertical eddy
+   !!                     diffusivity and viscosity coefficients 
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_qiao    ! routine called in zdf_ric
+
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: QBv, QBvu, QBvv
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE zdf_qiao( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE zdf_qiao ***
+      !!
+      !! ** Purpose :Compute the Qiao term (QBv) to be added to
+      !!             vertical viscosity and diffusivity coeffs.  
+      !!
+      !! ** Method  :QBv = alpha * A * Us(0) * exp (3 * k * z)
+      !!             
+      !! ** action  :Compute the Qiao wave dependent term 
+      !!             only if ln_wave=.true.
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt   ! ocean time step
+      INTEGER                ::  ierror   ! return error code
+      INTEGER                ::  jj,ji,jk
+      !!---------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+      !                                         ! -------------------- !
+         ALLOCATE(QBv(jpi,jpj,jpk))
+         ALLOCATE(QBvu(jpi,jpj,jpk))
+         ALLOCATE(QBvv(jpi,jpj,jpk))
+
+      ENDIF
+
+         QBv (:,:,:) = 0.0
+         QBvu(:,:,:) = 0.0
+         QBvv(:,:,:) = 0.0
+      !
+      ! Compute the Qiao term Bv (QBv) to be added to
+      ! vertical viscosity and diffusivity
+      ! QBv = alpha * A * Us(0) * exp (3 * k * z)
+      ! alpha here is set to 1
+      !---------------------------------------------------------------------------------
+      !
+      IF ( ln_wave ) THEN
+         DO jk = 1, jpk
+          DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               QBv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *     &
+               &              exp(3.0 * wnum(ji,jj) *                            &                     
+               &              (-MIN( gdept_0(ji  ,jj  ,jk) , gdept_0(ji+1,jj  ,jk),     & 
+               &               gdept_0(ji  ,jj+1,jk),gdept_0(ji+1,jj+1,jk))))
+
+            ENDDO
+          ENDDO
+         ENDDO
+
+         QBv(jpi,:,:)=QBv(jpim1,:,:)
+         QBv(:,jpj,:)=QBv(:,jpjm1,:)
+
+         !
+         ! Interpolate Qiao parameter QBv into the grid_U and grid_V
+         !-------------------------------------------------
+
+         DO jk = 1, jpk
+           DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               QBvu(ji,jj,jk)=0.5 * ( 2. - umask(ji,jj,jk) ) *                 &
+               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                   & 
+               &           + QBv(ji+1,jj,jk) * tmask(ji+1,jj,jk) )               
+               QBvv(ji,jj,jk)=0.5 * ( 2. - vmask(ji,jj,jk) ) *                 &
+               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                   &
+               &           + QBv(ji,jj+1,jk) * tmask(ji,jj+1,jk) )
+
+             ENDDO
+            ENDDO
+          ENDDO
+!
+         QBvu(jpi,:,:)=QBvu(jpim1,:,:)
+         QBvu(:,jpj,:)=QBvu(:,jpjm1,:)
+         QBvv(jpi,:,:)=QBvv(jpim1,:,:)
+         QBvv(:,jpj,:)=QBvv(:,jpjm1,:)
+!
+        ELSE
+        !
+        CALL ctl_stop( 'STOP', 'To use Qiao formulation you have to set: ln_wave=.true.')
+        !
+        ENDIF
+!
+   END SUBROUTINE zdf_qiao
+      
+   !!======================================================================
+END MODULE zdfqiao
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5983)
@@ -26,4 +26,5 @@
    !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs
    !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state
+   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
    !!            3.7  !  2014-10  (G. Madec)  LDF simplication 
    !!             -   !  2014-12  (G. Madec) remove KPP scheme
@@ -75,5 +76,5 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   ji,jj,jk ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
@@ -132,4 +133,17 @@
       IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
       IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
+      IF( ln_zdfqiao )   THEN
+        CALL zdf_qiao(kstp )                             ! Qiao vertical mixing 
+         DO jk = 1, jpkm1
+          DO jj = 1, jpj
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = (avmu(ji,jj,jk) + QBvu(ji,jj,jk)) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = (avmv(ji,jj,jk) + QBvv(ji,jj,jk)) * vmask(ji,jj,jk)
+               avt( ji,jj,jk) = (avt( ji,jj,jk) + QBv(ji,jj,jk))  * tmask(ji,jj,jk)
+            END DO
+          END DO
+         END DO
+      ENDIF
+      !
       IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
          avt (:,:,:) = rn_avt0 * wmask (:,:,:)
@@ -218,4 +232,5 @@
                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
+      IF( ln_stcor    )  CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg       ( kstp )  ! surface pressure gradient
Index: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5982)
+++ branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5983)
@@ -19,4 +19,5 @@
    USE sbcapr           ! surface boundary condition: atmospheric pressure
    USE sbctide          ! Tide initialisation
+   USE sbcwave          ! Wave intialisation
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
@@ -41,4 +42,5 @@
    USE dynzdf           ! vertical diffusion               (dyn_zdf routine)
    USE dynspg           ! surface pressure gradient        (dyn_spg routine)
+   USE dynstcor         ! simp. form of Stokes-Coriolis
 
    USE dynnxt           ! time-stepping                    (dyn_nxt routine)
@@ -71,4 +73,5 @@
    USE zdfric           ! Richardson vertical mixing       (zdf_ric routine)
    USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
+   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine)
 
    USE zpshde           ! partial step: hor. derivative     (zps_hde routine)
