Changeset 14021 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE

Timestamp:

2020-12-02T20:53:00+01:00 (3 years ago)

Author:

laurent

Message:

Caught up with trunk rev 14020...

Location:

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE

Files:

• C1D/step_c1d.F90 (modified) (1 diff)
• DYN/dynspg.F90 (modified) (11 diffs)
• DYN/dynvor.F90 (modified) (23 diffs)
• DYN/dynzad.F90 (modified) (6 diffs)
• LBC/lbc_lnk_nc_generic.h90 (added)
• LBC/mpp_nc_generic.h90 (added)
• SBC/cpl_oasis3.F90 (modified) (22 diffs)
• SBC/sbc_oce.F90 (modified) (5 diffs)
• SBC/sbcblk.F90 (modified) (2 diffs)
• SBC/sbcblk_algo_coare3p0.F90 (modified) (1 diff)
• SBC/sbcblk_algo_coare3p6.F90 (modified) (1 diff)
• SBC/sbcblk_algo_ecmwf.F90 (modified) (1 diff)
• SBC/sbcblk_algo_ice_an05.F90 (modified) (1 diff)
• SBC/sbcblk_algo_ice_cdn.F90 (modified) (1 diff)
• SBC/sbcblk_algo_ncar.F90 (modified) (1 diff)
• SBC/sbcblk_skin_coare.F90 (modified) (1 diff)
• SBC/sbcblk_skin_ecmwf.F90 (modified) (1 diff)
• SBC/sbccpl.F90 (modified) (17 diffs)
• SBC/sbcmod.F90 (modified) (11 diffs)
• SBC/sbcwave.F90 (modified) (12 diffs)
• TRA/eosbn2.F90 (modified) (2 diffs)
• TRA/traadv_fct_lf.F90 (added)
• TRA/traadv_mus_lf.F90 (added)
• TRA/trazdf.F90 (modified) (3 diffs)
• ZDF/zdf_oce.F90 (modified) (3 diffs)
• ZDF/zdfmfc.F90 (added)
• ZDF/zdfphy.F90 (modified) (16 diffs)
• ZDF/zdfsh2.F90 (modified) (6 diffs)
• ZDF/zdftke.F90 (modified) (14 diffs)
• step.F90 (modified) (1 diff)
• step_oce.F90 (modified) (1 diff)

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/C1D/step_c1d.F90

@@ -104 +104 @@
       IF( ln_tradmp )   CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs  )  ! internal damping trends- tracers
       IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection
+      IF( ln_zdfmfc  )  CALL tra_mfc( kstp, Nbb     , ts, Nrhs  )  ! Mass Flux Convection
       IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes
                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg.F90

@@ -6 +6 @@
    !! History :  1.0  ! 2005-12  (C. Talandier, G. Madec, V. Garnier)  Original code
    !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dyn_spg     : update the dynamics trend with surface pressure gradient 
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) add Bernoulli Head for
+   !!                            wave coupling
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dyn_spg     : update the dynamics trend with surface pressure gradient
    !!   dyn_spg_init: initialization, namelist read, and parameters control
    !!----------------------------------------------------------------------
@@ -19 +21 @@
    USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
    USE sbcapr         ! surface boundary condition: atmospheric pressure
+   USE sbcwave,  ONLY : bhd_wave
    USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
@@ -36 +39 @@
    PUBLIC   dyn_spg_init   ! routine called by opa module
 
-   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_... 
+   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...
 
    !                       ! Parameter to control the surface pressure gradient scheme
@@ -49 +52 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -58 +61 @@
       !!                  ***  ROUTINE dyn_spg  ***
       !!
-      !! ** Purpose :   compute surface pressure gradient including the 
+      !! ** Purpose :   compute surface pressure gradient including the
       !!              atmospheric pressure forcing (ln_apr_dyn=T).
       !!
@@ -65 +68 @@
       !!              - split-explicit : a time splitting technique is used
       !!
-      !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied 
+      !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied
       !!             as the gradient of the inverse barometer ssh:
       !!                apgu = - 1/rho0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
@@ -86 +89 @@
       !
       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
-         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
          ztrdu(:,:,:) = puu(:,:,:,Krhs)
          ztrdv(:,:,:) = pvv(:,:,:,Krhs)
@@ -140 +143 @@
                spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
             END_2D
-            DEALLOCATE( zpice )         
+            DEALLOCATE( zpice )
+         ENDIF
+         !
+         IF( ln_wave .and. ln_bern_srfc ) THEN          !== Add J terms: depth-independent Bernoulli head
+            DO_2D( 0, 0, 0, 0 )
+               spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
+               spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
+            END_2D
          ENDIF
          !
@@ -149 +159 @@
          !
 !!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
-         !    
+         !
       ENDIF
       !
@@ -156 +166 @@
       CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting
       END SELECT
-      !                    
+      !
       IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics
          ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
          ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt, Kmm )
-         DEALLOCATE( ztrdu , ztrdv ) 
+         DEALLOCATE( ztrdu , ztrdv )
       ENDIF
       !                                      ! print mean trends (used for debugging)
@@ -175 +185 @@
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_spg_init  ***
-      !!                
-      !! ** Purpose :   Control the consistency between namelist options for 
+      !!
+      !! ** Purpose :   Control the consistency between namelist options for
       !!              surface pressure gradient schemes
       !!----------------------------------------------------------------------

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynvor.F90

@@ -15 +15 @@
    !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
-   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity 
+   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity
    !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
    !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
@@ -21 +21 @@
    !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
    !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T)
    !!----------------------------------------------------------------------
 
@@ -37 +38 @@
    USE trddyn         ! trend manager: dynamics
    USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
-   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
+   USE sbc_oce,  ONLY : ln_stcor, ln_vortex_force   ! use Stoke-Coriolis force
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -70 +71 @@
    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
 
-   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
+   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
    !                               ! associated indices:
    INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
@@ -79 +80 @@
 
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2v)/(2*e1e2f)  used in F-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       - 
-   
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       -
+
    REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
    REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
    REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
-   
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -105 +106 @@
       !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) and send them to trd_dyn 
+      !!               and planetary vorticity trends) and send them to trd_dyn
       !!               for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
@@ -121 +122 @@
          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
          !
-         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend
          ztrdv(:,:,:) = pvv(:,:,:,Krhs)
          SELECT CASE( nvor_scheme )
          CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
-            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
-            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
-            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
-            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
-            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          END SELECT
          ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
@@ -161 +157 @@
          CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
                              CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
                              CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_ENE )                        !* energy conserving scheme
                              CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_ENS )                        !* enstrophy conserving scheme
                              CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
+
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          CASE( np_MIX )                        !* mixed ene-ens scheme
                              CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
                              CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
-            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor )        CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )        ! add the Stokes-Coriolis trend
+            IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add vortex force
          CASE( np_EEN )                        !* energy and enstrophy conserving scheme
                              CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
+                             CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
+                             CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
+            ENDIF
          END SELECT
          !
@@ -195 +213 @@
       !!                  ***  ROUTINE vor_enT  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** Purpose :   Compute the now total vorticity trend and add it to
       !!      the general trend of the momentum equation.
       !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
+      !! ** Method  :   Trend evaluated using now fields (centered in time)
       !!       and t-point evaluation of vorticity (planetary and relative).
       !!       conserves the horizontal kinetic energy.
-      !!         The general trend of momentum is increased due to the vorticity 
+      !!         The general trend of momentum is increased due to the vorticity
       !!       term which is given by:
       !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
@@ -235 +253 @@
                   &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
+            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
                DO_2D( 1, 0, 1, 0 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
@@ -250 +268 @@
                   &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
             END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
+            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
                DO_2D( 1, 0, 1, 0 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
@@ -304 +322 @@
                !
             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
-               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
-               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
+               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &
+               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )
          END_2D
          !                                             ! ===============
@@ -317 +335 @@
       !!                  ***  ROUTINE vor_ene  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** Purpose :   Compute the now total vorticity trend and add it to
       !!      the general trend of the momentum equation.
       !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
+      !! ** Method  :   Trend evaluated using now fields (centered in time)
       !!       and the Sadourny (1975) flux form formulation : conserves the
       !!       horizontal kinetic energy.
-      !!         The general trend of momentum is increased due to the vorticity 
+      !!         The general trend of momentum is increased due to the vorticity
       !!       term which is given by:
       !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
@@ -356 +374 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
+            zwz(:,:) = ff_f(:,:)
          CASE ( np_RVO )                           !* relative vorticity
             DO_2D( 1, 0, 1, 0 )
@@ -402 +420 @@
             zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
+            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
          END_2D
          !                                             ! ===============
@@ -452 +470 @@
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwz(:,:) = ff_f(:,:) 
+            zwz(:,:) = ff_f(:,:)
          CASE ( np_RVO )                           !* relative vorticity
             DO_2D( 1, 0, 1, 0 )
@@ -510 +528 @@
       !!                ***  ROUTINE vor_een  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** Purpose :   Compute the now total vorticity trend and add it to
       !!      the general trend of the momentum equation.
       !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves 
+      !! ** Method  :   Trend evaluated using now fields (centered in time)
+      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
       !!      both the horizontal kinetic energy and the potential enstrophy
       !!      when horizontal divergence is zero (see the NEMO documentation)
@@ -654 +672 @@
       !!                ***  ROUTINE vor_eeT  ***
       !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to 
+      !! ** Purpose :   Compute the now total vorticity trend and add it to
       !!      the general trend of the momentum equation.
       !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Arakawa and Lamb (1980) vector form formulation using 
+      !! ** Method  :   Trend evaluated using now fields (centered in time)
+      !!      and the Arakawa and Lamb (1980) vector form formulation using
       !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
-      !!      The change consists in 
+      !!      The change consists in
       !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
       !!
@@ -677 +695 @@
       REAL(wp) ::   zua, zva       ! local scalars
       REAL(wp) ::   zmsk, z1_e3t   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy 
+      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy
       REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
       REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
@@ -837 +855 @@
       !
       IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
-      !                      
+      !
       IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
       ncor = np_COR                       ! planetary vorticity
@@ -846 +864 @@
          ntot = np_COR        !     -         -
       CASE( np_VEC_c2  )
-         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
+         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'
          nrvm = np_RVO        ! relative vorticity
-         ntot = np_CRV        ! relative + planetary vorticity         
+         ntot = np_CRV        ! relative + planetary vorticity
       CASE( np_FLX_c2 , np_FLX_ubs  )
          IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
@@ -873 +891 @@
          !
       END SELECT
-      
+
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
@@ -883 +901 @@
          CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
          CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
-         END SELECT         
+         END SELECT
       ENDIF
       !

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynzad.F90

@@ -7 +7 @@
    !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
    !!----------------------------------------------------------------------
-   
+
    !!----------------------------------------------------------------------
    !!   dyn_zad       : vertical advection momentum trend
@@ -16 +16 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   USE sbcwave, ONLY: wsd   ! Surface Waves (add vertical Stokes-drift)
    !
    USE in_out_manager ! I/O manager
@@ -24 +25 @@
    IMPLICIT NONE
    PRIVATE
-   
+
    PUBLIC   dyn_zad       ! routine called by dynadv.F90
 
@@ -40 +41 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dynzad  ***
-      !! 
-      !! ** Purpose :   Compute the now vertical momentum advection trend and 
+      !!
+      !! ** Purpose :   Compute the now vertical momentum advection trend and
       !!      add it to the general trend of momentum equation.
       !!
@@ -72 +73 @@
 
       IF( l_trddyn )   THEN           ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends
-         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
-         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
-         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
+         ztrdu(:,:,:) = puu(:,:,:,Krhs)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
       ENDIF
-      
+
       DO jk = 2, jpkm1                ! Vertical momentum advection at level w and u- and v- vertical
          DO_2D( 0, 1, 0, 1 )              ! vertical fluxes
+          IF( ln_vortex_force ) THEN
+            zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
+          ELSE
             zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
+          ENDIF
          END_2D
          DO_2D( 0, 0, 0, 0 )              ! vertical momentum advection at w-point
@@ -106 +111 @@
          ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt, Kmm )
