Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 7807)
@@ -297,6 +297,4 @@
               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
               ln_tra_dmp(ib_bdy)=.false.
-           ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
-              CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7807)
@@ -767,6 +767,6 @@
 
       IF( lk_vvl ) THEN
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T  )   ! salt content
          CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 7807)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 7807)
@@ -0,0 +1,126 @@
+MODULE dynstcor
+!! PHILLIPPS PARAMETRIZATION
+   !!======================================================================
+   !!                       ***  MODULE  dynstcor  ***
+   !! Ocean dynamics: Stokes-Coriolis effect.
+   !!
+   !!======================================================================
+   !! History :  0.1  !  2012-10  (Oyvind Breivik)
+   !!            0.2  !  2016-06  (oyvind.breivik@met.no)
+   !!                 !           Changed Stokes profile to the Phillips approximation 
+   !!                        Ref: Breivik, O, J-R Bidlot, P A Janssen (2016). 
+   !!                             A Stokes drift approximation based on the Phillips spectrum, 
+   !!                             Ocean Model, 100, pp 49-56, doi:10.1016/j.ocemod.2016.01.005
+   !!            0.2.1! 2016-12   WAVE2NEMO: adapted to NEMO3.6, this is the PHILLIPS approximation
+   !!----------------------------------------------------------------------
+   !! dyn_stcor        : Add the Stokes-Coriolis forcing to the momentum equation
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+!  USE obc_oce          ! ocean lateral open boundary condition !! WAVE2NEMO: does not exist in NEMO36
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! distributed memory computing
+   USE prtctl           ! Print control
+   USE phycst
+   USE lbclnk
+   USE wrk_nemo         ! Memory Allocation
+   USE sbcmod           ! Access to ln_stcor (sbc_oce) and wave parameters (sbc_wave)
+   USE sbcwave          ! Wave module
+
+   IMPLICIT NONE
+
+   REAL(wp) :: rn_deptmaxstcor = 150.0_wp ! maximum depth [m] to be affected by Stokes-Coriolis effect
+   !PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dyn_stcor           ! routine called by step.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ustc, vstc  ! Stokes-Coriolis u and v
+
+   !! * Module variables
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
+   !!----------------------------------------------------------------------
+
+ CONTAINS
+
+   INTEGER FUNCTION dynstcor_alloc()
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dynstcor_alloc  ***
+      !!
+      !! History :  1.0  !  2012-10  (Oyvind Breivik)
+      !!----------------------------------------------------------------------
+      ALLOCATE( ustc(jpi,jpj,jpk) , vstc(jpi,jpj,jpk) , STAT=dynstcor_alloc )
+         !
+      IF( dynstcor_alloc /= 0 )   CALL ctl_warn('dynstcor_alloc: array allocate failed.')
+   END FUNCTION dynstcor_alloc
+
+   SUBROUTINE dyn_stcor( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_stcor ***
+      !!
+      !! ** Purpose:  Add Stokes-Coriolis forcing to horizontal momentum equation.
+      !!
+      !! ** History:  0.1  !   2012-10  oyvind.breivik@ecmwf.int
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt  ! ocean time-step index
+      !!
+      INTEGER  ::  ji, jj, jk        ! dummy loop indices
+      REAL(wp) ::  ztransp, zsp0, zk, zfac
+      REAL(wp) ::  zus, zvs
+      !
+
+      ! Allocation at first time step.
+
+      IF ( kt == nit000 ) THEN
+         IF( dynstcor_alloc() /= 0 ) CALL ctl_stop('dyn_stcor_phil: array allocate failed.')
+      ENDIF
+
+      !!----------------------------------------------------------------------
+      !
+      ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities ustc, vstc
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! Skip deep levels where the Stokes-Coriolis effect is negligible
+               IF (fsdept(ji,jj,jk) <= rn_deptmaxstcor) THEN
+                  ! Stokes transport speed estimated from Hs and Tmean
+                  ztransp = 2.0_wp*rpi*swh_wavepar(ji,jj)**2.0_wp/(16.0_wp*MAX(mwp_wavepar(ji,jj),0.0000001_wp))
+
+                  ! Stokes surface speed
+                  zsp0 = SQRT(ust_wavepar(ji,jj)**2 + vst_wavepar(ji,jj)**2)
+
+                  ! Wavenumber scale
+                  !zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
+                  zk = (1.0_wp-2.0_wp/3.0_wp)*zsp0/MAX(2.0_wp*ztransp,0.0000001_wp)
+
+                  ! Depth attenuation
+                  !zfac = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk))
+                  zfac = EXP(-2.0*zk*fsdept(ji,jj,jk)) - SQRT(2.0_wp*rpi*zk*fsdept(ji,jj,jk))*ERFC(SQRT(2.0_wp*zk*fsdept(ji,jj,jk)))
+
+                  ! The Stokes-Coriolis forcing
+                  zus =  ff(ji,jj)*vst_wavepar(ji,jj)*zfac
+                  zvs = -ff(ji,jj)*ust_wavepar(ji,jj)*zfac
+
+                  ! Store arrays of tendencies for diagnostics output
+                  ! This may be removed later for efficiency
+                  ustc(ji,jj,jk) = zus
+                  vstc(ji,jj,jk) = zvs
+
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + zus * umask(ji,jj,jk)
+                  va(ji,jj,jk) = va(ji,jj,jk) + zvs * vmask(ji,jj,jk)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+   END SUBROUTINE dyn_stcor
+
+
+END MODULE dynstcor
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 7807)
@@ -62,7 +62,14 @@
    !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
    !                                             !:  = 2 annual global mean of e-p-r set to zero
