Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7905)
@@ -266,9 +266,4 @@
                            !     =2 annual global mean of e-p-r set to zero
    ln_wave = .false.       !  Activate coupling with wave (T => fill namsbc_wave)  
-   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)   
-   ln_tauoc= .false.       !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)  
-   ln_phioc= .false.       !  Activate wave to ocean energy (T => ln_wave=.true. & fill namsbc_wave)  
-   ln_stcor= .false.       !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)  
    nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
                            !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
@@ -278,12 +273,4 @@
                            !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled
                            !  = 2  Redistribute a single flux over categories (coupled mode only)
-   nn_drag   = 0           !  formula to calculate momentum from the wind components
-                           !  = 0 UKMO SHELF formulation
-                           !  = 1 standard formulation with forced of coupled drag coefficient
-                           !  = 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,6 +381,5 @@
    sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
-   sn_rcv_sdrfx  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
-   sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
+   sn_rcv_sdrft  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
@@ -948,5 +934,4 @@
    ln_zdfexp   = .false.   !  time-stepping: split-explicit (T) or implicit (F) time stepping
    nn_zdfexp   =    3            !  number of sub-timestep for ln_zdfexp=T
-   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010) (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -1022,7 +1007,4 @@
    nn_stab_func  =     2   !  stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB)
    nn_clos       =     1   !  predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen)
-   nn_wmix       =     0   !  type of wave breaking mixing
-           !                             ! = 0 Craig and Banner formulation (original NEMO formulation) 
-           !                             ! = 1 Janssen formulation (no assumption of direct energy conversion) 
 /
 !-----------------------------------------------------------------------
@@ -1302,4 +1284,22 @@
 !              !  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      !
+   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model
+   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift
+   ln_stcor    = .false.   !  Activate Stokes Coriolis term
+   ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress
+   ln_phioc    = .false.   !  Activate wave to ocean energy
+   ln_rough    = .false.   !  Wave roughness equals the significant wave height
+   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010)
+   nn_drag   = 0      !  formula to calculate momentum from the wind components
+                           ! = 0 UKMO SHELF formulation
+                           ! = 1 standard formulation with forced of coupled drag coefficient
+                           ! = 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))]
+   nn_wmix     =  0   !  type of wave breaking mixing
+                           ! = 0 Craig and Banner formulation (original NEMO formulation) 
+                           ! = 1 Janssen formulation (no assumption of direct energy conversion) 
    sn_cdg      =  'sdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
@@ -1311,6 +1311,6 @@
    sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    sn_phioc    =  'sdw_wave' ,        1          , 'wave_energy',     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
+   cn_dir      = './'  !  root directory for the location of wave forcing files
 ! 
-   cn_dir      = './'  !  root directory for the location of drag coefficient files
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7905)
@@ -163,4 +163,5 @@
                              CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
                              CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
+            IF( ln_stcor )   CALL vor_ens( kt, nrvm, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
          ENDIF
@@ -182,5 +183,5 @@
          ELSE
                              CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity
-            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
+            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
          ENDIF
          !
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 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7905)
@@ -66,10 +66,12 @@
    LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model
    LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model
-   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
+   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d Stokes drift from wave model
+   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
    LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used 
    LOGICAL , PUBLIC ::   ln_phioc       !: true if wave energy to ocean is used 
-   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
+   LOGICAL , PUBLIC ::   ln_rough       !: true if wave roughness equals significant wave height
+   LOGICAL , PUBLIC ::   ln_zdfqiao     !: Enhanced wave vertical mixing Qiao(2010) formulation flag
    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
+   INTEGER , PUBLIC ::   nn_wmix        ! type of wave breaking mixing
    !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
@@ -107,8 +109,9 @@
 
    !!----------------------------------------------------------------------
-   !!           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))]
+   !!           Wave mixing vertical parameterization
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0  ! Craig and Banner formulation (original NEMO formulation -
+                                                       !    direct conversion of mechanical to turbulent energy)
+   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1  ! Janssen formulation - no assumption on direct energy conversion 
 
    !!----------------------------------------------------------------------
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 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7905)
@@ -168,5 +168,5 @@
    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, &
+   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrft,sn_rcv_wper, &
                     sn_rcv_wfreq,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_wdrag
    ! Other namelist parameters                        !
@@ -244,8 +244,7 @@
          &                  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_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
+         &                  sn_rcv_sdrft, 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
       !!---------------------------------------------------------------------
       !
@@ -290,6 +289,5 @@
          WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'  
          WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'  
-         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'  
-         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'  
+         WRITE(numout,*)'      Surface Stokes drift u,v        = ', TRIM(sn_rcv_sdrft%cldes ), ' (', TRIM(sn_rcv_sdrft%clcat ), ')'  
          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  ), ')'  
