Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/ARCH/INGV/arch-IBM_EKMAN_INGV.fcm
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/ARCH/INGV/arch-IBM_EKMAN_INGV.fcm	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/ARCH/INGV/arch-IBM_EKMAN_INGV.fcm	(revision 7383)
@@ -34,2 +34,4 @@
 %USER_INC            -I%XIOS_ROOT/inc %NCDF_INC %MPI_INTEL -I/srv/lib/zlib-last/include
 %USER_LIB            -L%XIOS_ROOT/lib -lxios %NCDF_LIB -L/srv/lib/zlib-last/lib -lz
+%CC                  icc
+%CFLAGS              -O0
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7383)
@@ -288,5 +288,9 @@
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
    ln_isf      = .false.   !  ice shelf                                 (T   => fill namsbc_isf)
-   ln_wave     = .false.   !  coupling with surface wave                (T => fill namsbc_wave)
+   ln_wave     = .false.   !  Activate coupling with wave  (T => fill namsbc_wave)
+   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)
+   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 
+   ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)
+   ln_stcor    = .false.   !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & 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)
@@ -380,4 +384,7 @@
    sn_snd_crt    =   'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'
    sn_snd_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_crtw   =       'none'                 ,    'no'    ,     ''      ,         ''           , 'U,V'
+   sn_snd_ifrac  =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_snd_wlev   =       'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
 ! receive
    sn_rcv_w10m   =   'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
@@ -391,4 +398,14 @@
    sn_rcv_cal    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_co2    =   'coupled'              ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_hsig   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_iceflx =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_sdrfx  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wstrf  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
 !
    nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -535,16 +552,17 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_wave   ! External fields from wave model                        (ln_wave=T)
-!-----------------------------------------------------------------------
-!              ! file name ! frequency (hours) ! variable    ! time interp.!  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-!              !           !  (if <0  months)  !   name      !  (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_cdg      = 'cdg_wave',        1          , 'drag_coeff',   .true.    , .false., 'daily'   ,    ''      , ''       , ''
-   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'   ,    ''      , ''       , ''
-!
-   cn_dir_cdg  = './'      !  root directory for the location of drag coefficient files
-   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model
-   ln_sdw      = .false.   !  Computation of 3D stokes drift               
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_cdg      =  'sdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   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_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+!
+   cn_dir  = './'  !  root directory for the location of drag coefficient files
 /
 !-----------------------------------------------------------------------
@@ -973,4 +991,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) (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 7383)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 7383)
@@ -0,0 +1,82 @@
+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"
+   !!----------------------------------------------------------------------
+   !! 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('dyn_stcor')
+      !
+   END SUBROUTINE dyn_stcor
+
+END MODULE dynstcor
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 7383)
@@ -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=50   ! Maximum number of coupling fields
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=55   ! 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: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7383)
@@ -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/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7383)
@@ -745,5 +745,5 @@
 
       !! Neutral coefficients at 10m:
-      IF( ln_cdgw ) THEN      ! wave drag case
+      IF( ln_wave .AND. ln_cdgw ) THEN      ! wave drag case
          cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
          ztmp0   (:,:) = cdn_wave(:,:)
@@ -791,5 +791,5 @@
          END IF
        
-         IF( ln_cdgw ) THEN      ! surface wave case
+         IF( ln_wave .AND. ln_cdgw ) THEN      ! surface wave case
             sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
             Cd      = sqrt_Cd * sqrt_Cd
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 7383)
@@ -17,5 +17,4 @@
    USE fldread        ! read input fields
    USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE sbcwave  ,ONLY :   cdn_wave !wave module
    !
    USE iom            ! I/O manager library
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7383)
@@ -18,10 +18,11 @@
    !!   sbc_cpl_snd     : send     fields to the atmosphere
    !!----------------------------------------------------------------------
-   USE dom_oce        ! ocean space and time domain
-   USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE sbc_ice        ! Surface boundary condition: ice fields
-   USE sbcapr         ! Stochastic param. : ???
-   USE sbcdcy         ! surface boundary condition: diurnal cycle
-   USE phycst         ! physical constants
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbc_ice         ! Surface boundary condition: ice fields
+   USE sbcapr
+   USE sbcdcy          ! surface boundary condition: diurnal cycle
+   USE sbcwave         ! surface boundary condition: waves
+   USE phycst          ! physical constants
 #if defined key_lim3
    USE ice            ! ice variables