-   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_wave        !: true if some coupling with wave model
+   LOGICAL , PUBLIC ::   ln_cdgw = .FALSE.        !: true if neutral drag coefficient from wave model
+  ! LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
+         ! WAVE2NEMO 11.12.2016
+      LOGICAL , PUBLIC ::   ln_wavetke  = .FALSE.   !: true if wave parametersare read
+      LOGICAL , PUBLIC ::   ln_stcor    = .FALSE.   !: true if Stokes-Coriolisforcing is included 
+      LOGICAL , PUBLIC ::   ln_tauoc    = .FALSE.   !: true if wave-modifiedwater-side stress is used
+      ! END  WAVE2NEMO
+   
+
    !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
@@ -125,5 +132,11 @@
 #endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
-
+   ! START WAVE2NEMO 26.11.2016
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tauoc_wavepar !:normalized stress into the ocean from wave model  CCC
+   REAL(wp), PUBLIC,ALLOCATABLE,SAVE,DIMENSION (:,:)       :: wspd_wavepar  !CCCmoved to sbc_oce   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_wave !:wave model drag coefficient [N/m2]
+  ! END 
+!WAVE2NEMO (14.12.2016) rn_crban moved to here, it needs to be a more global variable
+   REAL(wp),PUBLIC,ALLOCATABLE,SAVE, DIMENSION (:,:)  ::  rn_crban
    !!----------------------------------------------------------------------
    !!                     Sea Surface Mean fields
@@ -176,5 +189,12 @@
       ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
 #endif
-         !
+      ! Start WAVE2NEMO 11.12.2016
+      ALLOCATE( cdn_wave(jpi,jpj) )
+      ALLOCATE( wspd_wavepar(jpi,jpj) )
+      ! End
+   !
+          ALLOCATE( rn_crban(jpi,jpj) ) !WAVE2NEMO (13.12.2016)
+
+!
       sbc_oce_alloc = MAXVAL( ierr )
       IF( lk_mpp            )   CALL mpp_sum ( sbc_oce_alloc )
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 7807)
@@ -17,5 +17,6 @@
    !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk
    !!----------------------------------------------------------------------
-
+   !!            3.6     2016-06  (V.Alari) Added wave dependent drag
+   !coefficient and ocean side stress
    !!----------------------------------------------------------------------
    !!   sbc_blk_core    : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