-         DEALLOCATE( ztrdu, ztrdv ) 
+         DEALLOCATE( ztrdu, ztrdv )
       ENDIF
       !                               ! Control print

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/cpl_oasis3.F90

@@ -14 +14 @@
    !!            3.6  !  2014-11  (S. Masson) OASIS3-MCT
    !!----------------------------------------------------------------------
-   
+
    !!----------------------------------------------------------------------
    !!   'key_oasis3'                    coupled Ocean/Atmosphere via OASIS3-MCT
@@ -63 +63 @@
 #endif
 
-   INTEGER                    ::   nrcv         ! total number of fields received 
-   INTEGER                    ::   nsnd         ! total number of fields sent 
+   INTEGER                    ::   nrcv         ! total number of fields received
+   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=60   ! Maximum number of coupling fields
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=62   ! 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
-   
+
    TYPE, PUBLIC ::   FLD_CPL               !: Type for coupling field information
       LOGICAL               ::   laction   ! To be coupled or not
-      CHARACTER(len = 8)    ::   clname    ! Name of the coupling field   
-      CHARACTER(len = 1)    ::   clgrid    ! Grid type  
+      CHARACTER(len = 8)    ::   clname    ! Name of the coupling field
+      CHARACTER(len = 1)    ::   clgrid    ! Grid type
       REAL(wp)              ::   nsgn      ! Control of the sign change
       INTEGER, DIMENSION(nmaxcat,nmaxcpl) ::   nid   ! Id of the field (no more than 9 categories and 9 extrena models)
@@ -98 +98 @@
       !!    exchange between AGCM, OGCM and COUPLER. (OASIS3 software)
       !!
-      !! ** Method  :   OASIS3 MPI communication 
+      !! ** Method  :   OASIS3 MPI communication
       !!--------------------------------------------------------------------
       CHARACTER(len = *), INTENT(in   ) ::   cd_modname   ! model name as set in namcouple file
@@ -132 +132 @@
       !!    exchange between AGCM, OGCM and COUPLER. (OASIS3 software)
       !!
-      !! ** Method  :   OASIS3 MPI communication 
+      !! ** Method  :   OASIS3 MPI communication
       !!--------------------------------------------------------------------
       INTEGER, INTENT(in) ::   krcv, ksnd     ! Number of received and sent coupling fields
@@ -180 +180 @@
       !
       ! -----------------------------------------------------------------
-      ! ... Define the partition, excluding halos as we don't want them to be "seen" by oasis    
+      ! ... Define the partition, excluding halos as we don't want them to be "seen" by oasis
       ! -----------------------------------------------------------------
-      
+
       paral(1) = 2                                      ! box partitioning
-      paral(2) = Ni0glo * mjg0(nn_hls) + mig0(nn_hls)   ! NEMO lower left corner global offset, without halos 
+      paral(2) = Ni0glo * mjg0(nn_hls) + mig0(nn_hls)   ! NEMO lower left corner global offset, without halos
       paral(3) = Ni_0                                   ! local extent in i, excluding halos
       paral(4) = Nj_0                                   ! local extent in j, excluding halos
       paral(5) = Ni0glo                                 ! global extent in x, excluding halos
-      
+
       IF( sn_cfctl%l_oasout ) THEN
          WRITE(numout,*) ' multiexchg: paral (1:5)', paral
@@ -195 +195 @@
          WRITE(numout,*) ' multiexchg: Njs0, Nje0, njmpp =', Njs0, Nje0, njmpp
       ENDIF
-   
+
       CALL oasis_def_partition ( id_part, paral, nerror, Ni0glo*Nj0glo )   ! global number of points, excluding halos
       !
-      ! ... Announce send variables. 
+      ! ... Announce send variables.
       !
       ssnd(:)%ncplmodel = kcplmodel
@@ -210 +210 @@
                RETURN
             ENDIF
-            
+
             DO jc = 1, ssnd(ji)%nct
                DO jm = 1, kcplmodel
@@ -225 +225 @@
                   ENDIF
 #if defined key_agrif
-                  IF( agrif_fixed() /= 0 ) THEN 
+                  IF( agrif_fixed() /= 0 ) THEN
                      zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
                   ENDIF
@@ -243 +243 @@
       END DO
       !
-      ! ... Announce received variables. 
+      ! ... Announce received variables.
       !
       srcv(:)%ncplmodel = kcplmodel
       !
       DO ji = 1, krcv
-         IF( srcv(ji)%laction ) THEN 
-            
+         IF( srcv(ji)%laction ) THEN
+
             IF( srcv(ji)%nct > nmaxcat ) THEN
                CALL oasis_abort ( ncomp_id, 'cpl_define', 'Number of categories of '//   &
@@ -255 +255 @@
                RETURN
             ENDIF
-            
+
             DO jc = 1, srcv(ji)%nct
                DO jm = 1, kcplmodel
-                  
+
                   IF( srcv(ji)%nct .GT. 1 ) THEN
                      WRITE(cli2,'(i2.2)') jc
@@ -270 +270 @@
                   ENDIF
 #if defined key_agrif
-                  IF( agrif_fixed() /= 0 ) THEN 
+                  IF( agrif_fixed() /= 0 ) THEN
                      zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
                   ENDIF
@@ -288 +288 @@
          ENDIF
       END DO
-      
+
       !------------------------------------------------------------------
       ! End of definition phase
       !------------------------------------------------------------------
-      !     
+      !
 #if defined key_agrif
       IF( agrif_fixed() == Agrif_Nb_Fine_Grids() ) THEN
@@ -303 +303 @@
       !
    END SUBROUTINE cpl_define
-   
-   
+
+
    SUBROUTINE cpl_snd( kid, kstep, pdata, kinfo )
       !!---------------------------------------------------------------------
@@ -324 +324 @@
       DO jc = 1, ssnd(kid)%nct
          DO jm = 1, ssnd(kid)%ncplmodel
-        
+
             IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN   ! exclude halos from data sent to oasis
                CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0:Nie0, Njs0:Nje0,jc), kinfo )
-               
-               IF ( sn_cfctl%l_oasout ) THEN        
+
+               IF ( sn_cfctl%l_oasout ) THEN
                   IF ( kinfo == OASIS_Sent     .OR. kinfo == OASIS_ToRest .OR.   &
                      & kinfo == OASIS_SentOut  .OR. kinfo == OASIS_ToRestOut ) THEN
@@ -342 +342 @@
                   ENDIF
                ENDIF
-               
+
             ENDIF
-            
+
          ENDDO
       ENDDO
@@ -379 +379 @@
             IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
 
-               CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )         
-               
+               CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )
+
                llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   &
                   &        kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
-               
+
                IF ( sn_cfctl%l_oasout )   &
                   &  WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
-               
+
                IF( llaction ) THEN   ! data received from oasis do not include halos
-                  
+
                   kinfo = OASIS_Rcv
-                  IF( ll_1st ) THEN 
+                  IF( ll_1st ) THEN
                      pdata(Nis0:Nie0,Njs0:Nje0,jc) =   exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
                      ll_1st = .FALSE.
@@ -397 +397 @@
                         &                                + exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
                   ENDIF
-                  
-                  IF ( sn_cfctl%l_oasout ) THEN        
+
+                  IF ( sn_cfctl%l_oasout ) THEN
                      WRITE(numout,*) '****************'
                      WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname
@@ -409 +409 @@
                      WRITE(numout,*) '****************'
                   ENDIF
-                  
+
                ENDIF
-               
+
             ENDIF
-            
+
          ENDDO
 
          !--- we must call lbc_lnk to fill the halos that where not received.
          IF( .NOT. ll_1st ) THEN
-            CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )   
+            CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
          ENDIF
- 
+
       ENDDO
       !
@@ -426 +426 @@
 
 
-   INTEGER FUNCTION cpl_freq( cdfieldname )  
+   INTEGER FUNCTION cpl_freq( cdfieldname )
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_freq  ***
@@ -491 +491 @@
       DEALLOCATE( exfld )
       IF(nstop == 0) THEN
-         CALL oasis_terminate( nerror )         
+         CALL oasis_terminate( nerror )
       ELSE
          CALL oasis_abort( ncomp_id, "cpl_finalize", "NEMO ABORT STOP" )
-      ENDIF       
+      ENDIF
       !
    END SUBROUTINE cpl_finalize
@@ -544 +544 @@
       WRITE(numout,*) 'oasis_enddef: Error you sould not be there...'
    END SUBROUTINE oasis_enddef
-  
+
    SUBROUTINE oasis_put(k1,k2,p1,k3)
       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::  p1
@@ -574 +574 @@
       WRITE(numout,*) 'oasis_terminate: Error you sould not be there...'
    END SUBROUTINE oasis_terminate
-   
+
 #endif
 

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_oce.F90

@@ -12 +12 @@
    !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk)
    !!            4.0  ! 2019-03  (F. Lemarié, G. Samson) add compatibility with ABL mode
