Index: /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7808)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7809)
@@ -269,4 +269,5 @@
    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) ,
@@ -1006,5 +1007,5 @@
    rn_clim_galp  = 0.267   !  galperin limit
    ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case
-   rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux
+   rn_crban_default = 100. !  Craig and Banner 1994 constant for wb tke flux
    rn_charn      = 70000.  !  Charnock constant for wb induced roughness length
    rn_hsro       =  0.02   !  Minimum surface roughness
@@ -1300,4 +1301,5 @@
    sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    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 drag coefficient files
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 7808)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7809)
@@ -68,4 +68,5 @@
    LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
    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 
    !
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 7808)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7809)
@@ -1195,4 +1195,11 @@
       !                                                      ! ========================= !  
       IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1) 
+      
+      !                                                      ! ========================= !  
+      !                                                      !   Wave to ocean energy    ! 
+      !                                                      ! ========================= !  
+      IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN
+         rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1)
+      ENDIF
       
       !  Fields received by SAS when OASIS coupling
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 7808)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7809)
@@ -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_wavcpl , nn_drag
+         &             ln_phioc  , ln_wavcpl , nn_drag
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -136,4 +136,5 @@
          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
@@ -224,16 +225,17 @@
       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 ) )   THEN  
-             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F') 
+         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
       ELSE
-         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) &  
+         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 Stokes-Coriolis term (ln_stcori=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)
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 7808)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7809)
@@ -56,6 +56,8 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy 
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !:
@@ -222,4 +224,11 @@
       ENDIF
 
+      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
+         WHERE( rn_crban <  10.0 ) rn_crban =  10.0
+         WHERE( rn_crban > 300.0 ) rn_crban = 300.0
+      ENDIF
+
       IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==! 
          !
@@ -267,8 +276,8 @@
       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,  &
+      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
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc, sn_phioc
       !!---------------------------------------------------------------------
       !
@@ -285,5 +294,5 @@
          IF( .NOT. cpl_wdrag ) THEN
             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' )
             !
                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
@@ -297,5 +306,5 @@
          IF( .NOT. cpl_tauoc ) THEN
             ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' )
             !
                                     ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
@@ -304,4 +313,16 @@
          ENDIF
          ALLOCATE( tauoc_wave(jpi,jpj) )
+      ENDIF
+
+      IF( ln_phioc ) THEN
+         IF( .NOT. cpl_phioc ) THEN
+            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' )
+            !
+                                    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' )
+         ENDIF
+         ALLOCATE( rn_crban(jpi,jpj) )
       ENDIF
 
@@ -334,5 +355,5 @@
             IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
             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' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' )
             !
             DO ifpr= 1, jpfld
@@ -354,5 +375,5 @@
          IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
             ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' )
                                    ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
             IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
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 7808)
+++ /branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7809)
@@ -24,5 +24,5 @@
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer
-   USE sbcwave ,  ONLY: hsw   ! significant wave height 
+   USE sbcwave, ONLY: hsw,rn_crban
    ! 
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -48,4 +48,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_tke1
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_tke3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_psi1
 
    !                              !! ** Namelist  namzdf_gls  **
@@ -61,7 +65,7 @@
    REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
    REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
-   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
+   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
-   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1) 
 
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
@@ -94,6 +98,6 @@
    REAL(wp) ::   rtrans        =  0.1_wp
    REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
-   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
-   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
+   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        -
+   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        -
    REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
    REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
@@ -117,6 +121,8 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
-         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
+      ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) ,     &
+         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     &
+         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     &
+         &      rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
@@ -163,4 +169,12 @@
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
+
+      ! 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
+      ENDIF
+
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
          avt (:,:,:) = avt_k (:,:,:)
@@ -200,5 +214,9 @@
          zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
-         zhsro(:,:) = hsw(:,:)
+         WHERE( hsw == 0._wp ) ! surface roughness length according to Charnock formula when sign. wave height 0
+            zhsro = MAX(rn_charn / grav * ustars2, rn_hsro)
+         ELSEWHERE
+            zhsro = MAX(hsw, rn_hsro)
+         END WHERE
       END SELECT
 
@@ -317,40 +335,65 @@
       !
       CASE ( 0 )             ! Dirichlet case
-      ! First level
-      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
-      !
-      !
+         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
+         ELSE                 ! wave induced mixing case with default values
+            en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(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_default * ((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
       CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
-      !
-      ! 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_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
-           &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
-
-      en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
-      !
-      !
+         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)
+         ELSE                  ! Shear free case: d(e)/dz=Fw with default values
+            ! Dirichlet conditions at k=1
+            en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(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_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
+                 &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
+
+            en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
+            !
+            !
+         ENDIF
       END SELECT
 
@@ -539,45 +582,86 @@
       !
       CASE ( 0 )             ! Dirichlet boundary conditions
-      !
-      ! 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
-      ! 
-      !
+         IF( ln_phioc ) THEN   ! Wave induced mixing case
+                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
+                               ! 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
+            ! 
+            !
+         ELSE                  ! Wave induced mixing case with default values
+            ! 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
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
-      !
-      ! 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)
-
-      !   
-      !
+         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)
+            !
+         ELSE                 ! Wave induced mixing case with default values
+            ! 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_default * (rnn + rmm*ra_sf) * zdep(:,:)) * &
+                       (1._wp + rsbc_tke1_default * 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
       END SELECT
 
@@ -870,5 +954,5 @@
       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
          &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
-         &            rn_crban, rn_charn, rn_frac_hs,        &
+         &            rn_crban_default, rn_charn, rn_frac_hs,&
          &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
          &            nn_stab_func, nn_clos
@@ -898,5 +982,5 @@
          WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
          WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
-         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
+         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default
          WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
          WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
@@ -915,5 +999,5 @@
       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_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 
+      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_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' )
@@ -1089,5 +1173,5 @@
                &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
 
-      IF ( rn_crban==0._wp ) THEN
+      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN
          rl_sf = vkarmn
       ELSE
@@ -1128,8 +1212,8 @@
       rc03  = rc02 * rc0
       rc04  = rc03 * rc0
-      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
-      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
-      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf
+      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
       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 
