Index: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7168)
+++ branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7193)
@@ -41,5 +41,4 @@
    USE sbcapr          ! surface boundary condition: atmospheric pressure
    USE dynadv, ONLY: ln_dynadv_vec
-   USE sbcwave, ONLY: usd2d, vsd2d
 #if defined key_agrif
    USE agrif_opa_interp ! agrif
Index: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7168)
+++ branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 7193)
@@ -39,5 +39,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-   USE sbcwave,  ONLY:  zusd2dt, zvsd2dt,wsd3d
+   USE sbcwave,  ONLY:  usd3dt, vsd3dt,wsd3d
 
    IMPLICIT NONE
@@ -162,5 +162,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::  sshnu, sshnv, dsshnu, dsshnv
+      REAL(wp) ::  dsshnu, dsshnv
       !!----------------------------------------------------------------------
       !
@@ -210,37 +210,20 @@
       ENDIF
 !
-!     In case ln_wave and ln_sdw the surface vertical velocity is modified
+!     In case ln_wave and ln_sdw, the surface vertical velocity is modified
 !     accounting for Sokes Drift velocity
+!
       IF ( ln_wave .AND. ln_sdw )  THEN
-       ALLOCATE(sshnu(jpi,jpj),sshnv(jpi,jpj),dsshnu(jpi,jpj),dsshnv(jpi,jpj))
-       sshnu (:,:) = 0._wp
-       sshnv (:,:) = 0._wp
-       dsshnu(:,:) = 0._wp
-       dsshnv(:,:) = 0._wp
-       ! sshn interpolated on U-V grid
-       !---------------------------------
-       DO jj = 1 , jpjm1
-         DO ji = 1 , jpim1
-           sshnu(ji,jj) =  0.5 * ( 2. - umask(ji,jj,1) ) *               &
-                        &        ( sshn(ji  ,jj) * tmask(ji  ,jj,1)      &
-                        &        + sshn(ji+1,jj) * tmask(ji+1,jj,1) )
-           sshnv(ji,jj) =  0.5 * ( 2. - vmask(ji,jj,1) ) *               &
-                        &        ( sshn(ji,jj  ) * tmask(ji,jj  ,1)      &
-                        &        + sshn(ji,jj+1) * tmask(ji,jj+1,1) )
-         ENDDO
+         ! Compute d(ssh)/dx  and d(ssh)/dy   
+         ! Compute the surface vertical velocity accounting for the Stokes Drift 
+         !--------------------------------- 
+         DO jj = 2 , jpjm1 
+            DO ji = 2 , jpim1 
+               dsshnu = ( ssha(ji+1,jj) - ssha(ji,jj) ) / e1u(ji,jj) 
+               dsshnv = ( ssha(ji,jj+1) - ssha(ji,jj) ) / e2v(ji,jj) 
+               wn(ji,jj,1) = wn(ji,jj,1) +( usd3dt(ji,jj,1) * dsshnu     & 
+                     &     + vsd3dt(ji,jj,1) * dsshnv                    & 
+                     &     - ( wsd3d (ji,jj,1) )) * tmask(ji,jj,1) 
+            ENDDO 
        ENDDO
-       ! Compute d(ssh)/dx  and d(ssh)/dy  
-       !---------------------------------
-       DO jj = 1 , jpjm1
-         DO ji = 1 , jpim1
-           dsshnu(ji,jj) = ( sshnu(ji+1,jj) - sshnu(ji,jj) ) / e1u(ji,jj)
-           dsshnv(ji,jj) = ( sshnv(ji,jj+1) - sshnv(ji,jj) ) / e2v(ji,jj)
-         ENDDO
-       ENDDO
-       ! Compute the surface vertical velocity accounting for the Stokes Drift 
-       !---------------------------------------------------------------------
-       wn(:,:,1) = wn(:,:,1) + zusd2dt(:,:) * dsshnu(:,:)     &
-                 &           + zvsd2dt(:,:) * dsshnv(:,:)      &
-                 &           - ( wsd3d (:,:,1) ) * tmask(:,:,1)
       ENDIF
       !
