Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 4644)
@@ -188,2 +188,18 @@
 !-----------------------------------------------------------------------
 /
+!-----------------------------------------------------------------------
+&namsbc_wave  ! External fields from wave model
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_wave_ecmwf   ! ECWMF external wave model drag coefficient
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_wavepar ! namsbc_wavepar  parameters from wave model
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_waveparlim
+!-----------------------------------------------------------------------
+/
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4644)
@@ -240,4 +240,7 @@
    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_stcor    = .false.   !   Stokes drift read from wave model
+   ln_wavetke  = .false.   !   Wave parameters from wave model for the TKE BC
+   ln_tauoc    = .false.   !   Wave-modified stress from wave model
    cn_iceflx = 'linear'    !  redistribution of solar input into ice categories during coupling ice/atm.
 /
@@ -1164,4 +1167,12 @@
 /
 !-----------------------------------------------------------------------
+&namsbc_wave_ecmwf   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  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      =  'cdww' ,        6          , 'cdww'       ,     .true.   , .false. , 'monthly'  , ''      , ''       , '' !
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
+!-----------------------------------------------------------------------
 &namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
 !-----------------------------------------------------------------------
@@ -1178,2 +1189,30 @@
    rn_htrmax         =  200.0   ! max. depth of transition range
 /
+!-----------------------------------------------------------------------
+&namsbc_wavepar ! namsbc_wavepar  parameters from wave model (CCC not identical to namsbc_wave)
+!-----------------------------------------------------------------------
+!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!                                !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_ust      = 'ust'            ,         6         , 'ust'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_vst      = 'vst'            ,         6         , 'vst'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_swh      = 'swh'            ,         6         , 'swh'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_mwp      = 'mwp'            ,         6         , 'mwp'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_wspd     = 'wspd'           ,         6         , 'wspd'    ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_phioc    = 'physphioc'      ,         6         , 'physphioc',  .true.    , .false. , 'monthly'  , ''       , ''
+   sn_tauoc    = 'tauoc'          ,         6         , 'tauoc'   ,   .true.    , .false. , 'monthly'  , ''       , ''
+   cn_dir_wavepar  = './'      !  root directory for the location of the ECWAM files
+/
+!-----------------------------------------------------------------------
+&namsbc_waveparlim
+!-----------------------------------------------------------------------
+!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!                                !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_strn   = 'istr'            ,         6         , 'icestrain'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   cn_dir_waveparlim  = './'      !  root directory for the location of the ECWAM files
+/
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+sn_cdg     = 'cdww' ,        6          , 'cdww' , .true.   , .false. , 'monthly'  ,''         , ''  !
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
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 4644)
@@ -0,0 +1,116 @@
+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
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4644)
@@ -53,4 +53,8 @@
    LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
    !
+   LOGICAL , PUBLIC ::   ln_wavetke  = .FALSE.   !: true if wave parameters are read
+   LOGICAL , PUBLIC ::   ln_stcor    = .FALSE.   !: true if Stokes-Coriolis forcing is included 
+   LOGICAL , PUBLIC ::   ln_tauoc    = .FALSE.   !: true if wave-modified water-side stress is used
+
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
    !
@@ -69,4 +73,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_wave          !: wave model drag coefficient [N/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wspd_wavepar      !: wave model 10-m neutral wind [m/s] CCC
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tauoc_wavepar     !: normalized stress into the ocean from wave model  CCC
    !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
@@ -131,4 +138,11 @@
       ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
          &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )
+         ! Initialize these since it may not done elsewhere in the code.
+                rnf        (:,:) = 0.0_wp
+                sbc_tsc  (:,:,:) = 0.0_wp
+                qsr_hc   (:,:,:) = 0.0_wp
+                rnf_b      (:,:) = 0.0_wp
+                sbc_tsc_b(:,:,:) = 0.0_wp
+                qsr_hc_b (:,:,:) = 0.0_wp
          !
       ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
@@ -138,4 +152,6 @@
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       &
          &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