@@ -106,5 +107,14 @@
    INTEGER, PARAMETER ::   jpr_e3t1st = 41   ! first T level thickness 
    INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jprcv      = 42   ! total number of fields received
+   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 
+   INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period
+   INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber
+   INTEGER, PARAMETER ::   jpr_wstrf  = 50   ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient
+   INTEGER, PARAMETER ::   jprcv      = 51   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -136,5 +146,9 @@
    INTEGER, PARAMETER ::   jps_e3t1st = 27   ! first level depth (vvl)
    INTEGER, PARAMETER ::   jps_fraqsr = 28   ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jpsnd      = 28   ! total number of fields sended
+   INTEGER, PARAMETER ::   jps_ficet  = 29   ! total ice fraction  
+   INTEGER, PARAMETER ::   jps_ocxw   = 30   ! currents on grid 1  
+   INTEGER, PARAMETER ::   jps_ocyw   = 31   ! currents on grid 2
+   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level 
+   INTEGER, PARAMETER ::   jpsnd      = 32   ! total number of fields sent 
 
    !                                  !!** namelist namsbc_cpl **
@@ -150,5 +164,9 @@
    !                                   ! Received from the atmosphere
    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
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                        
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp                           
+   ! 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_wstrf,sn_rcv_wdrag
    !                                   ! Other namelist parameters
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -163,5 +181,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
 
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo          ! OASIS info argument
+   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
+   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0) 
+
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
 
    !! Substitution
@@ -178,5 +199,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(3)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -189,4 +210,6 @@
       ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
       !
+      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
+
       sbc_cpl_alloc = MAXVAL( ierr )
       IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
@@ -214,8 +237,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      &
-         &                  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_cal   , sn_rcv_iceflx,   &
-         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask
+      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      & 
+         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
+         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , 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_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &
+         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp 
       !!---------------------------------------------------------------------
       !
@@ -258,8 +283,17 @@
          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,*)'      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,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
          WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
          WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
+         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 
          WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
          WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
@@ -267,4 +301,10 @@
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
+         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')' 
+         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%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 
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
@@ -305,4 +345,5 @@
       ! 
       ! Vectors: change of sign at north fold ONLY if on the local grid
+      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled
       IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
       
@@ -372,5 +413,6 @@
          srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
       ENDIF
-      !
+      ENDIF
+
       !                                                      ! ------------------------- !
       !                                                      !    freshwater budget      !   E-P
@@ -468,4 +510,10 @@
       !                                                      ! ------------------------- !
       srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
+
+      !                                                      ! ------------------------- ! 
+      !                                                      ! Mean Sea Level Pressure   ! 
+      !                                                      ! ------------------------- ! 
+      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE. 
+
       !                                                      ! ------------------------- !
       !                                                      !   topmelt and botmelt     !   
@@ -481,4 +529,48 @@
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
+      !                                                      ! ------------------------- !
+      !                                                      !      Wave breaking        !    
+      !                                                      ! ------------------------- ! 
+      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height
+      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN
+         srcv(jpr_hsig)%laction = .TRUE.
+         cpl_hsig = .TRUE.
+      ENDIF
+      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy
+      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN
+         srcv(jpr_phioc)%laction = .TRUE.
+         cpl_phioc = .TRUE.
+      ENDIF
+      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction
+      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN
+         srcv(jpr_sdrftx)%laction = .TRUE.
+         cpl_sdrftx = .TRUE.
+      ENDIF
+      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction
+      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN
+         srcv(jpr_sdrfty)%laction = .TRUE.
+         cpl_sdrfty = .TRUE.
+      ENDIF
+      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period
+      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN
+         srcv(jpr_wper)%laction = .TRUE.
+         cpl_wper = .TRUE.
+      ENDIF
+      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number
+      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wnum)%laction = .TRUE.
+         cpl_wnum = .TRUE.
+      ENDIF
+      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
+      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN
+         srcv(jpr_wdrag)%laction = .TRUE.
+         cpl_wdrag = .TRUE.
+      ENDIF
+      ! 
       !                                                      ! ------------------------------- !
       !                                                      !   OPA-SAS coupling - rcv by opa !   
@@ -635,4 +727,5 @@
       !                                                      ! ------------------------- !
       ssnd(jps_fice)%clname = 'OIceFrc'
+      ssnd(jps_ficet)%clname = 'OIceFrcT' 
       ssnd(jps_hice)%clname = 'OIceTck'
       ssnd(jps_hsnw)%clname = 'OSnwTck'