@@ -41,5 +42,6 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE sbcwave, ONLY   :  cdn_wave ! wave module
+!   USE sbcwave, ONLY   :  cdn_wave ! wave module (V.Alari 27.06.2016) We get
+!   info from sbc_oce
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE lib_fortran     ! to use key_nosignedzero
@@ -183,4 +185,7 @@
          !
          lhftau = ln_taudif                        ! do we use HF tau information?
+
+
+
          jfld = jpfld - COUNT( (/.NOT. lhftau/) )
          !
@@ -192,4 +197,5 @@
          END DO
          !                                         ! fill sf with slf_i and control print
+
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
@@ -258,4 +264,6 @@
       REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height
       !!---------------------------------------------------------------------
+      REAL(wp) ::   ztheta               ! local variable, wind direction,V.Alari 27.06.2016
+
       !
       IF( nn_timing == 1 )  CALL timing_start('blk_oce_core')
@@ -285,4 +293,25 @@
       END DO
 #endif
+        ! Start 27.06.2016 (V.Alari)
+           ! Wind vector at T points
+      IF ( ln_cdgw ) THEN
+         ! Use neutral 10-m wind speed from wave model if available
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               ! Local wind direction from non-neutral 10-m wind vector
+               ! (relative to grid and no current)
+               ! since we do not have the directional information from the wave
+               ! model
+               ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1),sf(jp_wndj)%fnow(ji,jj,1))
+               ! Wind vector magnitude from 10-m neutral wind speed from wave
+               ! model
+               zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta)
+               zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta)
+               ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0
+               zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj  )+pu(ji,jj)))
+               zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji  ,jj-1)+pv(ji,jj)))
+            END DO
+         END DO
+      ELSE ! If no neutral wind from wave model
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
@@ -291,4 +320,6 @@
          END DO
       END DO
+    ENDIF ! End (V.Alari)
+
       CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
       CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
@@ -320,4 +351,18 @@
     
       ! ... tau module, i and j component
+! Start 27.06.2016 (V.Alari)
+     IF ( ln_cdgw ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj)
+               zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj)
+               !taum(ji,jj) = tauoc_wavepar(ji,jj)
+               !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj)
+               taum  (ji,jj) = zztmp * wndm  (ji,jj)
+               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+            ENDDO
+         ENDDO
+      ELSE
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -328,5 +373,5 @@
          END DO
       END DO
-
+     ENDIF ! End, (V.Alari)
       ! ... add the HF tau contribution to the wind stress module?
       IF( lhftau ) THEN
@@ -349,5 +394,13 @@
       CALL lbc_lnk( vtau(:,:), 'V', -1. )
 
-    
+        ! START: 27.06.2015 (V. Alari). Added tauoc
+
+         IF (ln_tauoc) THEN
+         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:)
+         ENDIF
+
+        ! END  (V.Alari)
       !  Turbulent fluxes over ocean
       ! -----------------------------
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90	(revision 7807)
@@ -24,5 +24,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE sbcwave,ONLY : cdn_wave !wave module
+!   USE sbcwave,ONLY : cdn_wave !wave module
 
    IMPLICIT NONE
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 7807)
@@ -22,4 +22,5 @@
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE sbcwave         ! Wave parameters
 
    IMPLICIT NONE
@@ -152,4 +153,9 @@
                zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
                zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+               IF (ln_cdgw) THEN
+                  zcdrag = 1.5e-3
+                  IF (cdn_wave(ji,jj)>=0.) zcdrag = cdn_wave(ji,jj)
+                  zcoef = 1./(zrhoa*zcdrag)
+               ENDIF
                taum(ji,jj) = zmod
                wndm(ji,jj) = SQRT( zmod * zcoef )
@@ -158,4 +164,10 @@
          taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
          CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
+
+         IF (ln_tauoc) THEN
+         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:)
+         ENDIF
 
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7807)
@@ -85,9 +85,11 @@
       !!----------------------------------------------------------------------
       INTEGER ::   icpt   ! local integer
