
 MODULE wave

#ifdef key_wave

   !&E=====================================================================
   !&E
   !&E                      ***  MODULE  waveread  ***
   !&E ** Purpose : Read, spatial and time interpolation of waves datas
   !&E
   !&E
   !&E ** Called by:
   !&E       step (in file step.F90) 
   !&E
   !&E ** HOWTO: add New Waves Quantity (-nwq-)
   !&E       !    Search -nwq- in present file and add your quantity as others
   !&E       !    Open comwave.F90, search -nwq- and add your quantity as others
   !&E       !    Open outputsave.F90 search -nwq- and add your quantity as others
   !&E
   !&E ** History :
   !&E       !  2010-02 (G. Simon) Built on the basis of sflxread.F90
   !&E       Rq val OMP+MPI not done yet
   !&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E       !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Introduction of casrip test case
   !&E                    and review of the wave-current coupling
   !&E       !  2015-06  (V. Garnier, C. Penard) merge wave / comwave
   !&E
   !&E=====================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,riosh,rlg,    &
                          kmax,valmanq,     &
                          imin,imax,        &
                          jmin,jmax,        &
                          limin,limax,      &
                          ljmin,ljmax,      &
                          limaxp1,ljmaxp1,  &
                          ljminp1,liminp1,  &
                          limaxp2,ljmaxp2,  &
                          liminm2,ljminm2,  &
                          liminm1,ljminm1,  &
                          pi,epsilon,lchain

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC wave_read             ! function called by step
!#ifdef key_wave_qeuler_dynamics
   PUBLIC wave_alloc,wave_dealloc      ! function called by step
!#endif
#ifdef key_wave_readww3txt
   PUBLIC wave_read_ww3txt             ! function called by step
