Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5899)
@@ -635,4 +635,5 @@
    rn_bfrien   =    50.    !  local multiplying factor of bfr (ln_bfr2d=T)
    ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true)
+   ln_zdfqiao  = .false.    !  Enhanced wave vertical mixing Qiao (2010)
 /
 !-----------------------------------------------------------------------
@@ -1177,4 +1178,6 @@
    sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_swh      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
 !
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/CONFIG/cfg.txt	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/CONFIG/cfg.txt	(revision 5899)
@@ -1,5 +1,4 @@
 GYRE_PISCES OPA_SRC TOP_SRC
 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
-GYRE OPA_SRC
 GYRE_XIOS OPA_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
@@ -11,2 +10,3 @@
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
+GYRE OPA_SRC
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5899)
@@ -405,4 +405,21 @@
          zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
       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
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 5899)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 5899)
@@ -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_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5899)
@@ -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 
    !!----------------------------------------------------------------------
 
@@ -42,4 +43,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE sbcwave,  ONLY: usd2dt, vsd2dt,wsd3d
 
    IMPLICIT NONE
@@ -172,4 +174,6 @@
       INTEGER             ::   ji, jj, jk   ! dummy loop indices
       REAL(wp)            ::   z1_2dt       ! local scalars
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::  sshnu, sshnv, dsshnu,dsshnv
       !!----------------------------------------------------------------------
       
@@ -220,4 +224,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
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 5899)
@@ -52,4 +52,5 @@
    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_icebergs    !: Icebergs
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5899)
@@ -84,5 +84,6 @@
       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
-         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, cn_iceflx
+         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, ln_tauoc,        &
+         &             nn_lsm, cn_iceflx
       INTEGER  ::   ios
       !!----------------------------------------------------------------------
@@ -337,5 +338,16 @@
                        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   !
       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
       !                                            !==  Misc. Options  ==!
       
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 5899)
@@ -4,15 +4,16 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3.1  !   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 
+   !!         :  3.3.1  !2011-09  (Adani M)  Original code: Drag Coefficient 
+   !!         :  3.4    !2012-10  (Adani M)                 Stokes Drift 
    !!----------------------------------------------------------------------
    USE iom             ! I/O manager library
    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 oce
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
+   USE sbc_oce         ! Surface boundary condition: ocean fields
    USE domvvl
-
+   USE phycst
    
    !!----------------------------------------------------------------------
@@ -25,13 +26,21 @@
    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 ::   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_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
+   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),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
+   REAL(wp),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
+   LOGICAL, PUBLIC :: ln_stcor = .FALSE.
 
    !! * Substitutions
@@ -54,5 +63,6 @@
       !!              - Read stokes drift 2d in netcdf files 
       !!              - Read wave number      in netcdf files 
-      !!              - Compute 3d stokes drift using monochromatic
+      !!              - Compute 3d stokes drift using Breivik et al.,2014
+      !!                formulation
       !! ** action  :   
       !!               
@@ -61,4 +71,5 @@
       USE divcur
       USE wrk_nemo
+      USE zdf_oce,  ONLY : ln_zdfqiao
 #if defined key_bdy
       USE bdy_oce, ONLY : bdytmask
@@ -69,10 +80,13 @@
       INTEGER                ::   ios     ! Local integer output status for namelist read
       REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
-      REAL                                          ::  z2dt,z1_2dt
+      REAL(wp)                                      ::  z2dt,z1_2dt
+      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
-      !!---------------------------------------------------------------------
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
+      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
+                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, ln_stcor, sn_tauoc
       !!---------------------------------------------------------------------
 
@@ -103,7 +117,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' )
             !
@@ -112,12 +139,25 @@
                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) )
-            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
-            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
+            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
+            usd2d(:,:)   = 0.0 ;  vsd2d(:,:)   = 0.0 ;
+            usd2dt(:,:)  = 0.0 ;  vsd2dt(:,:)  = 0.0 ;
+            usd3d(:,:,:) = 0.0 ;  vsd3d(:,:,:) = 0.0 ;
+            wsd3d(:,:,:) = 0.0 ;
+            swh(:,:)     = 0.0 ;  wmp(:,:)     = 0.0 ;
+            IF ( ln_zdfqiao ) THEN
+               ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fillsf_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.0; tsd2d(:,:) = 0.0
+            ENDIF
          ENDIF
       ENDIF
-         !
          !
       IF ( ln_cdgw ) THEN
@@ -125,47 +165,91 @@
          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
       ENDIF
+!
+      IF ( ln_tauoc ) THEN
+         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
-          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing
-
-         ! Interpolate wavenumber, stokes drift into the grid_V and grid_V
+          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
+          swh(:,:)    = sf_sd(jp_swh)%fnow(:,:,1)
+          wmp(:,:)    = sf_sd(jp_wmp)%fnow(:,:,1)
+          usd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1) 
+          vsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)
          !-------------------------------------------------
-
+         ! Interpolate stokes drift into the grid_V and grid_V
+         !-------------------------------------------------
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
-
-               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(2)%fnow(ji,jj+1,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) )
             END DO
          END DO
-
-          !Computation of the 3d Stokes Drift
+         !
+         CALL lbc_lnk( usd2d(:,:), 'U', -1. )
+         CALL lbc_lnk( vsd2d(:,:), 'V', -1. )
+
+          !Computation of the 3d Stokes Drift according to Breivik et al.,2014
+          !(DOI: 10.1175/JPO-D-14-0020.1)
           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, jpj-1
                 DO ji = 1, jpi-1