@@ -549,12 +547,9 @@
       ENDIF 
       srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction 
-      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN 
+      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction 
+      IF( TRIM(sn_rcv_sdrft%cldes ) == 'coupled' )  THEN 
          srcv(jpr_sdrftx)%laction = .TRUE. 
-         cpl_sdrftx = .TRUE. 
-      ENDIF 
-      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction 
-      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN 
          srcv(jpr_sdrfty)%laction = .TRUE. 
-         cpl_sdrfty = .TRUE. 
+         cpl_sdrft = .TRUE. 
       ENDIF 
       srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period 
@@ -964,5 +959,4 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      USE zdf_oce,  ONLY : ln_zdfqiao
       USE sbcflx ,  ONLY : ln_shelf_flx
 
@@ -1171,12 +1165,10 @@
       IF( ln_sdw ) THEN  ! Stokes Drift correction activated 
       !                                                      ! ========================= !  
-      !                                                      !       Stokes drift u      ! 
+      !                                                      !     Stokes drift u,v      ! 
       !                                                      ! ========================= !  
-         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
-      ! 
-      !                                                      ! ========================= !  
-      !                                                      !       Stokes drift v      ! 
-      !                                                      ! ========================= !  
-         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
+         IF( srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction ) THEN
+                                        ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
+                                        vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
+         ENDIF
       ! 
       !                                                      ! ========================= !  
@@ -1201,6 +1193,6 @@
       
          ! 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_wfreq)%laction) & 
+         IF( (srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction) .OR. srcv(jpr_wper)%laction & 
+                                        .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) & 
             CALL sbc_stokes() 
       ENDIF 
@@ -1215,4 +1207,6 @@
       IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN
          rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1)
+         WHERE( rn_crban <    0.0 ) rn_crban = 0.0
+         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0
       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 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7905)
@@ -88,7 +88,6 @@
       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
-         &             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, nn_sdrift
+         &             ln_ssr    , nn_isf    , nn_fwb, ln_wave    , nn_lsm     , nn_limflx,   &
+                       nn_components, ln_cpl , ln_wavcpl
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -134,9 +133,4 @@
          WRITE(numout,*) '              forced-coupled wav mixed formulation       ln_wavcpl   = ', ln_wavcpl
          WRITE(numout,*) '              wave physics                               ln_wave     = ', ln_wave 
-         WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw 
-         WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc    = ', ln_tauoc 
-         WRITE(numout,*) '                 wave to ocean energy - wave breaking    ln_phioc    = ', ln_phioc
-         WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor    = ', ln_stcor 
-         WRITE(numout,*) '                 neutral drag coefficient                ln_cdgw     = ', ln_cdgw
          WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
          WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
@@ -153,5 +147,4 @@
          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,*) '              momentum formulation                       nn_drag     = ', nn_drag
       ENDIF
 
@@ -222,29 +215,4 @@
       IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   &
          &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
-      
-      IF ( ln_wave ) THEN
-         !Activated wave module but neither drag nor stokes drift activated 
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc) )   THEN  
-             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F, ln_phioc=F') 
-         ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN  
-             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
-         ENDIF
-         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) & 
-             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3')
-         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) &
-             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3')
-         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 ) &  
-            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
-            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
-            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
-            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
-            &                  'or wave to ocean energy modification (ln_phioc=T) ',           & 
-            &                  'or Stokes-Coriolis term (ln_stcor=T)'  )
-      ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
       ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl .AND. .NOT. ln_wavcpl
@@ -414,6 +382,10 @@
             SELECT CASE( nsbc ) 
             CASE(  0,1,2,3,5,-1 )  ; 
-                IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. & 
-                        & If not requested select ln_tauoc=.false' 
+                IF(lwp .AND. kt == nit000 ) THEN
+                   write(numout,*)
+                   WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. & 
+                                                      & If not requested select ln_tauoc=.false' 
+                   write(numout,*)
+                END IF
             END SELECT 
       ! 
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 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7905)
@@ -18,5 +18,4 @@
    USE oce            ! ocean variables
    USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE zdf_oce,  ONLY : ln_zdfqiao
    USE bdy_oce        ! open boundary condition variables
    USE domvvl         ! domain: variable volume layers
@@ -39,6 +38,5 @@
    LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_sdrft  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
@@ -46,4 +44,9 @@
    LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
+
+   INTEGER ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
+   INTEGER, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
+   INTEGER, 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, PARAMETER ::   jp_peakph   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale
 
    INTEGER ::   jpfld    ! number of files to read for stokes drift
@@ -103,9 +106,8 @@
       !
       ! select parameterization for the calculation of vertical Stokes drift
-      SELECT CASE ( nn_sdrift )
-      !
-      CASE ( jp_breivik )
+      ! 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               ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) )
+         DO jj = 1, jpj               
             DO ji = 1, jpi
                ! Stokes drift velocity estimated from Hs and Tmean