#endif

   !! * Shared module variables

   !! * Public variables
   INTEGER,PUBLIC               :: idimt ! Number of sample in the wave file
   INTEGER,PUBLIC               :: iwavemin,iwavemax,jwavemin,jwavemax ! Extrem spatial size of wave array
   INTEGER,PUBLIC               :: ilec_wave ! For manage time sample
   INTEGER,PUBLIC               :: ilecmem_wave ! For manage time sample
   INTEGER,PUBLIC               :: iwmin0  ! shift to read the correct index in the netcdf file

   CHARACTER(LEN=lchain),PUBLIC :: file_wave ! Location of the NetCDF file
   CHARACTER(LEN=lchain),PUBLIC :: longitude_name_wave,latitude_name_wave ! name of the longitude / latitude variable in netCDF file
   REAL(KIND=rsh),PUBLIC        :: valfill_wave ! Fill value (ground, the same for each variable of waves file)
   INTEGER,PUBLIC               :: cdf_id_wave
   REAL(KIND=rlg),PUBLIC        :: tnc1_wave,tnc2_wave

   LOGICAL,PUBLIC                      :: l_swan=.false.

   !! For locate MARS grid from WAVE grid
   !!   Related to center of the MARS cell
   INTEGER,PUBLIC,       DIMENSION(:,:),ALLOCATABLE :: i1xe_wave,j1xe_wave
   REAL(KIND=rsh),PUBLIC,DIMENSION(:,:),ALLOCATABLE :: xdpxe_wave,ydpxe_wave
   !!   Related to zonal velocity (East face)
   INTEGER,PUBLIC,       DIMENSION(:,:),ALLOCATABLE :: i1u_wave,j1u_wave
   REAL(KIND=rsh),PUBLIC,DIMENSION(:,:),ALLOCATABLE :: xdpu_wave,ydpu_wave
   !!   Related to meridional velocity (North face)
   INTEGER,PUBLIC,       DIMENSION(:,:),ALLOCATABLE :: i1v_wave,j1v_wave
   REAL(KIND=rsh),PUBLIC,DIMENSION(:,:),ALLOCATABLE :: xdpv_wave,ydpv_wave

   !! * Waves quantities -neww-
   !!   - Significant height
   LOGICAL,PUBLIC               :: l_wave_rhs
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_hs
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_hs1,wave_hs2
   !!   - Bottom orbital velocity
   !!  - rms amplitude of orbital bottom velocity of the waves (uba in m/s) ! Bottom rms velocity
   LOGICAL,PUBLIC               :: l_wave_rubot
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_ubot
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_ubot1,wave_ubot2
   LOGICAL,PUBLIC               :: l_wave_ruubot
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_uubot
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_uubot1,wave_uubot2
   LOGICAL,PUBLIC               :: l_wave_rvubot
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_vubot
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_vubot1,wave_vubot2
   !! Radiation stress
   LOGICAL,PUBLIC               :: l_wave_rradstr
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_sx
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_sx1,wave_sx2
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_sy
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_sy1,wave_sy2
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_sxy
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_sxy1,wave_sxy2
   !!  - Wave-to-ocean TKE flux (phioc in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
   LOGICAL,PUBLIC               :: l_wave_rfoc
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_foc
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_foc1,wave_foc2
   !!  - Wave-to-ocean TKE flux due to wave bottom friction (phibbl in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
   LOGICAL,PUBLIC               :: l_wave_rfbb
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_fbb
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_fbb1,wave_fbb2
   !!  - Wave-to-ocean momentum due to wave bottom friction (taubbl in m2 s-2)
   LOGICAL,PUBLIC               :: l_wave_rtbb
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tbbx,name_wave_tbby
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tbbx1,wave_tbbx2
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tbby1,wave_tbby2
   !!  - Mean wave period (tmn in s) (m0,-1)
   LOGICAL,PUBLIC               :: l_wave_rt0m1
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_t0m1
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_t0m11,wave_t0m12
   !!  - Directionality
   LOGICAL,PUBLIC               :: l_wave_rdir
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_dir
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_dir1,wave_dir2
   !!  - Wave peak period
   LOGICAL,PUBLIC               :: l_wave_rtp
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tp
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tp1,wave_tp2
   !!  - Net wave-supported stress, x-component (tauwix in m2.s-2)
   LOGICAL,PUBLIC               :: l_wave_rtauwx
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tauwx
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tauwx1,wave_tauwx2
   !!  - Net wave-supported stress, y-component (tauwiy in m2.s-2)
   LOGICAL,PUBLIC               :: l_wave_rtauwy
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tauwy
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tauwy1,wave_tauwy2
   !!  - Wave-ocean momentum flux, x-component (tauox in m2.s-2)
   LOGICAL,PUBLIC               :: l_wave_rtauox
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tauox
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tauox1,wave_tauox2
   !!  - Wave-ocean momentum flux, y-component (tauoy in m2.s-2)
   LOGICAL,PUBLIC               :: l_wave_rtauoy
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tauoy
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tauoy1,wave_tauoy2
   !!  - eastward surface stokes drift (uuss in m2 s-1)
   LOGICAL,PUBLIC               :: l_wave_russ
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_ussx,name_wave_ussy
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_ussx1,wave_ussx2
   !!  - northward surface stokes drift (vuss in m2 s-1)
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_ussy1,wave_ussy2
   !!  - eastward stokes transport (utus in m2 s-1)
   LOGICAL,PUBLIC               :: l_wave_rtus
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_tusx,name_wave_tusy
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tusx1,wave_tusx2
   !!  - northward stokes transport (vtus in m2 s-1)
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_tusy1,wave_tusy2
   !!  - Mean wave length (wlm in m)
   LOGICAL,PUBLIC               :: l_wave_rwlm
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_wlm
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_wlm1,wave_wlm2
   !!  - pression induite par les vagues
   LOGICAL,PUBLIC               :: l_wave_rbhd
   CHARACTER(LEN=lchain),PUBLIC :: name_wave_bhd
   REAL(KIND=rsh),DIMENSION(:,:),ALLOCATABLE :: wave_bhd1,wave_bhd2

   !!  - If data have been interpolated on the MARS grid in advance
   LOGICAL,PUBLIC               :: l_wave_marsgrid

   ! Saving variables
   CHARACTER(len=lchain), PUBLIC :: name_out_wave_dir,name_out_wave_hs,name_out_wave_ubot, &
                                    name_out_wave_ubotx,name_out_wave_uboty,               &
                                    name_out_wave_foc,name_out_wave_tp,name_out_wave_t0m1, &
                                    name_out_wave_tauox,name_out_wave_tauoy,name_out_wave_wlm, &
                                    name_out_wave_tauwx,name_out_wave_tauwy,              &
                                    name_out_wave_bhd,name_out_wave_ussx,name_out_wave_ussy, &
                                    name_out_wave_tusx,name_out_wave_tusy,                &
                                    name_out_wave_usx,name_out_wave_usy,name_out_wave_usz,&
                                    name_out_wave_fbb,name_out_wave_tbbx,name_out_wave_tbby,&
                                    name_out_wave_sx,name_out_wave_sy

   REAL(KIND=riosh), PUBLIC      :: riog_valid_min_wave_dir,riog_valid_max_wave_dir,      &
                                    riog_valid_min_wave_ubot,riog_valid_max_wave_ubot,    &
                                    riog_valid_min_wave_hs,riog_valid_max_wave_hs,        &
                                    riog_valid_min_wave_foc,riog_valid_max_wave_foc,      &
                                    riog_valid_min_wave_tp,riog_valid_max_wave_tp,        &
                                    riog_valid_min_wave_t0m,riog_valid_max_wave_t0m,      &
                                    riog_valid_min_wave_two,riog_valid_max_wave_two,      &
                                    riog_valid_min_wave_wlm,riog_valid_max_wave_wlm,      &
                                    riog_valid_min_wave_bhd,riog_valid_max_wave_bhd,      &
                                    riog_valid_min_wave_stokes,riog_valid_max_wave_stokes

   REAL(KIND=rsh), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::  &
                   wave_hs,    & ! significant wave height (m)
                   wave_tp,    & ! Wave peak period (s)
                   wave_uubot, & ! x-component of rms amplitude of wave orbital bottom velocity (uba in m/s)
                   wave_vubot, & ! y-component of rms amplitude of wave orbital bottom velocity (uba in m/s)
                   wave_sx,    & ! x-component of radiation stress
                   wave_sy,    & ! y-component of radiation stress 
                   wave_sxy    ! xy-component of radiation stress 


#if defined key_wave_qeuler_dynamics || defined key_wave

   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE :: mdptc  ! water depth at T location
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hrms
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: wave_usz,wave_vsz,wave_wsz
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: duzdsig
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: dvzdsig
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: dsigdx_u,dsigdy_v
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: wave_botdiss_u,wave_botdiss_v
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: taubbl
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:,:),ALLOCATABLE  :: wave_tke,peps
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE  :: wave_freq
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE  :: tausurf
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE :: usz_zint,vsz_zint
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE :: wave_number
   REAL(KIND=rsh),PUBLIC, DIMENSION(:,:), ALLOCATABLE :: wave_km,wave_km_u,wave_km_v
   INTEGER,       PUBLIC, DIMENSION(:,:), ALLOCATABLE :: wave_km_iter
   REAL(KIND=rsh),PUBLIC, PARAMETER :: c1eps=1.44_rsh
   REAL(KIND=rsh),PUBLIC, PARAMETER :: cflin=0.05_rsh

   REAL(KIND=rsh), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::  &
                   wave_t0m1, &   ! mean period of the waves
                   wave_dir, &    ! (wave_thm) mean direction of the waves
                   wave_bhd, &    ! pression induite par les vagues
                   wave_tauox, &  ! tension liee au deferlement, composante en x
                   wave_tauoy, &  !tension liee au deferlement, composante en y
                   wave_tauwx, &  ! tension de vent supportee par les vagues, composante en x
                   wave_tauwy, &  ! tension de vent supportee par les vagues, composante en y
                   wave_ubot,  &  ! amplitude of orbital velocity of the waves
                   wave_phioc, &  ! flux d energie associee au deferlement (in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
                   wave_phibbl, &  ! tke flux associated to dissipation in the bottom boundary layer (positive)
                                   ! (in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
                   wave_taubbl_x, &  ! stress associated to dissipation in the bottom boundary layer (same direction of wave propagation)
                   wave_taubbl_y, &  ! stress associated to dissipation in the bottom boundary layer (same direction of wave propagation)
                   wave_wlm,  &   ! mean length wave
                   wave_ussx, &   ! eastward surface stokes drift
                   wave_ussy, &   ! northward surface stokes drift
                   wave_tusx, &   ! eastward stokes transport
                   wave_tusy      ! northward stokes transport

   ! longueur d onde resolue par WW3 (2015) 10 cm a 30 km
   REAL(KIND=rsh),PUBLIC,PARAMETER :: wave_wlm_min=0.1_rsh, wave_wlm_max=30000.0_rsh
   ! période des vaguesresolue par WW3 (2015) 0.2 a 400 s
!cval REVOIR
   REAL(KIND=rsh),PUBLIC,PARAMETER :: wave_t0m1_min=0.095_rsh, wave_t0m1_max=400.0_rsh
#endif

!#ifdef key_wave_soulsby_bottom_stress
   REAL(KIND=rsh), PUBLIC, DIMENSION(:,:), ALLOCATABLE  :: wave_uzk1, & ! quasi-eulerian velocity at bottom layer
                                                           wave_vzk1
   REAL(KIND=rsh), PUBLIC, DIMENSION(:,:), ALLOCATABLE  :: uorb   ! bottom wave orbital velocity
   REAL(KIND=rsh), PUBLIC, DIMENSION(:,:), ALLOCATABLE  :: fwst   ! wave friction factor
!#endif

   !! * Private variables
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: lonwave2d,latwave2d
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: lonwave2dtmp,latwave2dtmp
#ifdef key_sflx_interp_scrip
   LOGICAL,DIMENSION(:,:),ALLOCATABLE        :: maskdst,masksrc
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: matrix,matrix_u,matrix_v
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: tabvar2d,tabinterp2d
   INTEGER,DIMENSION(:),ALLOCATABLE          :: src_add,dst_add,src_add_u,dst_add_u,src_add_v,dst_add_v
   INTEGER :: nxfin,nyfin
#endif

 CONTAINS

   !!======================================================================

  SUBROUTINE wave_read

   !&E---------------------------------------------------------------------
   !&E                         ***  ROUTINE wave_read  ***
   !&E
   !&E ** Purpose : Read waves file data generated from differents waves models.
   !&E              These data are provided in a NetCDF file with an AO structure :
   !&E              the grid is oriented from South to North and from West to East
   !&E              and spatial steps are constant.
   !&E
   !&E ** Description : 
   !&E              Read the following variables :
   !&E                Mean direction of the total sea (dir=ATAN2(b1tot/a1tot)) (in [degrees])
   !&E                Significant wave height (in [m])
   !&E                Energy flux from waves to current in in ([W.m^{-2}])
   !&E
   !&E ** Map : 
   !&E             -onr- Open Netcdf, Read lat and lonfile (init: only do for nbouc=1)
   !&E             -faa- Find size of 2D Array to Allocate (init)
   !&E             -lmw- Locate Mars cell from Wave cell (init)
   !&E             -ctd- Checkout Time Dependences of inputed waves data (init)
   !&E             -ri2- Read and spatial Interpolate of 2 time step of input file (init)
   !&E             -cou- CheckOut Units (init)
   !&E             -fse- Find time Step in the file that Exceed computational time step
   !&E             -2i1- We transfert array 2 In array 1
   !&E             -rin- Read and spatial Interpolate Next time step of input file while computational time step is not exceed
   !&E             -tli- Do time Linear Interpolation
   !&E             (-mmg- Mpi Managing for Gradient computation)
   !&E
   !&E ** Called by : 
   !&E                 step
   !&E
   !&E ** External calls : tool_interp2d_sflx,ionc4_openr,ionc4_corres,ionc4_flimax,
   !&E                     ionc4_fljmax,ionc4_read_lat,ionc4_read_lon, ionc4_read_dimt,
   !&E                     ionc4_read_time,ionc4_read_subxyt
   !&E
   !&E ** Reference :
   !&E      SWAN: TODO
   !&E      WAVEWATCH III:
   !&E        TOLMAN Hendrik L., 2002, User manual and system documentation of WAVEWATCH-III version 2.22
   !&E        TOLMAN Hendrik L., 2007, The 2007 release of WAVEWATCH III
   !&E
   !&E ** History :
   !&E       !  2010-02 (G. Simon) Build on the basis of module sflx_read
   !&E       !  2010-07 (R. Verney) Introduction of orbital velocity and wave period variables
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used: (definitions of variable in module)
   USE parameters, ONLY : tref

   USE comvars2d, ONLY  : t,h0,                              &
                          nbouc,                             &
                          iscreenlog,ierrorlog,iwarnlog,     &
                          lon2d,lat2d,lon2d_g,lat2d_g,       &
                          lon2du,lon2dv,lat2du,lat2dv,       &
                          lonmin,latmin,lonmax,latmax
   USE ionc4, ONLY : ionc4_openr,ionc4_corres,ionc4_flimax,ionc4_fljmax, &
                     ionc4_read_lat,ionc4_read_lon,ionc4_read_dimt,ionc4_read_time, &
                     ionc4_read_subxyt,ionc4_read_xy,ionc4_read_var_dim
#if defined key_grid_rotated || defined key_sflx_interp_scrip
   USE comvars2d, ONLY  : grid_angle
#endif

   USE_MPI toolmpi,   ONLY : ex_i_rsh,ex_j_rsh,BCAST_MPI,MASTER
#ifdef key_sflx_interp_scrip
   USE bicubic_interp,ONLY : get_remap_bicub,make_bicubic_remap
#endif

   IMPLICIT NONE
   !! * Local declarations
   LOGICAL           :: l_pb,IERR_MPI
   CHARACTER(LEN=19) :: tool_sectodat,tdate
   REAL(KIND=rlg)    :: tool_datosec ! Function for convert time date->second
   INTEGER           :: ncerror ! Handle IONC-NetCDF error
   INTEGER           :: dimlon ! dimension of the longitude, latitude variables (1D or 2D)
   INTEGER           :: iwavemin_ref,iwavemax_ref,jwavemin_ref,jwavemax_ref
   INTEGER           :: ixwave,jxwave ! size of the lon and lat array
   INTEGER           :: i,j ! Loop indice
   REAL(KIND=rsh)    :: dgwave,dfiwave ! differentail element of longitude and latitude [degret]
   REAL(KIND=rlg),  ALLOCATABLE,DIMENSION(:)   :: lonwave,latwave ! longitude and latitude in wave file [degrees]
   REAL(KIND=riosh),ALLOCATABLE,DIMENSION(:,:) :: tab1,tab2 ! array for process wave quantity
   REAL(KIND=rsh)    :: tabmin,tabmax ! size maximum of process array
   REAL(KIND=rlg)    :: torigin ! Used for comute offset between Mars and wave time [s]
   REAL(KIND=rlg)    :: twave ! The same as t [s]
   REAL(KIND=rlg)    :: dt1,dt2 ! Used for linear interpolation in time [s]
   REAL(KIND=rsh)    :: radeg ! Identity for convert radien angle in  degret
   REAL(KIND=rsh)    :: valfillm ! Left limit to valfill_wave
   REAL(KIND=rsh)    :: valfillp ! Right limit to valfill_wave

   !! For locate MARS cell from Wave cell
   REAL(KIND=rsh)    :: xm,ym

   LOGICAL           :: dbug ! dbug flag
  
   REAL(KIND=rlg)   :: t_lag

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter SFLX_WAVE routine, nbouc=',nbouc

!!!!!!$OMP PARALLEL DEFAULT(none)                                          &
!!!!!!$OMP& PRIVATE (                                                      &
!!!!!!$OMP& IERR_MPI,l_pb,tdate,dimlon,                                    &
!!!!!!$OMP& iwavemin_ref,iwavemax_ref,jwavemin_ref,jwavemax_ref,               &
!!!!!!$OMP&ncerror,                                                        &
!!!!!!$OMP&ixwave,jxwave,                                                  &
!!!!!!$OMP&i,j,                                                            &
!!!!!!$OMP&dgwave,dfiwave,                                                 &
!!!!!!$OMP&lonwave,latwave,                                                &
!!!!!!$OMP&tab1,tab2,                                                      &
!!!!!!$OMP&tabmin,tabmax,                                                  &
!!!!!!$OMP&torigin,                                                        &
!!!!!!$OMP&twave,                                                          &
!!!!!!$OMP&dt1,dt2,                                                        &
!!!!!!$OMP&radeg,                                                          &
!!!!!!$OMP& xm,ym,                                                         &
!!!!!!$OMP& dbug,                                                          &
!!!!!!$OMP& valfillm,                                                      &
!!!!!!$OMP& valfillp                                                       &
!!!!!!$OMP& )                                                              &
!!!!!!$OMP& SHARED(                                                        &
#ifdef key_MPI_2D
!!!!!!$OMP& MASTER,                                                        &
#endif
!!!!!!$OMP& lon2d,lat2d,lon2du,lat2du,lon2dv,lat2dv,lon2d_g,lat2d_g,       &
!!!!!!$OMP& lonmin,latmin,lonmax,latmax,                                   &
#if defined key_grid_rotated || defined key_sflx_interp_scrip
!!!!!!$OMP& grid_angle,                                                    &
#endif
!!!!!!$OMP& lonwave2d,latwave2d,lonwave2dtmp,latwave2dtmp,                 &
#ifdef key_sflx_interp_scrip
!!!!!!$OMP& maskdst,masksrc,matrix,matrix_u,matrix_v,tabvar2d,tabinterp2d, &
#endif
#ifdef key_sflx_interp_scrip
!!!!!!$OMP& src_add,dst_add,src_add_u,dst_add_u,src_add_v,dst_add_v,       &
!!!!!!$OMP& nxfin,nyfin,                                                   &
#endif
!!!!!!$OMP& file_wave,valfill_wave,cdf_id_wave,tnc1_wave,tnc2_wave,        & 
!!!!!!$OMP& longitude_name_wave,latitude_name_wave,                        &
!!!!!!$OMP& idimt,h0,                                                      &
!!!!!!$OMP& ilecmem_wave,ilec_wave,iwmin0,                                 &
!!!!!!$OMP& iwavemin,iwavemax,jwavemin,jwavemax,l_wave_marsgrid,           &
!!!!!!$OMP& l_wave_rdir,wave_dir,name_wave_dir,wave_dir1,wave_dir2,        &
!!!!!!$OMP& l_wave_rfoc,wave_phioc,name_wave_foc,wave_foc1,wave_foc2,      &
!!!!!!$OMP& l_wave_rfbb,wave_phibbl,name_wave_fbb,wave_fbb1,wave_fbb2,     &
!!!!!!$OMP& l_wave_rtbb,wave_taubbl_x,wave_taubbl_y,name_wave_tbbx,        &
!!!!!!$OMP& name_wave_tbby,wave_tbbx1,wave_tbbx2,wave_tbby1,wave_tbby2,    &
!!!!!!$OMP& l_wave_rhs,wave_hs,name_wave_hs,wave_hs1,wave_hs2,             &
!!!!!!$OMP& l_wave_rubot,wave_ubot,name_wave_ubot,wave_ubot1,wave_ubot2,   &
!!!!!!$OMP& l_wave_ruubot,wave_uubot,name_wave_uubot,wave_uubot1,wave_uubot2,   &
!!!!!!$OMP& l_wave_rvubot,wave_vubot,name_wave_vubot,wave_vubot1,wave_vubot2,   &
!!!!!!$OMP& l_wave_rtp,wave_tp,name_wave_tp,wave_tp1,wave_tp2,             &
!!!!!!$OMP& i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,                     &
!!!!!!$OMP& i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,                         &
!!!!!!$OMP& i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,                         &
!!!!!!$OMP& t,nbouc,limin,limax,liminp1,limaxp1,ljmin,ljmax,               &
!!!!!!$OMP& ljminp1,limaxp2,ljmaxp2,iscreenlog,ierrorlog,iwarnlog,         &
!!!!!!$OMP& l_wave_rtauwx,l_wave_rtauwy,l_wave_rtauox,l_wave_rtauoy,l_wave_rwlm, &
!!!!!!$OMP& wave_tauwx1,wave_tauwx2,wave_tauwy1,wave_tauwy2,               &
!!!!!!$OMP& wave_tauox1,wave_tauox2,wave_tauoy1,wave_tauoy2,               &
!!!!!!$OMP& wave_tauwx,wave_tauwy,wave_tauox,wave_tauoy,wave_wlm,          &
!!!!!!$OMP& wave_wlm1,wave_wlm2,name_wave_tauwx,name_wave_tauwy,           &
!!!!!!$OMP& name_wave_tauox,name_wave_tauoy,name_wave_wlm,                 &
!!!!!!$OMP& l_wave_rt0m1,name_wave_t0m1,l_wave_rbhd,name_wave_bhd,         &
!!!!!!$OMP& wave_t0m1,wave_bhd,wave_bhd1,wave_bhd2,wave_t0m11,wave_t0m12,  &
!!!!!!$OMP& wave_tusx,wave_tusx1,wave_tusx2,l_wave_russ,l_wave_rtus,       &
!!!!!!$OMP& wave_tusy,wave_tusy1,wave_tusy2,l_wave_rradstr,                &
!!!!!!$OMP& name_wave_sx,name_wave_sy,name_wave_sxy,                       &
!!!!!!$OMP& name_wave_ussx,name_wave_ussy,name_wave_tusx,name_wave_tusy    &
!!!!!!$OMP )

   dbug = .true.
   l_pb = .false. ! No problem

   radeg = pi/180.0_rsh

   twave=t
   torigin=tref
   torigin=0.0_rlg
   t_lag=tool_datosec('01/01/1990 00:00:00')-tool_datosec('01/01/1900 00:00:00')
   ! include time lag if time origine is not 01/01/1900 in wave file
!#ifdef key_castest
!   torigin=tool_datosec('01/01/1900 00:00:00')-tool_datosec('01/01/1990 00:00:00')
!#endif

   IF (dbug) print *,'ENTER WAVEREAD'

   IF (nbouc==1) THEN ! START if (nbouc==1) then -init-
   !!----------------------------------------------------------------------
   !! * Open NetCDF file, read lat and lon -onr-
   !!----------------------------------------------------------------------

      iwmin0=1

!!!!!!$OMP single
     CALL ionc4_openr(file_wave,l_out_nc4par=.false.)
     CALL ionc4_corres(file_wave,cdf_id_wave)
     IF (dbug) print *,'cfd_id',cdf_id_wave
     ixwave=ionc4_flimax(cdf_id_wave)
     IF (dbug) print *,'ixwave',ixwave
     jxwave=ionc4_fljmax(cdf_id_wave)  
     IF (dbug) print *,'jxwave',jxwave


     ! Get the dimensions of the longitude or latitude variables
     CALL ionc4_read_var_dim(file_wave,TRIM(longitude_name_wave),dimlon)

     ALLOCATE(latwave2dtmp(1:ixwave,1:jxwave),lonwave2dtmp(1:ixwave,1:jxwave))
     IF (dimlon == 2) THEN
       CALL ionc4_read_xy(file_wave,TRIM(latitude_name_wave),latwave2dtmp(1:ixwave,1:jxwave),1,ixwave,1,jxwave)
       CALL ionc4_read_xy(file_wave,TRIM(longitude_name_wave),lonwave2dtmp(1:ixwave,1:jxwave),1,ixwave,1,jxwave)
     ELSE
       ALLOCATE(latwave(1:jxwave),lonwave(1:ixwave))
       CALL ionc4_read_lat(file_wave,latwave,1,jxwave)
       CALL ionc4_read_lon(file_wave,lonwave,1,ixwave)
       IF (dbug) print *,'latwave',latwave
       IF (dbug) print *,'lonwave',lonwave
       DO j=1,jxwave
         lonwave2dtmp(1:ixwave,j)=lonwave(1:ixwave)
       ENDDO
       DO i=1,ixwave
         latwave2dtmp(i,1:jxwave)=latwave(1:jxwave)
       ENDDO
     END IF

     dgwave =ABS(lonwave2dtmp(2,1)-lonwave2dtmp(1,1))
     dfiwave=ABS(latwave2dtmp(1,2)-latwave2dtmp(1,1))


   !!----------------------------------------------------------------------
   !! * Find size of 2D array to allocate -faa-
   !!----------------------------------------------------------------------
#ifdef key_grid_rotated
     !toseelater
     ! to estimate min and max indexes relative to the rotated grid
     ! to be used in atmospherical file
     ! Work indifferently if wave. model longitudes decrease or increase.
     ! dimensions of the wave file usefull for each cpu
     iwavemin=1               ! global dimensions of the wave file by default
     iwavemax=ixwave
     jwavemin=1
     jwavemax=jxwave
     DO i=1,ixwave-1
       IF ((lonmin-lonwave(i))*(lonmin-lonwave(i+1)) <= 0.0_rlg) iwavemin=i
       IF ((lonmax-lonwave(i)+epsilon)*(lonmax-lonwave(i+1)-epsilon) <= 0.0_rlg) iwavemax=i+1
     END DO
     ! Work indifferently if wave. model latitudes decrease or increase.
     DO j=1,jxwave-1
       IF ((latmin-latwave(j))*(latmin-latwave(j+1)) <= 0.0_rlg) jwavemin=j
       IF ((latmax-latwave(j)+epsilon)*(latmax-latwave(j+1)-epsilon) <= 0.0_rlg) jwavemax=j+1
     END DO

#else
     IF (.NOT. l_wave_marsgrid) THEN
     ! dimensions of the wave file usefull for each cpu
     iwavemin=ixwave          ! global dimensions of the wave file by default
     iwavemax=1
     jwavemin=jxwave
     jwavemax=1
     iwavemin_ref=iwavemin
     iwavemax_ref=iwavemax
     jwavemin_ref=jwavemin
     jwavemax_ref=jwavemax
     ! check the useful indexes of the wave file
     ! Works indifferently if wave. model longitudes decrease or increase.
     DO j=1,jxwave
     DO i=1,ixwave-1
       IF((lonmin-lonwave2dtmp(i,j)) * (lonmin-lonwave2dtmp(i+1,j)) <= 0.0_rlg) THEN
         IF (i<iwavemin .AND. MAXVAL(latwave2dtmp(i:ixwave,j))>latmin-dfiwave .AND. &
                              MINVAL(latwave2dtmp(i:ixwave,j))<latmax+dfiwave) THEN
           iwavemin=i
           IF (j<jwavemin) jwavemin=j
         END IF
       END IF
     ENDDO
     ENDDO
     DO i=iwavemin,ixwave-1
     DO j=jwavemin,jxwave-1
       IF((latmax-latwave2dtmp(i,j)  +epsilon) * (latmax-latwave2dtmp(i,j+1)-epsilon) <= 0.0_rlg) THEN
         IF (j+1>=jwavemax  .AND. MAXVAL(lonwave2dtmp(i+1,jwavemin:j+1))>lonmin-dgwave .AND. &
                                  MINVAL(lonwave2dtmp(i+1,jwavemin:j+1))<lonmax+dgwave ) THEN
           jwavemax=j+1
           IF (i+1>iwavemax) iwavemax=i+1
         END IF
       END IF
     ENDDO
     ENDDO

#endif

     IF (iwavemin==iwavemin_ref .AND. iwavemax==iwavemax_ref) THEN
       WRITE(iwarnlog,*) 'None index has been found'
       WRITE(iwarnlog,*) 'iwavemin =',iwavemin,' iwavemax =',iwavemax
       WRITE(iwarnlog,*) 'jwavemin =',jwavemin,' jwavemax =',jwavemax
       WRITE(iwarnlog,*) 'You use too many cpus for a so coarse wave grid'
       WRITE(iwarnlog,*) 'or the file is not large enough'
       WRITE(iwarnlog,*) 'iwavemin is set to 1 and iwavemax to ixwave'
!       WRITE(ierrorlog,*) 'Simulation stopped.'
!       CALL_MPI MPI_FINALIZE(IERR_MPI)
!       STOP
       iwavemin=1
       iwavemax=ixwave
     END IF

     IF (jwavemin==jwavemin_ref .AND. jwavemax==jwavemax_ref) THEN
       WRITE(iwarnlog,*) 'None index has been found'
       WRITE(iwarnlog,*) 'iwavemin =',iwavemin,' iwavemax =',iwavemax
       WRITE(iwarnlog,*) 'jwavemin =',jwavemin,' jwavemax =',jwavemax
       WRITE(iwarnlog,*) 'You use too many cpus for a so coarse wave grid'
       WRITE(iwarnlog,*) 'or the file is not large enough'
       WRITE(iwarnlog,*) 'jwavemin is set to 1 and jwavemax to jxwave'
!       WRITE(ierrorlog,*) 'Simulation stopped.'
!       CALL_MPI MPI_FINALIZE(IERR_MPI)
!       STOP
       jwavemin=1
       jwavemax=jxwave
     END IF

     IF (iwavemin>=iwavemax .OR. jwavemin>=jwavemax) THEN
       WRITE(ierrorlog,*) 'The size of the domain read from the wave field is not large enough'
       WRITE(ierrorlog,*) 'iwavemin =',iwavemin,' iwavemax =',iwavemax
       WRITE(ierrorlog,*) 'jwavemin =',jwavemin,' jwavemax =',jwavemax
       WRITE(ierrorlog,*) 'You use too many cpus for a so coarse wave grid'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF

     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'Waves file : spatial limits'
     WRITE(iscreenlog,*) 'iwavemin = ',iwavemin
     WRITE(iscreenlog,*) 'iwavemax = ',iwavemax
     WRITE(iscreenlog,*) 'jwavemin = ',jwavemin
     WRITE(iscreenlog,*) 'jwavemax = ',jwavemax
     WRITE(iscreenlog,*) ' '
     ENDIF_MPI

#ifdef key_sflx_interp_scrip
     IF (iwavemin <= 2) THEN
       WRITE(ierrorlog,*) 'Wave file too small'
       WRITE(ierrorlog,*) 'At the Eastern boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
     IF (jwavemin <= 2) THEN
       WRITE(ierrorlog,*) 'Wave file too small'
       WRITE(ierrorlog,*) 'At the Southern boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
     IF (iwavemax >= ixwave - 1) THEN
       WRITE(ierrorlog,*) 'Wave file too small'
       WRITE(ierrorlog,*) 'At the Western boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
     IF (jwavemax >= jxwave - 1) THEN
       WRITE(ierrorlog,*) 'Wave file too small'
       WRITE(ierrorlog,*) 'At the Northern boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
#endif

     iwavemin = MAX0(1 , iwavemin - 2)
     jwavemin = MAX0(1 , jwavemin - 2)
     iwavemax = MIN0(ixwave, iwavemax + 2)
     jwavemax = MIN0(jxwave, jwavemax + 2)
     ENDIF
     IF (l_wave_marsgrid) THEN  ! the grid of the wave file is equal to the MARS grid
                            ! but on 1:imax+1,1:jmax+1
       iwmin0=0
       iwavemin=limin
       iwavemax=limax
       jwavemin=ljmin
       jwavemax=ljmax
!       iwmin0=1
!       iwavemin=imin+1
 !      iwavemax=imax+1
  !     jwavemin=jmin+1
   !    jwavemax=jmax+1
     END IF

     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'Wave file : spatial limits'
     WRITE(iscreenlog,*) 'iwavemin = ',iwavemin
     WRITE(iscreenlog,*) 'iwavemax = ',iwavemax
     WRITE(iscreenlog,*) 'jwavemin = ',jwavemin
     WRITE(iscreenlog,*) 'jwavemax = ',jwavemax
     WRITE(iscreenlog,*) ' '
     ENDIF_MPI

     ! Save the longitude latitude in 2D variables
     ALLOCATE(lonwave2d(iwavemin:iwavemax,jwavemin:jwavemax))
     ALLOCATE(latwave2d(iwavemin:iwavemax,jwavemin:jwavemax))

     IF (l_wave_marsgrid) THEN
       DO j=jwavemin,jwavemax
         lonwave2d(iwavemin:iwavemax,j)=lonwave2dtmp(iwavemin+1:iwavemax+1,j+1)
       ENDDO
       DO i=iwavemin,iwavemax
         latwave2d(i,jwavemin:jwavemax)=latwave2dtmp(i+1,jwavemin+1:jwavemax+1)
       ENDDO
     ELSE
       DO j=jwavemin,jwavemax
         lonwave2d(iwavemin:iwavemax,j)=lonwave2dtmp(iwavemin:iwavemax,j)
       ENDDO
       DO i=iwavemin,iwavemax
         latwave2d(i,jwavemin:jwavemax)=latwave2dtmp(i,jwavemin:jwavemax)
       ENDDO
     END IF

     DO j=jwavemin,jwavemax
     DO i=iwavemin,iwavemax
       IF(lonwave2d(i,j) >= 180.0_rlg) lonwave2d(i,j)=lonwave2d(i,j)-360.0_rlg
     END DO
     END DO

     ALLOCATE(tab1(iwavemin:iwavemax,jwavemin:jwavemax))
     ALLOCATE(tab2(iwavemin:iwavemax,jwavemin:jwavemax))
     tab1(:,:)=0.0_riosh
     tab2(:,:)=0.0_riosh

#ifdef key_sflx_interp_scrip
     ALLOCATE(masksrc(iwavemin:iwavemax,jwavemin:jwavemax))
     ALLOCATE(tabvar2d(iwavemin:iwavemax,jwavemin:jwavemax))

     ALLOCATE(tabinterp2d(limin:limax,ljmin:ljmax))
     ALLOCATE(maskdst(limin:limax,ljmin:ljmax))

     masksrc(:,:)=.true.
     maskdst(:,:)=.true.
     tabvar2d(:,:)=0.0_rlg
#endif

!!!!!!$OMP end single

   ! -2- Ponderation weights for the spatial interpolation
   ! -----------------------------------------------------

#ifdef key_sflx_interp_scrip
    IF (.NOT. l_wave_marsgrid) THEN
     nxfin=limax-limin+1
     nyfin=ljmax-ljmin+1

!!!!!!$OMP single
     maskdst(:,:)=.true.
     print*,"valfill_wave",valfill_wave
     print*,"min max lon",minval(lonwave2d),maxval(lonwave2d)
     print*,"mask wave",count(maskdst .eqv. .false.)
     where(lon2d(limin:limax,ljmin:ljmax)==-9.e+33) maskdst=.False.
     print*,"mask wave2",count(maskdst .eqv. .false.)
     CALL get_remap_bicub(latwave2d,REAL(lat2d(limin:limax,ljmin:ljmax),rlg),lonwave2d,REAL(lon2d(limin:limax,ljmin:ljmax),rlg), &
                          masksrc,maskdst,matrix,src_add,dst_add)

     maskdst(:,:)=.true.
     print*,"mask wave",count(maskdst .eqv. .false.)
     where(lon2du(limin:limax,ljmin:ljmax)==-9.e+33) maskdst=.False.
     print*,"mask wave2",count(maskdst .eqv. .false.)
     CALL get_remap_bicub(latwave2d,REAL(lat2du(limin:limax,ljmin:ljmax),rlg),lonwave2d,REAL(lon2du(limin:limax,ljmin:ljmax),rlg),&
                          masksrc,maskdst,matrix_u,src_add_u,dst_add_u)

     maskdst(:,:)=.true.
     print*,"mask wave",count(maskdst .eqv. .false.)
     where(lon2dv(limin:limax,ljmin:ljmax)==-9.e+33) maskdst=.False.
     print*,"mask wave2",count(maskdst .eqv. .false.)
     CALL get_remap_bicub(latwave2d,REAL(lat2dv(limin:limax,ljmin:ljmax),rlg),lonwave2d,REAL(lon2dv(limin:limax,ljmin:ljmax),rlg),&
                          masksrc,maskdst,matrix_v,src_add_v,dst_add_v)
!!!!!!$OMP end single
    ENDIF
#else

   !!----------------------------------------------------------------------
   !! * Locate MARS cell from WAVE cell -lmw-
   !!----------------------------------------------------------------------
    IF (.NOT. l_wave_marsgrid) THEN
     dgwave =ABS(lonwave2dtmp(2,1)-lonwave2dtmp(1,1))
     dfiwave=ABS(latwave2dtmp(1,2)-latwave2dtmp(1,1))

!!!!!!$OMP do schedule(static)
     DO j=ljminp1,ljmax
       DO i=liminp1,limax
         xm = (lon2d(i,j)-lonwave2dtmp(1,1))/dgwave+1        !pour les points xe
         ym = (lat2d(i,j)-latwave2dtmp(1,1))/dfiwave+1
         xdpxe_wave(i,j) = xm - REAL(INT(xm),rsh)
         ydpxe_wave(i,j) = ym - REAL(INT(ym),rsh)
         i1xe_wave(i,j) = INT(xm)
         j1xe_wave(i,j) = INT(ym)

         xm = (lon2du(i,j)-lonwave2dtmp(1,1))/dgwave+1        !pour les points u 
         ym = (lat2du(i,j)-latwave2dtmp(1,1))/dfiwave+1
         xdpu_wave(i,j) = xm - REAL(INT(xm),rsh)
         ydpu_wave(i,j) = ym - REAL(INT(ym),rsh)
         i1u_wave(i,j) = INT(xm)
         j1u_wave(i,j) = INT(ym)

         xm = (lon2dv(i,j)-lonwave2dtmp(1,1))/dgwave+1        !pour les points v 
         ym = (lat2dv(i,j)-latwave2dtmp(1,1))/dfiwave+1
         xdpv_wave(i,j) = xm - REAL(INT(xm),rsh)
         ydpv_wave(i,j) = ym - REAL(INT(ym),rsh)
         i1v_wave(i,j) = INT(xm)
         j1v_wave(i,j) = INT(ym)
       END DO
     END DO
    ENDIF
!!!!!!$OMP end do

#endif

   ! Time index in the wave file
   ! ---------------------------

     IF_MPI (MASTER) THEN

!!!!!!$OMP single
     idimt=ionc4_read_dimt(file_wave)
     ilec_wave=0 
     tnc1_wave=0.0_rlg
     tnc2_wave=0.0_rlg
     DO WHILE((twave-(tnc1_wave-torigin))*(twave-(tnc2_wave-torigin)) > 0.0_rlg .AND. ilec_wave < idimt)
       ilec_wave=ilec_wave+1
       IF (ilec_wave == idimt) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'All records are read, dates are not compatible with simulation time.'
         WRITE(ierrorlog,*) 'Make sure your wave file is correct: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'Make sure the origin time is 01/01/1900 00:00:00 in the wave file.'
         WRITE(ierrorlog,*) 'If not, modify variable "torigin" in the routine WAVE.F90.'
         WRITE(ierrorlog,*) ''
         tdate = tool_sectodat(twave)
         WRITE(ierrorlog,*) 't in MARS model =',tdate
         tdate = tool_sectodat(tnc1_wave-torigin)
         WRITE(ierrorlog,*) 't1 in wave file =',tdate
         tdate = tool_sectodat(tnc2_wave-torigin)
         WRITE(ierrorlog,*) 't2 in wave file =',tdate
         WRITE(ierrorlog,*) 'Simulation stopped.'
         WRITE(ierrorlog,*) ''
         l_pb = .true.
       END IF
       IF (.NOT. l_pb) THEN
         CALL ionc4_read_time(file_wave,ilec_wave,tnc1_wave)
         tnc1_wave=tnc1_wave*86400.0_rsh+t_lag
         print*,tool_sectodat(tnc1_wave)
         CALL ionc4_read_time(file_wave,ilec_wave+1,tnc2_wave)
         tnc2_wave=tnc2_wave*86400.0_rsh+t_lag
         !print*,'time wave',tnc1_wave,tnc1_wave,twave
         print*,tool_sectodat(tnc2_wave)
       END IF
     END DO
#ifdef key_castest
     ilec_wave=10
#endif
     ilecmem_wave=ilec_wave
     WRITE(iscreenlog,*) 'ilec_wave',ilec_wave
!!!!!!$OMP end single

     ENDIF_MPI

     ! the MASTER sends informations to every cpus
OMPMPI barrier
OMPMPI master
     CALL_MPI BCAST_MPI(l_pb)
     CALL_MPI BCAST_MPI(idimt)
     CALL_MPI BCAST_MPI(ilec_wave)
     CALL_MPI BCAST_MPI(ilecmem_wave)
     CALL_MPI BCAST_MPI(tnc1_wave)
     CALL_MPI BCAST_MPI(tnc2_wave)
OMPMPI end master
OMPMPI barrier
OMPMPI flush(l_pb,idimt,ilec_wave,ilecmem_wave,tnc1_wave,tnc2_wave)

     IF (l_pb) THEN
       WRITE(ierrorlog,*) "PB in WAVE_READ"
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF

     ! -3- lecture des tableaux deux indices qui encadrent le pas de temps courant

   !!----------------------------------------------------------------------
   !! * Read and spatial interpolate of 2 time step of input file -ri2-
   !!----------------------------------------------------------------------
   !! -nwq-
     IF (l_wave_rdir) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_dir),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_dir),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_dir1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_dir2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_dir1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_dir2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_dir1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_dir2)