+   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) modified wave parameters in namelist
    !!----------------------------------------------------------------------
 
@@ -36 +37 @@
    LOGICAL , PUBLIC ::   ln_blk         !: bulk formulation
    LOGICAL , PUBLIC ::   ln_abl         !: Atmospheric boundary layer model
+   LOGICAL , PUBLIC ::   ln_wave        !: wave in the system (forced or coupled)
 #if defined key_oasis3
    LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used
@@ -56 +58 @@
    !                                             !:  = 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_tauwoc       !: 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
-   !
-   INTEGER , PUBLIC ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
-   !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
    !
@@ -71 +64 @@
    !                                   !!* namsbc_cpl namelist *
    INTEGER , PUBLIC ::   nn_cats_cpl    !: Number of sea ice categories over which the coupling is carried out
-
+   !
+   !                                   !!* namsbc_wave namelist *
+   LOGICAL , PUBLIC ::   ln_sdw         !: =T 3d stokes drift from wave model
+   LOGICAL , PUBLIC ::   ln_stcor       !: =T if Stokes-Coriolis and tracer advection terms are used
+   LOGICAL , PUBLIC ::   ln_cdgw        !: =T neutral drag coefficient from wave model
+   LOGICAL , PUBLIC ::   ln_tauoc       !: =T if normalized stress from wave is used
+   LOGICAL , PUBLIC ::   ln_wave_test   !: =T wave test case (constant Stokes drift)
+   LOGICAL , PUBLIC ::   ln_charn       !: =T Chranock coefficient from wave model
+   LOGICAL , PUBLIC ::   ln_taw         !: =T wind stress corrected by wave intake
+   LOGICAL , PUBLIC ::   ln_phioc       !: =T TKE surface BC from wave model
+   LOGICAL , PUBLIC ::   ln_bern_srfc   !: Bernoulli head, waves' inuced pressure
+   LOGICAL , PUBLIC ::   ln_breivikFV_2016 !: Breivik 2016 profile
+   LOGICAL , PUBLIC ::   ln_vortex_force !: vortex force activation
+   LOGICAL , PUBLIC ::   ln_stshear     !: Stoked Drift shear contribution in zdftke
+   !
    !!----------------------------------------------------------------------
    !!           switch definition (improve readability)
@@ -81 +88 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 5        !: Pure ocean-atmosphere Coupled formulation
    INTEGER , PUBLIC, PARAMETER ::   jp_none    = 6        !: for OPA when doing coupling via SAS module
-
-   !!----------------------------------------------------------------------
-   !!           Stokes drift parametrization definition
-   !!----------------------------------------------------------------------
-   INTEGER , PUBLIC, PARAMETER ::   jp_breivik_2014 = 0     !: Breivik  2014: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
-   INTEGER , PUBLIC, PARAMETER ::   jp_li_2017      = 1     !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016)
-   !  with depth averaged profile
-   INTEGER , PUBLIC, PARAMETER ::   jp_peakfr       = 2     !: Li et al 2017: using the peak wave number read from wave model instead
-   !  of the inverse depth scale
-   LOGICAL , PUBLIC            ::   ll_st_bv2014  = .FALSE. !  logical indicator, .true. if Breivik 2014 parameterisation is active.
-   LOGICAL , PUBLIC            ::   ll_st_li2017  = .FALSE. !  logical indicator, .true. if Li 2017 parameterisation is active.
-   LOGICAL , PUBLIC            ::   ll_st_bv_li   = .FALSE. !  logical indicator, .true. if either Breivik or Li parameterisation is active.
-   LOGICAL , PUBLIC            ::   ll_st_peakfr  = .FALSE. !  logical indicator, .true. if using Li 2017 with peak wave number
-
+   !
    !!----------------------------------------------------------------------
    !!           component definition

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90

@@ -19 +19 @@
    !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)
    !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
+   !!            4.2  !  2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
    !!----------------------------------------------------------------------
 
@@ -383 +384 @@
             !
             IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   &
-               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
-               &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )
+         &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
+         &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )
          ENDIF
       END DO
-      !
-      IF( ln_wave ) THEN
-         !Activated wave module but neither drag nor stokes drift activated
-         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
-            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
-            !drag coefficient read from wave model definable only with mfs bulk formulae and core
-         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN
-            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR bulk formulae')
-         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN
-            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
-         ENDIF
-      ELSE
-         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
-            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
-            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
-            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
-            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &
-            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
-      ENDIF
       !
       IF( ln_abl ) THEN       ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_coare3p0.F90

@@ -15 +15 @@
    !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
-   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !! History :  4.0  ! 2016-02  (L.Brodeau)   Original code
+   !!            4.2  ! 2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
    !!----------------------------------------------------------------------
 

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_coare3p6.F90

@@ -16 +16 @@
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!            4.2  !  2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
    !!----------------------------------------------------------------------
 

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90

@@ -17 +17 @@
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!            4.2  !  2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
    !!----------------------------------------------------------------------
 

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ice_an05.F90

@@ -19 +19 @@
    !!----------------------------------------------------------------------
    USE par_kind, ONLY: wp
-   USE par_oce,  ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls
+   USE par_oce,  ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls, ntsi, ntsj, ntei, ntej
    USE lib_mpp,  ONLY: ctl_stop         ! distribued memory computing library
    USE phycst          ! physical constants

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ice_cdn.F90

@@ -13 +13 @@
    !!====================================================================================
    USE par_kind, ONLY: wp
-   USE par_oce,  ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls
+   USE par_oce,  ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls, ntsi, ntsj, ntei, ntej
    USE phycst          ! physical constants
    USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ncar.F90

@@ -16 +16 @@
    !!=====================================================================
    !! History :  3.6  !  2016-02  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90
+   !!            4.2  !  2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
    !!----------------------------------------------------------------------
 

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_coare.F90

@@ -13 +13 @@
    !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
-   !! History :  4.x  !  2019-11  (L.Brodeau)   Original code
+   !! History :  4.0  !  2019-11  (L.Brodeau)   Original code
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_ecmwf.F90

@@ -28 +28 @@
    !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
-   !! History :  4.x  !  2019-11  (L.Brodeau)   Original code
+   !! History :  4.0  ! 2019-11  (L.Brodeau)   Original code
+   !!            4.2  ! 2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbccpl.F90

@@ -8 +8 @@
    !!            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-12  (G. Madec, E. Clementi)  wave coupling updates
    !!----------------------------------------------------------------------
 
@@ -108 +109 @@
    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
+   !**  surface wave coupling  **
    INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig
    INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux
@@ -114 +116 @@
    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
@@ -186 +193 @@
       &             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
@@ -276 +283 @@
          &                  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,  &
-         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  &
-         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
-         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
-         &                  sn_rcv_ts_ice
+         &                  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   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice
+
       !!---------------------------------------------------------------------
       !
@@ -321 +328 @@
          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 ), ')'
@@ -327 +336 @@
          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  ), ')'
@@ -353 +360 @@
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd
       ENDIF
-
+      IF( lwp .AND. ln_wave) THEN                        ! control print
+      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 ), ')'
+         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'
+         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'
+         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,*)'      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,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift grid u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift grid v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
+         WRITE(numout,*)'      Bernouilli pressure head        = ', TRIM(sn_rcv_bhd%cldes   ), ' (', TRIM(sn_rcv_bhd%clcat  ), ')'
+         WRITE(numout,*)'Wave to ocean momentum flux and Net wave-supported stress = ', TRIM(sn_rcv_taw%cldes ), ' (', TRIM(sn_rcv_taw%clcat ), ')'
+         WRITE(numout,*)'      Surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'
+         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref
+         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor
+         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd
+      ENDIF
       !                                   ! allocate sbccpl arrays
       IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
@@ -631 +657 @@
          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
@@ -641 +662 @@
          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
@@ -658 +672 @@
          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'     ! Chranock coefficient
+      IF( TRIM(sn_rcv_charn%cldes ) == 'coupled' )  THEN
+         srcv(jpr_charn)%laction = .TRUE.
+         cpl_charn = .TRUE.
+      ENDIF
+      srcv(jpr_bhd)%clname = 'O_Bhd'     ! Bernoulli head. waves' induced surface pressure
+      IF( TRIM(sn_rcv_bhd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_bhd)%laction = .TRUE.
+         cpl_bhd = .TRUE.
+      ENDIF
+      srcv(jpr_tusd)%clname = 'O_Tusd'     ! zonal stokes transport
+      IF( TRIM(sn_rcv_tusd%cldes ) == 'coupled' )  THEN
+         srcv(jpr_tusd)%laction = .TRUE.
+         cpl_tusd = .TRUE.
+      ENDIF
+      srcv(jpr_tvsd)%clname = 'O_Tvsd'     ! meridional stokes tranmport
+      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
       !
       !                                                      ! ------------------------------- !
@@ -1060 +1103 @@
       !   initialisation of the coupler  !
       ! ================================ !
-
       CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
 
@@ -1073 +1115 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
+      !
       !
    END SUBROUTINE sbc_cpl_init
@@ -1148 +1191 @@
          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
       !
@@ -1323 +1370 @@
       !
       !                                                      ! ========================= !
-      !                                                      !    Wave peak frequency    !
-      !                                                      ! ========================= !
-         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)
-      !
-      !                                                      ! ========================= !
       !                                                      !    Vertical mixing Qiao   !
       !                                                      ! ========================= !
@@ -1333 +1375 @@
 
          ! 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
@@ -1341 +1383 @@
       !                                                      ! 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)
-
+      !
+      !                                                      ! ========================= !
+      !                                                      !   Chranock coefficient    !
+      !                                                      ! ========================= !
+      IF( srcv(jpr_charn)%laction .AND. ln_charn )  charn(:,:) = frcv(jpr_charn)%z3(:,:,1)
+      !
+      !                                                      ! ========================= !
+      !                                                      ! net wave-supported stress !
+      !                                                      ! ========================= !
+      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)
+      !
+      !                                                      ! ========================= !
+      !                                                      !wave to ocean momentum flux!
+      !                                                      ! ========================= !
+      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)

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcmod.F90

@@ -16 +16 @@
    !!            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
    !!----------------------------------------------------------------------
 
@@ -55 +56 @@
    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)
@@ -71 +73 @@
 
    INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
-
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -100 +103 @@
          &             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
       !!----------------------------------------------------------------------
       !
@@ -134 +136 @@
          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
@@ -151 +154 @@
          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)