@@ -126,6 +128,18 @@
             END DO
          END DO
-         !
-         !                       !==  horizontal Stokes Drift 3D velocity  ==!
+      ELSE IF( nn_sdrift==jp_peakph ) 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  ==!
+      IF( nn_sdrift==jp_breivik ) THEN
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -145,16 +159,5 @@
             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  ==!
+      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -165,7 +168,7 @@
                   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))
+                  !                                ! 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)
@@ -174,5 +177,5 @@
             END DO
          END DO
-      END SELECT
+      ENDIF
 
       CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. )
@@ -265,9 +268,10 @@
          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 < -1000.0 ) rn_crban = 0.0
-         WHERE( rn_crban >  1000.0 ) rn_crban = 0.0
-      ENDIF
-
-      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==! 
+         WHERE( rn_crban > 1.e8   ) rn_crban = 0.0    !remove first mask mistmatch points, then cap values in case of low friction velocity
+         WHERE( rn_crban < 0.0    ) rn_crban = 0.0
+         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0
+      ENDIF
+
+      IF( ln_sdw .OR. ln_rough )  THEN             !==  Computation of the 3d Stokes Drift  ==! 
          !
          IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
@@ -285,6 +289,6 @@
             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
+               WHERE( wfreq <    0.0 ) wfreq = 0.0 
+               WHERE( wfreq >  100.0 ) wfreq = 0.0
             ENDIF
             IF( jp_usd > 0 ) THEN
@@ -299,5 +303,7 @@
             ENDIF
          ENDIF
-         !
+      ENDIF
+      !
+      IF( ln_sdw ) THEN
          ! Read also wave number if needed, so that it is available in coupling routines
          IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
@@ -308,6 +314,7 @@
          !                                         !==  Computation of the 3d Stokes Drift  ==! 
          !
-         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) ) &
+         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_peakph .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
@@ -341,5 +348,7 @@
                              &   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
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc, sn_phioc, &
+                             ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_zdfqiao, ln_rough,                 &
+                             nn_drag, nn_sdrift, nn_wmix
       !!---------------------------------------------------------------------
       !
@@ -353,4 +362,53 @@
       IF(lwm) WRITE ( numond, namsbc_wave )
       !
+      IF(lwp) THEN               !* Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave_init: wave physics'
+         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) '   Namelist namsbc_wave : set wave physics parameters'
+         WRITE(numout,*) '      Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw 
+         WRITE(numout,*) '        vertical parametrization              nn_sdrift   = ', nn_sdrift 
+         WRITE(numout,*) '      Stokes coriolis term                    ln_stcor    = ', ln_stcor 
+         WRITE(numout,*) '      wave modified ocean stress              ln_tauoc    = ', ln_tauoc 
+         WRITE(numout,*) '      wave to ocean energy                    ln_phioc    = ', ln_phioc
+         WRITE(numout,*) '        vertical mixing parametrization       nn_wmix     = ', nn_wmix 
+         WRITE(numout,*) '      neutral drag coefficient                ln_cdgw     = ', ln_cdgw
+         WRITE(numout,*) '        momentum formulation                  nn_drag     = ', nn_drag
+         WRITE(numout,*) '      wave roughness length modification      ln_rough    = ', ln_rough 
+         WRITE(numout,*) '      Qiao vertical mixing formulation        ln_zdfqiao  = ', ln_zdfqiao
+      ENDIF
+
+      IF ( ln_wave ) THEN
+         ! Activated wave physics but no wave physics components activated 
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao) )   THEN  
+             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F, ln_phioc=F ', &
+                                                      'ln_rough=F, ln_zdfqiao=F' ) 
+         ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN  
+             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+         ENDIF
+         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) & 
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3')
+         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) &
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3')
+         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_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) & 
+            CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' )
+         IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) & 
+            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
+         IF( ln_zdfqiao .AND. .NOT.ln_sdw ) & 
+            CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' )
+      ELSE
+         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) &  
+            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
+            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
+            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
+            &                  'or Stokes-Coriolis term (ln_stcor=T)',                         &
+            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
+            &                  'or wave to ocean energy modification (ln_phioc=T) ',           & 
+            &                  'or wave surface roughness (ln_rough=T) ',                      & 
+            &                  'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' )
+      ENDIF 
+      !
       IF( ln_cdgw ) THEN
          IF( .NOT. cpl_wdrag ) THEN
@@ -360,5 +418,5 @@
                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
             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_init', 'Wave module ', 'namsbc_wave' )
+            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( cdn_wave(jpi,jpj) )
@@ -372,5 +430,5 @@
                                     ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
             IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