#endif
       END IF

     END IF

     IF (l_wave_rhs) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_hs),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_hs),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_hs1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_hs2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_hs1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_hs2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_hs1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_hs2)
#endif
       END IF

     END IF


     IF (l_wave_rubot) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ubot),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ubot),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_ubot1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_ubot2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ubot1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ubot2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_ubot1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_ubot2)
#endif
       END IF

     END IF

     IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_uubot),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_uubot),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_uubot1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_uubot2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_uubot1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_uubot2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_uubot1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_uubot2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_vubot),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_vubot),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_vubot1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_vubot2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_vubot1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_vubot2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_vubot1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_vubot2)
#endif
       END IF

     END IF

     IF (l_wave_rfoc) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_foc),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_foc),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_foc1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_foc2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_foc1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_foc2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_foc1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_foc2)
#endif
       END IF

     END IF

     IF (l_wave_rfbb) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_fbb),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_fbb),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_fbb1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_fbb2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_fbb1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_fbb2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_fbb1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_fbb2)
#endif
       END IF

     END IF

     IF (l_wave_rtp) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tp),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tp),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tp1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tp2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tp1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tp2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_tp1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_tp2)
#endif
       END IF

     END IF

     IF (l_wave_rtauwx .AND. l_wave_rtauwy) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauwx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauwx),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauwx1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tauwx2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauwx1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauwx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tauwx1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_tauwx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauwy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauwy),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauwy1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tauwy2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauwy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauwy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tauwy1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_tauwy2)
