Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7877)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7878)
@@ -283,4 +283,7 @@
                            !  = 2 standard formulation with constant drag coefficient
                            !  = 3 momentum calculated from core forcing fields
+   nn_sdrift = 0           !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift
+                           !  = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
+                           !  = 1 Phillips:                      v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]
 /
 !-----------------------------------------------------------------------
@@ -394,4 +397,5 @@
    sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
+   sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    sn_rcv_tauoc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
@@ -1303,4 +1307,5 @@
    sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
+   sn_wfr      =  'sdw_wave' ,        1          , 'wave_freq'  ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7877)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7878)
@@ -71,4 +71,5 @@
    LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
    INTEGER , PUBLIC ::   nn_drag        ! type of formula to calculate wind stress from wind components
+   INTEGER , PUBLIC ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
    !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
@@ -104,4 +105,10 @@
    INTEGER, PUBLIC, PARAMETER ::   jp_const = 2        ! standard formulation with constant drag coefficient 
    INTEGER, PUBLIC, PARAMETER ::   jp_mcore = 3        ! momentum calculated from core forcing fields 
+
+   !!----------------------------------------------------------------------
+   !!           Stokes drift parameterization
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
+   INTEGER, PUBLIC, PARAMETER ::   jp_phillips = 1     ! Phillips:     v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]
 
    !!----------------------------------------------------------------------
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7877)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7878)
@@ -115,5 +115,6 @@
    INTEGER, PARAMETER ::   jpr_tauoc  = 50            ! Stress fraction adsorbed by waves 
    INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient 
-   INTEGER, PARAMETER ::   jprcv      = 51            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_wfreq  = 52            ! Wave peak frequency 
+   INTEGER, PARAMETER ::   jprcv      = 52            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -167,5 +168,6 @@
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev  
    ! Received from waves  
-   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_wdrag
+   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper, &
+                    sn_rcv_wfreq,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_wdrag
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -242,7 +244,8 @@
          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &  
          &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &  
-         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_tauoc ,   & 
-         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   & 
-         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp
+         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wfreq,    &
+         &                  sn_rcv_tauoc, sn_rcv_wdrag, sn_rcv_qns   , sn_rcv_emp   , sn_rcv_rnf,      &
+         &                  sn_rcv_cal , sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp  , nn_cplmodel,     &
+         &                  ln_usecplmask
       !!---------------------------------------------------------------------
       !
@@ -291,4 +294,5 @@
          WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'  
          WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'  
+         WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'  
          WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')'  
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
@@ -559,4 +563,9 @@
          cpl_wper = .TRUE. 
       ENDIF 
+      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency 
+      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN 
+         srcv(jpr_wfreq)%laction = .TRUE. 
+         cpl_wfreq = .TRUE. 
+      ENDIF 
       srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number 
       IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN 
@@ -1182,4 +1191,9 @@
       ! 
       !                                                      ! ========================= !  
+      !                                                      !    Wave peak frequency    ! 
+      !                                                      ! ========================= !  
+         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1) 
+      !
+      !                                                      ! ========================= !  
       !                                                      !    Vertical mixing Qiao   ! 
       !                                                      ! ========================= !  
@@ -1188,5 +1202,5 @@
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode 
          IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction & 
-                                                                    .OR. srcv(jpr_hsig)%laction ) & 
+                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) & 
             CALL sbc_stokes() 
       ENDIF 
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7877)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7878)
@@ -90,5 +90,5 @@
          &             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_phioc  , ln_wavcpl , nn_drag
+         &             ln_phioc  , ln_wavcpl , nn_drag, nn_sdrift
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -236,4 +236,6 @@
          IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) &
              CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2')
+         IF ( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) &
+             CALL ctl_stop( 'The chosen nn_sdrift for vertical Stokes drift must be 0, or 1')
       ELSE
          IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc ) &  
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7877)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7878)
@@ -42,4 +42,5 @@
    LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE.
@@ -51,4 +52,5 @@
    INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point
    INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
+   INTEGER ::   jp_wfr   ! index of wave peak frequency         (s^-1)   at T-point
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient
@@ -59,4 +61,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !: 
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:
@@ -99,43 +102,78 @@
       CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
       !
