Index: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
===================================================================
--- branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 8750)
+++ branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 8755)
@@ -66,5 +66,5 @@
    INTEGER                    ::   nsnd         ! total number of fields sent 
    INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
-   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=55   ! Maximum number of coupling fields
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=60   ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields
Index: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8750)
+++ branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8755)
@@ -65,4 +65,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_tauw        !: true if ocean stress components from wave is used
    LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
    !
Index: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8750)
+++ branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8755)
@@ -113,11 +113,13 @@
    INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period
    INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber
-   INTEGER, PARAMETER ::   jpr_wstrf  = 50   ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_tauoc  = 50   ! Stress fraction adsorbed by waves
    INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient
    INTEGER, PARAMETER ::   jpr_isf    = 52
    INTEGER, PARAMETER ::   jpr_icb    = 53
    INTEGER, PARAMETER ::   jpr_wfreq  = 54   ! Wave peak frequency
-
-   INTEGER, PARAMETER ::   jprcv      = 54   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_tauwx  = 55   ! x component of the ocean stress from waves
+   INTEGER, PARAMETER ::   jpr_tauwy  = 56   ! y component of the ocean stress from waves
+
+   INTEGER, PARAMETER ::   jprcv      = 56   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -166,10 +168,11 @@
    TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                        
    !                                   ! Received from the atmosphere
-   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
+   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr, &
+                                                            sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf                              
    ! Send to waves 
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
    ! Received from waves 
-   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag, &
+   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_wdrag, &
                     sn_rcv_wfreq
    !                                   ! Other namelist parameters
@@ -244,8 +247,8 @@
          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
          &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   & 
-         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   &
+         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_tauoc ,   &
          &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &
          &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp  ,   &
-         &                  sn_rcv_icb , sn_rcv_isf   , sn_rcv_wfreq
+         &                  sn_rcv_icb , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw
 
       !!---------------------------------------------------------------------
@@ -298,5 +301,6 @@
          WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')' 
          WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'
-         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')' 
+         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
          WRITE(numout,*)'  sent fields (multiple ice categories)'
@@ -591,8 +595,15 @@
          cpl_wnum = .TRUE.
       ENDIF
-      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave
-      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN
-         srcv(jpr_wstrf)%laction = .TRUE.
-         cpl_wstrf = .TRUE.
+      srcv(jpr_tauoc)%clname = 'O_TauOce'    ! stress fraction adsorbed by the wave
+      IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tauoc)%laction = .TRUE.
+         cpl_tauoc = .TRUE.
+      ENDIF
+      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction
+      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction
+      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tauwx)%laction = .TRUE.
+         srcv(jpr_tauwy)%laction = .TRUE.
+         cpl_tauw = .TRUE.
       ENDIF
       srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
@@ -602,4 +613,8 @@
       ENDIF
       ! 
+      IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) &
+            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
+                                     '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' )
+      !
       !                                                      ! ------------------------------- !
       !                                                      !   OPA-SAS coupling - rcv by opa !   
@@ -1193,5 +1208,13 @@
       !                                                      ! Stress adsorbed by waves  !
       !                                                      ! ========================= ! 
-      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)
+      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1)
+
+      !                                                      ! ========================= !  
+      !                                                      ! Stress component by waves ! 
+      !                                                      ! ========================= !  
+      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN
+         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1)
+         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1)
+      ENDIF
 
       !                                                      ! ========================= ! 
Index: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 8750)
+++ branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 8755)
@@ -29,10 +29,10 @@
    PUBLIC sbc_flx       ! routine called by step.F90
 
-   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read 
-   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
-   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
-   INTEGER , PARAMETER ::   jp_qtot = 3   ! index of total (non solar+solar) heat file
-   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
-   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
+   INTEGER             ::   jpfld         ! maximum number of files to read 
+   INTEGER             ::   jp_utau       ! index of wind stress (i-component) file
+   INTEGER             ::   jp_vtau       ! index of wind stress (j-component) file
+   INTEGER             ::   jp_qtot       ! index of total (non solar+solar) heat file
+   INTEGER             ::   jp_qsr        ! index of solar heat file
+   INTEGER             ::   jp_emp        ! index of evaporation-precipation file
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
 
@@ -82,4 +82,5 @@
       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
       REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
+      LOGICAL  ::   ln_readtau            ! Is it necessary to read utau, vtau from file?
       !!
       CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
@@ -90,4 +91,20 @@
       !
       IF( kt == nit000 ) THEN                ! First call kt=nit000  
+         ln_readtau = .NOT. (ln_wave .AND. ln_tauw )
+
+         ! prepare the index of the fields that have to be read
+         jpfld = 0
+         IF( ln_readtau ) THEN
+            jp_utau = jpfld+1
+            jp_vtau = jpfld+2
+            jpfld = jpfld+2
+         ELSE
+            jp_utau = 0   ;  jp_vtau = 0
+         ENDIF
+         jp_qtot = jpfld+1
+         jp_qsr = jpfld+2
+         jp_emp = jpfld+3
+         jpfld = jpfld+3
+
          ! set file information
          REWIND( numnam_ref )              ! Namelist namsbc_flx in reference namelist : Files for fluxes
@@ -105,5 +122,7 @@
          !
          !                                         ! store namelist information in an array
-         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
+         IF( ln_readtau ) THEN
+            slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
+         ENDIF
          slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
          slf_i(jp_emp ) = sn_emp
@@ -133,6 +152,8 @@
          DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
             DO ji = 1, jpi
-               utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
-               vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
+               IF( ln_readtau ) THEN
+                  utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
+                  vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
+               ENDIF
                qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
                emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
@@ -143,16 +164,18 @@
          !
          !                                                        ! module of wind stress and wind speed at T-point