#endif
       END IF

     END IF


     IF (l_wave_rtauox .AND. l_wave_rtauoy) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauox),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauox),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauox1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tauox2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauox1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauox2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tauox1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_tauox2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauoy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauoy),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauoy1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tauoy2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauoy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauoy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tauoy1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_tauoy2)
#endif
       END IF

     END IF

     IF (l_wave_rtbb) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tbbx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tbbx),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tbbx1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tbbx2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tbbx1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tbbx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tbbx1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_tbbx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tbby),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tbby),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tbby1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tbby2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tbby1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tbby2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tbby1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_tbby2)
#endif
       END IF

     END IF
     IF (l_wave_rtus) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tusx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tusx),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tusx1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tusx2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tusx1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tusx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tusx1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_tusx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tusy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tusy),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tusy1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_tusy2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tusy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tusy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tusy1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_tusy2)
#endif
       END IF

     END IF

     IF (l_wave_russ) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ussx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ussx),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_ussx1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_ussx2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ussx1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ussx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_ussx1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_ussx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ussy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ussy),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_ussy1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_ussy2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ussy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ussy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_ussy1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_ussy2)
#endif
       END IF

     END IF

     IF (l_wave_rwlm) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_wlm),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_wlm),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_wlm1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_wlm2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_wlm1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_wlm2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_wlm1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_wlm2)
#endif
       END IF

     END IF

     IF (l_wave_rt0m1) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_t0m1),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_t0m1),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_t0m11(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_t0m12(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_t0m11(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_t0m12(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_t0m11)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_t0m12)
#endif
       END IF

     END IF

     IF (l_wave_rbhd) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_bhd),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_bhd),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_bhd1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_bhd2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_bhd1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_bhd2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_bhd1)
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab2,wave_bhd2)
#endif
       END IF

     END IF

     IF (l_wave_rradstr) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sx),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_sx1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_sx2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_sx1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_sx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_sx1)
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab2,wave_sx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sy),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_sy1(limin:limax,ljmin:ljmax)=tab1(:,:)
         wave_sy2(limin:limax,ljmin:ljmax)=tab2(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_sy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
         tabvar2d(:,:)=tab2(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_sy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_sy1)
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_sy2)
#endif
       END IF

       IF (.not. l_swan) THEN ! if not swan we have to read sxy. in SWAN Fx and Fy are given and there is no need to compute it
!!!!!!!$OMP single
        CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sxy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave, &
                               iwmin0,iwmin0,fillvar=0.0_riosh)
        CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sxy),tab2,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                               iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!!$OMP end single
 
        IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!!$OMP single
          wave_sxy1(limin:limax,ljmin:ljmax)=tab1(:,:)
          wave_sxy2(limin:limax,ljmin:ljmax)=tab2(:,:)
 !!!!!!$OMP end single
        ELSE
