Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbccpl.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbccpl.F90	(revision 12991)
@@ -8,4 +8,5 @@
    !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
    !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
+   !!            4.2  ! 2020-05  (G. Madec, E. Clementi)  wave coupling updates
    !!----------------------------------------------------------------------
 
@@ -102,20 +103,26 @@
    INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
    INTEGER, PARAMETER ::   jpr_mslp   = 43   ! mean sea level pressure 
-   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig 
-   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux 
-   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1 
-   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2 
+   !**  surface wave coupling  **
+   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig
+   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux
+   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1
+   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2
    INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period
    INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber
-   INTEGER, PARAMETER ::   jpr_tauwoc = 50   ! Stress fraction adsorbed by waves
+   INTEGER, PARAMETER ::   jpr_wstrf  = 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 ::   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 ::   jpr_ts_ice = 57   ! Sea ice surface temp
-
-   INTEGER, PARAMETER ::   jprcv      = 57   ! total number of fields received  
+   INTEGER, PARAMETER ::   jpr_charn  = 52   ! Chranock coefficient
+   INTEGER, PARAMETER ::   jpr_twox   = 53   ! wave to ocean momentum flux
+   INTEGER, PARAMETER ::   jpr_twoy   = 54   ! wave to ocean momentum flux
+   INTEGER, PARAMETER ::   jpr_tawx   = 55   ! net wave-supported stress
+   INTEGER, PARAMETER ::   jpr_tawy   = 56   ! net wave-supported stress
+   INTEGER, PARAMETER ::   jpr_bhd    = 57   ! Bernoulli head. waves' induced surface pressure
+   INTEGER, PARAMETER ::   jpr_tusd   = 58   ! zonal stokes transport
+   INTEGER, PARAMETER ::   jpr_tvsd   = 59   ! meridional stokes tranmport
+   INTEGER, PARAMETER ::   jpr_isf    = 60
+   INTEGER, PARAMETER ::   jpr_icb    = 61
+   INTEGER, PARAMETER ::   jpr_ts_ice = 62   ! Sea ice surface temp
+
+   INTEGER, PARAMETER ::   jprcv      = 62   ! total number of fields received  
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
@@ -172,12 +179,12 @@
       &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
    !                                   ! Received from the atmosphere
-   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr,  &
+   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, sn_rcv_ts_ice
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
-   ! Send to waves 
+   !                                   ! 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_tauwoc, &
-                    sn_rcv_wdrag, sn_rcv_wfreq
+   !                                   ! 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, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
    !                                   ! Other namelist parameters
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -248,13 +255,14 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   & 
-         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  & 
-         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   & 
-         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   & 
-         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  &
+      NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   &
+         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  &
+         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   &
+         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   &
+         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_wstrf,   &
+         &                  sn_rcv_charn , sn_rcv_taw   , sn_rcv_bhd  , sn_rcv_tusd  , sn_rcv_tvsd,    &
          &                  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_tauw, nn_cats_cpl  ,   &
-         &                  sn_rcv_ts_ice
+         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_ts_ice, nn_cats_cpl
+
 
       !!---------------------------------------------------------------------
@@ -294,4 +302,6 @@
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
+         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
+         WRITE(numout,*)'      surface waves:'
          WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')' 
          WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' 
@@ -300,9 +310,7 @@
          WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')' 
          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_tauwoc%cldes), ' (', TRIM(sn_rcv_tauwoc%clcat ), ')' 