-      !!
+      !! V.Alari 26.06.2016 Added logicals for waves, deleted ln_wave ln_sdw, which came
+      !with UKMO
       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   ,   &
-         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
+         &             ln_ssr    , nn_isf    , nn_fwb,                                        &
+         &             nn_lsm    , nn_limflx , nn_components, ln_cpl,ln_stcor, ln_wavetke,    &
+         &             ln_tauoc, ln_cdgw
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
@@ -120,5 +122,5 @@
           nn_ice      =   0
       ENDIF
-
+     ! V.Alari added wave logicalc (26.06.2016)
       IF(lwp) THEN               ! Control print
          WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
@@ -146,5 +148,11 @@
          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
-      ENDIF
+         WRITE(numout,*) '         Stokes-Coriolis forcing from wavemodel ln_stcor       = ', ln_stcor
+         WRITE(numout,*) '              TKE wave breaking BC from wave model ln_wavetke        = ', ln_wavetke
+         WRITE(numout,*) '              Wave-modified stress from wave model ln_tauoc          = ', ln_tauoc
+         WRITE(numout,*) '              Drag coefficient from wave model ln_cdgw              = ', ln_cdgw
+ 
+
+     ENDIF
 
       ! LIM3 Multi-category heat flux formulation
@@ -179,5 +187,5 @@
 
       !                          ! Checks:
-      IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf 
+      IF( nn_isf .EQ. 0 ) THEN                      ! no specific treatment in vicinity of ice shelf 
          IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
          fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
@@ -214,18 +222,20 @@
       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
+   
+     ! V.Alari comment out UKMO wave stuff   
+    !  IF ( ln_wave ) THEN
       !Activated wave module but neither drag nor stokes drift activated
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
-            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
+     !    IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
+      !      CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
       !drag coefficient read from wave model definable only with mfs bulk formulae and core 
-         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
-             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
-         ENDIF
-      ELSE
-      IF ( ln_cdgw .OR. ln_sdw  )                                         & 
-         &   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) ')
-      ENDIF 
+       !  ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
+        !     CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+       !  ENDIF
+     ! ELSE
+     ! IF ( ln_cdgw .OR. ln_sdw  )                                         & 
+      !   &   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) ')
+     ! ENDIF 
+      ! End V.Alari
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
       ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
@@ -352,6 +362,17 @@
       IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
       !                                                        ! averaged over nf_sbc time-step
-
-      IF (ln_wave) CALL sbc_wave( kt )
+            ! START V.Alari 26.06.2016
+     IF (ln_cdgw) CALL sbc_wave( kt )
+      ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for
+      ! TKE are needed
+      ! Also read it if drag coefficient to ensure that we use the neutral 10m
+      ! from WAM.
+      IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc) THEN
+            CALL sbc_wavepar( kt )
+      ENDIF
+
+       ! END
+
+!      IF (ln_wave) CALL sbc_wave( kt )
                                                    !==  sbc formulation  ==!
                                                             
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7807)
@@ -4,16 +4,13 @@
    !! Wave module 
    !!======================================================================
-   !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient 
-   !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift 
+   !! History :  3.3.1  !   2011-09  (Adani M)  Original code
    !!----------------------------------------------------------------------
    USE iom             ! I/O manager library
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distribued memory computing library
-   USE fldread	       ! read input fields
-   USE oce
-   USE sbc_oce	       ! Surface boundary condition: ocean fields
-   USE domvvl
-
-   
+   USE fldread         ! read input fields
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
+   ! USE sbcblk_core, ONLY: Cd_n10   
    !!----------------------------------------------------------------------
    !!   sbc_wave       : read drag coefficient from wave model in netcdf files 
@@ -23,18 +20,29 @@
    PRIVATE
 