#ifdef key_sflx_interp_scrip
!!!!!!!$OMP single
          tabvar2d(:,:)=tab1(:,:)
          CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
          wave_sxy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
          tabvar2d(:,:)=tab2(:,:)
          CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
          wave_sxy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!!$OMP end single
#else
          CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_sxy1)
          CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab2,wave_sxy2)
#endif
        END IF
       END IF
     END IF

!!!!!!$OMP single
     DEALLOCATE(tab1,tab2)
     IF (dimlon == 1) DEALLOCATE(latwave,lonwave) 
     DEALLOCATE(lonwave2d,latwave2d,lonwave2dtmp,latwave2dtmp)
!!!!!!$OMP end single

   !!----------------------------------------------------------------------
   !! * Checkout units -cou-
   !!----------------------------------------------------------------------
!!!!!!$OMP single
     ncerror=0

     IF (l_wave_rdir) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_dir1(i,j),wave_dir2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_dir1(i,j),wave_dir2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<0_rsh).OR.(tabmax>360_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the direction of the wave is not expressed in [degree].'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_dir)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF

     IF (l_wave_rhs) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_hs1(i,j),wave_hs2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_hs1(i,j),wave_hs2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<0_rsh).OR.(tabmax>20_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'wave hight (Hs) is not expressed in [m].'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_hs)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF

     IF (l_wave_rfoc) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_foc1(i,j),wave_foc2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_foc1(i,j),wave_foc2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<-1.0_rsh).OR.(tabmax>100.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the coefficient of energy flux of wave to ocean (due to surface breaking) expressed in [W.m-2] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_foc)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF

     IF (l_wave_rfbb) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_fbb1(i,j),wave_fbb2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_fbb1(i,j),wave_fbb2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<-1.0_rsh).OR.(tabmax>100.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the coefficient of energy flux of wave to ocean (due to bottom dissipation) expressed in [W.m-2] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_fbb)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF

     IF (l_wave_rubot) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_ubot1(i,j),wave_ubot2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_ubot1(i,j),wave_ubot2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<-5.0_rsh).OR.(tabmax>5.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the orbital velocity expressed in [m/s] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_ubot)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF

     IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_uubot1(i,j),wave_uubot2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_uubot1(i,j),wave_uubot2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<-5.0_rsh).OR.(tabmax>5.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the Ubot (U) expressed in [m/s] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_uubot)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_vubot1(i,j),wave_vubot2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_vubot1(i,j),wave_vubot2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<-5.0_rsh).OR.(tabmax>5.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the Ubot (V) expressed in [m/s] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_vubot)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF



     IF (l_wave_rtp) THEN

       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_tp1(i,j),wave_tp2(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_tp1(i,j),wave_tp2(i,j))))
           ENDIF
         ENDDO
       ENDDO
       ! print*,'tabmin tp=SQRT',tabmin,tabmax
       IF ((tabmin<0_rsh).OR.(tabmax>50_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the peak wave period expressed in [s] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_tp)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF



     IF (l_wave_rwlm) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_t0m11(i,j),wave_t0m12(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_t0m11(i,j),wave_t0m12(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<0.0_rsh).OR.(tabmax>wave_wlm_max)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the wave length expressed in [m] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_tp)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF


     IF (l_wave_rt0m1) THEN
       tabmin=1000000.0_rsh
       tabmax=-10000000.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,(MIN(wave_t0m11(i,j),wave_t0m12(i,j))))
             tabmax=MAX(tabmax,(MAX(wave_t0m11(i,j),wave_t0m12(i,j))))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<0.0_rsh).OR.(tabmax>wave_t0m1_max)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
         WRITE(ierrorlog,*) 'In the waves file: ',TRIM(file_wave)
         WRITE(ierrorlog,*) 'the mean wave period expressed in [s] is not correct.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_wave_tp)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       END IF
     END IF


     IF (ncerror>0) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
!!!!!!$OMP end single

   END IF ! ENDOF if (nbouc==1) then -init-

#if ! defined key_castest
   ! deroulement normal

   !!----------------------------------------------------------------------
   !! * Find time step in the file that exceed computational time step -fse-
   !!----------------------------------------------------------------------
!!!!!!$OMP single
   DO WHILE((twave-(tnc1_wave-torigin))*(twave-(tnc2_wave-torigin)) > 0.0_rlg .AND. ilec_wave < idimt)
     ilec_wave=ilec_wave+1 
     WRITE(iscreenlog,*) __LINE__,'ilec_wave',ilec_wave
     IF (ilec_wave == idimt) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'ERROR: routine WAVE_READ'
       WRITE(ierrorlog,*) 'All records are already read.'
       WRITE(ierrorlog,*) 'The waves file is not long enough: ',TRIM(file_wave)
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
     CALL ionc4_read_time(file_wave,ilec_wave,tnc1_wave)
     CALL ionc4_read_time(file_wave,ilec_wave+1,tnc2_wave)
     tnc1_wave=tnc1_wave*86400.0_rsh+t_lag
     tnc2_wave=tnc2_wave*86400.0_rsh+t_lag
   END DO
!!!!!!$OMP end single

   IF (ilecmem_wave/=ilec_wave) THEN

     WRITE(iscreenlog,*) __LINE__,'ilec_wave',ilec_wave,ilecmem_wave
     IF (nbouc/=1) THEN
!!!!!!$OMP single
       ALLOCATE(tab1(iwavemin:iwavemax,jwavemin:jwavemax))
!!!!!!$OMP end single
     END IF

     ilecmem_wave=ilec_wave
   !!----------------------------------------------------------------------
   !! * We transfert array 2 in array 1 -2i1-
   !!----------------------------------------------------------------------
!!!!!!$OMP do schedule(static)
     DO j=ljmin,ljmax
   !! -nwq-
       IF (l_wave_rdir) wave_dir1(:,j)=wave_dir2(:,j)
       IF (l_wave_rhs)  wave_hs1(:,j) =wave_hs2(:,j)
       IF (l_wave_rubot) wave_ubot1(:,j)=wave_ubot2(:,j)
       IF (l_wave_ruubot) wave_uubot1(:,j)=wave_uubot2(:,j)
       IF (l_wave_rvubot) wave_vubot1(:,j)=wave_vubot2(:,j)
       IF (l_wave_rtp) wave_tp1(:,j)=wave_tp2(:,j)
       IF (l_wave_rfoc) wave_foc1(:,j)=wave_foc2(:,j)
       IF (l_wave_rfbb) wave_fbb1(:,j)=wave_fbb2(:,j)
       IF (l_wave_rtauwx) wave_tauwx1(:,j)=wave_tauwx2(:,j)
       IF (l_wave_rtauwy) wave_tauwy1(:,j)=wave_tauwy2(:,j)
       IF (l_wave_rtauox) wave_tauox1(:,j)=wave_tauox2(:,j)
       IF (l_wave_rtauoy) wave_tauoy1(:,j)=wave_tauoy2(:,j)
       IF (l_wave_rtbb) wave_tbbx1(:,j)=wave_tbbx2(:,j)
       IF (l_wave_rtbb) wave_tbby1(:,j)=wave_tbby2(:,j)
       IF (l_wave_rtus) wave_tusx1(:,j)=wave_tusx2(:,j)
       IF (l_wave_rtus) wave_tusy1(:,j)=wave_tusy2(:,j)
       IF (l_wave_russ) wave_ussx1(:,j)=wave_ussx2(:,j)
       IF (l_wave_russ) wave_ussy1(:,j)=wave_ussy2(:,j)
       IF (l_wave_rwlm) wave_wlm1(:,j)=wave_wlm2(:,j)
       IF (l_wave_rt0m1) wave_t0m11(:,j)=wave_t0m12(:,j)
       IF (l_wave_rbhd) wave_bhd1(:,j)=wave_bhd2(:,j)
       IF (l_wave_rradstr) THEN
         wave_sx1(:,j)=wave_sx2(:,j)
         wave_sy1(:,j)=wave_sy2(:,j)
         IF (.NOT. l_swan) wave_sxy1(:,j)=wave_sxy2(:,j)
       ENDIF
     END DO
!!!!!!$OMP end do

   !!----------------------------------------------------------------------
   !! * Read and spatial interpolate next time step of input file while computational time step is not exceed -rin-
   !!----------------------------------------------------------------------
   !! -nwq-
     IF (l_wave_rdir) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_dir),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_dir2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_dir2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_dir2)
#endif
       END IF

     END IF

     IF (l_wave_rhs) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_hs),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_hs2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_hs2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_hs2)
#endif
       END IF

     END IF

     IF (l_wave_rubot) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ubot),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_ubot2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ubot2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_ubot2)
#endif
       END IF

     END IF

     IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_uubot),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_uubot2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_uubot2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_uubot2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_vubot),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_vubot2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_vubot2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_vubot2)
#endif
       END IF

     END IF

     IF (l_wave_rfoc) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_foc),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_foc2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_foc2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_foc2)
#endif
       END IF

     END IF
     
     IF (l_wave_rfbb) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_fbb),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_fbb2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_fbb2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_fbb2)
#endif
       END IF

     END IF
     
     IF (l_wave_rtp) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tp),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tp2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tp2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_tp2)
#endif
       END IF

     END IF

     IF (l_wave_rtauwx .AND. l_wave_rtauwy) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauwx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauwx2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauwx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tauwx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauwy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
          wave_tauwy2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauwy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tauwy2)
#endif
       END IF

     END IF

     IF (l_wave_rtauox .AND. l_wave_rtauoy) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauox),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauox2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauox2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tauox2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tauoy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tauoy2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tauoy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tauoy2)
#endif
       END IF

     END IF


     IF (l_wave_rtbb) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tbbx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tbbx2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tbbx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tbbx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tbby),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tbby2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tbby2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tbby2)
#endif
       END IF

     END IF

     IF (l_wave_rtus) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tusx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tusx2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tusx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_tusx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_tusy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_tusy2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_tusy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_tusy2)
#endif
       END IF

     END IF

     IF (l_wave_russ) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ussx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_ussx2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ussx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_ussx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_ussy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_ussy2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_ussy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_ussy2)
#endif
       END IF

     END IF

     IF (l_wave_rwlm) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_wlm),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_wlm2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_wlm2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_wlm2)
#endif
       END IF

     END IF
 
     IF (l_wave_rt0m1) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_t0m1),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_t0m12(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_t0m12(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_t0m12)
#endif
       END IF

     END IF

     IF (l_wave_rbhd) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_bhd),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_bhd2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_bhd2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1xe_wave,j1xe_wave,xdpxe_wave,ydpxe_wave,tab1,wave_bhd2)
#endif
       END IF

     END IF

     IF (l_wave_rradstr) THEN
!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sx),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_sx2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_sx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1u_wave,j1u_wave,xdpu_wave,ydpu_wave,tab1,wave_sx2)
#endif
       END IF

!!!!!!$OMP single
       CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                              iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!$OMP end single

       IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!$OMP single
         wave_sy2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!$OMP end single
       ELSE
#ifdef key_sflx_interp_scrip
!!!!!!$OMP single
         tabvar2d(:,:)=tab1(:,:)
         CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
         wave_sy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!$OMP end single
#else
         CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_sy2)
#endif
       END IF
       IF (.NOT. l_swan) THEN
!!!!!!!$OMP single
        CALL ionc4_read_subxyt(file_wave,TRIM(name_wave_sxy),tab1,iwavemin,iwavemax,jwavemin,jwavemax,ilec_wave+1, &
                               iwmin0,iwmin0,fillvar=0.0_riosh)
!!!!!!!$OMP end single
 
        IF (l_wave_marsgrid) THEN  ! no interpolation : the grid of the wave file is equal to the MARS grid
!!!!!!!$OMP single
          wave_sxy2(limin:limax,ljmin:ljmax)=tab1(:,:)
!!!!!!!$OMP end single
        ELSE
#ifdef key_sflx_interp_scrip
!!!!!!!$OMP single
          tabvar2d(:,:)=tab1(:,:)
          CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
          wave_sxy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!!!!!!!$OMP end single
#else
          CALL tool_interp2d_wave(i1v_wave,j1v_wave,xdpv_wave,ydpv_wave,tab1,wave_sxy2)
#endif
        END IF
       END IF
     END IF


!!!!!!$OMP single
     DEALLOCATE(tab1)
!!!!!!$OMP end single

   END IF  ! test ilecmem_wave/=ilec_wave

