Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7852)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7853)
@@ -1017,4 +1017,7 @@
    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 standard formulation (original NEMO formulation) 
+           !                             ! = 1 Craig and Banner formulation 
 /
 !-----------------------------------------------------------------------
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 7852)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7853)
@@ -226,5 +226,5 @@
       IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==!
          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(rw/ra)sbc_wa
+         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
@@ -235,8 +235,24 @@
          IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
             CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
-            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_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
+            IF( jp_hsw > 0 ) THEN
+               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
+               WHERE( hsw > 100.0 ) hsw = 0.0
+               WHERE( hsw <   0.0 ) hsw = 0.0
+            ENDIF
+            IF( jp_wmp > 0 ) THEN
+               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+               WHERE( wmp > 100.0 ) wmp = 0.0
+               WHERE( wmp <   0.0 ) wmp = 0.0
+            ENDIF
+            IF( jp_usd > 0 ) THEN
+               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
+               WHERE( ut0sd < -100.0 ) ut0sd = 1.0
+               WHERE( ut0sd >  100.0 ) ut0sd = 1.0
+            ENDIF
+            IF( jp_vsd > 0 ) THEN
+               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
+               WHERE( vt0sd < -100.0 ) vt0sd = 1.0
+               WHERE( vt0sd >  100.0 ) vt0sd = 1.0
+            ENDIF
          ENDIF
          !
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 7852)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7853)
@@ -52,4 +52,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_tke3
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_psi1
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rtrans
 
    !                              !! ** Namelist  namzdf_gls  **
@@ -61,4 +62,8 @@
    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)
@@ -96,5 +101,5 @@
    REAL(wp) ::   rm7           =  0.0_wp
    REAL(wp) ::   rm8           =  0.318_wp
-   REAL(wp) ::   rtrans        =  0.1_wp
+   REAL(wp) ::   rtrans_default =  0.1_wp
    REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
    REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        -
@@ -124,5 +129,5 @@
          &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     &
          &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     &
-         &      rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc )
+         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
@@ -172,7 +177,13 @@
       ! variable initialization
       IF( ln_phioc ) THEN
-         rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking 
-         rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking
-         rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking
+         IF( nn_wmix==jp_janssen ) THEN
+            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking 
+            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking
+            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking
+         ELSE
+            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf
+            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf
+            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp )
+         ENDIF
       ENDIF
 
@@ -336,15 +347,32 @@
       CASE ( 0 )             ! Dirichlet case
          IF( ln_phioc ) THEN  ! wave induced mixing case with forced/coupled fields
-            ! First level
-            en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,1) = en(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-
-            ! One level below
-            en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
-            z_elem_a(:,:,2) = 0._wp
-            z_elem_c(:,:,2) = 0._wp
-            z_elem_b(:,:,2) = 1._wp
+            IF( nn_wmix==jp_janssen ) THEN
+               ! First level
+               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )
+               z_elem_a(:,:,1) = en(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+
+               ! One level below
+               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
+               z_elem_a(:,:,2) = 0._wp
+               z_elem_c(:,:,2) = 0._wp
+               z_elem_b(:,:,2) = 1._wp
+            ELSE
+               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)
+               en(:,:,1) = MAX(en(:,:,1), rn_emin) 
+               z_elem_a(:,:,1) = en(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               ! 
+               ! One level below
+               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) &
+                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
+               en(:,:,2) = MAX(en(:,:,2), rn_emin )
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp
+               z_elem_b(:,:,2) = 1._wp
+               !
+            ENDIF
          ELSE                 ! wave induced mixing case with default values
             en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
@@ -366,14 +394,34 @@
       CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
          IF( ln_phioc ) THEN   ! Shear free case: d(e)/dz=Fw with forced/coupled fields
-            ! Dirichlet conditions at k=1
-            en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,1) = en(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            ! at k=2, set de/dz=Fw
-            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-            z_elem_a(:,:,2) = 0._wp
-            zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)
-            en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+            IF( nn_wmix==jp_janssen ) THEN
+               ! Dirichlet conditions at k=1
+               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )
+               z_elem_a(:,:,1) = en(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               ! at k=2, set de/dz=Fw
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+               z_elem_a(:,:,2) = 0._wp
+               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)
+               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+            ELSE
+               ! Dirichlet conditions at k=1
+               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)
+               en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
+               z_elem_a(:,:,1) = en(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               !
+               ! at k=2, set de/dz=Fw
+               !cbr
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+               z_elem_a(:,:,2) = 0._wp
+               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) ))
+               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) &
+                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
+
+               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
+               !
+            ENDIF
          ELSE                  ! Shear free case: d(e)/dz=Fw with default values
             ! Dirichlet conditions at k=1