-   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
-   
-   INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift
-   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd	  ! structure of input fields (file informations, fields read) Drag Coefficient
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd	  ! structure of input fields (file informations, fields read) Stokes Drift
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
-   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
+   PUBLIC   sbc_wavepar    ! routine called in sbc_blk_core or sbc_blk_mfs
+    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
+
+   INTEGER , PARAMETER ::   jpfld_wavepar = 7    ! maximum number of fields to be read 
+   INTEGER , PARAMETER ::   jp_ust = 1           ! index of Stokes velocity (east component) (m/s)  at T-point
+   INTEGER , PARAMETER ::   jp_vst = 2           ! index of Stokes velocity (north component) (m/s) at T-point
+   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave height (m)
+   INTEGER , PARAMETER ::   jp_mwp = 4           ! index of mean wave period (s)
+   INTEGER , PARAMETER ::   jp_phioc = 5         ! index of normalized energy flux into the ocean (non-dim)
+   INTEGER , PARAMETER ::   jp_tauoc = 6         ! index of normalized wave stress into the ocean (non-dim)
+   !INTEGER , PARAMETER ::   jp_cdww = 7          ! index of wave drag coeff (non-dim)
+   INTEGER , PARAMETER ::   jp_wspd = 7          ! index of 10 m neutral wind speed (m/s)
+
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wave ! structure of input fields (file informations, fields read)
+   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave ! CCC moved to sbc_oce
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wavepar ! structure of input fields (file informations, fields read)
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: ust_wavepar, vst_wavepar, swh_wavepar, mwp_wavepar
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: phioc_wavepar  !WAVE2NEMO 14.12.2016 moved rn_crban to sbc_oce
+   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: tauoc_wavepar ! CCC moved to sbc_oce
+   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: wspd_wavepar  ! CCC moved to sbc_oce
+   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdww_wavepar
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
@@ -51,28 +59,14 @@
       !!
       !! ** Method  : - Read namelist namsbc_wave
-      !!              - Read Cd_n10 fields in netcdf files 
-      !!              - Read stokes drift 2d in netcdf files 
-      !!              - Read wave number      in netcdf files 
-      !!              - Compute 3d stokes drift using monochromatic
-      !! ** action  :   
-      !!               
-      !!---------------------------------------------------------------------
-      USE oce,  ONLY : un,vn,hdivn,rotn
-      USE divcur
-      USE wrk_nemo
-#if defined key_bdy
-      USE bdy_oce, ONLY : bdytmask
-#endif
-      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
+      !!              - Read Cd_n10 fields in netcdf files
+      !! ** action  :
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt   ! ocean time step
       INTEGER                ::  ierror   ! return error code
-      INTEGER                ::  ifpr, jj,ji,jk 
-      INTEGER                ::   ios     ! Local integer output status for namelist read
-      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
-      REAL                                          ::  z2dt,z1_2dt
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
-      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
-      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read
-      !!---------------------------------------------------------------------
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
+      CHARACTER(len=100)     ::  cn_dir_cdg                      ! Rootdirectory for location of drag coefficient files
+      TYPE(FLD_N)            ::  sn_cdg                          ! informationsabout the fields to be read
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir_cdg
       !!---------------------------------------------------------------------
 
@@ -83,112 +77,162 @@
       IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
          !                                      ! -------------------- !
