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_ecmwf    ! Wave module

   IMPLICIT NONE

   REAL(wp) :: rn_deptmaxstcor = 150.0_wp ! maximum depth [m] to be affected by Stokes-Coriolis effect
   !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

   !! * Substitutions
#  include "vectopt_loop_substitute.h90"
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
   !!----------------------------------------------------------------------

 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 )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_stcor ***
      !!
      !! ** Purpose:  Add Stokes-Coriolis forcing to horizontal momentum equation.
      !!
      !! ** History:  0.1  !   2012-10  oyvind.breivik@ecmwf.int
      !!----------------------------------------------------------------------
      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
      !

      ! Allocation at first time step.

      IF ( kt == nit000 ) THEN
         IF( dynstcor_alloc() /= 0 ) CALL ctl_stop('dyn_stcor: array allocate failed.')
      ENDIF

      !!----------------------------------------------------------------------
      !
      ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities ustc, vstc
      !
      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))

                  ! 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
            ENDDO
         ENDDO
      ENDDO
      !
   END SUBROUTINE dyn_stcor


END MODULE dynstcor