@@ -643,4 +736,6 @@
       ENDIF
       
+      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
+
       SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
       CASE( 'none'         )       ! nothing to do
@@ -663,4 +758,6 @@
       ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
       ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
+      ssnd(jps_ocxw)%clname = 'O_OCurxw' 
+      ssnd(jps_ocyw)%clname = 'O_OCuryw' 
       !
       ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
@@ -683,8 +780,29 @@
       END SELECT
 
+      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold 
+        
+      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN 
+         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V' 
+      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN 
+         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' ) 
+      ENDIF 
+      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1. 
+      SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
+         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE. 
+         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE. 
+         CASE( 'weighted oce and ice' )   !   nothing to do 
+         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE. 
+         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' ) 
+      END SELECT 
+
       !                                                      ! ------------------------- !
       !                                                      !          CO2 flux         !
       !                                                      ! ------------------------- !
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
+
+      !                                                      ! ------------------------- ! 
+      !                                                      !     Sea surface height    ! 
+      !                                                      ! ------------------------- ! 
+      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE. 
 
       !                                                      ! ------------------------------- !
@@ -781,5 +899,5 @@
       IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
          &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
-      ncpl_qsr_freq = 86400 / ncpl_qsr_freq
+      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
       CALL wrk_dealloc( jpi,jpj,   zacs, zaos )
@@ -835,8 +953,11 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
-      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation 
-      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
-
+      USE zdf_oce,  ONLY : ln_zdfqiao
+
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
+      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
+      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
       !!
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
@@ -990,4 +1111,62 @@
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
 #endif
+      ! 
+      !                                                      ! ========================= ! 
+      !                                                      ! Mean Sea Level Pressure   !   (taum) 
+      !                                                      ! ========================= ! 
+      ! 
+      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH 
+          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields 
+
+          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization 
+          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer) 
+          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                         !atmospheric pressure 
+    
+          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible) 
+      END IF 
+      !
+      IF( ln_sdw ) THEN  ! Stokes Drift correction activated
+      !                                                      ! ========================= ! 
+      !                                                      !       Stokes drift u      !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !       Stokes drift v      !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !      Wave mean period     !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !  Significant wave height  !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1)
+      !
+      !                                                      ! ========================= ! 
+      !                                                      !    Vertical mixing Qiao   !
+      !                                                      ! ========================= ! 
+         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)
+
+         ! 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 ) THEN
+            CALL sbc_stokes()
+            IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao()
+         ENDIF
+         IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao()
+      ENDIF
+      !                                                      ! ========================= ! 
+      !                                                      ! Stress adsorbed by waves  !
+      !                                                      ! ========================= ! 
+      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)
 
       !  Fields received by SAS when OASIS coupling
@@ -2063,4 +2242,148 @@
       ENDIF
       !