-         REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
-         READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
-
-         REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
-         READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp )
-         IF(lwm) WRITE ( numond, namsbc_wave )
-         !
-
-         IF ( ln_cdgw ) THEN
-            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
-            !
-                                   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', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( cdn_wave(jpi,jpj) )
-            cdn_wave(:,:) = 0.0
-        ENDIF
-         IF ( ln_sdw ) THEN
-            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
-            ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave 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', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) )
-            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
-            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
-            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
+         !                                            !* set file information
+         !                                            (default values)
+         ! ... default values (NB: frequency positive => hours, negative =>
+         ! months)
+         !              !   file   ! frequency !  variable  ! time intep !  clim
+         !              ! 'yearly' or ! weights  ! rotation !
+         !              !   name   !  (hours)  !   name     !   (T/F)    !
+         !              (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_cdg = FLD_N('cdg_wave'  ,    1     ,'drag_coeff',  .true.    ,.false. ,   'daily'   , ''       , '',''       )
+         cn_dir_cdg = './'          ! directory in which the Patm data are
+
+
+         REWIND( numnam_ref )                             !* read in namlist namsbc_wave
+         READ  ( numnam_ref, namsbc_wave )
+	 REWIND( numnam_ref )
+         REWIND( numnam_cfg )                             !* read in namlist namsbc_wave
+         READ  ( numnam_cfg, namsbc_wave )
+	 REWIND( numnam_cfg )
+         !
+
+         ALLOCATE( sf_wave(1), STAT=ierror )           !* allocate and fillsf_wave with sn_cdg
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocatesf_wave structure' )
+         !
+         CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave', 'Wavemodule ', 'namsbc_wave' )
+                                ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1)   )
+         IF( sn_cdg%ln_tint )   ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) )
+         !ALLOCATE( cdn_wave(jpi,jpj) )
+         ! Allocation done by sbc_oce
+         cdn_wave(:,:) = 0.0
+      ENDIF
+         !
+         !
+      CALL fld_read( kt, nn_fsbc, sf_wave )      !* read drag coefficient fromexternal forcing
+      cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1)
+   END SUBROUTINE sbc_wave
+
+
+   SUBROUTINE sbc_wavepar( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_wavepar  ***
+      !!
+      !! ** Purpose :   Provide at each time step wave model parameters, including the
+      !!                Stokes drift (east and north components), significant wave height and the
+      !!                mean wave period as well as the normalized stress and energy flux into the 
+      !!                ocean for TKE.
+      !!
+      !!
+      !! ** Method  : - Read namelist namsbc_wavepar
+      !!              - Read fields in NetCDF files 
+      !! ** action  :   
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt  ! ocean time step
+      !!
+      INTEGER  ::  ierror ! return error code
+      INTEGER  ::   ji,iii,jjj      ! dummy loop index, Victor added iii, jjj
+      INTEGER  ::   jfld    ! dummy loop arguments
+      !!
+      CHARACTER(len=100)     ::  cn_dir_wavepar ! Root directory for location of ECWAM wave parameter fields
+      TYPE(FLD_N), DIMENSION(jpfld_wavepar) :: slf_i ! array of namelist informations on the fields to read
+      TYPE(FLD_N)            ::  sn_ust, sn_vst, sn_swh, sn_mwp ! information about the fields to be read
+      TYPE(FLD_N)            ::  sn_phioc, sn_tauoc
+      TYPE(FLD_N)            ::  sn_wspd
+      !TYPE(FLD_N)            ::  sn_cdww
+      !!---------------------------------------------------------------------
+      !NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_phioc, sn_tauoc, cn_dir_wavepar
+      NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_wspd, sn_phioc, sn_tauoc, cn_dir_wavepar
+      !!---------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         !                                      ! -------------------- !
+         ! set file information (default values)
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !            !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !            !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_ust = FLD_N( 'ust'   ,     6     ,  'ust'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''       )
+         sn_vst = FLD_N( 'vst'   ,     6     ,  'vst'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''      )
+         sn_swh = FLD_N( 'swh'   ,     6     ,  'swh'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''      )
+         sn_mwp = FLD_N( 'mwp'   ,     6     ,  'mwp'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''      )
+         !sn_cdww = FLD_N( 'cdww'  ,    6     ,  'cdww'    ,  .true.    , .false. ,   'monthly' , ''       , ''       )
+         sn_wspd = FLD_N( 'wspd'   ,   6     ,  'wspd'    ,  .true.    , .false. ,   'monthly' , ''       , '',''       )
+         sn_phioc = FLD_N( 'phioc'   , 6     ,  'phioc'   ,  .true.    , .false. ,   'monthly' , ''       , '',''       )
+         sn_tauoc = FLD_N( 'tauoc'   , 6     ,  'tauoc'   ,  .true.    , .false. ,   'monthly' , ''       , '',''       )
+         cn_dir_wavepar = './'          ! directory in which the wave data are found
+
+         REWIND( numnam_ref )                             !* read in namlist namsbc_wave
+         READ  ( numnam_ref, namsbc_wavepar ) 
+         REWIND ( numnam_ref )
+
+         REWIND( numnam_cfg )                             !* read in namlist namsbc_wave
+         READ  ( numnam_cfg, namsbc_wavepar ) 
+         REWIND ( numnam_cfg )
+         !
+
+
+         slf_i(jp_ust) = sn_ust 
+         slf_i(jp_vst) = sn_vst
+         slf_i(jp_swh) = sn_swh
+         slf_i(jp_mwp) = sn_mwp
+         !slf_i(jp_cdww) = sn_cdww
+         slf_i(jp_wspd) = sn_wspd
+         slf_i(jp_phioc) = sn_phioc
+         slf_i(jp_tauoc) = sn_tauoc
+
+         ALLOCATE( sf_wavepar(jpfld_wavepar), STAT=ierror )          !* set sf_wavepar structure
+         IF ( ierror > 0 ) THEN
+            CALL ctl_stop( 'STOP', 'sbc_wavepar: unable to allocate sf_wavepar structure' ) ; RETURN
          ENDIF
+         !
+         jfld = jpfld_wavepar
+         DO ji = 1, jpfld_wavepar
+            ALLOCATE( sf_wavepar(ji)%fnow(jpi,jpj,1)   )
+            IF ( slf_i(ji)%ln_tint ) ALLOCATE( sf_wavepar(ji)%fdta(jpi,jpj,1,2) )
+         ENDDO
+         ALLOCATE( ust_wavepar(jpi,jpj) )
+         ALLOCATE( vst_wavepar(jpi,jpj) )
+         ALLOCATE( swh_wavepar(jpi,jpj) )
+         ALLOCATE( mwp_wavepar(jpi,jpj) )
+         !ALLOCATE( cdww_wavepar(jpi,jpj) )
+         !ALLOCATE( wspd_wavepar(jpi,jpj) )
+         ALLOCATE( phioc_wavepar(jpi,jpj) )
+         ALLOCATE( tauoc_wavepar(jpi,jpj) )
+         !
+         CALL fld_fill( sf_wavepar, slf_i, cn_dir_wavepar, 'sbc_wavepar', 'Wave module ', 'namsbc_wavepar' )
+         !
       ENDIF
          !
-         !
-      IF ( ln_cdgw ) THEN
-         CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing
-         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
-      ENDIF
-      IF ( ln_sdw )  THEN
-          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing
-
-         ! Interpolate wavenumber, stokes drift into the grid_V and grid_V
-         !-------------------------------------------------
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
-
-               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
-
-               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) &
-               &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
-            END DO
-         END DO
-
-          !Computation of the 3d Stokes Drift
-          DO jk = 1, jpk
-             DO jj = 1, jpj-1
-                DO ji = 1, jpi-1
-                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk))))
-                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk))))
-                END DO
-             END DO
-             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) )
-             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) )
-          END DO
-
-          CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-          
-          udummy(:,:,:)=un(:,:,:)
-          vdummy(:,:,:)=vn(:,:,:)
-          hdivdummy(:,:,:)=hdivn(:,:,:)
-          rotdummy(:,:,:)=rotn(:,:,:)
-          un(:,:,:)=usd3d(:,:,:)
-          vn(:,:,:)=vsd3d(:,:,:)
-          CALL div_cur(kt)
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-          z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
-
-          z1_2dt = 1.e0 / z2dt
-          DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-             ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-                &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-                &                         * tmask(:,:,jk) * z1_2dt
-#if defined key_bdy
-             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
-#endif
-          END DO
-          hdivn(:,:,:)=hdivdummy(:,:,:)
-          rotn(:,:,:)=rotdummy(:,:,:)
-          vn(:,:,:)=vdummy(:,:,:)
-          un(:,:,:)=udummy(:,:,:)
-          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
-      ENDIF
-   END SUBROUTINE sbc_wave
+      CALL fld_read( kt, nn_fsbc, sf_wavepar )      !* read wave parameters from ECWAM NetCDF file
+      ust_wavepar(:,:) = sf_wavepar(jp_ust)%fnow(:,:,1)
+      vst_wavepar(:,:) = sf_wavepar(jp_vst)%fnow(:,:,1)
+      swh_wavepar(:,:) = sf_wavepar(jp_swh)%fnow(:,:,1)
+      mwp_wavepar(:,:) = sf_wavepar(jp_mwp)%fnow(:,:,1)
+      !cdww_wavepar(:,:) = sf_wavepar(jp_cdww)%fnow(:,:,1)
+      wspd_wavepar(:,:) = sf_wavepar(jp_wspd)%fnow(:,:,1)
+      phioc_wavepar(:,:) = sf_wavepar(jp_phioc)%fnow(:,:,1)
+      tauoc_wavepar(:,:) = sf_wavepar(jp_tauoc)%fnow(:,:,1)
+ 
+     
+       rn_crban(:,:)=29*phioc_wavepar(:,:) ! Alfa is phioc*sqrt(rw/ra)sbc_wa
+       ! Limit cr_ban between 10 and 300
+       DO iii=1,jpi
+           DO jjj=1,jpj
+
+                IF (rn_crban(iii,jjj) < 10 ) THEN
+                rn_crban(iii,jjj)=10                
+                ELSEIF (rn_crban(iii,jjj) > 300) THEN               
+                rn_crban(iii,jjj)=300
+                ENDIF
+        
+           ENDDO
+       ENDDO
+ 
+     END SUBROUTINE sbc_wavepar
       