+      ALLOCATE( cdn_wave(jpi,jpj) )
+      ALLOCATE( wspd_wavepar(jpi,jpj) )
          !
 #if defined key_vvl
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4644)
@@ -38,5 +38,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE sbcwave,ONLY :  cdn_wave !wave module 
+   !USE sbcwave,ONLY :  cdn_wave !wave module  ! moved to sbc_oce 
 #if defined key_lim3 || defined key_cice
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -64,4 +64,5 @@
    
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Cd_n10 ! 10m neutral drag coefficient for output
          
    !                                             !!! CORE bulk parameters
@@ -83,4 +84,6 @@
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
+   LOGICAL  ::   ln_cdec   = .FALSE.   ! logical flag for using ec neutral wind drag coef 
+   ! References : P. Janssen, 2008, ECMWF Workshop on Ocean-Atmosphere Interactions, 10-12 November p47-60
 
    !! * Substitutions
@@ -137,6 +140,6 @@
       TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  &
-         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
-         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
+         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr , ln_cdec,           &
+         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,                    &
          &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu
       !!---------------------------------------------------------------------
@@ -163,4 +166,5 @@
             sn_qsr%ln_tint = .false.
          ENDIF
+         IF( ln_cdec .AND. lwp ) WRITE(numout,*)'Using Hans Hersbach formula for drag.'
          !                                         ! store namelist information in an array
          slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
@@ -182,4 +186,7 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
+         ! Drag coefficent so we can write to disk
+         ALLOCATE(Cd_n10(jpi,jpj))
+
          sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
@@ -239,4 +246,5 @@
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable
+      REAL(wp) ::   ztheta               ! local variable, wind direction
       REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
       REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst
@@ -283,10 +291,29 @@
 !CDIR COLLAPSE
 #endif
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
+      ! Wind vector at T points
+      IF ( ln_cdgw ) THEN
+         ! Use neutral 10-m wind speed from wave model if available
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               ! Local wind direction from non-neutral 10-m wind vector (relative to grid and no current)
+               ! since we do not have the directional information from the wave model
+               ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1), sf(jp_wndj)%fnow(ji,jj,1))
+               ! Wind vector magnitude from 10-m neutral wind speed from wave model
+               zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta)
+               zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta)
+               ! Correct for surface current, 0.0 <= rn_vfac <= 1.0
+               zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj  ) + pu(ji,jj)))
+               zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji  ,jj-1) + pv(ji,jj)))
+            END DO
          END DO
-      END DO
+      ELSE
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               ! Correct for surface current, 0.0 <= rn_vfac <= 1.0
+               zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - 0.5*rn_vfac*( pu(ji-1,jj  ) + pu(ji,jj) )  )
+               zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - 0.5*rn_vfac*( pv(ji  ,jj-1) + pv(ji,jj) )  )
+            END DO
+         END DO
+      ENDIF
       CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
       CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
@@ -353,12 +380,27 @@
 
       ! ... tau module, i and j component
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
-            taum  (ji,jj) = zztmp * wndm  (ji,jj)
-            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
-         END DO
-      END DO
+      ! Use cdn_wave instead of Cd if we have wave parameters
+      IF ( ln_cdgw ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj)
+               zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj)
+               !taum(ji,jj) = tauoc_wavepar(ji,jj)
+               !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj)
+               taum  (ji,jj) = zztmp * wndm  (ji,jj)
+               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+            ENDDO
+         ENDDO
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
+               taum  (ji,jj) = zztmp * wndm  (ji,jj)
+               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+            ENDDO
+         ENDDO
+      ENDIF
 
       ! ... add the HF tau contribution to the wind stress module?
@@ -379,4 +421,10 @@
       CALL lbc_lnk( utau(:,:), 'U', -1. )
       CALL lbc_lnk( vtau(:,:), 'V', -1. )
+
+      IF (ln_tauoc) THEN 
+         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:)
+      ENDIF
 
       !  Turbulent fluxes over ocean
@@ -755,5 +803,4 @@
       REAL(wp), DIMENSION(:,:), POINTER  ::   dT            ! air/sea temperature differeence      [K]
       REAL(wp), DIMENSION(:,:), POINTER  ::   dq            ! air/sea humidity difference          [K]
-      REAL(wp), DIMENSION(:,:), POINTER  ::   Cd_n10        ! 10m neutral drag coefficient
       REAL(wp), DIMENSION(:,:), POINTER  ::   Ce_n10        ! 10m neutral latent coefficient
       REAL(wp), DIMENSION(:,:), POINTER  ::   Ch_n10        ! 10m neutral sensible coefficient
