Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/SBC/sbcmod.F90	(revision 13992)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/SBC/sbcmod.F90	(revision 14009)
@@ -16,4 +16,5 @@
    !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation
    !!            4.0  ! 2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) modified wave forcing and coupling  
    !!----------------------------------------------------------------------
 
@@ -54,4 +55,5 @@
    USE usrdef_sbc     ! user defined: surface boundary condition
    USE closea         ! closed sea
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    !
    USE prtctl         ! Print control                    (prt_ctl routine)
@@ -70,5 +72,6 @@
 
    INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
-
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -99,7 +102,6 @@
          &             nn_ice   , ln_ice_embd,                                       &
          &             ln_traqsr, ln_dm2dc ,                                         &
-         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_apr_dyn,                  &
-         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc , ln_stcor  ,      &
-         &             ln_tauw  , nn_lsm, nn_sdrift
+         &             ln_rnf   , nn_fwb     , ln_ssr   , ln_apr_dyn,                &
+         &             ln_wave  , nn_lsm
       !!----------------------------------------------------------------------
       !
@@ -133,4 +135,5 @@
          WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk
          WRITE(numout,*) '         ABL          formulation                   ln_abl        = ', ln_abl
+         WRITE(numout,*) '         Surface wave (forced or coupled)           ln_wave       = ', ln_wave
          WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : '
          WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
@@ -150,32 +153,5 @@
          WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
          WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
-         WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
-         WRITE(numout,*) '               Stokes drift corr. to vert. velocity ln_sdw        = ', ln_sdw
-         WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift
-         WRITE(numout,*) '               wave modified ocean stress           ln_tauwoc     = ', ln_tauwoc
-         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,NCAR) ln_cdgw       = ', ln_cdgw
-      ENDIF
-      !
-      IF( .NOT.ln_wave ) THEN
-         ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false.
-      ENDIF 
-      IF( ln_sdw ) THEN
-         IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) &
-            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
-      ENDIF
-      ll_st_bv2014  = ( nn_sdrift==jp_breivik_2014 )
-      ll_st_li2017  = ( nn_sdrift==jp_li_2017 )
-      ll_st_bv_li   = ( ll_st_bv2014 .OR. ll_st_li2017 )
-      ll_st_peakfr  = ( nn_sdrift==jp_peakfr )
-      IF( ln_tauwoc .AND. ln_tauw ) &
-         CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
-                                  '(ln_tauwoc=.true. and ln_tauw=.true.)' )
-      IF( ln_tauwoc ) &
-         CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauwoc=.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' )
+      ENDIF
       !
       IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case)
@@ -357,5 +333,22 @@
       IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc, Kbb, Kmm )   ! CICE initialization
       !
-      IF( ln_wave     )   CALL sbc_wave_init                     ! surface wave initialisation
+      IF( ln_wave     ) THEN
+                          CALL sbc_wave_init                     ! surface wave initialisation
+      ELSE
+                          IF(lwp) WRITE(numout,*)
+                          IF(lwp) WRITE(numout,*) '   No surface waves : all wave related logical set to false'
+                          ln_sdw       = .false.
+                          ln_stcor     = .false.
+                          ln_cdgw      = .false.
+                          ln_tauoc     = .false.
+                          ln_wave_test = .false.
+                          ln_charn     = .false.
+                          ln_taw       = .false.
+                          ln_phioc     = .false.
+                          ln_bern_srfc = .false.
+                          ln_breivikFV_2016 = .false.
+                          ln_vortex_force = .false.
+                          ln_stshear  = .false.
+      ENDIF
       !
    END SUBROUTINE sbc_init
@@ -380,4 +373,5 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      INTEGER  ::   jj, ji          ! dummy loop argument
       !
       LOGICAL ::   ll_sas, ll_opa   ! local logical
@@ -412,9 +406,8 @@
       !
       IF( .NOT.ll_sas )   CALL sbc_ssm ( kt, Kbb, Kmm )  ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
-      IF( ln_wave     )   CALL sbc_wave( kt, Kmm )       ! surface waves
-
       !
       !                                            !==  sbc formulation  ==!
       !                                                   
+      !
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
       !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
@@ -423,4 +416,9 @@
       CASE( jp_blk     )
          IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OPA-SAS coupling: SAS receiving fields from OPA
+!!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+         IF( ln_wave )   THEN
+             IF ( lk_oasis )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OPA-wave coupling
+             CALL sbc_wave ( kt, Kmm )
+         ENDIF
                                CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean
                                !
@@ -436,5 +434,33 @@
       IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing
       !
-      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves 
+      IF( ln_wave .AND. ln_tauoc )  THEN            ! Wave stress reduction
+         DO_2D( 0, 0, 0, 0)
+            utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji-1,jj) ) * 0.5_wp
+            vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj-1) ) * 0.5_wp
+         END_2D
+         !
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         !
+         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+         !
+         IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   &
+            &                                'If not requested select ln_tauoc=.false.' )
+         !
+      ELSEIF( ln_wave .AND. ln_taw ) THEN                  ! Wave stress reduction
+         utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:)
+         vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:)
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         !
+         DO_2D( 0, 0, 0, 0)
+             taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
+         END_2D
+         !
+         IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   &
+            &                                'If not requested select ln_taw=.false.' )
+         !
+      ENDIF
+      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. )
       !
       !                                            !==  Misc. Options  ==!