@@ -358 +334 @@
       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
@@ -381 +374 @@
       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
@@ -413 +407 @@
       !
       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
@@ -424 +417 @@
       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
                                !
@@ -437 +435 @@
       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  ==!

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcwave.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  sbcwave  ***
-   !! Wave module 
+   !! Wave module
    !!======================================================================
-   !! History :  3.3  !  2011-09  (M. Adani)  Original code: Drag Coefficient 
-   !!         :  3.4  !  2012-10  (M. Adani)  Stokes Drift 
+   !! History :  3.3  !  2011-09  (M. Adani)  Original code: Drag Coefficient
+   !!         :  3.4  !  2012-10  (M. Adani)  Stokes Drift
    !!            3.6  !  2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation
    !!             -   !  2016-12  (G. Madec, E. Clementi) update Stoke drift computation
    !!                                                    + add sbc_wave_ini routine
+   !!            4.2  !  2020-12  (G. Madec, E. Clementi) updates, new Stoke drift computation
+   !!                                                    according to Couvelard et al.,2019
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   sbc_stokes    : calculate 3D Stokes-drift velocities
-   !!   sbc_wave      : wave data from wave model in netcdf files 
-   !!   sbc_wave_init : initialisation fo surface waves 
+   !!   sbc_wave      : wave data from wave model: forced (netcdf files) or coupled mode
+   !!   sbc_wave_init : initialisation fo surface waves
    !!----------------------------------------------------------------------
-   USE phycst         ! physical constants 
+   USE phycst         ! physical constants
    USE oce            ! ocean variables
-   USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE zdf_oce,  ONLY : ln_zdfswm
+   USE dom_oce        ! ocean domain variables
+   USE sbc_oce        ! Surface boundary condition: ocean fields
    USE bdy_oce        ! open boundary condition variables
    USE domvvl         ! domain: variable volume layers
@@ -26 +28 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! distribued memory computing library
-   USE fldread        ! read input fields
+   USE fldread        ! read input fields
 
    IMPLICIT NONE
@@ -32 +34 @@
 
    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
-   
+
    ! Variables checking if the wave parameters are coupled (if not, they are read from file)
-   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_tauwoc = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_hsig          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_phioc         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_sdrftx        = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_sdrfty        = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wper          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wnum          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wstrf         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wdrag         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_charn         = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_taw           = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_bhd           = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tusd          = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tvsd          = .FALSE.
 
    INTEGER ::   jpfld    ! number of files to read for stokes drift
@@ -53 +57 @@
    INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point
    INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
-   INTEGER ::   jp_wfr   ! index of wave peak frequency         (1/s)    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
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_wn      ! structure of input fields (file informations, fields read) wave number for Qiao
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauwoc  ! 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
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp.
-
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauoc   ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
+
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave        !: Neutral drag coefficient at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw             !: Significant Wave Height at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wmp             !: Wave Mean Period at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wnum            !: Wave Number at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave      !: stress reduction factor  at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d           !: Surface Stokes Drift module at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd          !: barotropic stokes drift divergence
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd    !: surface Stokes drift velocities at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd, vsd, wsd   !: Stokes drift velocities at u-, v- & w-points, resp.u
+!
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   charn           !: charnock coefficient at t-point
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tawx            !: Net wave-supported stress, u
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tawy            !: Net wave-supported stress, v
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   twox            !: wave-ocean momentum flux, u
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   twoy            !: wave-ocean momentum flux, v
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wavex     !: stress reduction factor  at, u component
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wavey     !: stress reduction factor  at, v component
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   phioc           !: tke flux from wave model
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   KZN2            !: Kz*N2
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   bhd_wave        !: Bernoulli head. wave induce pression
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tusd, tvsd      !: Stokes drift transport
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   ZMX             !: Kz*N2
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -88 +102 @@
       !!                2014 (DOI: 10.1175/JPO-D-14-0020.1)
       !!
-      !! ** Method  : - Calculate Stokes transport speed 
-      !!              - Calculate horizontal divergence 
-      !!              - Integrate the horizontal divergenze from the bottom 
-      !! ** action  
+      !! ** Method  : - Calculate the horizontal Stokes drift velocity (Breivik et al. 2014)
+      !!              - Calculate its horizontal divergence
+      !!              - Calculate the vertical Stokes drift velocity
+      !!              - Calculate the barotropic Stokes drift divergence
+      !!
+      !! ** action  : - tsd2d         : module of the surface Stokes drift velocity
+      !!              - usd, vsd, wsd : 3 components of the Stokes drift velocity
+      !!              - div_sd        : barotropic Stokes drift divergence
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: Kmm ! ocean time level index
       INTEGER  ::   jj, ji, jk   ! dummy loop argument
-      INTEGER  ::   ik           ! local integer 
-      REAL(wp) ::  ztransp, zfac, zsp0
-      REAL(wp) ::  zdepth, zsqrt_depth,  zexp_depth, z_two_thirds, zsqrtpi !sqrt of pi
-      REAL(wp) ::  zbot_u, zbot_v, zkb_u, zkb_v, zke3_u, zke3_v, zda_u, zda_v
-      REAL(wp) ::  zstokes_psi_u_bot, zstokes_psi_v_bot
-      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd     ! 2D workspace
-      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zstokes_psi_u_top, zstokes_psi_v_top ! 2D workspace
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh                              ! 3D workspace
-      !!---------------------------------------------------------------------
-      !
-      ALLOCATE( ze3divh(jpi,jpj,jpkm1) )   ! jpkm1 -> avoid lbc_lnk on jpk that is not defined
+      INTEGER  ::   ik           ! local integer
+      REAL(wp) ::  ztransp, zfac, ztemp, zsp0, zsqrt, zbreiv16_w
+      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v, sdtrp
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd ! 2D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh, zInt_w                  ! 3D workspace
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( ze3divh(jpi,jpj,jpkm1) ) ! jpkm1 -> avoid lbc_lnk on jpk that is not defined
+      ALLOCATE( zInt_w(jpi,jpj,jpk) )
       ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) )
+      zk_t    (:,:) = 0._wp
+      zk_u    (:,:) = 0._wp
+      zk_v    (:,:) = 0._wp
+      zu0_sd  (:,:) = 0._wp
+      zv0_sd  (:,:) = 0._wp
+      ze3divh (:,:,:) = 0._wp
+
       !
       ! select parameterization for the calculation of vertical Stokes drift
       ! exp. wave number at t-point
-      IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) )
+      IF( ln_breivikFV_2016 ) THEN
+      ! Assumptions :  ut0sd and vt0sd are surface Stokes drift at T-points
+      !                sdtrp is the norm of Stokes transport
+      !
+         zfac = 0.166666666667_wp
+         DO_2D( 1, 1, 1, 1 ) ! In the deep-water limit we have ke = ||ust0||/( 6 * ||transport|| )
+            zsp0          = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) !<-- norm of Surface Stokes drift
+            tsd2d(ji,jj)  = zsp0
+            IF( cpl_tusd .AND. cpl_tvsd ) THEN  !stokes transport is provided in coupled mode
+               sdtrp      = SQRT( tusd(ji,jj)*tusd(ji,jj) + tvsd(ji,jj)*tvsd(ji,jj) )  !<-- norm of Surface Stokes drift transport
+            ELSE
+               ! Stokes drift transport estimated from Hs and Tmean
+               sdtrp      = 2.0_wp * rpi / 16.0_wp *                             &
+                   &        hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
+            ENDIF
+            zk_t (ji,jj)  = zfac * zsp0 / MAX ( sdtrp, 0.0000001_wp ) !<-- ke = ||ust0||/( 6 * ||transport|| )
+         END_2D
+      !# define zInt_w ze3divh
+         DO_3D( 1, 1, 1, 1, 1, jpk ) ! Compute the primitive of Breivik 2016 function at W-points
+            zfac             = - 2._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm)  !<-- zfac should be negative definite
+            ztemp            = EXP ( zfac )
+            zsqrt            = SQRT( -zfac )
+            zbreiv16_w       = ztemp - SQRT(rpi)*zsqrt*ERFC(zsqrt) !Eq. 16 Breivik 2016
+            zInt_w(ji,jj,jk) = ztemp - 4._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) * zbreiv16_w
+         END_3D
+!
+         DO jk = 1, jpkm1
+            zfac = 0.166666666667_wp
+            DO_2D( 1, 1, 1, 1 ) !++ Compute the FV Breivik 2016 function at T-points
+               zsp0          = zfac / MAX(zk_t (ji,jj),0.0000001_wp)
+               ztemp         = zInt_w(ji,jj,jk) - zInt_w(ji,jj,jk+1)
+               zu0_sd(ji,jj) = ut0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk)
+               zv0_sd(ji,jj) = vt0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk)
+            END_2D
+            DO_2D( 1, 0, 1, 0 ) ! ++ Interpolate at U/V points
+               zfac          =  1.0_wp / e3u(ji  ,jj,jk,Kmm)
+               usd(ji,jj,jk) =  0.5_wp * zfac * ( zu0_sd(ji,jj)+zu0_sd(ji+1,jj) ) * umask(ji,jj,jk)
+               zfac          =  1.0_wp / e3v(ji  ,jj,jk,Kmm)
+               vsd(ji,jj,jk) =  0.5_wp * zfac * ( zv0_sd(ji,jj)+zv0_sd(ji,jj+1) ) * vmask(ji,jj,jk)
+            END_2D
+         ENDDO
+      !# undef zInt_w
+      !
+      ELSE
          zfac = 2.0_wp * rpi / 16.0_wp
          DO_2D( 1, 1, 1, 1 )
@@ -128 +193 @@
             zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
          END_2D
-      ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
-         DO_2D( 1, 1, 1, 1 )
-            zk_t(ji,jj) = ( 2.0_wp * rpi * wfreq(ji,jj) ) * ( 2.0_wp * rpi * wfreq(ji,jj) ) / grav
-         END_2D
-         DO_2D( 1, 0, 1, 0 )
-            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
-            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
-            !
-            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
-            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
-         END_2D
-      ENDIF
-      !
+
       !                       !==  horizontal Stokes Drift 3D velocity  ==!
-      IF( ll_st_bv2014 ) THEN
+
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zdep_u = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji+1,jj,jk,Kmm) )
             zdep_v = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji,jj+1,jk,Kmm) )
-            !                          
+            !
             zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
             zkh_v = zk_v(ji,jj) * zdep_v
@@ -156 +209 @@
             vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
          END_3D
