Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7471)
@@ -264,7 +264,9 @@
                            !     =1 global mean of e-p-r set to zero at each time step
                            !     =2 annual global mean of e-p-r set to zero
-   ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave)
-   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
-   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
+   ln_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)
@@ -363,4 +365,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   =       'none'                 ,    'no'    ,     ''      ,         ''           ,   '' 
 ! receive
    sn_rcv_w10m   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
@@ -374,4 +379,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
@@ -921,4 +936,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)
 /
 !-----------------------------------------------------------------------
@@ -1270,10 +1286,13 @@
 !              !  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_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_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
-!
-   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+   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
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/cfg.txt	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/cfg.txt	(revision 7471)
@@ -11,2 +11,3 @@
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
+GYRE_LONG OPA_SRC
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 7471)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 7471)
@@ -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('dyn_stcor')
+      !
+   END SUBROUTINE dyn_stcor
+
+END MODULE dynstcor
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 7471)
@@ -65,5 +65,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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7471)
@@ -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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7471)
@@ -737,5 +737,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(:,:)
@@ -783,5 +783,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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7471)
@@ -23,4 +23,5 @@
    USE sbcapr
    USE sbcdcy          ! surface boundary condition: diurnal cycle
+   USE sbcwave         ! surface boundary condition: waves
    USE phycst          ! physical constants
 #if defined key_lim3
@@ -105,5 +106,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
@@ -135,5 +145,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 **
@@ -149,5 +163,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
@@ -161,4 +179,7 @@
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+   
+   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
@@ -179,5 +200,5 @@
       !!             ***  FUNCTION sbc_cpl_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(3)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -190,4 +211,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 )
@@ -216,8 +239,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
       !!---------------------------------------------------------------------
       !
@@ -260,8 +285,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 
@@ -269,4 +303,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
@@ -307,4 +347,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.
       
@@ -374,4 +415,5 @@
          srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
       ENDIF
+      ENDIF
        
       !                                                      ! ------------------------- !
@@ -470,4 +512,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     !   
@@ -483,4 +531,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 !   
@@ -637,4 +729,5 @@
       !                                                      ! ------------------------- !
       ssnd(jps_fice)%clname = 'OIceFrc'
+      ssnd(jps_ficet)%clname = 'OIceFrcT'
       ssnd(jps_hice)%clname = 'OIceTck'
       ssnd(jps_hsnw)%clname = 'OSnwTck'
@@ -645,4 +738,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
@@ -665,4 +760,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
@@ -685,8 +782,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. 
 
       !                                                      ! ------------------------------- !
@@ -783,5 +901,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 )
@@ -837,4 +955,6 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
+      USE zdf_oce,  ONLY : ln_zdfqiao
+
       INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
       INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
@@ -996,4 +1116,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
@@ -2019,4 +2197,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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7471)
@@ -89,5 +89,5 @@
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
-         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
+         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx , nn_components, ln_cpl
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -132,4 +132,9 @@
          WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
          WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
+         WRITE(numout,*) '              wave physics                               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
          WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
          WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
@@ -216,15 +221,20 @@
       
       IF ( ln_wave ) THEN
-      !Activated wave module but neither drag nor stokes drift activated
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
-            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
-      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
+         !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  )                                         & 
-         &   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) ')
+         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)
@@ -382,4 +392,16 @@
       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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7471)
@@ -4,37 +4,61 @@
    !! Wave module 
    !!======================================================================
-   !! History :  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 oce
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
-   USE domvvl
-
-   
-   !!----------------------------------------------------------------------
-   !!   sbc_wave       : read drag coefficient from wave model in netcdf files 
-   !!----------------------------------------------------------------------
+   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient 
+   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift 
+   !!            3.6  !   2014-09  (Clementi E, Oddo P)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 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 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),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
+   ! 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
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -44,39 +68,143 @@
 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*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(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_e12t(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) - fse3t_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  :   
-      !!               
-      !!---------------------------------------------------------------------
-      USE oce,  ONLY : un,vn,hdivn,rotn
-      USE divcur
-      USE wrk_nemo
-#if defined key_bdy
-      USE bdy_oce, ONLY : bdytmask
-#endif
-      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
-      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
-      REAL                                          ::  z2dt,z1_2dt
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
+      !!              - Read wave number in netcdf files 
+      !!              - Compute 3d stokes drift using Breivik et al.,2014
+      !!                formulation
+      !! ** action  
+      !!---------------------------------------------------------------------
+      USE zdf_oce,  ONLY : ln_zdfqiao
+
+      IMPLICIT NONE
+
+      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
+      !
+      INTEGER                ::   ierror   ! return error code
+      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
-      !!---------------------------------------------------------------------
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
-      !!---------------------------------------------------------------------
-
-      !!----------------------------------------------------------------------
-      !
+      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
+      !!---------------------------------------------------------------------
       !
       !                                         ! -------------------- !
@@ -92,101 +220,126 @@
          IF(lwm) WRITE ( numond, namsbc_wave )
          !
-
-         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( ln_cdgw ) THEN
+            IF( .NOT. cpl_wdrag ) THEN
+               ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: 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
-         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) )
+         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
+            ! 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) )
-            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
-            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
-         ENDIF
-      ENDIF
+            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
+            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_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  ==! 
          !