-                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk))))
-                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk))))
+                 ! 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) )
+
+                 usd3d(ji,jj,jk) = usd2d(ji,jj)*zus
+                 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*zvs
                 END DO
              END DO
-             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) )
-             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) )
           END DO
+         !
+         CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
+         CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
 
           CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-          
-          udummy(:,:,:)=un(:,:,:)
-          vdummy(:,:,:)=vn(:,:,:)
-          hdivdummy(:,:,:)=hdivn(:,:,:)
-          rotdummy(:,:,:)=rotn(:,:,:)
-          un(:,:,:)=usd3d(:,:,:)
-          vn(:,:,:)=vsd3d(:,:,:)
+      !-------------------------------------------------------------------          
+      ! Store horiz. velocity divergence and Rot in TMP array 
+      ! -------------------------------------------------------------------
+
+          udummy(:,:,:)     = un(:,:,:)
+          vdummy(:,:,:)     = vn(:,:,:)
+          hdivdummy(:,:,:)  = hdivn(:,:,:)
+          rotdummy(:,:,:)   = rotn(:,:,:)
+          un(:,:,:)         = usd3d(:,:,:)
+          vn(:,:,:)         = vsd3d(:,:,:)
+
+      ! Compute divergence using 3d stokes drift to calculate vertical SD
+      ! velocity
+      ! -------------------------------------------------------------------
+
           CALL div_cur(kt)
+
+      ! -------------------------------------------------------------------
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
@@ -177,16 +261,33 @@
              ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
              wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-                &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-                &                         * tmask(:,:,jk) * z1_2dt
+                &                            - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
+                &                                * tmask(:,:,jk) * z1_2dt
 #if defined key_bdy
              wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
 #endif
           END DO
-          hdivn(:,:,:)=hdivdummy(:,:,:)
-          rotn(:,:,:)=rotdummy(:,:,:)
-          vn(:,:,:)=vdummy(:,:,:)
-          un(:,:,:)=udummy(:,:,:)
+      ! Back to state variables  
+      ! -------------------------------------------------------------------
+          hdivn(:,:,:)  = hdivdummy(:,:,:)
+          rotn(:,:,:)   = rotdummy(:,:,:)
+          vn(:,:,:)     = vdummy(:,:,:)
+          un(:,:,:)     = udummy(:,:,:)
+
           CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-      ENDIF
+!
+        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
+
+
    END SUBROUTINE sbc_wave
       
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5899)
@@ -7,4 +7,5 @@
    !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!            4.0  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation
+   !!            3.6  !  2015-06  (E. Clementi) Addition of Stokes drift in case of wave coupling
    !!----------------------------------------------------------------------
 
@@ -32,5 +33,6 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-
+   USE sbcwave         ! wave module
+   USE sbc_oce         ! surface boundary condition: ocean
 
    IMPLICIT NONE
@@ -80,4 +82,8 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
+      !
+      zun(:,:,:) = 0.0
+      zvn(:,:,:) = 0.0
+      zwn(:,:,:) = 0.0
       !                                          ! set time step
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
@@ -89,10 +95,20 @@
       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_traadv( kt )       !==  Cross Land Advection  ==! (hor. advection)
       !
-      !                                               !==  effective transport  ==!
-      DO jk = 1, jpkm1
-         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
-         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
-      END DO
+      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) = e1t(:,:) * e2t(:,:)      *      &
+                        &  ( wn(:,:,jk) + wsd3d(:,:,jk) )
+         END DO
+      ELSE
+         DO jk = 1, jpkm1
+            zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
+            zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+            zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
+         END DO
+      ENDIF
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 5899)
@@ -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_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 5899)
@@ -53,6 +53,7 @@
       INTEGER ::   ios
       !!
-      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
       !!----------------------------------------------------------------------
 
@@ -83,4 +84,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_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 5899)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 5899)
@@ -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_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5899)
@@ -25,4 +25,5 @@
    !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
    !!                 !  2012-07  (J. Simeon, G. Madec. C. Ethe) Online coarsening of outputs
+   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
    !!----------------------------------------------------------------------
 
@@ -70,5 +71,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)
@@ -116,4 +117,17 @@
       IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
       IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP 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 * tmask(:,:,:)
@@ -121,4 +135,5 @@
          avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
       ENDIF
+!
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
          DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
@@ -184,4 +199,6 @@
                                   CALL dyn_ldf      ( kstp )   ! lateral mixing
           IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
+          IF( ln_stcor    )       CALL dyn_stcor    ( kstp )   ! Stokes-Coriolis forcing
+
 #if defined key_agrif
           IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
@@ -286,4 +303,6 @@
                                CALL dyn_ldf( kstp )         ! lateral mixing
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
+        IF( ln_stcor    )       CALL dyn_stcor    ( kstp )   ! Stokes-Coriolis forcing (ln_stcor set in SBC/sbc_oce.F90)
+
 #if defined key_agrif
         IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge
Index: /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5898)
+++ /branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5899)
@@ -27,4 +27,5 @@
    USE cpl_oasis3, ONLY : lk_cpl
    USE sbctide          ! Tide initialisation
+   USE sbcwave          ! Wave intialisation
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
@@ -51,4 +52,5 @@
    USE dynspg           ! surface pressure gradient        (dyn_spg routine)
    USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine)
+   USE dynstcor         ! simp. form of Stokes-Coriolis
 
    USE dynnxt           ! time-stepping                    (dyn_nxt routine)
@@ -80,4 +82,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)