-   !!======================================================================
+
 END MODULE sbcwave
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7807)
@@ -24,4 +24,6 @@
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer
+   USE sbcwave
+   ! 
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP manager
@@ -47,5 +49,10 @@
    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  **
+   LOGICAL  ::   ln_crban      = .FALSE. !! WAVE2NEMO added ln_crban (=.TRUE. use Craig and Banner scheme)
    LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
    LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
@@ -59,7 +66,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
@@ -92,6 +99,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                !     -           -           -        -
@@ -115,6 +122,7 @@
       !!                ***  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 )
@@ -161,4 +169,12 @@
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
+
+      ! variable initialization
+      IF( ln_crban ) 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 (:,:,:)
@@ -197,5 +213,10 @@
          zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
          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) 
+         WHERE( swh_wavepar == 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(swh_wavepar, rn_hsro)
+         END WHERE
       END SELECT
 
@@ -314,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_crban ) 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_crban ) 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
 
@@ -536,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_crban ) 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_crban ) 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
 
@@ -866,6 +953,6 @@
       !!
       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_clim_galp, ln_crban, ln_sigpsi, rn_hsro,      &
+         &            rn_crban_default, rn_charn, rn_frac_hs,&
          &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
          &            nn_stab_func, nn_clos
@@ -895,5 +982,6 @@
          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,*) '      Craig and Banner scheme                       ln_crban       = ', ln_crban
          WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
          WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
@@ -909,9 +997,10 @@
 
       !                                !* Check of some namelist values
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' )
-      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
-      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
+      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_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_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' )
 
       SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
@@ -1085,5 +1174,5 @@
                &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
 
-      IF ( rn_crban==0._wp ) THEN
+      IF( .NOT. ln_crban .AND. rn_crban_default==0._wp ) THEN
          rl_sf = vkarmn
       ELSE
@@ -1124,8 +1213,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 
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7807)
@@ -323,4 +323,5 @@
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis
                                CALL dyn_ldf( kstp )         ! lateral mixing
+        IF ( ln_stcor )        CALL dyn_stcor ( kstp )      ! WAVE2NEMO add Stokes-Coriolis forcing
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
 #if defined key_agrif
Index: /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7806)
+++ /branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7807)
@@ -110,4 +110,5 @@
    USE timing           ! Timing
 
+   USE dynstcor         ! WAVE2NEMO Stokes-Coriolis forcing
 #if defined key_agrif
    USE agrif_opa_sponge ! Momemtum and tracers sponges