! fin test sur cas test (choix d un unique pas de temps)
#endif

   !!----------------------------------------------------------------------
   !! * Time linear interpolation -tli-
   !!----------------------------------------------------------------------
   dt1=(tnc2_wave-torigin-twave)/(tnc2_wave-tnc1_wave)
   dt2=1.0_rlg-dt1
#ifdef key_castest
   dt1=1.0_rlg
   dt2=0.0_rlg
#endif

   !! -nwq-
!!!!!!$OMP do schedule(static)
   DO j=ljmin,ljmax
     IF (l_wave_rdir) wave_dir(limin:limax,j) =  wave_dir1(limin:limax,j)*dt1   &
           + wave_dir2(limin:limax,j)*dt2
     IF (l_wave_rhs) wave_hs(limin:limax,j)   =  wave_hs1(limin:limax,j)*dt1   &
           + wave_hs2(limin:limax,j)*dt2
     IF (l_wave_rubot) wave_ubot(limin:limax,j) =  wave_ubot1(limin:limax,j)*dt1   &
           + wave_ubot2(limin:limax,j)*dt2
     IF (l_wave_ruubot) wave_uubot(limin:limax,j) =  wave_uubot1(limin:limax,j)*dt1   &
           + wave_uubot2(limin:limax,j)*dt2
     IF (l_wave_rvubot) wave_vubot(limin:limax,j) =  wave_vubot1(limin:limax,j)*dt1   &
           + wave_vubot2(limin:limax,j)*dt2
     IF (l_wave_rfoc) THEN
       wave_phioc(limin:limax,j) =  (wave_foc1(limin:limax,j)*dt1   &
           + wave_foc2(limin:limax,j)*dt2) / 1000.   ! (in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
       WHERE(wave_phioc(limin:limax,j) < 0.0_rsh) wave_phioc(limin:limax,j)=0.0_rsh
     END IF
     IF (l_wave_rfbb) THEN
       wave_phibbl(limin:limax,j) =  (wave_fbb1(limin:limax,j)*dt1   &
           + wave_fbb2(limin:limax,j)*dt2) / 1000.   ! (in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
       WHERE(wave_phibbl(limin:limax,j) < 0.0_rsh) wave_phibbl(limin:limax,j)=0.0_rsh
     END IF
     IF (l_wave_rtp) wave_tp(limin:limax,j) =  wave_tp1(limin:limax,j)*dt1   &
           + wave_tp2(limin:limax,j)*dt2
     IF (l_wave_rtauwx) wave_tauwx(limin:limax,j) = wave_tauwx1(limin:limax,j)*dt1   &
           + wave_tauwx2(limin:limax,j)*dt2
     IF (l_wave_rtauwy) wave_tauwy(limin:limax,j) = wave_tauwy1(limin:limax,j)*dt1   &
           + wave_tauwy2(limin:limax,j)*dt2
     IF (l_wave_rtauox) wave_tauox(limin:limax,j) = wave_tauox1(limin:limax,j)*dt1   &
           + wave_tauox2(limin:limax,j)*dt2
     IF (l_wave_rtauoy) wave_tauoy(limin:limax,j) = wave_tauoy1(limin:limax,j)*dt1   &
           + wave_tauoy2(limin:limax,j)*dt2
     IF (l_wave_rtbb) wave_taubbl_x(limin:limax,j) = wave_tbbx1(limin:limax,j)*dt1   &
           + wave_tbbx2(limin:limax,j)*dt2
     IF (l_wave_rtbb) wave_taubbl_y(limin:limax,j) = wave_tbby1(limin:limax,j)*dt1   &
           + wave_tbby2(limin:limax,j)*dt2
     IF (l_wave_rtus) wave_tusx(limin:limax,j) = wave_tusx1(limin:limax,j)*dt1   &
           + wave_tusx2(limin:limax,j)*dt2
     IF (l_wave_rtus) wave_tusy(limin:limax,j) = wave_tusy1(limin:limax,j)*dt1   &
           + wave_tusy2(limin:limax,j)*dt2
     IF (l_wave_russ) wave_ussx(limin:limax,j) = wave_ussx1(limin:limax,j)*dt1   &
           + wave_ussx2(limin:limax,j)*dt2
     IF (l_wave_russ) wave_ussy(limin:limax,j) = wave_ussy1(limin:limax,j)*dt1   &
           + wave_ussy2(limin:limax,j)*dt2
     IF (l_wave_rwlm) wave_wlm(limin:limax,j) = wave_wlm1(limin:limax,j)*dt1   &
           + wave_wlm2(limin:limax,j)*dt2
     IF (l_wave_rt0m1) wave_t0m1(limin:limax,j) = wave_t0m11(limin:limax,j)*dt1   &
           + wave_t0m12(limin:limax,j)*dt2
     IF (l_wave_rbhd) wave_bhd(limin:limax,j) = wave_bhd1(limin:limax,j)*dt1   &
           + wave_bhd2(limin:limax,j)*dt2
     IF (l_wave_rradstr) THEN
      wave_sx(limin:limax,j) =  wave_sx1(limin:limax,j)*dt1   &
           + wave_sx2(limin:limax,j)*dt2
      wave_sy(limin:limax,j) =  wave_sy1(limin:limax,j)*dt1   &
           + wave_sy2(limin:limax,j)*dt2
      IF (.NOT. l_swan) wave_sxy(limin:limax,j) =  wave_sxy1(limin:limax,j)*dt1   &
           + wave_sxy2(limin:limax,j)*dt2
     ENDIF

   END DO
!!!!!!$OMP end do

   IF ((l_wave_ruubot .AND. l_wave_rvubot) .AND. (.NOT. l_wave_rubot)) THEN 
!!!!!!$OMP do schedule(static)
     DO j=ljmin,ljmax
        DO i=limin,limax
           wave_ubot(i,j)=sqrt(wave_uubot(i,j)**2+wave_vubot(i,j)**2)
        END DO
     END DO
!!!!!!$OMP end do
   END IF

!!!!!!$OMP end parallel

   !     !!----------------------------------------------------------------------
   !     !! * MPI managing for gradient computation -mmg-
   !     !!----------------------------------------------------------------------
   !     !! -nwq-
   !toseelater on ne fait les echanges que s il y a besoin. Donc voir quels indices
   !vont etre utiles: a garder en commentaire car ca peux etre utils.
OMPMPI barrier
OMPMPI master
   IF_MPI (l_wave_rdir) THEN
     CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,wave_dir(limin:limaxp1,ljmin:ljmaxp1))
     CALL_MPI ex_j_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,wave_dir(limin:limaxp1,ljmin:ljmaxp1))
   ENDIF_MPI
   IF_MPI (l_wave_rbhd) THEN
     CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,wave_bhd(liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,wave_bhd(liminm2:limaxp2,ljminm2:ljmaxp2))
   ENDIF_MPI
   IF_MPI (l_wave_rradstr) THEN
     CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmax,wave_sx(limin:limaxp1,ljmin:ljmax))
     CALL_MPI ex_j_rsh(0,1,1,limin,limax,ljmin,ljmaxp1,wave_sy(limin:limax,ljmin:ljmaxp1))
     CALL_MPI ex_i_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,wave_sxy(liminm1:limaxp1,ljminm1:ljmaxp1))
     CALL_MPI ex_j_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,wave_sxy(liminm1:limaxp1,ljminm1:ljmaxp1))
   ENDIF_MPI
   IF_MPI (l_wave_rtauox) CALL ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmax,wave_tauox(limin:limaxp1,ljmin:ljmax))
   IF_MPI (l_wave_rtauoy) CALL ex_j_rsh(0,1,1,limin,limax,ljmin,ljmaxp1,wave_tauoy(limin:limax,ljmin:ljmaxp1))
   IF_MPI (l_wave_rtauwx) CALL ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmax,wave_tauwx(limin:limaxp1,ljmin:ljmax))
   IF_MPI (l_wave_rtauwy) CALL ex_j_rsh(0,1,1,limin,limax,ljmin,ljmaxp1,wave_tauwy(limin:limax,ljmin:ljmaxp1))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(wave_tauwy,wave_tauwx,wave_tauoy,wave_tauox,wave_bhd)


   IF_MPI (MASTER) THEN
   IF (dbug.and.(mod(nbouc,2)==0)) THEN
   !! -nwq-
     print *,'dir=',sum(sum(wave_dir,1),1)/(size(wave_dir,1)*size(wave_dir,2))
     print *,'hs=',sum(sum(wave_hs,1),1)/(size(wave_hs,1)*size(wave_hs,2))
     print *,'ubot=',sum(sum(wave_ubot,1),1)/(size(wave_ubot,1)*size(wave_ubot,2))
     print *,'tp=',sum(sum(wave_tp,1),1)/(size(wave_tp,1)*size(wave_tp,2))
     print *,'foc=',sum(sum(wave_phioc,1),1)/(size(wave_phioc,1)*size(wave_phioc,2))
     print *,'fbb=',sum(sum(wave_phibbl,1),1)/(size(wave_phibbl,1)*size(wave_phibbl,2))