-         zcoef = 1. / ( zrhoa * zcdrag )
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vect. opt.
-               ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-               zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-               zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-               taum(ji,jj) = zmod
-               wndm(ji,jj) = SQRT( zmod * zcoef )
+         IF( ln_readtau ) THEN
+            zcoef = 1. / ( zrhoa * zcdrag )
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vect. opt.
+                  ztx = utau(ji-1,jj  ) + utau(ji,jj) 
+                  zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
+                  zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+                  taum(ji,jj) = zmod
+                  wndm(ji,jj) = SQRT( zmod * zcoef )
+               END DO
             END DO
-         END DO
-         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
-         CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
+            taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
+            CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
+         ENDIF
 
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
Index: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8750)
+++ branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8755)
@@ -96,5 +96,5 @@
          &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     &
          &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauoc  , ln_stcor   ,     &
-         &             nn_lsm, nn_sdrift
+         &             ln_tauw  , nn_lsm, nn_sdrift
       !!----------------------------------------------------------------------
       !
@@ -159,4 +159,5 @@
          WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift
          WRITE(numout,*) '               wave modified ocean stress           ln_tauoc      = ', ln_tauoc
+         WRITE(numout,*) '               wave modified ocean stress component ln_tauw       = ', ln_tauw
          WRITE(numout,*) '               Stokes coriolis term                 ln_stcor      = ', ln_stcor
          WRITE(numout,*) '               neutral drag coefficient (CORE, MFS) ln_cdgw       = ', ln_cdgw
@@ -167,4 +168,12 @@
             CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
       ENDIF
+      IF( ln_tauoc .AND. ln_tauw ) &
+         CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
+                                  '(ln_tauoc=.true. and ln_tauw=.true.)' )
+      IF( ln_tauoc ) &
+         CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauoc=.true.)' )
+      IF( ln_tauw ) &
+         CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', &
+                              'This will override any other specification of the ocean stress' )
       !
       IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case)
@@ -416,18 +425,8 @@
          IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
       END SELECT
-      IF ( ln_wave .AND. ln_tauoc) THEN                                 ! Wave stress subctracted
-            utau(:,:) = utau(:,:)*tauoc_wave(:,:)
-            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
-            taum(:,:) = taum(:,:)*tauoc_wave(:,:)
-      !
-            SELECT CASE( nsbc )
-            CASE(  0,1,2,3,5,-1 )  ;
-                IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. &
-                        & If not requested select ln_tauoc=.false'
-            END SELECT
-      !
-      END IF
+      !
       IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
-
+      !
+      IF ( ln_wave .AND. (ln_tauoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves 
       !
       !                                            !==  Misc. Options  ==!
Index: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 8750)
+++ branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 8755)
@@ -33,4 +33,5 @@
 
    PUBLIC   sbc_stokes      ! routine called in sbccpl
+   PUBLIC   sbc_wstress     ! routine called in sbcmod 
    PUBLIC   sbc_wave        ! routine called in sbcmod
    PUBLIC   sbc_wave_init   ! routine called in sbcmod
@@ -44,5 +45,6 @@
    LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wstrf  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
 
@@ -58,8 +60,11 @@
    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_tauw    ! structure of input fields (file informations, fields read) ocean stress components from wave model
+
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !: 
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !: 
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:  
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x, tauw_y      !:  
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence
@@ -225,4 +230,41 @@
 
 
+   SUBROUTINE sbc_wstress( )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_wstress  ***
+      !!
+      !! ** Purpose :   Updates the ocean momentum modified by waves
+      !!
+      !! ** Method  : - Calculate u,v components of stress depending on stress
+      !!                model 
+      !!              - Calculate the stress module
+      !!              - The wind module is not modified by waves 
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER  ::   jj, ji   ! dummy loop argument
+      !
+      IF( ln_tauoc ) THEN
+         utau(:,:) = utau(:,:)*tauoc_wave(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+      ENDIF
+      !
+      IF( ln_tauw ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! Stress components at u- & v-points
+               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
+               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
+               !
+               ! Stress module at t points
+               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
+            END DO
+         END DO
+
+      ENDIF
+      !
+   END SUBROUTINE sbc_wstress
+
+
    SUBROUTINE sbc_wave( kt )
       !!---------------------------------------------------------------------
@@ -247,7 +289,13 @@
       ENDIF
 
-      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==!
+      IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN    !==  Wave induced stress  ==!
          CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing
          tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1)
+      ENDIF
+
+      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid)
+         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1)
+         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1)
       ENDIF
 
@@ -300,9 +348,11 @@
       !!
       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), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read
       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
-                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc      ! informations about the fields to be read
-      !
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc
+                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, &
+                             &   sn_tauoc, sn_tauwx, sn_tauwy      ! informations about the fields to be read
+      !
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, &
+                             sn_wnum, sn_tauoc, sn_tauwx, sn_tauwy
       !!---------------------------------------------------------------------
       !
@@ -329,5 +379,5 @@
 
       IF( ln_tauoc ) THEN
-         IF( .NOT. cpl_wstrf ) THEN
+         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' )
@@ -338,4 +388,22 @@
          ENDIF
          ALLOCATE( tauoc_wave(jpi,jpj) )
+      ENDIF
+
+      IF( ln_tauw ) THEN
+         IF( .NOT. cpl_tauw ) THEN
+            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' )
+            !
+            ALLOCATE( slf_j(2) )
+            slf_j(1) = sn_tauwx
+            slf_j(2) = sn_tauwy
+                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   )
+                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   )
+            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) )
+            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( tauw_x(jpi,jpj) )
+         ALLOCATE( tauw_y(jpi,jpj) )
       ENDIF
 