-      ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN
-         ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) )
-         DO_2D( 1, 0, 1, 0 )
-            zstokes_psi_u_top(ji,jj) = 0._wp
-            zstokes_psi_v_top(ji,jj) = 0._wp
-         END_2D
-         zsqrtpi = SQRT(rpi)
-         z_two_thirds = 2.0_wp / 3.0_wp
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! exp. wave number & Stokes drift velocity at u- & v-points
-            zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth
-            zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth
-            zkb_u  = zk_u(ji,jj) * zbot_u                             ! 2 * k * bottom depth
-            zkb_v  = zk_v(ji,jj) * zbot_v                             ! 2 * k * bottom depth
-            !
-            zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u(ji,jj,jk,Kmm))     ! 2k * thickness
-            zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v(ji,jj,jk,Kmm))     ! 2k * thickness
-
-            ! Depth attenuation .... do u component first..
-            zdepth      = zkb_u
-            zsqrt_depth = SQRT(zdepth)
-            zexp_depth  = EXP(-zdepth)
-            zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
-                 &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
-                 &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
-            zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u
-            zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot
-
-            !         ... and then v component
-            zdepth      =zkb_v
-            zsqrt_depth = SQRT(zdepth)
-            zexp_depth  = EXP(-zdepth)
-            zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
-                 &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
-                 &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
-            zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v
-            zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot
-            !
-            usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
-            vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
-         END_3D
-         DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top )
       ENDIF
 
@@ -228 +240 @@
       !                       !==  Horizontal divergence of barotropic Stokes transport  ==!
       div_sd(:,:) = 0._wp
-      DO jk = 1, jpkm1                                 ! 
+      DO jk = 1, jpkm1                                 !
         div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk)
       END DO
@@ -235 +247 @@
       CALL iom_put( "vstokes",  vsd  )
       CALL iom_put( "wstokes",  wsd  )
-      !
-      DEALLOCATE( ze3divh )
+!      !
+      DEALLOCATE( ze3divh, zInt_w )
       DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
       !
    END SUBROUTINE sbc_stokes
-
-
-   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_tauwoc ) THEN
-         utau(:,:) = utau(:,:)*tauoc_wave(:,:)
-         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
-         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
-      ENDIF
-      !
-      IF( ln_tauw ) THEN
-         DO_2D( 1, 0, 1, 0 )
-            ! 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_2D
-         CALL lbc_lnk_multi( 'sbcwave', utau(:,:), 'U', -1.0_wp , vtau(:,:), 'V', -1.0_wp , taum(:,:) , 'T', -1.0_wp )
-      ENDIF
-      !
-   END SUBROUTINE sbc_wstress
-
-
+!
+!
    SUBROUTINE sbc_wave( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE sbc_wave  ***
       !!
-      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
-      !!
-      !! ** 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 Breivik et al.,2014
-      !!                formulation
-      !! ** action  
+      !! ** Purpose :   read wave parameters from wave model in netcdf files
+      !!                or from a coupled wave mdoel
+      !!
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt   ! ocean time step
       INTEGER, INTENT(in   ) ::   Kmm  ! ocean time index
       !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave : update the read waves fields'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
       !
       IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==!
@@ -300 +277 @@
       ENDIF
 
-      IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN  !==  Wave induced stress  ==!
-         CALL fld_read( kt, nn_fsbc, sf_tauwoc )         ! read wave norm stress from external forcing
-         tauoc_wave(:,:) = sf_tauwoc(1)%fnow(:,:,1) * tmask(:,:,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) * tmask(:,:,1)
-         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) * tmask(:,:,1)
-      ENDIF
-
-      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==! 
+      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read stress reduction factor due to wave from external forcing
+         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) * tmask(:,:,1)
+      ELSEIF ( ln_taw .AND. cpl_taw ) THEN
+         IF (kt < 1) THEN ! The first fields gave by OASIS have very high erroneous values ....
+            twox(:,:)=0._wp
+            twoy(:,:)=0._wp
+            tawx(:,:)=0._wp
+            tawy(:,:)=0._wp
+            tauoc_wavex(:,:) = 1._wp
+            tauoc_wavey(:,:) = 1._wp
+         ELSE
+            tauoc_wavex(:,:) = abs(twox(:,:)/tawx(:,:))
+            tauoc_wavey(:,:) = abs(twoy(:,:)/tawy(:,:))
+         ENDIF
+      ENDIF
+
+      IF ( ln_phioc .and. cpl_phioc .and.  kt == nit000 ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave : PHIOC from wave model'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+
+      IF( ln_sdw .AND. .NOT. cpl_sdrftx)  THEN       !==  Computation of the 3d Stokes Drift  ==!
          !
          IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
             CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
+            !                                            ! NB: test case mode, not read as jpfld=0
             IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) * tmask(:,:,1)  ! significant wave height
             IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) * tmask(:,:,1)  ! wave mean period
-            IF( jp_wfr > 0 )   wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1) * tmask(:,:,1)  ! Peak wave frequency
             IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) * tmask(:,:,1)  ! 2D zonal Stokes Drift at T point
             IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) * tmask(:,:,1)  ! 2D meridional Stokes Drift at T point
          ENDIF
          !
-         ! Read also wave number if needed, so that it is available in coupling routines
-         IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN
-            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing
-            wnum(:,:) = sf_wn(1)%fnow(:,:,1) * tmask(:,:,1)
-         ENDIF
-           
-         ! Calculate only if required fields have been read
-         ! In coupled wave model-NEMO case the call is done after coupling
+         IF( jpfld == 4 .OR. ln_wave_test )   &
+            &      CALL sbc_stokes( Kmm )                 ! Calculate only if all required fields are read
+            !                                            ! or in wave test case
+         !  !                                            ! In coupled case the call is done after (in sbc_cpl)
+      ENDIF
          !
-         IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &
-           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes( Kmm )
-         !
-      ENDIF
-      !
    END SUBROUTINE sbc_wave
 
@@ -343 +324 @@
       !!                     ***  ROUTINE sbc_wave_init  ***
       !!
-      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
+      !! ** Purpose :   Initialisation fo surface waves
       !!
       !! ** 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 Breivik et al.,2014
-      !!                formulation
-      !! ** action  
+      !!              - create the structure used to read required wave fields
+      !!                (its size depends on namelist options)
+      !! ** action
       !!---------------------------------------------------------------------
       INTEGER ::   ierror, ios   ! local integer
@@ -357 +335 @@
       !!
       CHARACTER(len=100)     ::  cn_dir                            ! Root directory for location of drag coefficient files
-      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i            ! array of namelist informations on the fields to read
       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
-                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, &
-                             &   sn_tauwoc, 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_tauwoc, sn_tauwx, sn_tauwy
-      !!---------------------------------------------------------------------
+                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc    ! informations about the fields to be read
+      !
+      NAMELIST/namsbc_wave/ cn_dir, sn_cdg, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc,   &
+         &                  ln_cdgw, ln_sdw, ln_tauoc, ln_stcor, ln_charn, ln_taw, ln_phioc,     &
+         &                  ln_wave_test, ln_bern_srfc, ln_breivikFV_2016, ln_vortex_force, ln_stshear
+      !!---------------------------------------------------------------------
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave_init : surface waves in the system'
+         WRITE(numout,*) '~~~~~~~~~~~~~ '
+      ENDIF
       !
       READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_wave in reference namelist' )
-         
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_wave in reference namelist')
+
       READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist' )
       IF(lwm) WRITE ( numond, namsbc_wave )
       !