@@ -388,5 +436,5 @@
             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
             z_elem_a(:,:,2) = 0._wp
-            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) ))
+            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) ))
             zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
                  &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
@@ -586,21 +634,40 @@
                                ! balance between the production and the
                                ! dissipation terms including the wave effect
-            ! Surface value
-            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
-            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,1) = psi(:,:,1)      
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! One level below
-            zex1 = (rmm*ra_sf+rnn)
-            zex2 = (rmm*ra_sf)
-            zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2
-            psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)
-            z_elem_a(:,:,2) = 0._wp
-            z_elem_c(:,:,2) = 0._wp
-            z_elem_b(:,:,2) = 1._wp
-            ! 
-            !
+            IF( nn_wmix==jp_janssen ) THEN
+               ! Surface value
+               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+               z_elem_a(:,:,1) = psi(:,:,1)      
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               !
+               ! One level below
+               zex1 = (rmm*ra_sf+rnn)
+               zex2 = (rmm*ra_sf)
+               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2
+               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)
+               z_elem_a(:,:,2) = 0._wp
+               z_elem_c(:,:,2) = 0._wp
+               z_elem_b(:,:,2) = 1._wp
+               ! 
+               !
+            ELSE
+               ! Surface value
+               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+               z_elem_a(:,:,1) = psi(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               !
+               ! One level below
+               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) )))
+               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
+               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+               z_elem_a(:,:,2) = 0._wp
+               z_elem_c(:,:,2) = 0._wp
+               z_elem_b(:,:,2) = 1._wp
+               ! 
+               !
+            ENDIF
          ELSE                  ! Wave induced mixing case with default values
             ! Surface value
@@ -612,5 +679,5 @@
             !
             ! One level below
-            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) )))
+            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) )))
             zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
             psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
@@ -623,21 +690,45 @@
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
          IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields
-            !
-            zdep(:,:) = rl_sf * zhsro(:,:)
-            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,1) = psi(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! Neumann condition at k=2
-            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-            z_elem_a(:,:,2) = 0._wp
-            !
-            ! Set psi vertical flux at the surface:
-            zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)
-            zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &
-               &                   * en(:,:,1)**rmm * zdep
-            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-            !
+            IF( nn_wmix==jp_janssen ) THEN
+               !
+               zdep(:,:) = rl_sf * zhsro(:,:)
+               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+               z_elem_a(:,:,1) = psi(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               !
+               ! Neumann condition at k=2
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+               z_elem_a(:,:,2) = 0._wp
+               !
+               ! Set psi vertical flux at the surface:
+               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)
+               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &
+                  &                   * en(:,:,1)**rmm * zdep
+               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+               !
+            ELSE
+               ! Surface value: Dirichlet
+               zdep(:,:)       = zhsro(:,:) * rl_sf
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+               z_elem_a(:,:,1) = psi(:,:,1)
+               z_elem_c(:,:,1) = 0._wp
+               z_elem_b(:,:,1) = 1._wp
+               !
+               ! Neumann condition at k=2
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+               z_elem_a(:,:,2) = 0._wp
+               !
+               ! Set psi vertical flux at the surface:
+               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
+               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
+               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * &
+                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
+               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
+                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)
+               zflxs(:,:) = zdep(:,:) * zflxs(:,:)
+               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+               !   
+            ENDIF
          ELSE                 ! Wave induced mixing case with default values
             ! Surface value: Dirichlet
@@ -653,5 +744,5 @@
             !
             ! Set psi vertical flux at the surface:
-            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
+            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
             zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
             zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * &
@@ -956,5 +1047,5 @@
          &            rn_crban_default, rn_charn, rn_frac_hs,&
          &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
-         &            nn_stab_func, nn_clos
+         &            nn_stab_func, nn_clos, nn_wmix
       !!----------------------------------------------------------
       !
@@ -1215,5 +1306,5 @@
       rsbc_tke2         = rdt * rn_crban_default / rl_sf
       zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans            = 0.2_wp / zcr
+      rtrans_default    = 0.2_wp / zcr
       rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
       rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