+            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( tauoc_wave(jpi,jpj) )
@@ -384,53 +442,60 @@
                                     ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   )
             IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
+            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( rn_crban(jpi,jpj) )
       ENDIF
 
-      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_wfr=0
-         IF( .NOT. cpl_sdrftx ) 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_wfr=0
+      IF( ln_sdw ) THEN
+         IF( .NOT. cpl_sdrft ) THEN
             jpfld  = jpfld + 1
             jp_usd = jpfld
-         ENDIF
-         IF( .NOT. cpl_sdrfty ) THEN
             jpfld  = jpfld + 1
             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 ) THEN
+         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN
             jpfld  = jpfld + 1
             jp_wfr = jpfld
          ENDIF
-
-         ! Read from file only the non-coupled fields 
-         IF( jpfld > 0 ) THEN
-            ALLOCATE( slf_i(jpfld) )
-            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
-            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
-            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' )
-            !
-            DO ifpr= 1, jpfld
-               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
-               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
-            END DO
-            !
-            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
-         ENDIF
+      ENDIF
+
+      IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN
+         jpfld  = jpfld + 1
+         jp_hsw = jpfld
+      ENDIF
+
+      ! Read from file only the non-coupled fields 
+      IF( jpfld > 0 ) THEN
+         ALLOCATE( slf_i(jpfld) )
+         IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
+         IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
+         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' )
+         !
+         DO ifpr= 1, jpfld
+            ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
+            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
+         END DO
+         !
+         CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+      ENDIF
+
+      IF( ln_sdw ) THEN
          ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
-         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     )
+         ALLOCATE( wmp  (jpi,jpj)  )
          ALLOCATE( wfreq (jpi,jpj) )
          ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
@@ -446,7 +511,11 @@
                                    ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
             IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( wnum(jpi,jpj) )
+      ENDIF
+
+      IF( ln_sdw .OR. ln_rough ) THEN
+         ALLOCATE( hsw  (jpi,jpj) )
       ENDIF
       !
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7905)
@@ -10,4 +10,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
+   USE sbc_oce, ONLY : ln_zdfqiao
 
    IMPLICIT NONE
@@ -35,5 +36,4 @@
    INTEGER , PUBLIC ::   nn_npc      !: non penetrative convective scheme call  frequency
    INTEGER , PUBLIC ::   nn_npcp     !: non penetrative convective scheme print frequency
-   LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag
 
 
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7905)
@@ -24,5 +24,5 @@
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer
-   USE sbcwave, ONLY: hsw,rn_crban
+   USE sbcwave
    ! 
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -62,8 +62,4 @@
    INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
    INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
-   INTEGER  ::   nn_wmix           ! type of wave breaking mixing
-   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0   ! Craig and Banner formulation (original NEMO formulation -
-                                                        !    direct conversion of mechanical to turbulent energy)
-   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1   ! Janssen formulation - no assumption on direct energy conversion 
    REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267
    REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3)
@@ -181,5 +177,5 @@
             rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking
             rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking
-         ELSE
+         ELSE IF( nn_wmix==jp_craigbanner ) THEN
             rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf
             rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf
@@ -359,5 +355,5 @@
                z_elem_c(:,:,2) = 0._wp
                z_elem_b(:,:,2) = 1._wp
-            ELSE
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN
                en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)
                en(:,:,1) = MAX(en(:,:,1), rn_emin) 
@@ -405,5 +401,5 @@
                zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)
                en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-            ELSE
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN
                ! Dirichlet conditions at k=1
                en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)
@@ -652,5 +648,5 @@
                ! 
                !
-            ELSE
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN
                ! Surface value
                zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
@@ -708,5 +704,5 @@
                psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
                !
-            ELSE
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN
                ! Surface value: Dirichlet
                zdep(:,:)       = zhsro(:,:) * rl_sf
@@ -1047,5 +1043,5 @@
          &            rn_crban_default, rn_charn, rn_frac_hs,&
          &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
-         &            nn_stab_func, nn_clos, nn_wmix
+         &            nn_stab_func, nn_clos
       !!----------------------------------------------------------
       !
@@ -1090,5 +1086,9 @@
       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
       IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 
-      IF( nn_z0_met == 3 .AND. .NOT.ln_sdw ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_sdw=T' ) 
+      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' ) 
+      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN
+         CALL ctl_warn('W A R N I N G:  ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3')
+         nn_z0_met = 3
+      ENDIF
       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7904)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7905)
@@ -54,6 +54,5 @@
       !!
       NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  & 
-         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp,       & 
-         &        ln_zdfqiao 
+         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp
       !!----------------------------------------------------------------------
 
@@ -84,5 +83,4 @@
          WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc
          WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp
-         WRITE(numout,*) '      Qiao formulation flag               ln_zdfqiao=', ln_zdfqiao
       ENDIF
 