-      IF( ln_cdgw ) THEN
-         IF( .NOT. cpl_wdrag ) THEN
-            ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist namsbc_wave'
+         WRITE(numout,*) '      Stokes drift                                  ln_sdw = ', ln_sdw
+         WRITE(numout,*) '      Breivik 2016                       ln_breivikFV_2016 = ', ln_breivikFV_2016
+         WRITE(numout,*) '      Stokes Coriolis & tracer advection terms    ln_stcor = ', ln_stcor
+         WRITE(numout,*) '      Vortex Force                         ln_vortex_force = ', ln_vortex_force
+         WRITE(numout,*) '      Bernouilli Head Pressure                ln_bern_srfc = ', ln_bern_srfc
+         WRITE(numout,*) '      wave modified ocean stress                  ln_tauoc = ', ln_tauoc
+         WRITE(numout,*) '      neutral drag coefficient (CORE bulk only)    ln_cdgw = ', ln_cdgw
+         WRITE(numout,*) '      charnock coefficient                        ln_charn = ', ln_charn
+         WRITE(numout,*) '      Stress modificated by wave                    ln_taw = ', ln_taw
+         WRITE(numout,*) '      TKE flux from wave                          ln_phioc = ', ln_phioc
+         WRITE(numout,*) '      Surface shear with Stokes drift           ln_stshear = ', ln_stshear
+         WRITE(numout,*) '      Test with constant wave fields          ln_wave_test = ', ln_wave_test
+      ENDIF
+
+      !                                ! option check
+      IF( .NOT.( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_charn) )   &
+         &     CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')
+      IF( ln_cdgw .AND. ln_blk )   &
+         &     CALL ctl_stop( 'drag coefficient read from wave model NOT available yet with aerobulk package')
+      IF( ln_stcor .AND. .NOT.ln_sdw )   &
+         &     CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+
+      !                             !==  Allocate wave arrays  ==!
+      ALLOCATE( ut0sd (jpi,jpj)    , vt0sd (jpi,jpj) )
+      ALLOCATE( hsw   (jpi,jpj)    , wmp   (jpi,jpj) )
+      ALLOCATE( wnum  (jpi,jpj) )
+      ALLOCATE( tsd2d (jpi,jpj)    , div_sd(jpi,jpj)    , bhd_wave(jpi,jpj)     )
+      ALLOCATE( usd   (jpi,jpj,jpk), vsd   (jpi,jpj,jpk), wsd     (jpi,jpj,jpk) )
+      ALLOCATE( tusd  (jpi,jpj)    , tvsd  (jpi,jpj)    , ZMX     (jpi,jpj,jpk) )
+      usd   (:,:,:) = 0._wp
+      vsd   (:,:,:) = 0._wp
+      wsd   (:,:,:) = 0._wp
+      hsw     (:,:) = 0._wp
+      wmp     (:,:) = 0._wp
+      ut0sd   (:,:) = 0._wp
+      vt0sd   (:,:) = 0._wp
+      tusd    (:,:) = 0._wp
+      tvsd    (:,:) = 0._wp
+      bhd_wave(:,:) = 0._wp
+      ZMX   (:,:,:) = 0._wp
+!
+      IF( ln_wave_test ) THEN       !==  Wave TEST case  ==!   set uniform waves fields
+         jpfld    = 0                   ! No field read
+         ln_cdgw  = .FALSE.             ! No neutral wave drag input
+         ln_tauoc = .FALSE.             ! No wave induced drag reduction factor
+         ut0sd(:,:) = 0.13_wp * tmask(:,:,1)   ! m/s
+         vt0sd(:,:) = 0.00_wp                  ! m/s
+         hsw  (:,:) = 2.80_wp                  ! meters
+         wmp  (:,:) = 8.00_wp                  ! seconds
+         !
+      ELSE                          !==  create the structure associated with fields to be read  ==!
+         IF( ln_cdgw ) THEN                       ! wave drag
+            IF( .NOT. cpl_wdrag ) THEN
+               ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+               !
+                                      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_init', 'Wave module ', 'namsbc_wave' )
+            ENDIF
+            ALLOCATE( cdn_wave(jpi,jpj) )
+            cdn_wave(:,:) = 0._wp
+         ENDIF
+         IF( ln_charn ) THEN                     ! wave drag
+            IF( .NOT. cpl_charn ) THEN
+               CALL ctl_stop( 'STOP', 'Charnock based wind stress can be used in coupled mode only' )
+            ENDIF
+            ALLOCATE( charn(jpi,jpj) )
+            charn(:,:) = 0._wp
+         ENDIF
+         IF( ln_taw ) THEN                     ! wind stress
+            IF( .NOT. cpl_taw ) THEN
+               CALL ctl_stop( 'STOP', 'wind stress from wave model can be used in coupled mode only, use ln_cdgw instead' )
+            ENDIF
+            ALLOCATE( tawx(jpi,jpj) )
+            ALLOCATE( tawy(jpi,jpj) )
+            ALLOCATE( twox(jpi,jpj) )
+            ALLOCATE( twoy(jpi,jpj) )
+            ALLOCATE( tauoc_wavex(jpi,jpj) )
+            ALLOCATE( tauoc_wavey(jpi,jpj) )
+            tawx(:,:) = 0._wp
+            tawy(:,:) = 0._wp
+            twox(:,:) = 0._wp
+            twoy(:,:) = 0._wp
+            tauoc_wavex(:,:) = 1._wp
+            tauoc_wavey(:,:) = 1._wp
+         ENDIF
+
+         IF( ln_phioc ) THEN                     ! TKE flux
+            IF( .NOT. cpl_phioc ) THEN
+                CALL ctl_stop( 'STOP', 'phioc can be used in coupled mode only' )
+            ENDIF
+            ALLOCATE( phioc(jpi,jpj) )
+            phioc(:,:) = 0._wp
+         ENDIF
+
+         IF( ln_tauoc ) THEN                    ! normalized wave stress into the ocean
+            IF( .NOT. cpl_wstrf ) 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_tauoc structure' )
+               !
+                                       ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
+               IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
+            ENDIF
+            ALLOCATE( tauoc_wave(jpi,jpj) )
+            tauoc_wave(:,:) = 0._wp
+         ENDIF
+
+         IF( ln_sdw ) THEN                      ! Stokes drift
+            ! 1. Find out how many fields have to be read from file if not coupled
+            jpfld=0
+            jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0
+            IF( .NOT. cpl_sdrftx ) THEN
+               jpfld  = jpfld + 1
+               jp_usd = jpfld
+            ENDIF
+            IF( .NOT. cpl_sdrfty ) THEN
+               jpfld  = jpfld + 1
+               jp_vsd = jpfld
+            ENDIF
+            IF( .NOT. cpl_hsig ) THEN
+               jpfld  = jpfld + 1
+               jp_hsw = jpfld
+            ENDIF
+            IF( .NOT. cpl_wper ) THEN
+               jpfld  = jpfld + 1
+               jp_wmp = jpfld
+            ENDIF
+            ! 2. Read from file only the non-coupled fields
+            IF( jpfld > 0 ) THEN
+               ALLOCATE( slf_i(jpfld) )
+               IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
+               IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
+               IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
+               IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
+               ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+               !
+               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_init', 'Wave module ', 'namsbc_wave' )
+            ENDIF
             !
