Index: /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8749)
+++ /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8750)
@@ -216,4 +216,8 @@
    ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)
    ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 
+   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))]
+                           !   = 2 Phillips as (1) but using the wave frequency from a wave model
    ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)
    ln_stcor    = .false.   !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)
@@ -298,4 +302,5 @@
    sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_wstrf  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
@@ -456,4 +461,5 @@
    sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_wfr      =  'sdw_wave' ,        1          , 'wfr'        ,     .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/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8749)
+++ /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8750)
@@ -67,4 +67,6 @@
    LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
    !
+   INTEGER , PUBLIC ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
+   !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
    !
@@ -79,4 +81,11 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_none    = 5        !: for OPA when doing coupling via SAS module
    
+   !!----------------------------------------------------------------------
+   !!           Stokes drift parametrization definition 
+   !!----------------------------------------------------------------------
+   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))]
+   INTEGER , PUBLIC, PARAMETER ::   jp_peakfr   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale
+
    !!----------------------------------------------------------------------
    !!           component definition
Index: /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8749)
+++ /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8750)
@@ -117,6 +117,7 @@
    INTEGER, PARAMETER ::   jpr_isf    = 52
    INTEGER, PARAMETER ::   jpr_icb    = 53
-
-   INTEGER, PARAMETER ::   jprcv      = 53   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_wfreq  = 54   ! Wave peak frequency
+
+   INTEGER, PARAMETER ::   jprcv      = 54   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -170,5 +171,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_wstrf,sn_rcv_wdrag
+   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag, &
+                    sn_rcv_wfreq
    !                                   ! Other namelist parameters
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -245,5 +247,5 @@
          &                  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_icb , sn_rcv_isf
+         &                  sn_rcv_icb , sn_rcv_isf   , sn_rcv_wfreq
 
       !!---------------------------------------------------------------------
@@ -295,4 +297,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_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
@@ -578,4 +581,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
@@ -1165,4 +1173,9 @@
       !                                                      ! ========================= ! 
          IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
+      ! 
+      !                                                      ! ========================= !  
+      !                                                      !    Wave peak frequency    ! 
+      !                                                      ! ========================= !  
+         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)
       !
       !                                                      ! ========================= ! 
@@ -1173,5 +1186,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 ) THEN
+                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction ) THEN
             CALL sbc_stokes()
          ENDIF
Index: /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8749)
+++ /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8750)
@@ -96,5 +96,5 @@
          &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     &
          &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauoc  , ln_stcor   ,     &
-         &             nn_lsm
+         &             nn_lsm, nn_sdrift
       !!----------------------------------------------------------------------
       !
@@ -157,7 +157,13 @@
          WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
          WRITE(numout,*) '               Stokes drift corr. to vert. velocity ln_sdw        = ', ln_sdw
+         WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift
          WRITE(numout,*) '               wave modified ocean stress           ln_tauoc      = ', ln_tauoc
          WRITE(numout,*) '               Stokes coriolis term                 ln_stcor      = ', ln_stcor
          WRITE(numout,*) '               neutral drag coefficient (CORE, MFS) ln_cdgw       = ', ln_cdgw
+      ENDIF
+      !
+      IF( ln_sdw ) THEN
+         IF( .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr) ) &
+            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
       ENDIF
       !
Index: /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 8749)
+++ /branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 8750)
@@ -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_wstrf  = .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         (1/s)    at T-point
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_cd      ! structure of input fields (file informations, fields read) Drag Coefficient
@@ -58,4 +60,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(:,:)   ::   tauoc_wave          !:  
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
@@ -96,45 +99,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
+      ! select parameterization for the calculation of vertical Stokes drift
+      ! exp. wave number at t-point
+      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) )
+         zfac = 2.0_wp * rpi / 16.0_wp
+         DO jj = 1, jpj
+            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 )
+               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
                ! Stokes surface speed
-               zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) )
-               tsd2d(ji,jj) = zsp0
+               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
                ! Wavenumber scale
-               zk_t(ji,jj) = ABS( zsp0 ) / 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
+               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
+      ELSE IF( nn_sdrift==jp_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
+         DO jj = 1, jpjm1
+            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
+      ENDIF
       !
       !                       !==  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   
+      IF( nn_sdrift==jp_breivik ) THEN
+         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
+      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr ) THEN
+         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 ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u))
+                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.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
+      ENDIF
+
       CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. )
       !
@@ -222,4 +258,5 @@
             IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
             IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+            IF( jp_wfr > 0 )   wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency
             IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
             IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
@@ -234,6 +271,9 @@
          !                                         !==  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 .OR. nn_sdrift==jp_phillips) .AND. &
+                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. &
+             (nn_sdrift==jp_peakfr .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
@@ -262,7 +302,7 @@
       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_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_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
       !!---------------------------------------------------------------------
       !
@@ -302,5 +342,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
@@ -311,11 +351,15 @@
             jp_vsd = jpfld
          ENDIF
-         IF( .NOT. cpl_hsig ) THEN
+         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
             jpfld  = jpfld + 1
             jp_hsw = jpfld
          ENDIF
-         IF( .NOT. cpl_wper ) THEN
+         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
             jpfld  = jpfld + 1
             jp_wmp = jpfld
+         ENDIF
+         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakfr ) THEN
+            jpfld  = jpfld + 1
+            jp_wfr = jpfld
          ENDIF
 
@@ -327,4 +371,6 @@
             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_wave structure' )
@@ -339,4 +385,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) )
