New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7562 for branches/2016 – NEMO

Changeset 7562 for branches/2016


Ignore:
Timestamp:
2017-01-16T15:52:00+01:00 (7 years ago)
Author:
emanuelaclementi
Message:

#1811 Merge with dev_INGV_UKMO_2016: improve Stokes drift (including dynspg_ts , Stokes-Coriolis force , and GLS surface roughness)

Location:
branches/2016/dev_merge_2016/NEMOGCM
Files:
1 deleted
13 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_merge_2016/NEMOGCM/CONFIG/SHARED/field_def.xml

    r7528 r7562  
    415415         <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_U_3D" /> 
    416416         <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_U_3D" /> 
     417         <field id="ustokes"      long_name="Stokes Drift Velocity i-axis"                           standard_name="StokesDrift_x_velocity"      unit="m/s"        grid_ref="grid_U_3D" /> 
     418         <field id="ustokes_e3u"  long_name="Stokes Drift Velocity i-axis  (thickness weighted)"                                                 unit="m/s"        grid_ref="grid_U_3D"  > ustokes * e3u </field> 
    417419 
    418420         <!-- u-eddy coefficients (ldftra) --> 
     
    463465         <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_V_3D" /> 
    464466         <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_V_3D" /> 
     467         <field id="vstokes"      long_name="Stokes Drift Velocity j-axis"                           standard_name="StokesDrift_y_velocity"      unit="m/s"        grid_ref="grid_V_3D" /> 
     468         <field id="vstokes_e3v"  long_name="Stokes Drift Velocity j-axis  (thickness weighted)"                                                 unit="m/s"        grid_ref="grid_V_3D"  > vstokes * e3v </field> 
    465469 
    466470         <!-- v-eddy coefficients (ldftra, ldfdyn) --> 
     
    503507        <field id="woce"         long_name="ocean vertical velocity"              standard_name="upward_sea_water_velocity"   unit="m/s"  /> 
    504508        <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                               unit="m3/s" /> 
     509        <field id="wstokes"      long_name="Stokes Drift vertical velocity"       standard_name="upward_StokesDrift_velocity" unit="m/s" /> 
    505510 
    506511        <!-- woce_eiv: available with EIV --> 
  • branches/2016/dev_merge_2016/NEMOGCM/CONFIG/SHARED/namelist_ref

    r7533 r7562  
    447447/ 
    448448!----------------------------------------------------------------------- 
    449 &namsbc_wave   ! External fields from wave model 
     449&namsbc_wave   ! External fields from wave model                        (ln_wave=T) 
    450450!----------------------------------------------------------------------- 
    451451!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     
    454454   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    455455   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    456    sn_swh      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
     456   sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    457457   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    458458   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
     
    948948   rn_hsro       =  0.02   !  Minimum surface roughness 
    949949   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2) 
    950    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2) 
     950   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3) 
     951   !                             ! =3 requires ln_wave=T 
    951952   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum) 
    952953   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum) 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r7514 r7562  
    4242   USE diahth          ! thermocline diagnostics 
    4343   USE wet_dry         ! wetting and drying 
     44   USE sbcwave         ! wave parameters 
    4445   ! 
    4546   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    718719         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un 
    719720            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
     721         IF( ln_wave .AND. ln_sdw) THEN 
     722            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd 
     723               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
     724         ENDIF 
    720725         !                                                                                      !!! nid_U : 2D 
    721726         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau 
     
    727732         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn 
    728733            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
     734         IF( ln_wave .AND. ln_sdw) THEN 
     735            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd 
     736               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
     737         ENDIF 
    729738         !                                                                                      !!! nid_V : 2D 
    730739         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau 
     
    743752         IF( lk_zdfddm ) THEN 
    744753            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs 
     754               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
     755         ENDIF 
     756          
     757         IF( ln_wave .AND. ln_sdw) THEN 
     758            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd 
    745759               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    746760         ENDIF 
     
    865879      ENDIF 
    866880 
     881      IF( ln_wave .AND. ln_sdw ) THEN 
     882         CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current 
     883         CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current 
     884         CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current 
     885      ENDIF 
     886 
    867887      ! 3. Close all files 
    868888      ! --------------------------------------- 
     
    979999            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
    9801000      ENDIF 
     1001      ! 
     1002      IF( ln_wave .AND. ln_sdw ) THEN 
     1003         CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal"    , "m/s"    , &   ! StokesDrift zonal current 
     1004            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
     1005         CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid"    , "m/s"    , &   ! StokesDrift meridonal current 
     1006            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
     1007         CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert"     , "m/s"    , &   ! StokesDrift vertical current 
     1008            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
     1009      ENDIF 
    9811010 
    9821011#if defined key_lim2 
     
    10231052         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:) , jpi*jpj*jpk, idex )!  T-cell thickness   
    10241053      END IF  
     1054  
     1055      IF( ln_wave .AND. ln_sdw ) THEN 
     1056         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity 
     1057         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity 
     1058         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity 
     1059      ENDIF 
     1060 
    10251061      ! 3. Close the file 
    10261062      ! ----------------- 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7527 r7562  
    1515   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1616   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification 
     17   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence 
    1718   !!--------------------------------------------------------------------- 
    1819 
     
    3839   USE sbctide         ! tides 
    3940   USE updtide         ! tide potential 
     41   USE sbcwave         ! surface wave 
    4042   ! 
    4143   USE in_out_manager  ! I/O manager 
     
    343345         END DO 
    344346      END DO 
     347       
     348!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 
     349!!gm  Is it correct to do so ?   I think so... 
     350       
     351       
    345352      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    346353      !                                   ! -------------------------------------------------------- 
     
    492499      ! 
    493500      !                                         ! Add top stress contribution from baroclinic velocities:       
    494       IF (ln_bt_fw) THEN 
     501      IF( ln_bt_fw ) THEN 
    495502         DO jj = 2, jpjm1 
    496503            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    556563                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    557564      ENDIF 
     565      ! 
     566      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
     567         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
     568      ENDIF 
     569      ! 
    558570#if defined key_asminc 
    559571      !                                         ! Include the IAU weighted SSH increment 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r7421 r7562  
    1717   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity  
    1818   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory 
     19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    3233   USE trd_oce        ! trends: ocean variables 
    3334   USE trddyn         ! trend manager: dynamics 
     35   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force) 
     36   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force 
    3437   ! 
    3538   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    7780#  include "vectopt_loop_substitute.h90" 
    7881   !!---------------------------------------------------------------------- 
    79    !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
     82   !! NEMO/OPA 3.7 , NEMO Consortium (2016) 
    8083   !! $Id$ 
    8184   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    108111            ztrdu(:,:,:) = ua(:,:,:) 
    109112            ztrdv(:,:,:) = va(:,:,:) 
    110             CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     113            CALL vor_ene( kt, nrvm, un , vn , ua, va )                    ! relative vorticity or metric trend 
    111114            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    112115            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    114117            ztrdu(:,:,:) = ua(:,:,:) 
    115118            ztrdv(:,:,:) = va(:,:,:) 
    116             CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend 
     119            CALL vor_ene( kt, ncor, un , vn , ua, va )                    ! planetary vorticity trend 
    117120            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    118121            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    119122            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    120          ELSE 
    121             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend 
     123         ELSE                                               ! total vorticity trend 
     124                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     125            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    122126         ENDIF 
    123127         ! 
     
    126130            ztrdu(:,:,:) = ua(:,:,:) 
    127131            ztrdv(:,:,:) = va(:,:,:) 
    128             CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     132            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend 
    129133            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    130134            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    132136            ztrdu(:,:,:) = ua(:,:,:) 
    133137            ztrdv(:,:,:) = va(:,:,:) 
    134             CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend 
     138            CALL vor_ens( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend 
    135139            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    136140            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    137141            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    138          ELSE 
    139             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend 
     142         ELSE                                               ! total vorticity trend 
     143                             CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend 
     144            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    140145         ENDIF 
    141146         ! 
     
    144149            ztrdu(:,:,:) = ua(:,:,:) 
    145150            ztrdv(:,:,:) = va(:,:,:) 
    146             CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens) 
     151            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend (ens) 
    147152            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    148153            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    150155            ztrdu(:,:,:) = ua(:,:,:) 
    151156            ztrdv(:,:,:) = va(:,:,:) 
    152             CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene) 
     157            CALL vor_ene( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend (ene) 
    153158            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    154159            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    155160            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    156          ELSE 
    157             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens) 
    158             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene) 
     161         ELSE                                               ! total vorticity trend 
     162                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens) 
     163                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene) 
     164            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    159165        ENDIF 
    160166         ! 
     
    163169            ztrdu(:,:,:) = ua(:,:,:) 
    164170            ztrdv(:,:,:) = va(:,:,:) 
    165             CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     171            CALL vor_een( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend 
    166172            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    167173            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    169175            ztrdu(:,:,:) = ua(:,:,:) 
    170176            ztrdv(:,:,:) = va(:,:,:) 
    171             CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend 
     177            CALL vor_een( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend 
    172178            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    173179            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    174180            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    175          ELSE 
    176             CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend 
     181         ELSE                                               ! total vorticity trend 
     182                             CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     183            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    177184         ENDIF 
    178185         ! 
     
    190197 
    191198 
    192    SUBROUTINE vor_ene( kt, kvor, pua, pva ) 
     199   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 
    193200      !!---------------------------------------------------------------------- 
    194201      !!                  ***  ROUTINE vor_ene  *** 
     
    210217      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    211218      !!---------------------------------------------------------------------- 
    212       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    213       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    214       !                                                           ! =nrvm (relative vorticity or metric) 
    215       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    216       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     219      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
     220      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
     221      !                                                                ! =nrvm (relative vorticity or metric) 
     222      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
     223      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend 
    217224      ! 
    218225      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     
    223230      IF( nn_timing == 1 )  CALL timing_start('vor_ene') 
    224231      ! 
    225       CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )  
     232      CALL wrk_alloc( jpi,jpj,  zwx, zwy, zwz )  
    226233      ! 
    227234      IF( kt == nit000 ) THEN 
     
    241248            DO jj = 1, jpjm1 
    242249               DO ji = 1, fs_jpim1   ! vector opt. 
    243                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    244                      &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     250                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     251                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    245252               END DO 
    246253            END DO 
     
    248255            DO jj = 1, jpjm1 
    249256               DO ji = 1, fs_jpim1   ! vector opt. 
    250                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    251                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     257                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     258                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    252259                       &     * 0.5 * r1_e1e2f(ji,jj) 
    253260               END DO 
     
    256263            DO jj = 1, jpjm1 
    257264               DO ji = 1, fs_jpim1   ! vector opt. 
    258                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    259                      &                        - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
    260                      &                     * r1_e1e2f(ji,jj) 
     265                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     266                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     267                     &                   * r1_e1e2f(ji,jj) 
    261268               END DO 
    262269            END DO 
     
    264271            DO jj = 1, jpjm1 
    265272               DO ji = 1, fs_jpim1   ! vector opt. 
    266                   zwz(ji,jj) = ff_f(ji,jj)                                                                      & 
    267                        &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    268                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     273                  zwz(ji,jj) = ff_f(ji,jj)                                                                        & 
     274                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     275                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    269276                       &     * 0.5 * r1_e1e2f(ji,jj) 
    270277               END DO 
     
    284291         IF( ln_sco ) THEN 
    285292            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    286             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    287             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     293            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     294            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
    288295         ELSE 
    289             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    290             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
     296            zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
     297            zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
    291298         ENDIF 
    292299         !                                   !==  compute and add the vorticity term trend  =! 
     
    311318 
    312319 
    313    SUBROUTINE vor_ens( kt, kvor, pua, pva ) 
     320   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 
    314321      !!---------------------------------------------------------------------- 
    315322      !!                ***  ROUTINE vor_ens  *** 
     
    331338      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    332339      !!---------------------------------------------------------------------- 
    333       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    334       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    335          !                                                        ! =nrvm (relative vorticity or metric) 
    336       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    337       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     340      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
     341      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
     342         !                                                             ! =nrvm (relative vorticity or metric) 
     343      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
     344      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend 
    338345      ! 
    339346      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    344351      IF( nn_timing == 1 )  CALL timing_start('vor_ens') 
    345352      ! 
    346       CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )  
     353      CALL wrk_alloc( jpi,jpj,  zwx, zwy, zwz )  
    347354      ! 
    348355      IF( kt == nit000 ) THEN 
     
    361368            DO jj = 1, jpjm1 
    362369               DO ji = 1, fs_jpim1   ! vector opt. 
    363                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    364                      &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     370                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     371                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    365372               END DO 
    366373            END DO 
     
    368375            DO jj = 1, jpjm1 
    369376               DO ji = 1, fs_jpim1   ! vector opt. 
    370                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    371                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     377                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     378                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    372379                       &     * 0.5 * r1_e1e2f(ji,jj) 
    373380               END DO 
     
    376383            DO jj = 1, jpjm1 
    377384               DO ji = 1, fs_jpim1   ! vector opt. 
    378                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    379                      &                        - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
    380                      &                     * r1_e1e2f(ji,jj) 
     385                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     386                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     387                     &                   * r1_e1e2f(ji,jj) 
    381388               END DO 
    382389            END DO 
     
    384391            DO jj = 1, jpjm1 
    385392               DO ji = 1, fs_jpim1   ! vector opt. 
    386                   zwz(ji,jj) = ff_f(ji,jj)                                                                      & 
    387                        &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    388                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     393                  zwz(ji,jj) = ff_f(ji,jj)                                                                       & 
     394                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     395                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    389396                       &     * 0.5 * r1_e1e2f(ji,jj) 
    390397               END DO 
     
    404411         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    405412            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    406             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    407             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     413            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     414            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
    408415         ELSE 
    409             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    410             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
     416            zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
     417            zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
    411418         ENDIF 
    412419         !                                   !==  compute and add the vorticity term trend  =! 
     
    431438 
    432439 
    433    SUBROUTINE vor_een( kt, kvor, pua, pva ) 
     440   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 
    434441      !!---------------------------------------------------------------------- 
    435442      !!                ***  ROUTINE vor_een  *** 
     
    448455      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    449456      !!---------------------------------------------------------------------- 
    450       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    451       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 
    452       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    453       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     457      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
     458      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
     459         !                                                             ! =nrvm (relative vorticity or metric) 
     460      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
     461      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend 
    454462      ! 
    455463      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    512520            DO jj = 1, jpjm1 
    513521               DO ji = 1, fs_jpim1   ! vector opt. 
    514                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    515                      &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     522                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     523                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    516524                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
    517525               END DO 
     
    520528            DO jj = 1, jpjm1 
    521529               DO ji = 1, fs_jpim1   ! vector opt. 
    522                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    523                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     530                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     531                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    524532                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
    525533               END DO 
     
    528536            DO jj = 1, jpjm1 
    529537               DO ji = 1, fs_jpim1   ! vector opt. 
    530                   zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    531                      &                           - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
    532                      &                        * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj) 
     538                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     539                     &                           - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     540                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj) 
    533541               END DO 
    534542            END DO 
     
    536544            DO jj = 1, jpjm1 
    537545               DO ji = 1, fs_jpim1   ! vector opt. 
    538                   zwz(ji,jj) = (  ff_f(ji,jj)                                                                      & 
    539                        &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    540                        &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     546                  zwz(ji,jj) = (  ff_f(ji,jj)                                                                        & 
     547                       &        + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     548                       &            - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    541549                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    542550               END DO 
     
    557565         ! 
    558566         !                                   !==  horizontal fluxes  ==! 
    559          zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    560          zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     567         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     568         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
    561569 
    562570         !                                   !==  compute and add the vorticity term trend  =! 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7421 r7562  
    11331133      !                                                      !       Stokes drift u      ! 
    11341134      !                                                      ! ========================= !  
    1135          IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
     1135         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
    11361136      ! 
    11371137      !                                                      ! ========================= !  
    11381138      !                                                      !       Stokes drift v      ! 
    11391139      !                                                      ! ========================= !  
    1140          IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
     1140         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
    11411141      ! 
    11421142      !                                                      ! ========================= !  
     
    11481148      !                                                      !  Significant wave height  ! 
    11491149      !                                                      ! ========================= !  
    1150          IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
     1150         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
    11511151      ! 
    11521152      !                                                      ! ========================= !  
     
    11591159                                                                    .OR. srcv(jpr_hsig)%laction ) THEN 
    11601160            CALL sbc_stokes() 
    1161             IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 
    1162          ENDIF 
    1163          IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 
     1161         ENDIF 
    11641162      ENDIF 
    11651163      !                                                      ! ========================= !  
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7431 r7562  
    334334      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )   ! CICE initialization 
    335335      ! 
     336      IF( ln_wave     )   CALL sbc_wave_init              ! surface wave initialisation 
     337      ! 
    336338   END SUBROUTINE sbc_init 
    337339 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7421 r7562  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3  !   2011-09  (M. Adani)  Original code: Drag Coefficient  
    7    !!         :  3.4  !   2012-10  (M. Adani)  Stokes Drift  
    8    !!            3.6  !   2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation 
    9    !!---------------------------------------------------------------------- 
    10  
    11    !!---------------------------------------------------------------------- 
     6   !! History :  3.3  !  2011-09  (M. Adani)  Original code: Drag Coefficient  
     7   !!         :  3.4  !  2012-10  (M. Adani)  Stokes Drift  
     8   !!            3.6  !  2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation 
     9   !!             -   !  2016-12  (G. Madec, E. Clementi) update Stoke drift computation 
     10   !!                                                    + add sbc_wave_ini routine 
     11   !!---------------------------------------------------------------------- 
     12 
     13   !!---------------------------------------------------------------------- 
     14   !!   sbc_stokes    : calculate 3D Stokes-drift velocities 
    1215   !!   sbc_wave      : wave data from wave model in netcdf files  
    13    !!---------------------------------------------------------------------- 
    14    USE oce            !  
    15    USE sbc_oce        ! Surface boundary condition: ocean fields 
    16    USE bdy_oce        ! 
    17    USE domvvl         ! 
     16   !!   sbc_wave_init : initialisation fo surface waves  
     17   !!---------------------------------------------------------------------- 
     18   USE phycst         ! physical constants  
     19   USE oce            ! ocean variables 
     20   USE sbc_oce        ! Surface boundary condition: ocean fields 
     21   USE zdf_oce,  ONLY : ln_zdfqiao 
     22   USE bdy_oce        ! open boundary condition variables 
     23   USE domvvl         ! domain: variable volume layers 
     24   ! 
    1825   USE iom            ! I/O manager library 
    1926   USE in_out_manager ! I/O manager 
    2027   USE lib_mpp        ! distribued memory computing library 
    21    USE fldread        ! read input fields 
     28   USE fldread        ! read input fields 
    2229   USE wrk_nemo       ! 
    23    USE phycst         ! physical constants  
    2430 
    2531   IMPLICIT NONE 
    2632   PRIVATE 
    2733 
    28    PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
    29    PUBLIC   sbc_wave    ! routine called in sbcmod 
     34   PUBLIC   sbc_stokes      ! routine called in sbccpl 
     35   PUBLIC   sbc_wave        ! routine called in sbcmod 
     36   PUBLIC   sbc_wave_init   ! routine called in sbcmod 
    3037    
    3138   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
    32    LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE. 
    33    LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE. 
    34    LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE. 
    35    LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE. 
    36    LOGICAL, PUBLIC     ::   cpl_wper=.FALSE. 
    37    LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE. 
    38    LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE. 
    39    LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE. 
    40  
    41    INTEGER ::   jpfld                ! number of files to read for stokes drift 
    42    INTEGER ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point 
    43    INTEGER ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point 
    44    INTEGER ::   jp_swh               ! index of significant wave hight      (m)      at T-point 
    45    INTEGER ::   jp_wmp               ! index of mean wave period            (s)      at T-point 
    46  
    47    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    48    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    49    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
    50    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
    51    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave  
    52    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum 
    53    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave 
    54    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d 
    55    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt 
    56    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d  
    57    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt 
     39   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE. 
     40   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE. 
     41   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE. 
     42   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE. 
     43   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE. 
     44   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
     45   LOGICAL, PUBLIC ::   cpl_wstrf  = .FALSE. 
     46   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE. 
     47 
     48   INTEGER ::   jpfld    ! number of files to read for stokes drift 
     49   INTEGER ::   jp_usd   ! index of stokes drift  (i-component) (m/s)    at T-point 
     50   INTEGER ::   jp_vsd   ! index of stokes drift  (j-component) (m/s)    at T-point 
     51   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point 
     52   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point 
     53 
     54   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_cd      ! structure of input fields (file informations, fields read) Drag Coefficient 
     55   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sd      ! structure of input fields (file informations, fields read) Stokes Drift 
     56   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_wn      ! structure of input fields (file informations, fields read) wave number for Qiao 
     57   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauoc  ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     58   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
     59   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:  
     60   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:   
     61   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !:  
     62   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence 
     63   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point 
     64   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp. 
    5865 
    5966   !! * Substitutions 
     
    7885      !! ** action   
    7986      !!--------------------------------------------------------------------- 
    80       INTEGER                ::   jj,ji,jk  
    81       REAL(wp)                       ::  ztransp, zfac, zsp0, zk, zus, zvs 
    82       REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
    83       !!--------------------------------------------------------------------- 
    84       ! 
    85  
    86       CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
    87       DO jk = 1, jpk 
    88          DO jj = 1, jpj 
    89             DO ji = 1, jpi 
    90                ! On T grid 
    91                ! Stokes transport speed estimated from Hs and Tmean 
    92                ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
     87      INTEGER  ::   jj, ji, jk   ! dummy loop argument 
     88      INTEGER  ::   ik           ! local integer  
     89      REAL(wp) ::  ztransp, zfac, ztemp, zsp0 
     90      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v 
     91      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace 
     92      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace 
     93      !!--------------------------------------------------------------------- 
     94      ! 
     95      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh ) 
     96      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
     97      ! 
     98      ! 
     99      zfac =  2.0_wp * rpi / 16.0_wp 
     100      DO jj = 1, jpj                ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) ) 
     101         DO ji = 1, jpi 
     102               ! Stokes drift velocity estimated from Hs and Tmean 
     103               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj) , 0.0000001_wp ) 
    93104               ! Stokes surface speed 
    94                zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2) 
     105               zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) 
     106               tsd2d(ji,jj) = zsp0 
    95107               ! Wavenumber scale 
    96                zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
    97                ! Depth attenuation 
    98                zfac = EXP(-2.0_wp*zk*gdept_n(ji,jj,jk))/(1.0_wp+8.0_wp*zk*gdept_n(ji,jj,jk)) 
     108               zk_t(ji,jj) = ABS( zsp0 ) / MAX( ABS( 5.97_wp*ztransp ) , 0.0000001_wp ) 
     109         END DO 
     110      END DO       
     111      DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
     112         DO ji = 1, jpim1 
     113            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
     114            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
     115            ! 
     116            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     117            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     118         END DO 
     119      END DO 
     120      ! 
     121      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
     122      DO jk = 1, jpkm1 
     123         DO jj = 2, jpjm1 
     124            DO ji = 2, jpim1 
     125               zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     126               zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     127               !                           
     128               zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     129               zkh_v = zk_v(ji,jj) * zdep_v 
     130               !                                ! Depth attenuation 
     131               zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 
     132               zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 
    99133               ! 
    100                usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
    101                vsd3dt(ji,jj,jk) = zfac * zvsd2dt(ji,jj) * tmask(ji,jj,jk) 
     134               usd(ji,jj,jk) = zda_u * zk_u(ji,jj) * umask(ji,jj,jk) 
     135               vsd(ji,jj,jk) = zda_v * zk_v(ji,jj) * vmask(ji,jj,jk) 
    102136            END DO 
    103137         END DO 
    104       END DO  
    105       ! Into the U and V Grid 
    106       DO jk = 1, jpkm1 
    107          DO jj = 1, jpjm1 
    108             DO ji = 1, fs_jpim1 
    109                usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   & 
    110                                &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) ) 
    111                vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   & 
    112                                &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) ) 
    113             END DO 
    114          END DO 
    115       END DO 
    116       ! 
    117       CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
    118       CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
    119       ! 
    120       DO jk = 1, jpkm1               ! Horizontal divergence 
     138      END DO    
     139      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 
     140      ! 
     141      !                       !==  vertical Stokes Drift 3D velocity  ==! 
     142      ! 
     143      DO jk = 1, jpkm1               ! Horizontal e3*divergence 
    121144         DO jj = 2, jpj 
    122145            DO ji = fs_2, jpi 
    123                ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     & 
    124                   &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     & 
    125                   &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     & 
    126                   &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj) 
     146               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji  ,jj,jk)    & 
     147                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    & 
     148                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    & 
     149                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj) 
    127150            END DO 
    128151         END DO 
     
    130153      ! 
    131154      IF( .NOT. AGRIF_Root() ) THEN 
    132          IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east 
    133          IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west 
    134          IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north 
    135          IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south 
    136       ENDIF 
    137       ! 
    138       CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
    139       ! 
    140       DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence 
    141          wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - e3t_n(:,:,jk) * ze3hdiv(:,:,jk) 
     155         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh(nlci-1,   :  ,:) = 0._wp      ! east 
     156         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh(  2   ,   :  ,:) = 0._wp      ! west 
     157         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh(  :   ,nlcj-1,:) = 0._wp      ! north 
     158         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh(  :   ,  2   ,:) = 0._wp      ! south 
     159      ENDIF 
     160      ! 
     161      CALL lbc_lnk( ze3divh, 'T', 1. ) 
     162      ! 
     163      IF( ln_linssh ) THEN   ;   ik = 1   ! none zero velocity through the sea surface 
     164      ELSE                   ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init) 
     165      ENDIF 
     166      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero) 
     167         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk) 
    142168      END DO 
    143169      ! 
    144170      IF( ln_bdy ) THEN 
    145171         DO jk = 1, jpkm1 
    146             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    147          END DO 
    148       ENDIF 
    149  
    150       CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 
     172            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:) 
     173         END DO 
     174      ENDIF 
     175      !                       !==  Horizontal divergence of barotropic Stokes transport  ==! 
     176      div_sd(:,:) = 0._wp 
     177      DO jk = 1, jpkm1                                 !  
     178        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk) 
     179      END DO 
     180      ! 
     181      CALL iom_put( "ustokes",  usd  ) 
     182      CALL iom_put( "vstokes",  vsd  ) 
     183      CALL iom_put( "wstokes",  wsd  ) 
     184      ! 
     185      CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh ) 
     186      CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
    151187      ! 
    152188   END SUBROUTINE sbc_stokes 
    153189 
    154    SUBROUTINE sbc_qiao 
    155       !!--------------------------------------------------------------------- 
    156       !!                     ***  ROUTINE sbc_qiao  *** 
    157       !! 
    158       !! ** Purpose :   Qiao formulation for wave enhanced turbulence 
    159       !!                2010 (DOI: 10.1007/s10236-010-0326)  
    160       !! 
    161       !! ** Method  : -  
    162       !! ** action   
    163       !!--------------------------------------------------------------------- 
    164       INTEGER :: jj, ji 
    165  
    166       ! Calculate the module of the stokes drift on T grid 
    167       !------------------------------------------------- 
    168       DO jj = 1, jpj 
    169          DO ji = 1, jpi 
    170             tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj) * zusd2dt(ji,jj) + zvsd2dt(ji,jj) * zvsd2dt(ji,jj) ) 
    171          END DO 
    172       END DO 
    173       ! 
    174    END SUBROUTINE sbc_qiao 
    175190 
    176191   SUBROUTINE sbc_wave( kt ) 
     
    188203      !! ** action   
    189204      !!--------------------------------------------------------------------- 
    190       USE zdf_oce,  ONLY : ln_zdfqiao 
    191  
    192       INTEGER, INTENT( in  ) :: kt       ! ocean time step 
    193       ! 
    194       INTEGER                ::   ierror   ! return error code 
    195       INTEGER                ::   ifpr 
    196       INTEGER                ::   ios      ! Local integer output status for namelist read 
    197       ! 
     205      INTEGER, INTENT(in   ) ::   kt   ! ocean time step 
     206      !!--------------------------------------------------------------------- 
     207      ! 
     208      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==! 
     209         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing 
     210         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
     211      ENDIF 
     212 
     213      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==! 
     214         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing 
     215         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     216      ENDIF 
     217 
     218      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==!  
     219         ! 
     220         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled 
     221            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing 
     222            IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
     223            IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     224            IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     225            IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     226         ENDIF 
     227         ! 
     228         ! Read also wave number if needed, so that it is available in coupling routines 
     229         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     230            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing 
     231            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     232         ENDIF 
     233            
     234         !                                         !==  Computation of the 3d Stokes Drift  ==!  
     235         ! 
     236         IF( jpfld == 4 )   CALL sbc_stokes()            ! Calculate only if required fields are read 
     237         !                                               ! In coupled wave model-NEMO case the call is done after coupling 
     238         ! 
     239      ENDIF 
     240      ! 
     241   END SUBROUTINE sbc_wave 
     242 
     243 
     244   SUBROUTINE sbc_wave_init 
     245      !!--------------------------------------------------------------------- 
     246      !!                     ***  ROUTINE sbc_wave_init  *** 
     247      !! 
     248      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
     249      !! 
     250      !! ** Method  : - Read namelist namsbc_wave 
     251      !!              - Read Cd_n10 fields in netcdf files  
     252      !!              - Read stokes drift 2d in netcdf files  
     253      !!              - Read wave number in netcdf files  
     254      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     255      !!                formulation 
     256      !! ** action   
     257      !!--------------------------------------------------------------------- 
     258      INTEGER ::   ierror, ios   ! local integer 
     259      INTEGER ::   ifpr 
     260      !! 
    198261      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    199262      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    200263      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    201                              &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
    202       !! 
    203       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc 
    204       !!--------------------------------------------------------------------- 
    205       ! 
    206       !                                         ! -------------------- ! 
    207       IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
    208          !                                      ! -------------------- ! 
    209          REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
    210          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    211 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
     264                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     265      ! 
     266      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 
     267      !!--------------------------------------------------------------------- 
     268      ! 
     269      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
     270      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
     271901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    212272          
    213          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    214          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    215 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
    216          IF(lwm) WRITE ( numond, namsbc_wave ) 
    217          ! 
    218          IF( ln_cdgw ) THEN 
    219             IF( .NOT. cpl_wdrag ) THEN 
    220                ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    221                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    222                ! 
    223                                       ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    224                IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    225                CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    226             ENDIF 
    227             ALLOCATE( cdn_wave(jpi,jpj) ) 
    228          ENDIF 
    229  
    230          IF( ln_tauoc ) THEN 
    231             IF( .NOT. cpl_wstrf ) THEN 
    232                ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
    233                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    234                ! 
    235                                        ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
    236                IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
    237                CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
    238             ENDIF 
    239             ALLOCATE( tauoc_wave(jpi,jpj) ) 
    240          ENDIF 
    241  
    242          IF( ln_sdw ) THEN 
    243             ! Find out how many fields have to be read from file if not coupled 
    244             jpfld=0 
    245             jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 
    246             IF( .NOT. cpl_sdrftx ) THEN 
    247                jpfld=jpfld+1 
    248                jp_usd=jpfld 
    249             ENDIF 
    250             IF( .NOT. cpl_sdrfty ) THEN 
    251                jpfld=jpfld+1 
    252                jp_vsd=jpfld 
    253             ENDIF 
    254             IF( .NOT. cpl_hsig ) THEN 
    255                jpfld=jpfld+1 
    256                jp_swh=jpfld 
    257             ENDIF 
    258             IF( .NOT. cpl_wper ) THEN 
    259                jpfld=jpfld+1 
    260                jp_wmp=jpfld 
    261             ENDIF 
    262  
    263             ! Read from file only the non-coupled fields  
    264             IF( jpfld > 0 ) THEN 
    265                ALLOCATE( slf_i(jpfld) ) 
    266                IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 
    267                IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 
    268                IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 
    269                IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 
    270                ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
    271                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    272                ! 
    273                DO ifpr= 1, jpfld 
    274                   ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    275                   IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    276                END DO 
    277  
    278                CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    279             ENDIF 
    280             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    281             ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) ) 
    282             ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 
    283             ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) ) 
    284             usd3d(:,:,:) = 0._wp 
    285             vsd3d(:,:,:) = 0._wp 
    286             wsd3d(:,:,:) = 0._wp 
    287             IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
    288                IF( .NOT. cpl_wnum ) THEN 
    289                   ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
    290                   IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
    291                                          ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
    292                   IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
    293                   CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
    294                ENDIF 
    295                ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
    296             ENDIF 
    297          ENDIF 
    298       ENDIF 
    299       ! 
    300       IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
    301          CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    302          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    303       ENDIF 
    304  
    305       IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
    306          CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
    307          tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
    308       ENDIF 
    309  
    310       IF( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    311          ! 
    312          ! Read from file only if the field is not coupled 
     273      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
     274      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
     275902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
     276      IF(lwm) WRITE ( numond, namsbc_wave ) 
     277      ! 
     278      IF( ln_cdgw ) THEN 
     279         IF( .NOT. cpl_wdrag ) THEN 
     280            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     281            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     282            ! 
     283                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
     284            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
     285            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 
     286         ENDIF 
     287         ALLOCATE( cdn_wave(jpi,jpj) ) 
     288      ENDIF 
     289 
     290      IF( ln_tauoc ) THEN 
     291         IF( .NOT. cpl_wstrf ) THEN 
     292            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     293            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     294            ! 
     295                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     296            IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     297            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 
     298         ENDIF 
     299         ALLOCATE( tauoc_wave(jpi,jpj) ) 
     300      ENDIF 
     301 
     302      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled 
     303         jpfld=0 
     304         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0 
     305         IF( .NOT. cpl_sdrftx ) THEN 
     306            jpfld  = jpfld + 1 
     307            jp_usd = jpfld 
     308         ENDIF 
     309         IF( .NOT. cpl_sdrfty ) THEN 
     310            jpfld  = jpfld + 1 
     311            jp_vsd = jpfld 
     312         ENDIF 
     313         IF( .NOT. cpl_hsig ) THEN 
     314            jpfld  = jpfld + 1 
     315            jp_hsw = jpfld 
     316         ENDIF 
     317         IF( .NOT. cpl_wper ) THEN 
     318            jpfld  = jpfld + 1 
     319            jp_wmp = jpfld 
     320         ENDIF 
     321 
     322         ! Read from file only the non-coupled fields  
    313323         IF( jpfld > 0 ) THEN 
    314             CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
    315             IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
    316             IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
    317             IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    318             IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
    319          ENDIF 
    320          ! 
    321          ! Read also wave number if needed, so that it is available in coupling routines 
    322          IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
    323             CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing 
    324             wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
    325          ENDIF 
    326             
    327          !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
    328          !(DOI: 10.1175/JPO-D-14-0020.1)==!  
    329          ! 
    330          ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 
    331          IF( jpfld == 4 ) THEN 
    332             CALL sbc_stokes() 
    333             IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
    334                CALL sbc_qiao() 
    335             ENDIF 
    336          ENDIF 
    337       ENDIF 
    338       ! 
    339    END SUBROUTINE sbc_wave 
    340        
     324            ALLOCATE( slf_i(jpfld) ) 
     325            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd 
     326            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd 
     327            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw 
     328            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
     329            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
     330            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     331            ! 
     332            DO ifpr= 1, jpfld 
     333               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     334               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     335            END DO 
     336            ! 
     337            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 
     338         ENDIF 
     339         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
     340         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     ) 
     341         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
     342         ALLOCATE( div_sd(jpi,jpj) ) 
     343         ALLOCATE( tsd2d (jpi,jpj) ) 
     344         usd(:,:,:) = 0._wp 
     345         vsd(:,:,:) = 0._wp 
     346         wsd(:,:,:) = 0._wp 
     347         ! Wave number needed only if ln_zdfqiao=T 
     348         IF( .NOT. cpl_wnum ) THEN 
     349            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     350            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' ) 
     351                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     352            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     353            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     354         ENDIF 
     355         ALLOCATE( wnum(jpi,jpj) ) 
     356      ENDIF 
     357      ! 
     358   END SUBROUTINE sbc_wave_init 
     359 
    341360   !!====================================================================== 
    342361END MODULE sbcwave 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r7421 r7562  
    110110      !                                         !==  effective transport  ==! 
    111111      IF( ln_wave .AND. ln_sdw )  THEN 
    112          DO jk = 1, jpkm1 
    113             zun(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) *      & 
    114                         &  ( un(:,:,jk) + usd3d(:,:,jk) )                       ! eulerian transport + Stokes Drift 
    115             zvn(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) *      & 
    116                         &  ( vn(:,:,jk) + vsd3d(:,:,jk) ) 
    117             zwn(:,:,jk) = e1e2t(:,:) *                    & 
    118                         &  ( wn(:,:,jk) + wsd3d(:,:,jk) ) 
     112         DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift 
     113            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) ) 
     114            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) ) 
     115            zwn(:,:,jk) = e1e2t(:,:)                 * ( wn(:,:,jk) + wsd(:,:,jk) ) 
    119116         END DO 
    120117      ELSE 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r6140 r7562  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
     26   USE sbcwave ,  ONLY: hsw   ! significant wave height 
     27   ! 
    2628   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2729   USE lib_mpp        ! MPP manager 
     
    194196         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    195197         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    196       ! 
     198      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
     199         zhsro(:,:) = hsw(:,:) 
    197200      END SELECT 
    198201 
     
    896899 
    897900      !                                !* Check of some namelist values 
    898       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    899       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    900       IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 
    901       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    902       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     901      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
     902      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
     903      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 
     904      IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 
     905      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
     906      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    903907 
    904908      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90

    r7403 r7562  
    7070         DO jj = 1, jpjm1 
    7171            DO ji = 1, fs_jpim1 
    72                qbv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *             & 
     72               qbv(ji,jj,jk) = 1.0 * 0.353553 * hsw(ji,jj) * tsd2d(ji,jj) *             & 
    7373            &                  EXP(3.0 * wnum(ji,jj) *                                  &                      
    7474            &                  (-MIN( gdepw_n(ji  ,jj  ,jk), gdepw_n(ji+1,jj  ,jk),     & 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7421 r7562  
    210210                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
    211211                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
    212       IF( ln_wave .AND. ln_sdw .AND. ln_stcor)           & 
    213                &         CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    214212                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
    215213                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r7421 r7562  
    4242   USE dynzdf           ! vertical diffusion               (dyn_zdf routine) 
    4343   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
    44    USE dynstcor         ! simp. form of Stokes-Coriolis 
    4544 
    4645   USE dynnxt           ! time-stepping                    (dyn_nxt routine) 
Note: See TracChangeset for help on using the changeset viewer.