+      !                                                      ! ------------------------- ! 
+      !                                                      !  Surface current to waves ! 
+      !                                                      ! ------------------------- ! 
+      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN 
+          !     
+          !                                                  j+1  j     -----V---F 
+          ! surface velocity always sent from T point                    !       | 
+          !                                                       j      |   T   U 
+          !                                                              |       | 
+          !                                                   j   j-1   -I-------| 
+          !                                               (for I)        |       | 
+          !                                                             i-1  i   i 
+          !                                                              i      i+1 (for I) 
+          SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
+          CASE( 'oce only'             )      ! C-grid ==> T 
+             DO jj = 2, jpjm1 
+                DO ji = fs_2, fs_jpim1   ! vector opt. 
+                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
+                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) )  
+                END DO 
+             END DO 
+          CASE( 'weighted oce and ice' )    
+             SELECT CASE ( cp_ice_msh ) 
+             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = fs_2, fs_jpim1   ! vector opt. 
+                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   
+                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj) 
+                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
+                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
+                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
+                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
+                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
+                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
+                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
+                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             END SELECT 
+             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. ) 
+          CASE( 'mixed oce-ice'        ) 
+             SELECT CASE ( cp_ice_msh ) 
+             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = fs_2, fs_jpim1   ! vector opt. 
+                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   & 
+                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
+                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
+                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &    
+                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
+                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &  
+                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
+                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T 
+                DO jj = 2, jpjm1 
+                   DO ji = 2, jpim1   ! NO vector opt. 
+                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &    
+                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
+                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &  
+                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
+                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
+                   END DO 
+                END DO 
+             END SELECT 
+          END SELECT 
+         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 
+         ! 
+         ! 
+         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components 
+         !                                                                        ! Ocean component 
+            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component  
+            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component  
+            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components  
+            zoty1(:,:) = ztmp2(:,:)  
+            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component 
+               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component  
+               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component  
+               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components  
+               zity1(:,:) = ztmp2(:,:) 
+            ENDIF 
+         ENDIF 
+         ! 
+!         ! spherical coordinates to cartesian -> 2 components to 3 components 
+!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN 
+!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents 
+!            ztmp2(:,:) = zoty1(:,:) 
+!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 ) 
+!            ! 
+!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities 
+!               ztmp1(:,:) = zitx1(:,:) 
+!               ztmp1(:,:) = zity1(:,:) 
+!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 ) 
+!            ENDIF 
+!         ENDIF 
+         ! 
+         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid 
+         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid 
+         !  
+      ENDIF 
+      ! 
+      IF( ssnd(jps_ficet)%laction ) THEN 
+         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 
+      END IF 
+      !                                                      ! ------------------------- ! 
+      !                                                      !   Water levels to waves   ! 
+      !                                                      ! ------------------------- ! 
+      IF( ssnd(jps_wlev)%laction ) THEN 
+         IF( ln_apr_dyn ) THEN  
+            IF( kt /= nit000 ) THEN  
+               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
+            ELSE  
+               ztmp1(:,:) = sshb(:,:)  
+            ENDIF  
+         ELSE  
+            ztmp1(:,:) = sshn(:,:)  
+         ENDIF  
+         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
+      END IF 
       !
       !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7383)
@@ -89,9 +89,9 @@
       NAMELIST/namsbc/ nn_fsbc  , ln_ana   , ln_flx, ln_blk_clio, ln_blk_core, ln_blk_mfs,   &
          &             ln_cpl   , ln_mixcpl, nn_components      , nn_limflx  ,               &
-         &             ln_traqsr, ln_dm2dc ,                                                 &  
+         &             ln_traqsr, ln_dm2dc ,                                                 &
          &             nn_ice   , nn_ice_embd,                                               &
          &             ln_rnf   , ln_ssr   , ln_isf   , nn_fwb    , ln_apr_dyn,              &
-         &             ln_wave  ,                                                            &
-         &             nn_lsm   
+         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauoc  , ln_stcor  ,              &
+         &             nn_lsm
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -153,5 +153,9 @@
          WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea
          WRITE(numout,*) '              nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
-         WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave  
+         WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave
+         WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw        = ', ln_sdw
+         WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc      = ', ln_tauoc
+         WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor      = ', ln_stcor
+         WRITE(numout,*) '                 neutral drag coefficient (CORE, MFS)    ln_cdgw       = ', ln_cdgw
       ENDIF
       !
@@ -220,4 +224,22 @@
          &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
       
+      IF ( ln_wave ) THEN
+      !Activated wave module but neither drag nor stokes drift activated
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN
+            CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')
+      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
+         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
+             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+         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_tauoc .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_tauoc=T) ',      &  
+         &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
+      ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
       ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
@@ -357,5 +379,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                                 ! Wave stress subctracted
+            utau(:,:) = utau(:,:)*tauoc_wave(:,:)
+            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
+            taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+      !
+            SELECT CASE( nsbc )
+            CASE(  0,1,2,3,5,-1 )  ;
+                IF(lwp .AND. kt == nit000 ) 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/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7383)
@@ -4,38 +4,56 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
-   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   sbc_wave      : read drag coefficient from wave model in netcdf files 
+   !! History :  3.3  !   2011-09  (M. Adani)  Original code: Drag Coefficient 
+   !!         :  3.4  !   2012-10  (M. Adani)  Stokes Drift 
+   !!            3.6  !   2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_wave      : wave data from wave model in netcdf files 
    !!----------------------------------------------------------------------
    USE oce            ! 
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
+   USE sbc_oce	      ! Surface boundary condition: ocean fields
    USE bdy_oce        !
    USE domvvl         !
-   !
    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 wrk_nemo       !
+   USE phycst         ! physical constants 
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
+   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl
+   PUBLIC   sbc_wave    ! routine called in sbcmod
    
