Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 4644)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 4960)
@@ -18,18 +18,12 @@
    USE wrk_nemo         ! Memory Allocation
    USE sbcmod           ! Access to ln_stcor (sbc_oce) and wave parameters (sbc_wave)
-   USE sbcwave_ecmwf    ! Wave module
+   USE sbcwave          ! Wave module
+   USE timing          ! Timing
 
    IMPLICIT NONE
-
-   REAL(wp) :: rn_deptmaxstcor = 150.0_wp ! maximum depth [m] to be affected by Stokes-Coriolis effect
-   !PRIVATE
+   PRIVATE
 
    !! * Routine accessibility
-   PUBLIC dyn_stcor           ! routine called by step.F90
-
-   !! * Shared module variables
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ustc, vstc  ! Stokes-Coriolis u and v
-
-   !! * Module variables
+   PUBLIC dyn_stcor      ! routine called by step.F90
 
    !! * Substitutions
@@ -37,19 +31,9 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
+   !! NEMO/OPA 3.6 , NEMO Consortium (2010)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
  CONTAINS
-
-   INTEGER FUNCTION dynstcor_alloc()
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE dynstcor_alloc  ***
-      !!
-      !! History :  1.0  !  2012-10  (Oyvind Breivik)
-      !!----------------------------------------------------------------------
-      ALLOCATE( ustc(jpi,jpj,jpk) , vstc(jpi,jpj,jpk) , STAT=dynstcor_alloc )
-         !
-      IF( dynstcor_alloc /= 0 )   CALL ctl_warn('dynstcor_alloc: array allocate failed.')
-   END FUNCTION dynstcor_alloc
 
    SUBROUTINE dyn_stcor( kt )
@@ -59,58 +43,41 @@
       !! ** Purpose:  Add Stokes-Coriolis forcing to horizontal momentum equation.
       !!
-      !! ** History:  0.1  !   2012-10  oyvind.breivik@ecmwf.int
+      !! ** 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) ::  ztransp, zsp0, zk, zfac
-      REAL(wp) ::  zus, zvs
+      REAL(wp) :: zx1, zx2, zy1, zy2
+      !!----------------------------------------------------------------------
       !
-
-      ! Allocation at first time step.
-
-      IF ( kt == nit000 ) THEN
-         IF( dynstcor_alloc() /= 0 ) CALL ctl_stop('dyn_stcor: array allocate failed.')
+      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 ustc, vstc
+      ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities 
       !
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! Skip deep levels where the Stokes-Coriolis effect is negligible
-               IF (fsdept(ji,jj,jk) <= rn_deptmaxstcor) THEN
-                  ! Stokes transport speed estimated from Hs and Tmean
-                  ztransp = 2.0_wp*rpi*swh_wavepar(ji,jj)**2.0_wp/(16.0_wp*MAX(mwp_wavepar(ji,jj),0.0000001_wp))
+      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) )
 
-                  ! Stokes surface speed
-                  zsp0 = SQRT(ust_wavepar(ji,jj)**2 + vst_wavepar(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))
-
-                  ! The Stokes-Coriolis forcing
-                  zus =  ff(ji,jj)*vst_wavepar(ji,jj)*zfac
-                  zvs = -ff(ji,jj)*ust_wavepar(ji,jj)*zfac
-
-                  ! Store arrays of tendencies for diagnostics output
-                  ! This may be removed later for efficiency
-                  ustc(ji,jj,jk) = zus
-                  vstc(ji,jj,jk) = zvs
-
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + zus * umask(ji,jj,jk)
-                  va(ji,jj,jk) = va(ji,jj,jk) + zvs * vmask(ji,jj,jk)
-               ENDIF
+               ua(ji,jj,jk) = ua(ji,jj,jk) + 0.25 * ( zy1 + zy2 ) *umask(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( zx1 + zx2 ) *vmask(ji,jj,jk)
             ENDDO
          ENDDO
       ENDDO
       !
+      IF( nn_timing == 1 )  CALL timing_stop('dynst_cor')
+      !
    END SUBROUTINE dyn_stcor
 
-
 END MODULE dynstcor