+         ! 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
          !
-      IF ( ln_cdgw ) THEN
-         CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing
-         cdn_wave(:,:) = sf_cd(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
-         !-------------------------------------------------
-
-         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) )
-            END DO
-         END DO
-
-          !Computation of the 3d Stokes Drift
-          DO jk = 1, jpk
-             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))))
-                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 wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-          
-          udummy(:,:,:)=un(:,:,:)
-          vdummy(:,:,:)=vn(:,:,:)
-          hdivdummy(:,:,:)=hdivn(:,:,:)
-          rotdummy(:,:,:)=rotn(:,:,:)
-          un(:,:,:)=usd3d(:,:,:)
-          vn(:,:,:)=vsd3d(:,:,:)
-          CALL div_cur(kt)
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-          z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
-
-          z1_2dt = 1.e0 / z2dt
-          DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-             ! - 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
-#if defined key_bdy
-             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
-#endif
-          END DO
-          hdivn(:,:,:)=hdivdummy(:,:,:)
-          rotn(:,:,:)=rotdummy(:,:,:)
-          vn(:,:,:)=vdummy(:,:,:)
-          un(:,:,:)=udummy(:,:,:)
-          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-      ENDIF
+         ! 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)==! 
+         !
+         ! 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
+      !
    END SUBROUTINE sbc_wave
       
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7471)
@@ -6,4 +6,5 @@
    !! History :  2.0  !  2005-11  (G. Madec)  Original code
    !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
+   !!            3.6  !  2015-06  (E. Clementi) Addition of Stokes drift in case of wave coupling
    !!            4.0  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation
    !!----------------------------------------------------------------------
@@ -34,6 +35,7 @@
    USE timing          ! Timing
    USE sbc_oce
-   USE diaptr          ! Poleward heat transport 
-
+   USE sbcwave        ! wave module 
+   USE sbc_oce        ! surface boundary condition: ocean 
+   USE diaptr         ! Poleward heat transport
 
    IMPLICIT NONE
@@ -85,4 +87,8 @@
       CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
       !                                          ! set time step
+      zun(:,:,:) = 0.0 
+      zvn(:,:,:) = 0.0 
+      zwn(:,:,:) = 0.0 
+      !     
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
          r2dtra(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
@@ -94,9 +100,20 @@
       !
       !                                               !==  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) = e1e2t(:,:) *                    & 
+                        &  ( 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) = e1e2t(:,:)                 * wn(:,:,jk) 
+         END DO 
+      ENDIF 
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7471)
@@ -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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7471)
@@ -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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7471)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7471)
@@ -0,0 +1,129 @@
+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 "domzgr_substitute.h90"
+#  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( fsdepw(ji  ,jj  ,jk), fsdepw(ji+1,jj  ,jk),     &
+            &                         fsdepw(ji  ,jj+1,jk), fsdepw(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/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7471)
@@ -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
    !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state
    !!----------------------------------------------------------------------
@@ -72,5 +73,5 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   ji,jj,jk ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
@@ -132,4 +133,6 @@
       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_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
       IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
@@ -209,4 +212,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
           IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7471)
@@ -28,4 +28,5 @@
    USE sbctide          ! Tide initialisation
    USE sbcapr           ! surface boundary condition: ssh_ib required by bdydta 
+   USE sbcwave          ! Wave intialisation
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
@@ -51,4 +52,5 @@
    USE dynspg_oce       ! surface pressure gradient        (dyn_spg routine)
    USE dynspg           ! surface pressure gradient        (dyn_spg routine)
+   USE dynstcor         ! simp. form of Stokes-Coriolis
    USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine)
 
@@ -84,4 +86,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)
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/SETTE/sette.sh
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/SETTE/sette.sh	(revision 7470)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/SETTE/sette.sh	(revision 7471)
@@ -88,6 +88,7 @@
 #
 # Compiler among those in NEMOGCM/ARCH
-COMPILER=X64_ADA
-export BATCH_COMMAND_PAR="llsubmit"
+module load cray-netcdf-hdf5parallel 
+COMPILER=XC40_METO 
+export BATCH_COMMAND_PAR="qsub"
 export BATCH_COMMAND_SEQ=$BATCH_COMMAND_PAR
 export INTERACT_FLAG="no"