-   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 
+   ! 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_wnum=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE.
+   LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE.
+
+   INTEGER ::   jpfld                ! number of files to read for stokes drift
+   INTEGER ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point
+   INTEGER ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point
+   INTEGER ::   jp_swh               ! index of significant wave hight      (m)      at T-point
+   INTEGER ::   jp_wmp               ! 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(:,:)       :: swh,wmp, wnum
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt
 
    !! * Substitutions
@@ -48,28 +66,140 @@
 CONTAINS
 
+   SUBROUTINE sbc_stokes( )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_stokes  ***
+      !!
+      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al.,
+      !!                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  
+      !!---------------------------------------------------------------------
+      INTEGER                ::   jj,ji,jk 
+      REAL(wp)                       ::  ztransp, zfac, zsp0, zk, zus, zvs
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace
+      !!---------------------------------------------------------------------
+      !
+
+      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( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
+               ! Wavenumber scale
+               zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
+               ! Depth attenuation
+               zfac = EXP(-2.0_wp*zk*gdept_n(ji,jj,jk))/(1.0_wp+8.0_wp*zk*gdept_n(ji,jj,jk))
+               !
+               usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk)
+               vsd3dt(ji,jj,jk) = zfac * zvsd2dt(ji,jj) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO 
+      ! Into the U and V Grid
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   &
+                               &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) )
+               vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   &
+                               &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
+      CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
+      !
+      DO jk = 1, jpkm1               ! Horizontal divergence
+         DO jj = 2, jpj
+            DO ji = fs_2, jpi
+               ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     &
+                  &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     &
+                  &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     &
+                  &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
+            END DO
+         END DO
+      END DO
+      !
+      IF( .NOT. AGRIF_Root() ) THEN
+         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east
+         IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west
+         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north
+         IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south
+      ENDIF
+      !
+      CALL lbc_lnk( ze3hdiv, 'T', 1. )
+      !
+      DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence
+         wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - e3t_n(:,:,jk) * ze3hdiv(:,:,jk)
+      END DO
+#if defined key_bdy
+      IF( lk_bdy ) THEN
+         DO jk = 1, jpkm1
+            wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
+         END DO
+      ENDIF
+#endif
+      CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv )
+      !
+   END SUBROUTINE sbc_stokes
+
+   SUBROUTINE sbc_qiao
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_qiao  ***
+      !!
+      !! ** Purpose :   Qiao formulation for wave enhanced turbulence
+      !!                2010 (DOI: 10.1007/s10236-010-0326) 
+      !!
+      !! ** Method  : - 
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER :: jj, ji
+
+      ! Calculate the module of the stokes drift on T grid
+      !-------------------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj) * zusd2dt(ji,jj) + zvsd2dt(ji,jj) * zvsd2dt(ji,jj) )
+         END DO
+      END DO
+      !
+   END SUBROUTINE sbc_qiao
+
    SUBROUTINE sbc_wave( kt )
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE sbc_apr  ***
-      !!
-      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
+      !!                     ***  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 monochromatic
-      !! ** action  :   
-      !!---------------------------------------------------------------------
-      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
+      !!              - 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
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
+      INTEGER                ::   ifpr
+      INTEGER                ::   ios      ! Local integer output status for namelist 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
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace
-      !!
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn, ln_cdgw , ln_sdw
+      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_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, sn_tauoc
       !!---------------------------------------------------------------------
       !
@@ -80,5 +210,5 @@
          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
-         !
+         
          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
@@ -86,99 +216,123 @@
          IF(lwm) WRITE ( numond, namsbc_wave )
          !
-         IF(lwp) THEN               ! Control print
-            WRITE(numout,*) '        Namelist namsbc_wave : surface wave setting' 
-            WRITE(numout,*) '           wave drag coefficient                      ln_cdgw  = ', ln_cdgw  
-            WRITE(numout,*) '           wave stokes drift                          ln_sdw   = ', ln_sdw
-         ENDIF
-         !
-         IF( .NOT.( ln_cdgw .OR. ln_sdw ) )    &
-            &  CALL ctl_warn( 'ln_sbcwave=T but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
-         IF( ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )   &       
-            &  CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
-         !
          IF( ln_cdgw ) THEN
-            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: 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', 'Wave module ', 'namsbc_wave' )
+            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: 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', 'Wave module ', 'namsbc_wave' )
+            ENDIF
             ALLOCATE( cdn_wave(jpi,jpj) )
-            cdn_wave(:,:) = 0.0
-         ENDIF
+         ENDIF
+
+         IF( ln_tauoc ) THEN
+            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: unable to allocate sf_wave 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', 'Wave module', 'namsbc_wave' )
+            ENDIF
+            ALLOCATE( tauoc_wave(jpi,jpj) )
+         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
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: 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', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj) , vsd2d(jpi,jpj) , uwavenum(jpi,jpj) , vwavenum(jpi,jpj) )
+            ! Find out how many fields have to be read from file if not coupled
+            jpfld=0
+            jp_usd=0; jp_vsd=0; jp_swh=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_swh=jpfld
+            ENDIF
+            IF( .NOT. cpl_wper ) THEN
+               jpfld=jpfld+1
+               jp_wmp=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_swh > 0 ) slf_i(jp_swh) = sn_swh
+               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: 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', 'Wave module ', 'namsbc_wave' )
+            ENDIF
             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