!     print *,'rtauwx=',sum(sum(wave_tauwx,1),1)/(size(wave_tauwx,1)*size(wave_tauwx,2))
!     print *,'rtauwy=',sum(sum(wave_tauwy,1),1)/(size(wave_tauwy,1)*size(wave_tauwy,2))
     print *,'rtauox=',sum(sum(wave_tauox,1),1)/(size(wave_tauox,1)*size(wave_tauox,2))
     print *,'rtauoy=',sum(sum(wave_tauoy,1),1)/(size(wave_tauoy,1)*size(wave_tauoy,2))
     print *,'rtbbx=',sum(sum(wave_taubbl_x,1),1)/(size(wave_taubbl_x,1)*size(wave_taubbl_x,2))
     print *,'rtbby=',sum(sum(wave_taubbl_y,1),1)/(size(wave_taubbl_y,1)*size(wave_taubbl_y,2))
     print *,'rtusx=',sum(sum(wave_tusx,1),1)/(size(wave_tusx,1)*size(wave_tusx,2))
     print *,'rtusy=',sum(sum(wave_tusy,1),1)/(size(wave_tusy,1)*size(wave_tusy,2))
     print *,'russx=',sum(sum(wave_ussx,1),1)/(size(wave_ussx,1)*size(wave_ussx,2))
     print *,'russy=',sum(sum(wave_ussy,1),1)/(size(wave_ussy,1)*size(wave_ussy,2))
     print *,'rwlm=',sum(sum(wave_wlm,1),1)/(size(wave_wlm,1)*size(wave_wlm,2))
     print *,'rt0m1=',sum(sum(wave_t0m1,1),1)/(size(wave_t0m1,1)*size(wave_t0m1,2))
     print *,'rbhd=',sum(sum(wave_bhd,1),1)/(size(wave_bhd,1)*size(wave_bhd,2))
   END IF
   IF (dbug) THEN
     print *,'OUT OF WAVEREAD'
   END IF
   ENDIF_MPI

   PRINT_DBG*, 'exit WAVE_READ routine'

  END SUBROUTINE wave_read

   !!======================================================================

  SUBROUTINE wave_alloc

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE wave_alloc  ***
   !&E
   !&E ** Purpose : allocation des variables lues depuis WW3
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E       ! module alloc (file alloc.F90)
   !&E
   !&E ** History :
   !&E       !  2010-02 (G. Simon) Creation  on the basis of alloc_cometdyn
   !&E       !  2013-08 (AC. Bennis) Introduction of arrays for 3D radiation stress tensor
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments

   !! * Local declarations


   !!----------------------------------------------------------------------
   !! * Executable part

   !! * Locate MARS cell from Wave cell
   ALLOCATE(i1v_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(j1v_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(i1u_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(j1u_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(i1xe_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(j1xe_wave(limin:limax,ljmin:ljmax))

   ALLOCATE(xdpxe_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(xdpu_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(xdpv_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(ydpxe_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(ydpu_wave(limin:limax,ljmin:ljmax))
   ALLOCATE(ydpv_wave(limin:limax,ljmin:ljmax))

#if ! defined key_wave_readww3txt && ! defined key_oasis_mars_ww3
   !! * Waves quantities to be read
   !!   - Mean wave direction (degrees)
!   IF (l_wave_rdir) THEN
   ALLOCATE(wave_dir1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_dir2(limin:limax,ljmin:ljmax))
   wave_dir1(:,:)=0.0_rsh ; wave_dir2(:,:)=0.0_rsh
!   END IF
   !!   - Significant wave height (hs in m)
!   IF (l_wave_rhs) THEN
   ALLOCATE(wave_hs1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_hs2(limin:limax,ljmin:ljmax))
   wave_hs1(:,:) =0.0_rsh ; wave_hs2(:,:) =0.0_rsh
!   END IF
   !!   - Wave-to-ocean TKE flux (phioc in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
!   IF (l_wave_rfoc) THEN
   ALLOCATE(wave_foc1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_foc2(limin:limax,ljmin:ljmax))
   wave_foc1(:,:)=0.0_rsh ; wave_foc2(:,:)=0.0_rsh
!   END IF
   !!   - Wave-to-ocean TKE flux due to bottom friction (phibbl in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
!   IF (l_wave_rfbb) THEN
   ALLOCATE(wave_fbb1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_fbb2(limin:limax,ljmin:ljmax))
   wave_fbb1(:,:)=0.0_rsh ; wave_fbb2(:,:)=0.0_rsh
!   END IF
   !!   - Wave-to-ocean momentim due to bottom friction (taubbl in m2 s-2)
!   IF (l_wave_rtbb) THEN
   ALLOCATE(wave_tbbx1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tbbx2(limin:limax,ljmin:ljmax))
   wave_tbbx1(:,:)=0.0_rsh ; wave_tbbx2(:,:)=0.0_rsh
   ALLOCATE(wave_tbby1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tbby2(limin:limax,ljmin:ljmax))
   wave_tbby1(:,:)=0.0_rsh ; wave_tbby2(:,:)=0.0_rsh
!   END IF
!   END IF
   !!   - orbital velocity
   !!   - rms amplitude of orbital bottom velocity of the waves (uba in m/s) (Bottom rms velocity (m/s))
!   IF (l_wave_rubot) THEN
   ALLOCATE(wave_ubot1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_ubot2(limin:limax,ljmin:ljmax))
!   END IF
!   IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
   wave_ubot1(:,:)=0.0_rsh ; wave_ubot2(:,:)=0.0_rsh
   ALLOCATE(wave_uubot1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_uubot2(limin:limax,ljmin:ljmax))
   wave_uubot1(:,:)=0.0_rsh ; wave_uubot2(:,:)=0.0_rsh
   ALLOCATE(wave_vubot1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_vubot2(limin:limax,ljmin:ljmax))
   wave_vubot1(:,:)=0.0_rsh ; wave_vubot2(:,:)=0.0_rsh
!   END IF
   !!   - Mean wave period (tmn in s) (m0,-1)
!   IF (l_wave_rt0m1) THEN
   ALLOCATE(wave_t0m11(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_t0m12(limin:limax,ljmin:ljmax))
   wave_t0m11(:,:)=0.0_rsh ; wave_t0m12(:,:)=0.0_rsh
!   END IF
   !!   - wave period
!   IF (l_wave_rtp) THEN
   ALLOCATE(wave_tp1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tp2(limin:limax,ljmin:ljmax))
   wave_tp1(:,:)=0.0_rsh ; wave_tp2(:,:)=0.0_rsh
!   END IF
   !!   - Wave-ocean momentum flux (m2.s-2)
!   IF (l_wave_rtauox .AND. l_wave_rtauoy) THEN
   ALLOCATE(wave_tauox1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tauox2(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tauoy1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tauoy2(limin:limax,ljmin:ljmax))
   wave_tauox1(:,:)=0.0_rsh
   wave_tauox2(:,:)=0.0_rsh
   wave_tauoy1(:,:)=0.0_rsh
   wave_tauoy2(:,:)=0.0_rsh
!   END IF
   !!   - Net wave-supported stress (m2.s-2)
!   IF (l_wave_rtauwx .AND. l_wave_rtauwy) THEN
   ALLOCATE(wave_tauwx1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tauwx2(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tauwy1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tauwy2(limin:limax,ljmin:ljmax))
   wave_tauwx1(:,:)=0.0_rsh
   wave_tauwx2(:,:)=0.0_rsh
   wave_tauwy1(:,:)=0.0_rsh
   wave_tauwy2(:,:)=0.0_rsh
!   END IF
   !!   - surface stokes drift (m2 s-1)
   ALLOCATE(wave_ussx1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_ussx2(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_ussy1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_ussy2(limin:limax,ljmin:ljmax))
   wave_ussx1(:,:)=0.0_rsh
   wave_ussx2(:,:)=0.0_rsh
   wave_ussy1(:,:)=0.0_rsh
   wave_ussy2(:,:)=0.0_rsh
   !!   - stokes transport (m2 s-1)
   ALLOCATE(wave_tusx1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tusx2(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tusy1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tusy2(limin:limax,ljmin:ljmax))
   wave_tusx1(:,:)=0.0_rsh
   wave_tusx2(:,:)=0.0_rsh
   wave_tusy1(:,:)=0.0_rsh
   wave_tusy2(:,:)=0.0_rsh
   !!   - Mean wave length (wlm in m)
!   IF (l_wave_rwlm) THEN
   ALLOCATE(wave_wlm1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_wlm2(limin:limax,ljmin:ljmax))
   wave_wlm1(:,:)=0.0_rsh
   wave_wlm2(:,:)=0.0_rsh
!   END IF
   !!   - Wave-induced Bernoulli head pressure (bhd in N.m-1) (J term, Smith JPO 2006)
!   IF (l_wave_rbhd) THEN
   ALLOCATE(wave_bhd1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_bhd2(limin:limax,ljmin:ljmax))
   wave_bhd1(:,:)=0.0_rsh
   wave_bhd2(:,:)=0.0_rsh
!   END IF
   !! radiation stress
   ALLOCATE(wave_sx1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_sx2(limin:limax,ljmin:ljmax))
   wave_sx1(:,:)=0.0_rsh ; wave_sx2(:,:)=0.0_rsh
   ALLOCATE(wave_sy1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_sy2(limin:limax,ljmin:ljmax))
   wave_sy1(:,:)=0.0_rsh ; wave_sy2(:,:)=0.0_rsh
   ALLOCATE(wave_sxy1(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_sxy2(limin:limax,ljmin:ljmax))
   wave_sxy1(:,:)=0.0_rsh ; wave_sxy2(:,:)=0.0_rsh

#endif

   !!   - Significant wave height (hs in m)
   ALLOCATE(wave_hs(limin:limax,ljmin:ljmax))
   wave_hs(:,:)=0.0_rsh
   !!   - Period of the peak wave (tp in s)
   ALLOCATE(wave_tp(limin:limax,ljmin:ljmax))
   wave_tp(:,:)=0.0_rsh
   !!   - rms amplitude of orbital wave bottom velocity (uba in m/s)
   ALLOCATE(wave_ubot(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_uubot(limin:limax,ljmin:ljmax))  ! x component
   ALLOCATE(wave_vubot(limin:limax,ljmin:ljmax))  ! y component
   wave_ubot(:,:)=0.0_rsh ; wave_uubot(:,:)=0.0_rsh ; wave_vubot(:,:)=0.0_rsh
   !!   - Wave-ocean momentum flux (m2.s-2) (due to breaking)
   ALLOCATE(wave_tauox(limin:limaxp1,ljmin:ljmax))
   ALLOCATE(wave_tauoy(limin:limax,ljmin:ljmaxp1))
   wave_tauox(:,:)=0.0_rsh
   wave_tauoy(:,:)=0.0_rsh
   !!   - Wave-ocean momentum flux (m2.s-2) (due to bottom friction)
   ALLOCATE(wave_taubbl_x(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_taubbl_y(limin:limax,ljmin:ljmax))
   wave_taubbl_x(:,:)=0.0_rsh
   wave_taubbl_y(:,:)=0.0_rsh
   !!   - Net wave-supported stress (m2.s-2)
   ALLOCATE(wave_tauwx(limin:limaxp1,ljmin:ljmax))
   ALLOCATE(wave_tauwy(limin:limax,ljmin:ljmaxp1))
   wave_tauwx(:,:)=0.0_rsh
   wave_tauwy(:,:)=0.0_rsh
   !!   - surface stokes drift (m2.s-1)
   ALLOCATE(wave_ussx(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_ussy(limin:limax,ljmin:ljmax))
   wave_ussx(:,:)=0.0_rsh
   wave_ussy(:,:)=0.0_rsh
   !!   - stokes transport (m2.s-1)
   ALLOCATE(wave_tusx(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tusy(limin:limax,ljmin:ljmax))
   wave_tusx(:,:)=0.0_rsh
   wave_tusy(:,:)=0.0_rsh
   !!   - Mean wave length (wlm in m)
   ALLOCATE(wave_wlm(limin:limax,ljmin:ljmax))
   wave_wlm(:,:)=0.0_rsh
   !!   - Mean wave period (tmn in s) (m0,-1)
   ALLOCATE(wave_t0m1(limin:limax,ljmin:ljmax))
   wave_t0m1(:,:)=0.0_rsh
   !!   - Mean wave direction (degrees)
   ALLOCATE(wave_dir(limin:limaxp1,ljmin:ljmaxp1))
   wave_dir(:,:)=0.0_rsh
   !!   - Wave-induced Bernoulli head pressure (bhd in N.m-1) (J term, Smith JPO 2006)
   ALLOCATE(wave_bhd(liminm2:limaxp2,ljminm2:ljmaxp2))
   wave_bhd(:,:)=0.0_rsh
   !!   - Wave-to-ocean TKE flux (phioc in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
   ALLOCATE(wave_phioc(limin:limax,ljmin:ljmax))
   wave_phioc(:,:)=0.0_rsh
   !!   - Wave-to-ocean TKE flux due to bottom friction (phibbl in m3 s-3) (i.e. W.m-2 / DWAT DWAT=1000)
   ALLOCATE(wave_phibbl(limin:limax,ljmin:ljmax))
   wave_phibbl(:,:)=0.0_rsh
   !!   - radiation stress 
   ALLOCATE(wave_sx(limin:limaxp1,ljmin:ljmax))
   ALLOCATE(wave_sy(limin:limax,ljmin:ljmaxp1))  ! x component
   ALLOCATE(wave_sxy(liminm1:limaxp1,ljminm1:ljmaxp1))  ! y component
   wave_sx(:,:)=0.0_rsh ; wave_sy(:,:)=0.0_rsh ; wave_sxy(:,:)=0.0_rsh

#ifdef key_wave_qeuler_dynamics
   ! 3D arrays
   ALLOCATE(wave_usz(1:kmax,liminm1:limax  ,ljmin  :ljmaxp1))
   ALLOCATE(wave_vsz(1:kmax,limin  :limaxp1,ljminm1:ljmax  ))
   ALLOCATE(wave_wsz(1:kmax,limin:limaxp1,ljmin:ljmaxp1))
   ALLOCATE(wave_botdiss_u(1:kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(wave_botdiss_v(1:kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(wave_tke(1:kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(peps(1:kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(taubbl(1:kmax,limin:limaxp1,ljmin:ljmaxp1))
   ALLOCATE(duzdsig(kmax,liminm1:limax,ljmin:ljmaxp1))
   ALLOCATE(dvzdsig(kmax,limin:limaxp1,ljminm1:ljmax))
   ALLOCATE(dsigdx_u(kmax,liminm1:limax,ljmin:ljmaxp1))
   ALLOCATE(dsigdy_v(kmax,limin:limaxp1,ljminm1:ljmax))
   wave_usz(:,:,:)=0.0_rsh
   wave_vsz(:,:,:)=0.0_rsh
   wave_wsz(:,:,:)=0.0_rsh
   wave_botdiss_u(:,:,:)=0.0_rsh
   wave_botdiss_v(:,:,:)=0.0_rsh
   wave_tke(:,:,:)=0.0_rsh
   peps(:,:,:)=0.0_rsh
   taubbl(:,:,:)=0.0_rsh
   duzdsig(:,:,:)=0.0_rsh
   dvzdsig(:,:,:)=0.0_rsh
   dsigdx_u(:,:,:)=0.0_rsh
   dsigdy_v(:,:,:)=0.0_rsh

   ! 2D arrays
   ALLOCATE(mdptc(liminm2:limaxp2,ljminm2:ljmaxp2))
   ALLOCATE(hrms(limin:limaxp1,ljmin:ljmaxp1))
   ALLOCATE(usz_zint(liminm1:limax,ljmin:ljmax))
   ALLOCATE(vsz_zint(limin:limax,ljminm1:ljmax))
   ALLOCATE(wave_number(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_km(limin:limax,ljmin:ljmax),   &
            wave_km_u(limin:limax,ljmin:ljmax), &
            wave_km_v(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_km_iter(limin:limax,ljmin:ljmax))
   ALLOCATE(wave_freq(limin:limax,ljmin:ljmax))
   ALLOCATE(tausurf(limin:limax,ljmin:ljmax))
   mdptc(:,:)=0.0_rsh
   hrms(:,:)=0.0_rsh
   wave_freq(:,:)=0.0_rsh
   tausurf(:,:)=0.0_rsh
   usz_zint(:,:)=0.0_rsh
   vsz_zint(:,:)=0.0_rsh
   wave_number(:,:)=0.0_rsh
   wave_km(:,:)=0.0_rsh
   wave_km_u(:,:)=0.0_rsh
   wave_km_v(:,:)=0.0_rsh
   wave_km_iter(:,:)=0

!#ifdef key_wave_soulsby_bottom_stress
   ALLOCATE(uorb(limin:limaxp1,ljmin:ljmaxp1))
   ALLOCATE(fwst(limin:limaxp1,ljmin:ljmaxp1))
   ALLOCATE(wave_uzk1(liminm1:limax  ,ljmin  :ljmaxp1))
   ALLOCATE(wave_vzk1(limin  :limaxp1,ljminm1:ljmax  ))
   uorb(:,:)=0.0_rsh
   fwst(:,:)=0.0_rsh
   wave_uzk1(:,:)=0.0_rsh
   wave_vzk1(:,:)=0.0_rsh
!#endif

#endif

  END SUBROUTINE wave_alloc

   !!======================================================================

  SUBROUTINE wave_dealloc

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dealloc_comwave  ***
   !&E
   !&E ** Purpose : Deallocate array use wave
   !&E
   !&E ** Called by :
   !&E       ! module dealloc (file dealloc.F90)
   !&E
   !&E ** History :
   !&E       !  2010-02 (G. Simon) Creation  on the basis of dealloc_cometdyn
   !&E       !  2013-08 (AC. Bennis) Introduction of arrays for 3D radiation stress tensor
   !&E       !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Review of the implementation of wave-current quasi-eulerian formalism
   !&E
   !&E---------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !! * Executable part

   !! * Locate MARS cell from Wave cell
   DEALLOCATE(xdpxe_wave,ydpxe_wave,xdpu_wave,ydpu_wave,xdpv_wave,ydpv_wave)
   DEALLOCATE(i1v_wave,j1v_wave,i1u_wave,j1u_wave,i1xe_wave,j1xe_wave)

#if ! defined key_wave_readww3txt && ! defined key_oasis_mars_ww3
   !! * Waves quantities -neww-
   !!   - directionality
   DEALLOCATE(wave_dir1,wave_dir2)
   !!   - significant height
   DEALLOCATE(wave_hs1,wave_hs2)
   !!   - turbulent kinetic energy flux from waves to ocean (breaking and bottom friction)
   DEALLOCATE(wave_foc1,wave_foc2,wave_fbb1,wave_fbb2)
   !!   - orbital velocity
   DEALLOCATE(wave_ubot1,wave_ubot2)
   DEALLOCATE(wave_uubot1,wave_uubot2)
   DEALLOCATE(wave_vubot1,wave_vubot2)
   !!   - wave period
   DEALLOCATE(wave_tp1,wave_tp2)
   DEALLOCATE(wave_tauox1,wave_tauoy1,wave_tauox2,wave_tauoy2)
   DEALLOCATE(wave_tusx1,wave_tusy1,wave_tusx2,wave_tusy2)
   DEALLOCATE(wave_ussx1,wave_ussy1,wave_ussx2,wave_ussy2)
   DEALLOCATE(wave_tauwx1,wave_tauwy1,wave_tauwx2,wave_tauwy2)
   DEALLOCATE(wave_t0m11,wave_t0m12)
   DEALLOCATE(wave_wlm1,wave_wlm2)
   DEALLOCATE(wave_bhd1,wave_bhd2)
   !!  - radiation stress
   DEALLOCATE(wave_sx1,wave_sx2)
   DEALLOCATE(wave_sy1,wave_sy2)
   DEALLOCATE(wave_sxy1,wave_sxy2)
#endif

   DEALLOCATE(wave_t0m1,wave_dir,wave_bhd,wave_tauox,wave_tauoy,wave_taubbl_x,wave_taubbl_y,wave_tauwx,wave_tauwy)
   DEALLOCATE(wave_hs,wave_ubot,wave_phioc,wave_phibbl,wave_wlm,wave_tusx,wave_tusy,wave_ussx,wave_ussy)
   DEALLOCATE(wave_uubot,wave_vubot,wave_tp)
   DEALLOCATE(wave_sx,wave_sy,wave_sxy)   

!#ifdef key_wave_soulsby_bottom_stress
   DEALLOCATE(uorb)
   DEALLOCATE(fwst)
   DEALLOCATE(wave_uzk1,wave_vzk1)
!#endif
#ifdef key_wave_qeuler_dynamics
   DEALLOCATE(duzdsig,dvzdsig)
   DEALLOCATE(dsigdx_u,dsigdy_v)
   DEALLOCATE(wave_botdiss_u,wave_botdiss_v)
   DEALLOCATE(wave_tke,peps,taubbl)
   DEALLOCATE(wave_freq)
   DEALLOCATE(hrms)
   DEALLOCATE(mdptc)
   DEALLOCATE(wave_usz,wave_vsz,wave_wsz)
   DEALLOCATE(wave_number,wave_km_u,wave_km_v,wave_km,wave_km_iter)
   DEALLOCATE(usz_zint,vsz_zint)
#endif

  END SUBROUTINE wave_dealloc

   !!======================================================================

#ifdef key_wave_readww3txt

  SUBROUTINE wave_read_ww3txt

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE wave_read_ww3txt  ***
   !&E
   !&E ** Purpose : lecture des champs de forcage venant de WW3
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  1987-02  (J.-C. Salomon, P. Lazure)  Original code
   !&E       !  1990-10  (F. Dumas ) new formulation
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters, ONLY : riosh,jshift_obc_north,jshift_obc_south

   !! * Arguments

   !! * Local declarations
   INTEGER ::  i,j
   REAL(KIND=riosh), DIMENSION(imin:imax,jmin:jmax) :: tmp_riosh

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter WAVE_READ_WW3TXT routine'

   tmp_riosh(:,:)=0.0_riosh

!variable: TMN (.t)
   OPEN(101,file='./../inputs/INPUTS_WW3/tmn.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(101,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(101)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_t0m1(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: HS (.hs)
   OPEN(102,file='./../inputs/INPUTS_WW3/hs.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(102,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(102)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_hs(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: THM (.dir)
   OPEN(103,file='./../inputs/INPUTS_WW3/thm.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(103,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(103)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_dir(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO
   CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,wave_dir(limin:limaxp1,ljmin:ljmaxp1))
   CALL_MPI ex_j_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,wave_dir(limin:limaxp1,ljmin:ljmaxp1))

!variable: J (.J)
   OPEN(104,file='./../inputs/INPUTS_WW3/jsmith.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(104,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(104)
   DO i=limin,limaxp1
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmaxp1,jmax-jshift_obc_north)
       wave_bhd(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: TAOUX  (.two) deferlement
   OPEN(105,file='./../inputs/INPUTS_WW3/tauox.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(105,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(105)
   DO i=limin,limaxp1
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_tauox(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: TAUOY (.two)  (champ nul dans cas test)
   OPEN(106,file='./../inputs/INPUTS_WW3/tauoy.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(106,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(106)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmaxp1,jmax-jshift_obc_north)
       wave_tauoy(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: TAUWX (voir avec F.Ardhuin si TAUWNX serait pas mieux...) vent
! (champ nul dans cas test)
   OPEN(107,file='./../inputs/INPUTS_WW3/tauwix.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(107,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(107)
   DO i=limin,limaxp1
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_tauwx(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: TAUWY (meme remarque) (champ nul dans cas test)
   OPEN(108,file='./../inputs/INPUTS_WW3/tauwiy.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(108,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(108)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmaxp1,jmax-jshift_obc_north)
       wave_tauwy(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: UBA (.ubr) (amplitude)
   OPEN(109,file='./../inputs/INPUTS_WW3/uba.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(109,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(109)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_ubot(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: PHIOC (.foc)
!dissipation liee au frott sur le fond
   OPEN(111,file='./../inputs/INPUTS_WW3/phioc.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(111,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(111)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_phioc(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

!variable: WLM (.l)
   OPEN(114,file='./../inputs/INPUTS_WW3/wlm.txt',status='old',form='formatted')
   DO i=imin+1,imax
     READ(114,*)(tmp_riosh(i,j),j=jmin,jmax-jshift_obc_north-jshift_obc_south)
   END DO
   CLOSE(114)
   DO i=limin,limax
     DO j=MAX0(ljmin,jmin+jshift_obc_south),MIN0(ljmax,jmax-jshift_obc_north)
       wave_wlm(i,j)=tmp_riosh(i,j-jshift_obc_south)
     END DO
   END DO

   PRINT_DBG*, 'exit WAVE_READ_WW3TXT routine'

  END SUBROUTINE wave_read_ww3txt

#endif

   !!======================================================================

  SUBROUTINE wavenumb(A,B)

   !&E-------------------------------------------------------------------------
   !&E                 ***  ROUTINE wavenum  ***
   !&E
   !&E ** Description : wave number computation
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E      - subroutine from WW3 (Tolman, 2009)
   !&E      - newton method (Mignot and Crouzeix, 1992)
   !&E
   !&E ** History :
   !&E       !  2013-08  (A-C. Bennis, F. Dumas ) creation - postdoc EPIGRAM
   !&E
   !&E-------------------------------------------------------------------------
   !! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh), dimension(limin:limax,ljmin:ljmax), intent(in) :: A,B
       ! A la periode moyenne des vagues
       ! B la profondeur moyenne

#if ! defined key_palm_ww3 && ! defined key_oasis_ww3

   !! * Local declarations
   REAL(KIND=rsh), dimension(limin:limax,ljmin:ljmax) :: KOLD,F,FD,DIF,RDIF,err,W
   REAL(KIND=rsh), PARAMETER :: EPS=1.0e-5
   INTEGER, PARAMETER :: mmax=1e7
   INTEGER :: m,i,j

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter WAVENUMB routine'
print*,'ATTENTION LA PARALLELISATION N A APS ETE NI TESTEE (bon a priori toutefois car tres simple) NI LA ROUTINE D AILLEURS'
   wave_number(:,:)=0.183_rsh
   KOLD(:,:)=0.0_rsh
   F(:,:)=0.0_rsh
   FD(:,:)=0.0_rsh
   DIF(:,:)=0.0_rsh
   RDIF(:,:)=0.0_rsh
   err(:,:)=0.0_rsh

   DO i=liminp1,limax
     DO j=ljminp1,ljmax

       IF ( A(i,j)*B(i,j) /= 0.0_rsh ) THEN
         W(i,j)=abs(2.0_rsh*pi/A(i,j))
         DO m=1,mmax
           DIF(i,j)=abs(wave_number(i,j)-KOLD(i,j))
           RDIF(i,j)=0.0_rsh
           IF (wave_number(i,j) /= 0.0_rsh) RDIF(i,j)=DIF(i,j)/wave_number(i,j)

           IF (DIF(i,j) < EPS .AND. RDIF(i,j) > EPS) THEN
             err(i,j)=abs(grav*wave_number(i,j)*TANH(wave_number(i,j)*B(i,j))-W(i,j)**2)
             exit
           ELSE
             KOLD(i,j)=wave_number(i,j)
             F(i,j)=grav*KOLD(i,j)*TANH(KOLD(i,j)*B(i,j))-W(i,j)**2
             IF (KOLD(i,j)*B(i,j) > 25.0_rsh) THEN
               FD(i,j)=grav*TANH(KOLD(i,j)*B(i,j))
             ELSE
               FD(i,j)=grav*TANH(KOLD(i,j)*B(i,j))  &
               +grav*KOLD(i,j)*B(i,j)/((COSH(KOLD(i,j)*B(i,j)))**2)
             END IF
             wave_number(i,j)=KOLD(i,j)-F(i,j)/FD(i,j)
           END IF

         END DO

       ELSE

         wave_number(i,j)=0.0_rsh

       END IF
     END DO
   END DO

#endif

  END SUBROUTINE wavenumb

   !!======================================================================

#endif

END MODULE wave