@@ -775,10 +822,10 @@
       !
       CALL wrk_alloc( jpi,jpj, stab )   ! integer
-      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
 
       !! * Start
       !! Air/sea differences
-      dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s
+      dU10 = MAX(0.5, dU)     ! we don't want to fall under 0.5 m/s
       dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu
       dq = q_a - q_sat
@@ -788,10 +835,17 @@
       !!
       !! Neutral Drag Coefficient
-      stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0 
-      IF  ( ln_cdgw ) THEN
-        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
-        Cd_n10(:,:) =   cdn_wave
+      stab   = 0.5 + SIGN(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0 
+
+      ! Use drag coefficient from a wave model ...
+      IF (ln_cdgw) THEN
+         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
+         Cd_n10(:,:) =  cdn_wave(:,:)
       ELSE
-        Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
+      ! ... or choose a drag law
+         IF (ln_cdec) THEN
+            Cd_n10  = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21)    ! PJ (6)
+         ELSE
+            Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )              ! L & Y eq. (6a)
+         ENDIF
       ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
@@ -824,5 +878,9 @@
 
            !! Updating the neutral 10m transfer coefficients :
-           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
+           IF (ln_cdec) THEN
+             Cd_n10  = 1E-3 * ( 1.03 + 0.04 * U_n10**1.48 ) / (U_n10**0.21)    ! PJ (6)
+           ELSE
+             Cd_n10  = 1E-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 )              ! L & Y eq. (6a)
+           ENDIF
            sqrt_Cd_n10 = sqrt(Cd_n10)
            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
@@ -847,5 +905,5 @@
       !!
       CALL wrk_dealloc( jpi,jpj, stab )   ! integer
-      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
       !
@@ -891,5 +949,4 @@
       REAL(wp), DIMENSION(:,:), POINTER ::   dT            ! air/sea temperature differeence   [K]
       REAL(wp), DIMENSION(:,:), POINTER ::   dq            ! air/sea humidity difference       [K]
-      REAL(wp), DIMENSION(:,:), POINTER ::   Cd_n10        ! 10m neutral drag coefficient
       REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient
       REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient
@@ -911,5 +968,5 @@
       IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_2Z')
       !
-      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
       CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
@@ -923,9 +980,16 @@
       !! Neutral Drag Coefficient :
       stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE
-      IF( ln_cdgw ) THEN
-        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
-        Cd_n10(:,:) =   cdn_wave
+
+      ! Use drag coefficient from a wave model ...
+      IF (ln_cdgw) THEN
+         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
+         Cd_n10(:,:) =  cdn_wave
       ELSE
-        Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
+      ! ... or choose a drag law
+         IF (ln_cdec) THEN
+            Cd_n10  = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21)    ! PJ (6)
+         ELSE
+            Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )              ! L & Y eq. (6a)
+         ENDIF
       ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
@@ -996,5 +1060,5 @@
       END DO
       !!
-      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
       CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4644)
@@ -50,4 +50,5 @@
    USE timing           ! Timing
    USE sbcwave          ! Wave module
+   USE sbcwave_ecmwf    ! ECMWF wave module
 
    IMPLICIT NONE
@@ -82,7 +83,8 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
-         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, cn_iceflx
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx ,  ln_blk_clio, ln_blk_core, ln_cpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice ,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
+         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave    , ln_sdw     , nn_lsm,   &
+         &             cn_iceflx , ln_stcor  , ln_tauoc, ln_wavetke 
       INTEGER  ::   ios
       !!----------------------------------------------------------------------
@@ -135,4 +137,8 @@
          WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
          WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm
+         WRITE(numout,*) '              Drag coefficient from wave model           ln_cdgw     = ', ln_cdgw  
+         WRITE(numout,*) '              Stokes-Coriolis forcing from wave model    ln_stcor    = ', ln_stcor  
+         WRITE(numout,*) '              TKE wave breaking BC from wave model       ln_wavetke  = ', ln_wavetke
+         WRITE(numout,*) '              Wave-modified stress from wave model       ln_tauoc    = ', ln_tauoc
       ENDIF
 