Index: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7168)
+++ branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7193)
@@ -53,4 +53,5 @@
    USE sbcwave          ! Wave module
    USE bdy_par          ! Require lk_bdy
+   USE zdf_oce,  ONLY : ln_zdfqiao
 
    IMPLICIT NONE
@@ -88,5 +89,6 @@
          &             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   ,   &
-         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx  , nn_components, ln_cpl 
+         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx  , nn_components, ln_cpl ,   &
+         &             ln_zdfqiao
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
Index: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7168)
+++ branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7193)
@@ -56,4 +56,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: usd2d, vsd2d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt
 
@@ -94,5 +95,5 @@
             ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp))
             ! Stokes surface speed
-            zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 + sf_sd(jp_vsd)%fnow(ji,jj,1)**2)
+            zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
             ! Wavenumber scale
             zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
@@ -101,18 +102,20 @@
             END DO
          END DO
-
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               usd3dt(ji,jj,jk) = zfac(ji,jj) * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
+               vsd3dt(ji,jj,jk) = zfac(ji,jj) * 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, jpim1
-               ! Into the U and V Grid 
-               zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
-                  &                                + zfac(ji+1,jj) * tmask(ji+1,jj,1) )
-               zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) &
-                  &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) )
-               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) &
-                  &                                         + zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) )
-               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) &
-                  &                                         + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) )
-               usd3d(ji,jj,jk) = usd2d(ji,jj) * zus
-               vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs
+               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
@@ -170,6 +173,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            tsd2d(ji,jj) = ((zusd2dt(ji,jj) * tmask(ji,jj,1))**2.0  + &
-            &               (zvsd2dt(ji,jj) * tmask(ji,jj,1))**2.0)**0.5
+            tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
          END DO
       END DO
@@ -287,4 +289,5 @@
             ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) )
             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
+            ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) )
             ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
             ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) )
@@ -292,4 +295,5 @@
             vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ; 
             wsd3d(:,:,:) = 0._wp   ;
+            usd3dt(:,:,:) = 0._wp  ;   vsd3dt(:,:,:) = 0._wp   ;
             swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ;
             IF ( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==!
Index: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7168)
+++ branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7193)
@@ -31,4 +31,6 @@
    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: QBv, QBvu, QBvv
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -82,6 +84,6 @@
                   QBv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *       &
                &              exp(3.0 * wnum(ji,jj) *                                &                     
-               &              (-MIN( gdept_0(ji  ,jj  ,jk) , gdept_0(ji+1,jj  ,jk),  & 
-               &               gdept_0(ji  ,jj+1,jk),gdept_0(ji+1,jj+1,jk))))
+               &              (-MIN( fsdept(ji  ,jj  ,jk) , fsdept(ji+1,jj  ,jk),    &
+               &                     fsdept(ji  ,jj+1,jk) , fsdept(ji+1,jj+1,jk))))
                END DO
             END DO
@@ -98,11 +100,10 @@
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  QBvu(ji,jj,jk)=0.5 * ( 2. - umask(ji,jj,jk) ) *                 &
-               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                      & 
-               &           + QBv(ji+1,jj,jk) * tmask(ji+1,jj,jk) )               
-                  QBvv(ji,jj,jk)=0.5 * ( 2. - vmask(ji,jj,jk) ) *                 &
-               &           ( QBv(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                  QBvu(ji,jj,jk) = 0.5 *  umask(ji,jj,jk)  *               &
+               &           ( QBv(ji  ,jj,jk) * tmask(ji  ,jj,jk)           &
+               &           + QBv(ji+1,jj,jk) * tmask(ji+1,jj,jk) )
+                  QBvv(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk)  *               &
+               &           ( QBv(ji,jj  ,jk) * tmask(ji,jj  ,jk)           &
                &           + QBv(ji,jj+1,jk) * tmask(ji,jj+1,jk) )
-
                END DO
             END DO