-                                   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_init', 'Wave module ', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( cdn_wave(jpi,jpj) )
-      ENDIF
-
-      IF( ln_tauwoc ) THEN
-         IF( .NOT. cpl_tauwoc ) THEN
-            ALLOCATE( sf_tauwoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwoc
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+            ! 3. Wave number (only needed for Qiao parametrisation, ln_zdfqiao=T)
+            IF( .NOT. cpl_wnum ) THEN
+               ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
+               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' )
+                                      ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
+               IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            ENDIF
             !
-                                     ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   )
-            IF( sn_tauwoc%ln_tint )  ALLOCATE( sf_tauwoc(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_tauwoc, (/ sn_tauwoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
-         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
-
-      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled
-         jpfld=0
-         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
-         IF( .NOT. cpl_sdrftx ) THEN
-            jpfld  = jpfld + 1
-            jp_usd = jpfld
-         ENDIF
-         IF( .NOT. cpl_sdrfty ) THEN
-            jpfld  = jpfld + 1
-            jp_vsd = jpfld
-         ENDIF
-         IF( .NOT. cpl_hsig  .AND. ll_st_bv_li  ) THEN
-            jpfld  = jpfld + 1
-            jp_hsw = jpfld
-         ENDIF
-         IF( .NOT. cpl_wper  .AND. ll_st_bv_li  ) THEN
-            jpfld  = jpfld + 1
-            jp_wmp = jpfld
-         ENDIF
-         IF( .NOT. cpl_wfreq .AND. ll_st_peakfr ) THEN
-            jpfld  = jpfld + 1
-            jp_wfr = jpfld
-         ENDIF
-
-         ! Read from file only the non-coupled fields 
-         IF( jpfld > 0 ) THEN
-            ALLOCATE( slf_i(jpfld) )
-            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
-            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
-            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
-            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
-            IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
-
-            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
-            !
-            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_init', 'Wave module ', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
-         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     )
-         ALLOCATE( wfreq(jpi,jpj) )
-         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
-         ALLOCATE( div_sd(jpi,jpj) )
-         ALLOCATE( tsd2d (jpi,jpj) )
-
-         ut0sd(:,:) = 0._wp
-         vt0sd(:,:) = 0._wp
-         hsw(:,:) = 0._wp
-         wmp(:,:) = 0._wp
-
-         usd(:,:,:) = 0._wp
-         vsd(:,:,:) = 0._wp
-         wsd(:,:,:) = 0._wp
-         ! Wave number needed only if ln_zdfswm=T
-         IF( .NOT. cpl_wnum ) THEN
-            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' )
-                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
-            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
-         ENDIF
-         ALLOCATE( wnum(jpi,jpj) )
+         ENDIF
+         !
       ENDIF
       !

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/TRA/eosbn2.F90

@@ -56 +56 @@
    !                  !! * Interface
    INTERFACE eos
-      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
+      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
    END INTERFACE
    !
@@ -574 +574 @@
       !
    END SUBROUTINE eos_insitu_2d_t
+
+
+   SUBROUTINE eos_insitu_pot_2d( pts, prhop )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_pot  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
+      !!      potential volumic mass (Kg/m3) from potential temperature and
+      !!      salinity fields using an equation of state selected in the
+      !!     namelist.
+      !!
+      !! ** Action  :
+      !!              - prhop, the potential volumic mass (Kg/m3)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                                ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      !
+      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
+      INTEGER  ::   jdof
+      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-pot')
+      !
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+            DO_2D( 1, 1, 1, 1 )
+               !
+               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+               ztm = tmask(ji,jj,1)                                         ! tmask
+               !
+               zn0 = (((((EOS060*zt   &
+                  &   + EOS150*zs+EOS050)*zt   &
+                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                  !
+               !
+               prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
+               !
+            END_2D
+
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_2D( 1, 1, 1, 1 )
+            zt  = pts  (ji,jj,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jp_sal) - 35._wp
+            ztm = tmask(ji,jj,1)
+            !                                                     ! potential density referenced at the surface
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
+               &  - rn_nu * zt * zs
+            prhop(ji,jj) = ( rho0 + zn ) * ztm
+            !
+         END_2D
+         !
+      END SELECT
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prhop, clinfo1=' pot: ', kdim=1 )
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prhop, clinfo1=' eos-pot: ' )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-pot')
+      !
+   END SUBROUTINE eos_insitu_pot_2d
 
 

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/TRA/trazdf.F90

@@ -17 +17 @@
    USE phycst         ! physical constant
    USE zdf_oce        ! ocean vertical physics variables
+   USE zdfmfc         ! Mass FLux Convection
    USE sbc_oce        ! surface boundary condition: ocean
    USE ldftra         ! lateral diffusion: eddy diffusivity
@@ -198 +199 @@
             ENDIF
             !
+            ! Modification of diagonal to add MF scheme
+            IF ( ln_zdfmfc ) THEN
+               CALL diag_mfc( zwi, zwd, zws, p2dt, Kaa )
+            END IF
+            !
             !! Matrix inversion from the first level
             !!----------------------------------------------------------------------
@@ -225 +231 @@
             !
          ENDIF
+         !
+         ! Modification of rhs to add MF scheme
+         IF ( ln_zdfmfc ) THEN
+            CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
+         END IF
          !
          DO_2D( 0, 0, 0, 0 )         !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ZDF/zdf_oce.F90

@@ -40 +40 @@
    LOGICAL , PUBLIC ::   ln_zdfswm   !: surface  wave-induced mixing flag
    LOGICAL , PUBLIC ::   ln_zdfiwm   !: internal wave-induced mixing flag
-   !                             ! coefficients 
+   LOGICAL , PUBLIC ::   ln_zdfmfc   !: convection: eddy diffusivity Mass Flux Convection
+   !                             ! coefficients
    REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s)
    REAL(wp), PUBLIC ::   rn_avt0     !: vertical eddy diffusivity (m2/s)
@@ -55 +56 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -66 +67 @@
       !
       ALLOCATE( avm (jpi,jpj,jpk) , avm_k(jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,   &
-         &      avt (jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   & 
+         &      avt (jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   &
          &      avmb(jpk)         , avtb(jpk)          , avtb_2d(jpi,jpj) , STAT = zdf_oce_alloc )
          !

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ZDF/zdfphy.F90

@@ -9 +9 @@
    !!----------------------------------------------------------------------
    !!   zdf_phy_init  : initialization of all vertical physics packages
-   !!   zdf_phy       : upadate at each time-step the vertical mixing coeff. 
+   !!   zdf_phy       : upadate at each time-step the vertical mixing coeff.
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
-   USE zdf_oce        ! vertical physics: shared variables         
+   USE zdf_oce        ! vertical physics: shared variables
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
    USE zdfsh2         ! vertical physics: shear production term of TKE
-   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing   
+   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing
    USE zdftke         ! vertical physics: TKE vertical mixing
    USE zdfgls         ! vertical physics: GLS vertical mixing
    USE zdfosm         ! vertical physics: OSMOSIS vertical mixing
-   USE zdfddm         ! vertical physics: double diffusion mixing      
-   USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion  
-   USE zdfiwm         ! vertical physics: internal wave-induced mixing  
+   USE zdfddm         ! vertical physics: double diffusion mixing
+   USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion
+   USE zdfmfc         ! vertical physics: Mass Flux Convection
+   USE zdfiwm         ! vertical physics: internal wave-induced mixing
    USE zdfswm         ! vertical physics: surface  wave-induced mixing
    USE zdfmxl         ! vertical physics: mixed layer
    USE tranpc         ! convection: non penetrative adjustment
-   USE trc_oce        ! variables shared between passive tracer & ocean           
+   USE trc_oce        ! variables shared between passive tracer & ocean
    USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test)
    USE sbcrnf         ! surface boundary condition: runoff variables
@@ -45 +46 @@
    PUBLIC   zdf_phy       ! called by step.F90
 
-   INTEGER ::   nzdf_phy   ! type of vertical closure used 
+   INTEGER ::   nzdf_phy   ! type of vertical closure used
    !                       ! associated indicators
    INTEGER, PARAMETER ::   np_CST = 1   ! Constant Kz
@@ -65 +66 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_phy_init  ***
-      !! 
+      !!
       !! ** Purpose :   initializations of the vertical ocean physics
       !!
-      !! ** Method  :   Read namelist namzdf, control logicals 
+      !! ** Method  :   Read namelist namzdf, control logicals
       !!                set horizontal shape and vertical profile of background mixing coef.
       !!----------------------------------------------------------------------
@@ -78 +79 @@
       NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
          &             ln_zdfosm,                                    &     ! type of closure scheme
+         &             ln_zdfmfc,                                    &     ! convection : mass flux
          &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
          &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
@@ -112 +114 @@
          WRITE(numout,*) '         OSMOSIS-OBL closure (OSM)               ln_zdfosm = ', ln_zdfosm
          WRITE(numout,*) '      convection: '
+         WRITE(numout,*) '         convection mass flux (mfc)              ln_zdfmfc = ', ln_zdfmfc
          WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
          WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
@@ -140 +143 @@
       IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter
          avmb(:) = rn_avm0
-         avtb(:) = rn_avt0                     
+         avtb(:) = rn_avt0
       ELSE                               ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
          avmb(:) = rn_avm0
@@ -147 +150 @@
       ENDIF
       !                                  ! 2D shape of the avtb
-      avtb_2d(:,:) = 1._wp                   ! uniform 
+      avtb_2d(:,:) = 1._wp                   ! uniform
       !
       IF( nn_havtb == 1 ) THEN               ! decrease avtb by a factor of ten in the equatorial band
@@ -172 +175 @@
       IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
       IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' )
+      IF( ln_zdfmfc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfevd' )
+      IF( ln_zdfmfc .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfnpc' )
+      IF( ln_zdfmfc .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfosm' )
       IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
       IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )
+      IF( lk_top    .AND. ln_zdfmfc )   CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' )
       IF(lwp) THEN
          WRITE(numout,*)
          IF    ( ln_zdfnpc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use non penetrative convective scheme'
          ELSEIF( ln_zdfevd ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use enhanced vertical diffusion scheme'
+         ELSEIF( ln_zdfmfc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use Mass Flux scheme'
          ELSE                      ;   WRITE(numout,*) '   ==>>>   convection: no specific scheme used'
          ENDIF
@@ -190 +198 @@
 
       !                          !==  type of vertical turbulent closure  ==!    (set nzdf_phy)
-      ioptio = 0 
+      ioptio = 0
       IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
       IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init          ;   ENDIF
@@ -205 +213 @@
       ELSE                   ;   l_zdfsh2 = .TRUE.
       ENDIF
-
+      !                          !== Mass Flux Convectiive algorithm  ==!
+      IF( ln_zdfmfc )   CALL zdf_mfc_init       ! Convection computed with eddy diffusivity mass flux
+      !
       !                          !== gravity wave-driven mixing  ==!
       IF( ln_zdfiwm )   CALL zdf_iwm_init       ! internal wave-driven mixing
@@ -226 +236 @@
       !! ** Purpose :  Update ocean physics at each time-step
       !!
-      !! ** Method  : 
+      !! ** Method  :
       !!
       !! ** Action  :   avm, avt vertical eddy viscosity and diffusivity at w-points
@@ -244 +254 @@
          !
          !                       !* bottom drag
-         CALL zdf_drg( kt, Kmm, mbkt , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in 
+         CALL zdf_drg( kt, Kmm, mbkt , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in
             &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   &
             &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s]
          IF( ln_isfcav ) THEN    !* top drag   (ocean cavities)
-            CALL zdf_drg( kt, Kmm, mikt , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in 
+            CALL zdf_drg( kt, Kmm, mikt , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in
                &              r_z0_top,   r_ke0_top,    rCd0_top,   &
                &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s]
@@ -263 +273 @@
       ENDIF
 #endif
-      ! 
+      !
       !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
       !
@@ -280 +290 @@
 !!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
       END SELECT
-      !  
+      !
       !                          !==  ocean Kz  ==!   (avt, avs, avm)
       !
@@ -302 +312 @@
       ENDIF
       !
-      !                                         !* wave-induced mixing 
-      IF( ln_zdfswm )   CALL zdf_swm( kt, Kmm, avm, avt, avs )   ! surface  wave (Qiao et al. 2004) 
+      !                                         !* wave-induced mixing
+      IF( ln_zdfswm )   CALL zdf_swm( kt, Kmm, avm, avt, avs )   ! surface  wave (Qiao et al. 2004)
       IF( ln_zdfiwm )   CALL zdf_iwm( kt, Kmm, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
 
-#if defined key_agrif 
+#if defined key_agrif
       ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
       IF( l_zdfsh2 )   CALL Agrif_avm
@@ -330 +340 @@
          IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
          IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
-         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' ) 
+         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' )
          ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated
       ENDIF

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ZDF/zdfsh2.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdfsh2  ***
-   !! Ocean physics:  shear production term of TKE 
+   !! Ocean physics:  shear production term of TKE
    !!=====================================================================
    !! History :   -   !  2014-10  (A. Barthelemy, G. Madec)  original code
    !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm
+   !!   NEMO     4.2  !  2020-12  (G. Madec, E. Clementi) add Stokes Drift Shear
+   !                  !           for wave coupling
    !!----------------------------------------------------------------------
 
@@ -13 +15 @@
    USE oce
    USE dom_oce        ! domain: ocean
+   USE sbcwave        ! Surface Waves (add Stokes shear)
+   USE sbc_oce , ONLY: ln_stshear  !Stoked Drift shear contribution
    !
    USE in_out_manager ! I/O manager
@@ -21 +25 @@
 
    PUBLIC   zdf_sh2        ! called by zdftke, zdfglf, and zdfric
-   
+
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -32 +36 @@
 CONTAINS
 
-   SUBROUTINE zdf_sh2( Kbb, Kmm, p_avm, p_sh2  ) 
+   SUBROUTINE zdf_sh2( Kbb, Kmm, p_avm, p_sh2  )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_sh2  ***
@@ -40 +44 @@
       !! ** Method  : - a stable discretization of this term is linked to the
       !!                time-space discretization of the vertical diffusion
-      !!                of the OGCM. NEMO uses C-grid, a leap-frog environment 
+      !!                of the OGCM. NEMO uses C-grid, a leap-frog environment
       !!                and an implicit computation of vertical mixing term,
       !!                so the shear production at w-point is given by:
-      !!                   sh2 = mi[   mi(avm) * dk[ub]/e3ub * dk[un]/e3un   ] 
-      !!                       + mj[   mj(avm) * dk[vb]/e3vb * dk[vn]/e3vn   ] 
+      !!                   sh2 = mi[   mi(avm) * dk[ub]/e3ub * dk[un]/e3un   ]
+      !!                       + mj[   mj(avm) * dk[vb]/e3vb * dk[vn]/e3vn   ]
       !!                NB: wet-point only horizontal averaging of shear
       !!
@@ -59 +63 @@
       !!--------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1
-         DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
-            zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
-               &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
-               &         * (   uu(ji,jj,jk-1,Kbb) -   uu(ji,jj,jk,Kbb) ) & 
-               &         / ( e3uw(ji,jj,jk  ,Kmm) * e3uw(ji,jj,jk,Kbb) ) &
-               &         * wumask(ji,jj,jk)
-            zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
-               &         * (   vv(ji,jj,jk-1,Kmm) -   vv(ji,jj,jk,Kmm) ) &
-               &         * (   vv(ji,jj,jk-1,Kbb) -   vv(ji,jj,jk,Kbb) ) &
-               &         / ( e3vw(ji,jj,jk  ,Kmm) * e3vw(ji,jj,jk,Kbb) ) &
-               &         * wvmask(ji,jj,jk)
-         END_2D
+      DO jk = 2, jpkm1                 !* Shear production at uw- and vw-points (energy conserving form)
+         IF ( cpl_sdrftx .AND. ln_stshear )  THEN       ! Surface Stokes Drift available  ===>>>  shear + stokes drift contibution
+            DO_2D( 1, 0, 1, 0 )
+               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) )        &
+                  &         * ( uu (ji,jj,jk-1,Kmm) -   uu (ji,jj,jk,Kmm)    &
+                  &           + usd(ji,jj,jk-1) -   usd(ji,jj,jk) )  &
+                  &         * ( uu (ji,jj,jk-1,Kbb) -   uu (ji,jj,jk,Kbb) )  &
+                  &         / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk)
+               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) )         &
+                  &         * ( vv (ji,jj,jk-1,Kmm) -   vv (ji,jj,jk,Kmm)     &
+                  &           + vsd(ji,jj,jk-1) -   vsd(ji,jj,jk) )   &
+                  &         * ( vv (ji,jj,jk-1,Kbb) -   vv (ji,jj,jk,Kbb) )   &
+                  &/ ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk)
+            END_2D
+         ELSE
+            DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
+               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
+                  &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
+                  &         * (   uu(ji,jj,jk-1,Kbb) -   uu(ji,jj,jk,Kbb) ) &
+                  &         / ( e3uw(ji,jj,jk  ,Kmm) * e3uw(ji,jj,jk,Kbb) ) &
+                  &         * wumask(ji,jj,jk)
+               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
+                  &         * (   vv(ji,jj,jk-1,Kmm) -   vv(ji,jj,jk,Kmm) ) &
+                  &         * (   vv(ji,jj,jk-1,Kbb) -   vv(ji,jj,jk,Kbb) ) &
+                  &         / ( e3vw(ji,jj,jk  ,Kmm) * e3vw(ji,jj,jk,Kbb) ) &
+                  &         * wvmask(ji,jj,jk)
+            END_2D
+         ENDIF
          DO_2D( 0, 0, 0, 0 )     !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
             p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
                &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
          END_2D
-      END DO 
+      END DO
       !
    END SUBROUTINE zdf_sh2

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ZDF/zdftke.F90

@@ -29 +29 @@
    !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
    !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition
+   !!            4.2  !  2020-12  (G. Madec, E. Clementi) add wave coupling
+   !                  !           following Couvelard et al., 2019
    !!----------------------------------------------------------------------
 
@@ -58 +60 @@
    USE prtctl         ! Print control
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE sbcwave        ! Surface boundary waves
 
    IMPLICIT NONE
@@ -68 +71 @@
    !                      !!** Namelist  namzdf_tke  **
    LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
+   LOGICAL  ::   ln_mxhsw  ! mixing length scale surface value as a fonction of wave height
    INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3)
    REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice
@@ -81 +85 @@
    INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
    INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
+   INTEGER  ::   nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
+   INTEGER  ::   nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
    REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
@@ -209 +215 @@
       REAL(wp) ::   zus   , zwlc  , zind       !   -      -
       REAL(wp) ::   zzd_up, zzd_lw             !   -      -
+      REAL(wp) ::   ztaui, ztauj, z1_norm
       INTEGER , DIMENSION(jpi,jpj)     ::   imlc
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3, zWlc2
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
@@ -219 +226 @@
       zfact2  = 1.5_wp * rn_Dt * rn_ediss
       zfact3  = 0.5_wp         * rn_ediss
+      !
+      zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used
       !
       ! ice fraction considered for attenuation of langmuir & wave breaking
@@ -232 +241 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      DO_2D( 0, 0, 0, 0 )         ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-!! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly
-!!       one way around would be to increase zbbirau
-!!          en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + &
-!!             &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1)
+      DO_2D( 0, 0, 0, 0 )
          en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+         zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
+         zd_lw(ji,jj,1) = 1._wp
+         zd_up(ji,jj,1) = 0._wp
       END_2D
       !
@@ -274 +282 @@
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke (Axell JGR 2002)
          !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
          !
-         !                        !* total energy produce by LC : cumulative sum over jk
+         !                       !* Langmuir velocity scale
+         !
+         IF ( cpl_sdrftx )  THEN       ! Surface Stokes Drift available
+            !                                ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2    with u* = (taum/rho0)^1/2
+            !                                ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s
+            !                                ! more precisely, it is the dot product that must be used :
+            !                                !     1/2  (W_lc)^2 = MAX( u* u_s + v* v_s , 0 )   only the positive part
+!!gm  ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
+!!gm  ! so we will overestimate the LC velocity....   !!gm I will do the work if !LC have an effect !
+            DO_2D( 0, 0, 0, 0 )
+!!XC                  zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )  )
+                  zWlc2(ji,jj) = 0.5_wp *  ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )
+            END_2D
+!
+!  Projection of Stokes drift in the wind stress direction
+!
+            DO_2D( 0, 0, 0, 0 )
+                  ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
+                  ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
+                  z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1)
+                  zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
+            END_2D
+         CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. )
+!
+         ELSE                          ! Surface Stokes drift deduced from surface stress
+            !                                ! Wlc = u_s   with u_s = 0.016*U_10m, the surface stokes drift  (Axell 2002, Eq.44)
+            !                                ! using |tau| = rho_air Cd |U_10m|^2 , it comes:
+            !                                ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2   and thus:
+            !                                ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag )
+            zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag
+            DO_2D( 1, 1, 1, 1 )
+               zWlc2(ji,jj) = zcof * taum(ji,jj)
+            END_2D
+            !
+         ENDIF
+         !
+         !                       !* Depth of the LC circulation  (Axell 2002, Eq.47)
+         !                             !- LHS of Eq.47
          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm)
          DO jk = 2, jpk