+            ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) )
+            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
+            ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) )
+            usd3d(:,:,:) = 0._wp
+            vsd3d(:,:,:) = 0._wp
             wsd3d(:,:,:) = 0._wp
-         ENDIF
-      ENDIF
-      !
-      IF( ln_cdgw ) THEN               !==  Neutral drag coefficient  ==!
+            IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==!
+               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: 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),tsd2d(jpi,jpj) )
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==!
          CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing
          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,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 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 wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
+         ! Read from file only if the field is not coupled
+         IF( jpfld > 0 ) THEN
+            CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
+            IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height
+            IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+            IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
+            IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
+         ENDIF
          !
-         CALL fld_read( kt, nn_fsbc, sf_sd )    !* read drag coefficient from external forcing
+         ! Read also wave number if needed, so that it is available in coupling routines
+         IF( ln_zdfqiao .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)
+         ENDIF
+           
+         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014
+         !(DOI: 10.1175/JPO-D-14-0020.1)==! 
          !
-         DO jk = 1, jpkm1                       !* distribute it on the vertical
-            zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) )
-            zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) )
-         END DO
-         DO jk = 1, jpkm1                       !* interpolate the stokes drift from t-point to u- and v-points
-            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
-         CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
-         CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
-         !
-         DO jk = 1, jpkm1                       !* e3t * Horizontal divergence  ==!
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     &
-                     &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     &
-                     &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     &
-                     &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
-               END DO  
-            END DO  
-            IF( .NOT. AGRIF_Root() ) THEN
-               IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,jk) = 0._wp      ! east
-               IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,jk) = 0._wp      ! west
-               IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,jk) = 0._wp      ! north
-               IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,jk) = 0._wp      ! south
-            ENDIF
-         END DO
-         CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
-         !
-         DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence
-            wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk)
-         END DO
-#if defined key_bdy
-         IF( lk_bdy ) THEN
-            DO jk = 1, jpkm1
-               wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
-            END DO
-         ENDIF
-#endif
-         CALL wrk_dealloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv )
-         ! 
+         ! Calculate only if no necessary fields are coupled, if not calculate later after coupling
+         IF( jpfld == 4 ) THEN
+            CALL sbc_stokes()
+            IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN
+               CALL sbc_qiao()
+            ENDIF
+         ENDIF
       ENDIF
       !
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7383)
@@ -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
    !!----------------------------------------------------------------------
 
@@ -35,6 +36,7 @@
    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 
 
    IMPLICIT NONE
@@ -96,4 +98,8 @@
       !
       !                                          ! set time step
+      zun(:,:,:) = 0.0
+      zvn(:,:,:) = 0.0
+      zwn(:,:,:) = 0.0
+      !    
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
          r2dt = rdt                                 ! = rdt (restarting with Euler time stepping)
@@ -103,9 +109,20 @@
       !
       !                                         !==  effective transport  ==!