-      zfac = 2.0_wp * rpi / 16.0_wp
-      DO jj = 1, jpj               ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) )
-         DO ji = 1, jpi
-            ! Stokes drift velocity estimated from Hs and Tmean
-            ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
-            ! Stokes surface speed
-            tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
-            ! Wavenumber scale
-            zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
-         END DO
-      END DO
-      DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
-         DO ji = 1, jpim1
-            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
-            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
-            !
-            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
-            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
-         END DO
-      END DO
-      !
-      !                       !==  horizontal Stokes Drift 3D velocity  ==!
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
-               zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
-               !                          
-               zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
-               zkh_v = zk_v(ji,jj) * zdep_v
-               !                                ! Depth attenuation
-               zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
-               zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
+      ! select parameterization for the calculation of vertical Stokes drift
+      SELECT CASE ( nn_sdrift )
+      !
+      CASE ( jp_breivik )
+         zfac = 2.0_wp * rpi / 16.0_wp
+         DO jj = 1, jpj               ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) )
+            DO ji = 1, jpi
+               ! Stokes drift velocity estimated from Hs and Tmean
+               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
+               ! Stokes surface speed
+               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
+               ! Wavenumber scale
+               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
+            END DO
+         END DO
+         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
+            DO ji = 1, jpim1
+               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
+               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
                !
-               usd(ji,jj,jk) = zda_u * zk_u(ji,jj) * umask(ji,jj,jk)
-               vsd(ji,jj,jk) = zda_v * zk_v(ji,jj) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+            END DO
+         END DO
+         !
+         !                       !==  horizontal Stokes Drift 3D velocity  ==!
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
+                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
+                  !                          
+                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
+                  zkh_v = zk_v(ji,jj) * zdep_v
+                  !                                ! Depth attenuation
+                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
+                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
+                  !
+                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      CASE ( jp_phillips )
+         DO jj = 1, jpjm1              ! Peak wavenumber & Stokes drift velocity at u- & v-points
+            DO ji = 1, jpim1
+               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav
+               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav
+               !
+               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+            END DO
+         END DO
+         !
+         !                       !==  horizontal Stokes Drift 3D velocity  ==!
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
+                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
+                  !                          
+                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
+                  zkh_v = zk_v(ji,jj) * zdep_v
+                  !                                ! Depth attenuation: beta=1 for Phillips
+                  zda_u = EXP( -2.0_wp*zkh_u ) - 1.0*SQRT(2.0*rpi*zkh_u) * ERFC(SQRT(2.0*zkh_u))
+                  zda_v = EXP( -2.0_wp*zkh_v ) - 1.0*SQRT(2.0*rpi*zkh_v) * ERFC(SQRT(2.0*zkh_v))
+                  !
+                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      END SELECT
+
       CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. )
       !
@@ -227,6 +265,6 @@
          CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing
          rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density
-         WHERE( rn_crban <  10.0 ) rn_crban =  10.0
-         WHERE( rn_crban > 300.0 ) rn_crban = 300.0
+         WHERE( rn_crban < -1000.0 ) rn_crban = 0.0
+         WHERE( rn_crban >  1000.0 ) rn_crban = 0.0
       ENDIF
 
@@ -245,4 +283,9 @@
                WHERE( wmp <   0.0 ) wmp = 0.0
             ENDIF
+            IF( jp_wfr > 0 ) THEN
+               wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency 
+               WHERE( wfreq <    0.0 ) wfreq = 0.001 
+               WHERE( wfreq >  100.0 ) wfreq = 0.001
+            ENDIF
             IF( jp_usd > 0 ) THEN
                ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
@@ -265,6 +308,8 @@
          !                                         !==  Computation of the 3d Stokes Drift  ==! 
          !
-         IF( jpfld == 4 )   CALL sbc_stokes()            ! Calculate only if required fields are read
-         !                                               ! In coupled wave model-NEMO case the call is done after coupling
+         IF( (nn_sdrift==jp_breivik  .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. &
+             (nn_sdrift==jp_phillips .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) &
+            CALL sbc_stokes()            ! Calculate only if required fields are read
+         !                               ! In coupled wave model-NEMO case the call is done after coupling
          !
       ENDIF
@@ -292,8 +337,9 @@
       CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read
-      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  sn_phioc, &
-                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
-      !
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc, sn_phioc
+      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_phioc, &
+                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum , &
+                             &   sn_tauoc      ! informations about the fields to be read
+      !
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc, sn_phioc
       !!---------------------------------------------------------------------
       !
@@ -345,5 +391,5 @@
       IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled
          jpfld=0
-         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0
+         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
          IF( .NOT. cpl_sdrftx ) THEN
             jpfld  = jpfld + 1
@@ -361,4 +407,8 @@
             jpfld  = jpfld + 1
             jp_wmp = jpfld
+         ENDIF
+         IF( .NOT. cpl_wfreq ) THEN
+            jpfld  = jpfld + 1
+            jp_wfr = jpfld
          ENDIF
 
@@ -370,4 +420,5 @@
             IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
             IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
+            IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
             ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' )
@@ -382,4 +433,5 @@
          ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
          ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     )
+         ALLOCATE( wfreq (jpi,jpj) )
          ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
          ALLOCATE( div_sd(jpi,jpj) )