@@ -283 +328 @@
                &        MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
          END DO
-         !                        !* finite Langmuir Circulation depth
-         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
+         !
+         !                             !- compare LHS to RHS of Eq.47
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
-         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! Last w-level at which zpelc>=0.5*us*us
-            zus = zcof * taum(ji,jj)          !      with us=0.016*wind(starting from jpk-1)
-            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
+         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
+            IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk
          END_3D
          !                               ! finite LC depth
@@ -294 +338 @@
             zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
          END_2D
+         !
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
          DO_2D( 0, 0, 0, 0 )
-            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
+            zus = SQRT( 2. * zWlc2(ji,jj) )             ! Stokes drift
             zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
          END_2D
@@ -351 +396 @@
             &                                ) * wmask(ji,jj,jk)
       END_3D
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface boundary condition on tke if
+      !                     !  coupling with waves
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF ( cpl_phioc .and. ln_phioc )  THEN
+         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
+
+         CASE ( 0 ) ! Dirichlet BC
+            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
+               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
+            END_2D
+
+         CASE ( 1 ) ! Neumann BC
+            DO_2D( 0, 0, 0, 0 )
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
+               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
+               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
+               zdiag(ji,jj,1) = 1._wp
+               zd_lw(ji,jj,2) = 0._wp
+            END_2D
+
+         END SELECT
+
+      ENDIF
+      !
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO_3D( 0, 0, 0, 0, 3, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_2D( 0, 0, 0, 0 )                          ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
-      END_2D
-      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
+!XC : commented to allow for neumann boundary condition
+!      DO_2D( 0, 0, 0, 0 )
+!         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
+!      END_2D
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
       END_3D
@@ -461 +537 @@
       zmxld(:,:,:)  = rmxl_min
       !
-     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
-         !
-         zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
+      IF(ln_sdw .AND. ln_mxhsw) THEN
+         zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
+         ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6
+         zcoef       = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff
+         zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
+      ELSE
+      !
+         IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
+         !
+            zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
 #if ! defined key_si3 && ! defined key_cice
-         DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
-            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
-         END_2D
+            DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
+               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
+            END_2D
 #else
-         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
-         !
-         CASE( 0 )                      ! No scaling under sea-ice
+            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
+            !
+            CASE( 0 )                      ! No scaling under sea-ice
+               DO_2D( 0, 0, 0, 0 )
+                  zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
+               END_2D
+               !
+            CASE( 1 )                      ! scaling with constant sea-ice thickness
+               DO_2D( 0, 0, 0, 0 )
+                  zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
+               END_2D
+               !
+            CASE( 2 )                      ! scaling with mean sea-ice thickness
+               DO_2D( 0, 0, 0, 0 )
+#if defined key_si3
+                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
+#elif defined key_cice
+                  zmaxice = MAXVAL( h_i(ji,jj,:) )
+                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
+#endif
+               END_2D
+               !
+            CASE( 3 )                      ! scaling with max sea-ice thickness
+               DO_2D( 0, 0, 0, 0 )
+                  zmaxice = MAXVAL( h_i(ji,jj,:) )
+                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
+                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
+               END_2D
+               !
+            END SELECT
+#endif
+            !
             DO_2D( 0, 0, 0, 0 )
-               zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
+               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
             END_2D
             !
-         CASE( 1 )                      ! scaling with constant sea-ice thickness
-            DO_2D( 0, 0, 0, 0 )
-               zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
-                  &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
-            END_2D
-            !
-         CASE( 2 )                      ! scaling with mean sea-ice thickness
-            DO_2D( 0, 0, 0, 0 )
-#if defined key_si3
-               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
-                  &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
-#elif defined key_cice
-               zmaxice = MAXVAL( h_i(ji,jj,:) )
-               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
-                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
-#endif
-            END_2D
-            !
-         CASE( 3 )                      ! scaling with max sea-ice thickness
-            DO_2D( 0, 0, 0, 0 )
-               zmaxice = MAXVAL( h_i(ji,jj,:) )
-               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
-                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
-            END_2D
-            !
-         END SELECT
-#endif
-         !
-         DO_2D( 0, 0, 0, 0 )
-            zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
-         END_2D
-         !
-      ELSE
-         zmxlm(:,:,1) = rn_mxl0
-      ENDIF
-
+         ELSE
+            zmxlm(:,:,1) = rn_mxl0
+         ENDIF
+      ENDIF
       !
       DO_3D( 0, 0, 0, 0, 2, jpkm1 )
@@ -624 +707 @@
          &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
          &                 nn_pdl  , ln_lc    , rn_lc    ,             &
-         &                 nn_etau , nn_htau  , rn_efr   , nn_eice
+         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  ,  &
+         &                 nn_bc_surf, nn_bc_bot, ln_mxhsw
       !!----------------------------------------------------------------------
       !
@@ -666 +750 @@
          WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
          WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
+         IF ( cpl_phioc .and. ln_phioc )  THEN
+            SELECT CASE( nn_bc_surf)             ! Type of scaling under sea-ice
+            CASE( 0 )   ;   WRITE(numout,*) '  nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
+            CASE( 1 )   ;   WRITE(numout,*) '  nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
+            END SELECT
+         ENDIF
          WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
          WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/step.F90

@@ -296 +296 @@
 
                             CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection ==> RHS
+         IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection
          IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
          IF( lrst_oce .AND. ln_zdfosm ) &

NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/step_oce.F90

@@ -70 +70 @@
    USE zdfphy          ! vertical physics manager      (zdf_phy_init routine)
    USE zdfosm  , ONLY : osm_rst, dyn_osm, tra_osm      ! OSMOSIS routines used in step.F90
+   USE zdfmfc          ! Mass FLux Convection routine used in step.F90
 
    USE diu_layers      ! diurnal SST bulk and coolskin routines