@@ -230,5 +236,5 @@
       ELSE
       IF ( ln_cdgw .OR. ln_sdw  )                                         & 
-         &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     &
+         &   CALL ctl_warn('Not Activated Wave Module (ln_wave=F) but     &
          & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
       ENDIF 
@@ -314,4 +320,11 @@
       CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
       !                                                  ! averaged over nf_sbc time-step
+
+      IF (ln_cdgw) CALL sbc_wave_ecmwf( kt )
+      ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for TKE are needed
+      ! Also read it if drag coefficient to ensure that we use the neutral 10m from WAM.
+      IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc .OR. ln_cdgw) THEN
+         CALL sbc_wavepar( kt )
+      ENDIF
 
       IF (ln_wave) CALL sbc_wave( kt )
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 4644)
@@ -31,5 +31,4 @@
    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 
@@ -61,4 +60,5 @@
       USE divcur
       USE wrk_nemo
+      USE sbc_oce, ONLY : cdn_wave
 #if defined key_bdy
       USE bdy_oce, ONLY : bdytmask
@@ -100,6 +100,4 @@
             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' )
-            ALLOCATE( cdn_wave(jpi,jpj) )
-            cdn_wave(:,:) = 0.0
         ENDIF
          IF ( ln_sdw ) THEN
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave_ecmwf.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave_ecmwf.F90	(revision 4644)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave_ecmwf.F90	(revision 4644)
@@ -0,0 +1,217 @@
+MODULE sbcwave_ecmwf
+   !!======================================================================
+   !!                       ***  MODULE  sbcwave  ***
+   !! Wave module 
+   !!======================================================================
+   !! History :  3.3.1  !   2011-09  (Adani M)  Original code
+   !!----------------------------------------------------------------------
+   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 sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
+   USE sbcblk_core, ONLY: Cd_n10
+   
+   !!----------------------------------------------------------------------
+   !!   sbc_wave_ecmwf : read drag coefficient from wave model in netcdf files 
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_wavepar       ! 
+   PUBLIC   sbc_wave_ecmwf    ! routine called in sbc_blk_core or sbc_blk_mfs
+   
+   INTEGER , PARAMETER ::   jpfld_wavepar = 7    ! maximum number of fields to be read 
+   INTEGER , PARAMETER ::   jp_ust = 1           ! index of Stokes velocity (east component) (m/s)  at T-point
+   INTEGER , PARAMETER ::   jp_vst = 2           ! index of Stokes velocity (north component) (m/s) at T-point
+   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave height (m)
+   INTEGER , PARAMETER ::   jp_mwp = 4           ! index of mean wave period (s)
+   INTEGER , PARAMETER ::   jp_phioc = 5         ! index of normalized energy flux into the ocean (non-dim)
+   INTEGER , PARAMETER ::   jp_tauoc = 6         ! index of normalized wave stress into the ocean (non-dim)
+   INTEGER , PARAMETER ::   jp_wspd = 7          ! index of 10 m neutral wind speed (m/s)
+
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wave    ! structure of input fields (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wavepar ! structure of input fields (file informations, fields read)
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: ust_wavepar, vst_wavepar, swh_wavepar, mwp_wavepar
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: phioc_wavepar
+
+   !! * 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 sbc_wave_ecmwf( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_apr  ***
+      !!
+      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
+      !!
+      !! ** Method  : - Read namelist namsbc_wave
+      !!              - Read Cd_n10 fields in netcdf files 
+      !! ** action  :   
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
+      INTEGER                ::  ierror   ! return error code
+      INTEGER                ::  ios      ! io return code
+      CHARACTER(len=100)     ::  cn_dir_cdg                      ! Root directory for location of drag coefficient files
+      TYPE(FLD_N)            ::  sn_cdg                          ! informations about the fields to be read
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_wave_ecmwf/  sn_cdg, cn_dir_cdg
+      !!---------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         !                                      ! -------------------- !
+         !                                            !* set file information (default values)
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !              !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !              !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_cdg = FLD_N( 'cdg_wave'  ,    1.    ,'drag_coeff',  .true.    , .false. ,   'daily'   , ''       , ''      , '' )
+         cn_dir_cdg = './'          ! directory in which the Patm data are 
+         
+
+         REWIND( numnam_ref )              !  Namelist namsbc_wavepar in reference namelist 
+         READ  ( numnam_ref, namsbc_wave_ecmwf, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave_ecmwf in reference namelist', lwp )
+
+         REWIND( numnam_cfg )              !  Namelist namsbc_wave_ecmwf in configuration namelist
+         READ  ( numnam_cfg, namsbc_wave_ecmwf, IOSTAT = ios, ERR = 902 )
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave_ecmwf in configuration namelist', lwp )
+         IF(lwm) WRITE ( numond, namsbc_wave_ecmwf )
+         !
+
+         ALLOCATE( sf_wave(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_ecmwf: unable to allocate sf_wave structure' )
+         !
+         CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave_ecmwf', 'Wave module ', 'namsbc_wave_ecmwf' )
+                                ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1)   )
+         IF( sn_cdg%ln_tint )   ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) )
+         ! Allocation done by sbc_oce 
+         cdn_wave(:,:) = 0.0
+      ENDIF
+         !
+         !
+      CALL fld_read( kt, nn_fsbc, sf_wave )      !* read drag coefficient from external forcing
+      cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1)
+
+   END SUBROUTINE sbc_wave_ecmwf
+
+
+   SUBROUTINE sbc_wavepar( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_wavepar  ***
+      !!
+      !! ** Purpose :   Provide at each time step wave model parameters, including the
+      !!                Stokes drift (east and north components), significant wave height and the
+      !!                mean wave period as well as the normalized stress and energy flux into the 
+      !!                ocean for TKE.
+      !!
+      !!
+      !! ** Method  : - Read namelist namsbc_wavepar
+      !!              - Read fields in NetCDF files 
+      !! ** action  :   
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt  ! ocean time step
+      !!
+      INTEGER  ::  ierror ! return error code
+      INTEGER  ::   ji      ! dummy loop index
+      INTEGER  ::   jfld    ! dummy loop arguments
+      INTEGER  ::   ios     ! io return code
+      !!
+      CHARACTER(len=100)     ::  cn_dir_wavepar ! Root directory for location of ECWAM wave parameter fields
+      TYPE(FLD_N), DIMENSION(jpfld_wavepar) :: slf_i ! array of namelist informations on the fields to read
+      TYPE(FLD_N)            ::  sn_ust, sn_vst, sn_swh, sn_mwp ! information about the fields to be read
+      TYPE(FLD_N)            ::  sn_phioc, sn_tauoc
+      TYPE(FLD_N)            ::  sn_wspd
+      !TYPE(FLD_N)            ::  sn_cdww
+      !!---------------------------------------------------------------------
+      !NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_phioc, sn_tauoc, cn_dir_wavepar
+      NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_wspd, sn_phioc, sn_tauoc, cn_dir_wavepar
+      !!---------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         !                                      ! -------------------- !
+         ! set file information (default values)
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !            !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !            !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_ust = FLD_N( 'ust'   ,     6     ,  'ust'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_vst = FLD_N( 'vst'   ,     6     ,  'vst'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_swh = FLD_N( 'swh'   ,     6     ,  'swh'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_mwp = FLD_N( 'mwp'   ,     6     ,  'mwp'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         !sn_cdww = FLD_N( 'cdww'  ,    6     ,  'cdww'    ,  .true.    , .false. ,   'monthly' , ''       , ''      ,  '' )
+         sn_wspd = FLD_N( 'wspd'   ,   6     ,  'wspd'    ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_phioc = FLD_N( 'phioc'   , 6     ,  'phioc'   ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_tauoc = FLD_N( 'tauoc'   , 6     ,  'tauoc'   ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         cn_dir_wavepar = './'          ! directory in which the wave data are found
+
+         REWIND( numnam_ref )              !  Namelist namsbc_wavepar in reference namelist 
+         READ  ( numnam_ref, namsbc_wavepar, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wavepar in reference namelist', lwp )
+
+         REWIND( numnam_cfg )              !  Namelist namsbc_wavepar in configuration namelist
+         READ  ( numnam_cfg, namsbc_wavepar, IOSTAT = ios, ERR = 902 )
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wavepar in configuration namelist', lwp )
+         IF(lwm) WRITE ( numond, namsbc_wavepar )
+         !
+         slf_i(jp_ust) = sn_ust 
+         slf_i(jp_vst) = sn_vst
+         slf_i(jp_swh) = sn_swh
+         slf_i(jp_mwp) = sn_mwp
+         !slf_i(jp_cdww) = sn_cdww
+         slf_i(jp_wspd) = sn_wspd
+         slf_i(jp_phioc) = sn_phioc
+         slf_i(jp_tauoc) = sn_tauoc
+
+         ALLOCATE( sf_wavepar(jpfld_wavepar), STAT=ierror )          !* set sf_wavepar structure
+         IF ( ierror > 0 ) THEN
+            CALL ctl_stop( 'STOP', 'sbc_wavepar: unable to allocate sf_wavepar structure' ) ; RETURN
+         ENDIF
+         !
+         jfld = jpfld_wavepar
+         DO ji = 1, jpfld_wavepar
+            ALLOCATE( sf_wavepar(ji)%fnow(jpi,jpj,1)   )
+            IF ( slf_i(ji)%ln_tint ) ALLOCATE( sf_wavepar(ji)%fdta(jpi,jpj,1,2) )
+         ENDDO
+         ALLOCATE( ust_wavepar(jpi,jpj) )
+         ALLOCATE( vst_wavepar(jpi,jpj) )
+         ALLOCATE( swh_wavepar(jpi,jpj) )
+         ALLOCATE( mwp_wavepar(jpi,jpj) )
+         !ALLOCATE( wspd_wavepar(jpi,jpj) ) ! allocation done by sbc_oce
+         ALLOCATE( phioc_wavepar(jpi,jpj) )
+         ALLOCATE( tauoc_wavepar(jpi,jpj) )
+
+         !
+         CALL fld_fill( sf_wavepar, slf_i, cn_dir_wavepar, 'sbc_wavepar', 'Wave module ', 'namsbc_wavepar' )
+         !
+      ENDIF
+         !
+      CALL fld_read( kt, nn_fsbc, sf_wavepar )      !* read wave parameters from ECWAM NetCDF file
+      ust_wavepar(:,:) = sf_wavepar(jp_ust)%fnow(:,:,1)
+      vst_wavepar(:,:) = sf_wavepar(jp_vst)%fnow(:,:,1)
+      swh_wavepar(:,:) = sf_wavepar(jp_swh)%fnow(:,:,1)
+      mwp_wavepar(:,:) = sf_wavepar(jp_mwp)%fnow(:,:,1)
+      !cdww_wavepar(:,:) = sf_wavepar(jp_cdww)%fnow(:,:,1)
+      wspd_wavepar(:,:) = sf_wavepar(jp_wspd)%fnow(:,:,1)
+      phioc_wavepar(:,:) = sf_wavepar(jp_phioc)%fnow(:,:,1)
+      tauoc_wavepar(:,:) = sf_wavepar(jp_tauoc)%fnow(:,:,1)
+
+   END SUBROUTINE sbc_wavepar
+
+END MODULE sbcwave_ecmwf
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 4644)
@@ -26,4 +26,5 @@
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!                 !  2012-12  (oyvind.breivik@ecwmf.int) Changed to wave breaking source term
    !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
@@ -42,4 +43,6 @@
    USE domvvl         ! domain: variable volume layer
    USE sbc_oce        ! surface boundary condition: ocean
+   USE sbcwave        ! wave parameters
+   USE sbcwave_ecmwf  ! ECMWF wave parameters
    USE zdf_oce        ! vertical physics: ocean variables
    USE zdfmxl         ! vertical physics: mixed layer
@@ -78,4 +81,9 @@
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
    REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
+   REAL(wp) ::   rn_alpavg = 100.0_wp      ! Average alpha (normalized, water-side energy flux, only used in
+                                           ! case we use the source term formulation with no wave model information
+                                           ! NB! The option of running with source term and no wave model has
+                                           ! not been properly tested yet but is included for consistency (Oyvind Breivik, 2013-07-16)
+   REAL(wp) ::   rn_deptmaxtkew = 50.0_wp   ! maximum depth [m] to be affected by TKE from breaking waves 
 
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
@@ -85,4 +93,6 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   e0             !: top level turbulent kinetic energy   [m2/s2] from waves
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dfac           !: wave depth scaling factor   
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
@@ -109,4 +119,5 @@
       !!                ***  FUNCTION zdf_tke_alloc  ***
       !!----------------------------------------------------------------------
+      INTEGER :: iexalloc
       ALLOCATE(                                                                    &
 #if defined key_c1d
@@ -118,4 +129,12 @@
          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      )
          !
+      IF ( ln_wavetke ) THEN
+         ALLOCATE( e0(jpi,jpj) , dfac(jpi,jpj), STAT= iexalloc )
+         ! Setting arrays to zero initially.
+         e0(:,:) = 0.0_wp
+         dfac(:,:) = 0.0_wp
+         !
+         zdf_tke_alloc = zdf_tke_alloc + iexalloc
+      ENDIF
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
       IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
@@ -219,4 +238,8 @@
       REAL(wp) ::   zus   , zwlc  , zind            !    -         -
       REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
+      REAL(wp) ::   lambda !    -         -
+      REAL(wp) ::   zphio_flux              !    -         -
+      REAL(wp) ::   sbr, z0, fb                     !    -         -
+
 !!bfr      REAL(wp) ::   zebot                           !    -         -
       INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
@@ -239,10 +262,32 @@
       !                     !  Surface boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      
+      ! Wave parameters read?
+      IF ( ln_wavetke ) THEN
+         ! Breaking-wave TKE injection as flux
+         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Roughness length proportional to wave height, capped at Hs = 0.1 m
+               z0 = 0.5_wp*MAX(swh_wavepar(ji,jj), 0.1_wp)
+               lambda = 2.38_wp/z0
+               ! The energy flux from the wave model in physical units
+               zphio_flux = MAX(phioc_wavepar(ji,jj), 0.0_wp)
+               ! TKE at surface (Mellor and Blumberg, 2004)
+               e0(ji,jj) = 0.5_wp*( 15.8_wp*zphio_flux/rau0 )**0.67_wp
+               ! Reducing the surface TKE by scaling with T-depth of first vertical level
+               dfac(ji,jj) = MAX( (2.0_wp*lambda*fsdept(ji,jj,1) )/3.0_wp, 0.00001_wp )
+               ! TKE surface boundary condition weighted by the thickness of the first level
+               en(ji,jj,1) = MAX( rn_emin0, e0(ji,jj) * ( 1.0_wp - EXP(-dfac(ji,jj)) ) / dfac(ji,jj) ) * tmask(ji,jj,1)
+            ENDDO
+         ENDDO
+         ! If wave parameters are not read, use rn_ebb, corresponding to alpha = 100.0
+      ELSE
+         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Estimate average energy flux from std value found in rn_ebb (corresponds to alpha=100.0)
+               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            ENDDO
+         ENDDO
+      ENDIF
+ 
 !!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -443,4 +488,8 @@
       ENDIF
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      IF ( ln_wavetke ) THEN
+         CALL lbc_lnk( e0,     'W', 1. ) ! Lateral boundary conditions (sign unchanged)
+         CALL lbc_lnk( dfac,   'W', 1. ) ! Lateral boundary conditions (sign unchanged)
+      ENDIF
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4644)
@@ -286,4 +286,5 @@
                                CALL dyn_ldf( kstp )         ! lateral mixing
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
+        IF( ln_stcor )         CALL dyn_stcor   ( kstp )    ! Stokes-Coriolis forcing (ln_stcor set in SBC/sbc_oce.F90)
 #if defined key_agrif
         IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge
Index: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 4644)
@@ -23,4 +23,5 @@
 
    USE sbcmod           ! surface boundary condition       (sbc     routine)
+   USE sbc_oce , ONLY: ln_stcor      ! seems to be necessary but should be included in sbcmod????
    USE sbcrnf           ! surface boundary condition: runoff variables
    USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step)
@@ -51,4 +52,5 @@
    USE dynspg           ! surface pressure gradient        (dyn_spg routine)
    USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine)
+   USE dynstcor         ! Stokes-Coriolis forcing          (dyn_stcor routine)
 
    USE dynnxt           ! time-stepping                    (dyn_nxt routine)