-         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
-         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 
+         WRITE(numout,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
@@ -605,9 +613,4 @@
          cpl_wper = .TRUE.
       ENDIF
-      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency 
-      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN
-         srcv(jpr_wfreq)%laction = .TRUE.
-         cpl_wfreq = .TRUE.
-      ENDIF
       srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number
       IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN
@@ -615,15 +618,8 @@
          cpl_wnum = .TRUE.
       ENDIF
-      srcv(jpr_tauwoc)%clname = 'O_TauOce'   ! stress fraction adsorbed by the wave
-      IF( TRIM(sn_rcv_tauwoc%cldes ) == 'coupled' )  THEN
-         srcv(jpr_tauwoc)%laction = .TRUE.
-         cpl_tauwoc = .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.
+      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.
       ENDIF
       srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
@@ -632,7 +628,36 @@
          cpl_wdrag = .TRUE.
       ENDIF
-      IF( srcv(jpr_tauwoc)%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_tauwoc=coupled and sn_rcv_tauw=coupled)' )
+      srcv(jpr_charn)%clname = 'O_Charn'     ! neutral surface drag coefficient
+      IF( TRIM(sn_rcv_charn%cldes ) == 'coupled' )  THEN
+         srcv(jpr_charn)%laction = .TRUE.
+         cpl_charn = .TRUE.
+      ENDIF
+      srcv(jpr_bhd)%clname = 'O_Bhd'     ! neutral surface drag coefficient
+      IF( TRIM(sn_rcv_bhd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_bhd)%laction = .TRUE.
+         cpl_bhd = .TRUE.
+      ENDIF
+      srcv(jpr_tusd)%clname = 'O_Tusd'     ! neutral surface drag coefficient
+      IF( TRIM(sn_rcv_tusd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tusd)%laction = .TRUE.
+         cpl_tusd = .TRUE.
+      ENDIF
+      srcv(jpr_tvsd)%clname = 'O_Tvsd'     ! neutral surface drag coefficient
+      IF( TRIM(sn_rcv_tvsd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tvsd)%laction = .TRUE.
+         cpl_tvsd = .TRUE.
+      ENDIF
+
+      srcv(jpr_twox)%clname = 'O_Twox'     ! wave to ocean momentum flux in the u direction
+      srcv(jpr_twoy)%clname = 'O_Twoy'     ! wave to ocean momentum flux in the v direction
+      srcv(jpr_tawx)%clname = 'O_Tawx'     ! Net wave-supported stress in the u direction
+      srcv(jpr_tawy)%clname = 'O_Tawy'     ! Net wave-supported stress in the v direction
+      IF( TRIM(sn_rcv_taw%cldes ) == 'coupled' )  THEN
+         srcv(jpr_twox)%laction = .TRUE.
+         srcv(jpr_twoy)%laction = .TRUE.
+         srcv(jpr_tawx)%laction = .TRUE.
+         srcv(jpr_tawy)%laction = .TRUE.
+         cpl_taw = .TRUE.
+      ENDIF
       !
       !                                                      ! ------------------------------- !
@@ -1043,4 +1068,5 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
+      !
       !
    END SUBROUTINE sbc_cpl_init
@@ -1118,4 +1144,8 @@
          IF( ncpl_qsr_freq /= 0) ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top
          
+         IF ( ln_wave .AND. nn_components == 0 ) THEN
+            ncpl_qsr_freq = 1;
+            WRITE(numout,*) 'ncpl_qsr_freq is set to 1 when coupling NEMO with wave (without SAS) '
+         ENDIF
       ENDIF
       !
@@ -1280,9 +1310,4 @@
          IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
       ! 
-      !                                                      ! ========================= !  
-      !                                                      !    Wave peak frequency    ! 
-      !                                                      ! ========================= !  
-         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)
-      !
       !                                                      ! ========================= ! 
       !                                                      !    Vertical mixing Qiao   !
@@ -1291,6 +1316,6 @@
 
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
-         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
-                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) THEN
+         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. &
+             srcv(jpr_wper)%laction .OR. srcv(jpr_hsig)%laction )   THEN
             CALL sbc_stokes( Kmm )
          ENDIF
@@ -1299,19 +1324,42 @@
       !                                                      ! Stress adsorbed by waves  !
       !                                                      ! ========================= ! 
-      IF( srcv(jpr_tauwoc)%laction .AND. ln_tauwoc ) tauoc_wave(:,:) = frcv(jpr_tauwoc)%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
-
+      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc )  tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)
+      !
       !                                                      ! ========================= ! 
       !                                                      !   Wave drag coefficient   !
       !                                                      ! ========================= ! 
       IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw )   cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
-
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Wave drag coefficient   !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_charn)%laction .AND. ln_charn )  charn(:,:) = frcv(jpr_charn)%z3(:,:,1)
+      !
+      !
+      IF( srcv(jpr_tawx)%laction .AND. ln_taw )     tawx(:,:) = frcv(jpr_tawx)%z3(:,:,1)
+      IF( srcv(jpr_tawy)%laction .AND. ln_taw )     tawy(:,:) = frcv(jpr_tawy)%z3(:,:,1)
+      IF( srcv(jpr_twox)%laction .AND. ln_taw )     twox(:,:) = frcv(jpr_twox)%z3(:,:,1)
+      IF( srcv(jpr_twoy)%laction .AND. ln_taw )     twoy(:,:) = frcv(jpr_twoy)%z3(:,:,1)
+      !                                                      
+      !                                                      ! ========================= !
+      !                                                      !    wave TKE flux at sfc   !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_phioc)%laction .AND. ln_phioc )     phioc(:,:) = frcv(jpr_phioc)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !      Bernoulli head       !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_bhd)%laction .AND. ln_bern_srfc )   bhd_wave(:,:) = frcv(jpr_bhd)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Stokes transport u dir  !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tusd)%laction .AND. ln_breivikFV_2016 )    tusd(:,:) = frcv(jpr_tusd)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Stokes transport v dir  !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_tvsd)%laction .AND. ln_breivikFV_2016 )     tvsd(:,:) = frcv(jpr_tvsd)%z3(:,:,1)
+      !
       !  Fields received by SAS when OASIS coupling
       !  (arrays no more filled at sbcssm stage)
