Index: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7481)
+++ branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7600)
@@ -41,4 +41,5 @@
    USE zdfddm          ! vertical  physics: double diffusion
    USE diahth          ! thermocline diagnostics
+   USE sbcwave         ! wave parameters
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
@@ -682,4 +683,8 @@
          CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
             &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
+         IF( ln_wave .AND. ln_sdw) THEN 
+            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd 
+               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
+         ENDIF 
          IF( ln_traldf_gdia ) THEN
             CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
@@ -700,4 +705,8 @@
          CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
             &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
+         IF( ln_wave .AND. ln_sdw) THEN 
+            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd 
+               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
+         ENDIF 
          IF( ln_traldf_gdia ) THEN
             CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
@@ -736,4 +745,9 @@
                &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
          ENDIF
+          
+         IF( ln_wave .AND. ln_sdw) THEN 
+            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd 
+               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
+         ENDIF 
          !                                                                                      !!! nid_W : 2D
 #if defined key_traldf_c2d
@@ -916,4 +930,10 @@
 #endif
 
+      IF( ln_wave .AND. ln_sdw ) THEN 
+         CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current 
+         CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current 
+         CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current 
+      ENDIF 
+      
       ! 3. Close all files
       ! ---------------------------------------
@@ -1021,4 +1041,13 @@
             &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
       END IF
+      ! 
+      IF( ln_wave .AND. ln_sdw ) THEN 
+         CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal", "m/s"    , &   ! StokesDrift zonal current 
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
+         CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid", "m/s"    , &   ! StokesDrift meridonal current 
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
+         CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert", "m/s"    , &   ! StokesDrift vertical current 
+            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
+      ENDIF 
 
 #if defined key_lim2
@@ -1054,4 +1083,10 @@
          CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpi*jpj*jpk, idex )!  T-cell thickness  
       END IF
+      
+      IF( ln_wave .AND. ln_sdw ) THEN 
+         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity 
+         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity 
+         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity 
+      ENDIF 
 
       ! 3. Close the file
Index: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7481)
+++ branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 7600)
@@ -34,5 +34,4 @@
    USE sbcwave         ! surface wave
    !
-   USE sbcwave         ! surface wave
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -334,4 +333,9 @@
          END DO
       END DO
+             
+      !!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 
+      !!gm  Is it correct to do so ?   I think so... 
+             
+             
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
@@ -463,4 +467,9 @@
                 &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
       ENDIF
+      ! 
+      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
+         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
+      ENDIF 
+      ! 
       ! 
       IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
Index: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7481)
+++ branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7600)
@@ -16,4 +16,5 @@
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
+   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    !!----------------------------------------------------------------------
 
@@ -222,9 +223,9 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
+      !                                                                ! =nrvm (relative vorticity or metric)
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
       !
@@ -259,5 +260,5 @@
                DO ji = 1, fs_jpim1   ! vector opt. 
                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
-                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e12f(ji,jj) 
+                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
                END DO 
             END DO 
@@ -267,5 +268,5 @@
                   zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
                        &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 * r1_e12f(ji,jj)
+                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO
             END DO
@@ -275,5 +276,5 @@
                   zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
                      &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
-                     &                   * r1_e12f(ji,jj) 
+                     &       / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO 
             END DO
@@ -284,12 +285,17 @@
                        &     + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
                        &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
-                       &     * 0.5 * r1_e12f(ji,jj) 
+                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO 
             END DO
          END SELECT
 
-         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
+         IF( ln_sco ) THEN
+            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
+            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk)
+         ELSE
+            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
+         ENDIF
 
          ! Compute and add the vorticity term trend
@@ -459,9 +465,9 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-         !                                                        ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
+         !                                                             ! =nrvm (relative vorticity or metric)
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
       !
@@ -496,5 +502,5 @@
                DO ji = 1, fs_jpim1   ! vector opt. 
                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
-                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e12f(ji,jj) 
+                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
                END DO 
             END DO
@@ -504,5 +510,5 @@
                   zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
                        &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
-                       &     * 0.5 * r1_e12f(ji,jj) 
+                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO 
             END DO 
@@ -512,5 +518,5 @@
                   zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
                      &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
-                     &                   * r1_e12f(ji,jj) 
+                     &       / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO 
             END DO
@@ -521,12 +527,26 @@
                        &     + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
                        &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
-                       &     * 0.5 * r1_e12f(ji,jj) 
+                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
                END DO 
             END DO
          END SELECT
-         !                                 !==  horizontal fluxes  ==! 
-         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
+         !
+         IF( ln_sco ) THEN
+            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
+               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
+                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
+                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
+                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
+               END DO
+            END DO
+         ELSE
+            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
+               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
+                  zwx(ji,jj) = e2u(ji,jj) * pun(ji,jj,jk)
+                  zwy(ji,jj) = e1v(ji,jj) * pvn(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+
          ! 
          ! Compute and add the vorticity term trend
@@ -570,9 +590,9 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
+      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
+      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
+      !                                                                ! =nrvm (relative vorticity or metric)
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
       !!
@@ -664,21 +684,21 @@
                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
                      &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
-                     &       * r1_e12f(ji,jj) * ze3f(ji,jj) 
-               END DO 
-            END DO
-        CASE ( 3 )                                                ! metric term
+                     &       / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj) 
+               END DO 
+            END DO
+         CASE ( 3 )                                                ! metric term
             DO jj = 1, jpjm1 
                DO ji = 1, fs_jpim1   ! vector opt. 
                   zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
                        &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
-                       &     * 0.5 * r1_e12f(ji,jj) * ze3f(ji,jj) 
-               END DO 
-            END DO
-        CASE ( 4 )                                                ! total (relative + planetary vorticity)
+                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj) 
+               END DO 
+            END DO
+         CASE ( 4 )                                                ! total (relative + planetary vorticity)
             DO jj = 1, jpjm1 
                DO ji = 1, fs_jpim1   ! vector opt. 
                   zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
                      &                         - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
-                     &                      * r1_e12f(ji,jj) ) * ze3f(ji,jj) 
+                     &       / ( e1f(ji,jj) * e2f(ji,jj) ) ) * ze3f(ji,jj) 
                END DO 
             END DO
@@ -689,5 +709,5 @@
                        &        + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
                        &            - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
-                       &        * 0.5 * r1_e12f(ji,jj)   ) * ze3f(ji,jj) 
+                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) ) * ze3f(ji,jj) 
                END DO 
             END DO 
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 7481)
+++ branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7600)
@@ -180,4 +180,8 @@
       END DO
       !
+      CALL iom_put( "ustokes",  usd  )
+      CALL iom_put( "vstokes",  vsd  )
+      CALL iom_put( "wstokes",  wsd  )
+      !
       CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh )
       CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