-      DO jk = 1, jpkm1
-         zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
-         zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
-         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
-      END DO
+      IF( ln_wave .AND. ln_sdw )  THEN
+         DO jk = 1, jpkm1
+            zun(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) *      &
+                        &  ( un(:,:,jk) + usd3d(:,:,jk) )                       ! eulerian transport + Stokes Drift
+            zvn(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) *      &
+                        &  ( vn(:,:,jk) + vsd3d(:,:,jk) )
+            zwn(:,:,jk) = e1e2t(:,:) *                    &
+                        &  ( wn(:,:,jk) + wsd3d(:,:,jk) )
+         END DO
+      ELSE
+         DO jk = 1, jpkm1
+            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)               ! eulerian transport only
+            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+            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/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7383)
@@ -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/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7383)
@@ -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/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7383)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7383)
@@ -0,0 +1,128 @@
+MODULE zdfqiao
+   !!======================================================================
+   !!                       ***  MODULE  zdfqiao  ***
+   !! Qiao module      : vertical mixing enhancement due to surface waves 
+   !!======================================================================
+   !! History :  3.6  !  2014-10  (E. Clementi)  Original code
+   !!----------------------------------------------------------------------
+   !!   zdf_qiao        : compute Qiao parameters
+   !!----------------------------------------------------------------------
+
+   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
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)  
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC zdf_qiao    ! routine called in step
+
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! 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 :: jj, ji, jk   ! dummy loop indices
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         IF( .NOT. ( ln_wave .AND. ln_sdw ) )   &
+            &   CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave   &
+            &                    and ln_sdw have to be activated')
+         IF( zdf_qiao_alloc() /= 0 )   &
+            &   CALL ctl_stop( 'STOP', 'zdf_qiao : unable to allocate arrays' )
+      ENDIF
+
+      !
+      ! 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
+      !---------------------------------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *             &
+            &                  EXP(3.0 * wnum(ji,jj) *                                  &                     
+            &                  (-MIN( gdepw_n(ji  ,jj  ,jk), gdepw_n(ji+1,jj  ,jk),     &
+            &                         gdepw_n(ji  ,jj+1,jk), gdepw_n(ji+1,jj+1,jk))))   &
+            &                          * wmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( qbv, 'W', 1. )   ! Lateral boundary conditions
+         
+      !
+      ! Interpolate Qiao parameter qbv into the grid_U and grid_V
+      !----------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               qbvu(ji,jj,jk) = 0.5 * wumask(ji,jj,jk)  *              &  
+            &                  ( qbv(ji,jj,jk) + qbv(ji+1,jj  ,jk) )
+               qbvv(ji,jj,jk) = 0.5 * wvmask(ji,jj,jk)  *              &
+            &                  ( qbv(ji,jj,jk) + qbv(ji  ,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      ! 
+      CALL lbc_lnk( qbvu, 'U', 1. ) ; CALL lbc_lnk( qbvv, 'V', 1. )   ! Lateral boundary conditions
+
+      ! Enhance vertical mixing coeff.         
+      !-------------------------------
+      !
+      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
+      !
+   END SUBROUTINE zdf_qiao
+
+   INTEGER FUNCTION zdf_qiao_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_qiao_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( qbv(jpi,jpj,jpk), qbvu(jpi,jpj,jpk), qbvv(jpi,jpj,jpk),   &
+         &      STAT = zdf_qiao_alloc )
+      !
+      IF( lk_mpp             )  CALL mpp_sum ( zdf_qiao_alloc )
+      IF( zdf_qiao_alloc > 0 )  CALL ctl_warn('zdf_qiao_alloc: allocation of arrays failed.')
+      !
+   END FUNCTION zdf_qiao_alloc
+      
+   !!======================================================================
+END MODULE zdfqiao
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7383)
@@ -205,5 +205,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               IF( fsdept(ji,jj,jk) < ekm_dep(ji,jj) ) THEN
+               IF( gdepw_n(ji,jj,jk) < ekm_dep(ji,jj) ) THEN
                   avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix )
                   avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix )
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7383)
@@ -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
@@ -73,5 +74,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)
@@ -128,8 +129,10 @@
                          CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
       !                                               ! Vertical eddy viscosity and diffusivity coefficients
-      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
-      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( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
+      IF( lk_zdfric  )   CALL zdf_ric ( kstp )             ! Richardson number dependent Kz
+      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 )   CALL zdf_qiao( kstp )             ! Qiao vertical mixing 
+      !
+      IF( lk_zdfcst  ) THEN                                ! Constant Kz (reset avt, avm[uv] to the background value)
          avt (:,:,:) = rn_avt0 * wmask (:,:,:)
          avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
@@ -207,4 +210,6 @@
                          CALL dyn_adv       ( kstp )  ! advection (vector or flux form)
                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
+      IF( ln_wave .AND. ln_sdw .AND. ln_stcor)           &
+               &         CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
Index: branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7382)
+++ branches/2016/dev_INGV_METO_merge_2016/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7383)
@@ -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 step_diu        ! Time stepping for diurnal sst
