!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#if defined  WAVE_OFFLINE && (defined MUSTANG || defined BBL || defined MRL_WCI)

                              ! Read point or gridded wind wave
      subroutine get_wwave    ! height, direction, period at approp.
                              ! time from forcing NetCDF file.
# ifdef OPENACC
      use exchange_device_host
# endif
      implicit none
# include "param.h"
# include "scalars.h"
# include "ncscrum.h"
# include "forces.h"
      real cff
      integer i,ierr, lstr,lvar, lenstr, nf_fread, advance_cycle
# include "netcdf.inc"
!
! Initialization: Inquire about the contents of forcing NetCDF file:
!================  variables and dimensions. Check for consistency.
!


# ifndef MUSTANG
      ncidwave=ncidfrc
      wave_file=frcname
# endif

      if (may_day_flag.ne.0) return      !-->  EXIT
      if (itww.eq.0 .or. iic.eq.0) then
        lstr=lenstr(wave_file)
        if (may_day_flag.ne.0) return    !-->  EXIT
!
! If not opened yet, open forcing NetCDF file for reading.
! Find and save IDs for relevant variables, determine whether
! wind wave data are field or scalar values.
!

        if (ncidwave.eq.-1) then
          ierr=nf_open(wave_file(1:lstr), nf_nowrite, ncidwave)
          if (ierr .ne. nf_noerr) goto 4                !--> ERROR
        endif

        ierr=nf_inq_varid (ncidwave, 'wwv_time', ww_tid)
        if (ierr .ne. nf_noerr) then
          write(stdout,3) 'wwv_time', wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif

# ifndef MUSTANG
        ncidfrc=ncidwave
# endif

! Check if wave amplitude present, and if so, if field or scalar series

        lvar=lenstr(vname(1,indxWWA))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWWA)(1:lvar),wwa_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wwa_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwagrd=1.
            else
              wwagrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWWA)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif

! Check if wave direction present, and if so, if field or scalar series

        lvar=lenstr(vname(1,indxWWD))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWWD)(1:lvar),wwd_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wwd_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwdgrd=1.
            else
              wwdgrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWWD)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif


! Check if wave period present, and if so, if field or scalar series

        lvar=lenstr(vname(1,indxWWP))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWWP)(1:lvar),wwp_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wwp_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwpgrd=1.
            else
              wwpgrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWWP)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif

# if defined MUSTANG
! Check if wave orbital velocity present, and if so, if field or scalar series

        lvar=lenstr(vname(1,indxWWU))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWWU)(1:lvar),wwu_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wwu_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwugrd=1.
            else
              wwugrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWWU)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif

# endif

! Check if breaking dissipation present, and if so, if field or scalar series

# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
        lvar=lenstr(vname(1,indxWEB))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWEB)(1:lvar),wweb_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wweb_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwebgrd=1.
            else
              wwebgrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWEB)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif
# endif

! Check if frictional dissipation present, and if so, if field or scalar series

# if defined MRL_WCI && defined WAVE_OFFLINE_FRICTION
        lvar=lenstr(vname(1,indxWED))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWED)(1:lvar),wwed_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wwed_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwedgrd=1.
            else
              wwedgrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWED)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif
# endif

! Check if roller dissipation present, and if so, if field or scalar series

# if defined MRL_WCI && defined WAVE_OFFLINE_ROLLER
        lvar=lenstr(vname(1,indxWER))
        ierr=nf_inq_varid (ncidwave,vname(1,indxWER)(1:lvar),wwer_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidwave, wwer_id, i)
          if (ierr. eq. nf_noerr) then
            if (i.gt.1) then
              wwergrd=1.
            else
              wwergrd=0.
            endif
          endif
        endif
        if (ierr .ne. nf_noerr) then
          write(stdout,3) vname(1,indxWER)(1:lvar),wave_file(1:lstr)
          goto 99                                         !--> ERROR
        endif
# endif

!
! Determine whether cycling is defined to reuse the input data
! and find cycling period "ww_cycle", set initial cycling
! index "ww_ncycle" and record index "ww_rec".
! Set initial value for time index "itww" and both time record
! bounds to large negative artificial values, so that it will
! trigger the logic in the reading part below.
! Also set scale factor to convert input data to model units:
! convert wave direction from degrees to radians;
! no conversion for amplitude and period.
!
        call set_cycle (ncidwave,   ww_tid,    ntww,
     &                  ww_cycle, ww_ncycle, ww_rec)
        if (may_day_flag.ne.0) return    !-->  EXIT

        itww=2
        wwv_time(1)=-1.E+20
        wwv_time(2)=-1.E+20


        wwa_scale=1.
        wwd_scale=deg2rad
        wwp_scale=1.
        wweb_scale=1.
        wwed_scale=1.
        wwer_scale=1.
# ifdef MUSTANG
        wwu_scale=1.
# endif

      endif                   !<-- itww.eq.0 .or. iic.eq.0

!
! Reading data from the forcing file:  Get out, if model time is
!======== ==== ==== === ======= =====  already within the interval
! set by the past and future data times. Otherwise flip the time
! index, increment record and cyclin indices and read a new portion
! of data. Repeat it until model time is between the two times from
! data.
!
  1   i=3-itww
      cff=time+0.5*dt
      if (wwv_time(i).le.cff .and. cff.lt.wwv_time(itww)) return
      ierr=advance_cycle (ww_cycle, ntww, ww_ncycle, ww_rec)
      if (ierr.ne.0) goto 7                           !--> ERROR

      ierr=nf_get_var1_FTYPE(ncidwave, ww_tid, ww_rec, cff)
      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wwv_time', ww_rec
        goto 99                                         !--> ERROR
      endif
# ifdef USE_CALENDAR
        call tool_origindate(ncidwave,ww_tid,
     &                       ww_origin_date_in_sec)
        cff=cff+ww_origin_date_in_sec*sec2day
# endif
      wwv_time(i)=cff*day2sec+ww_cycle*ww_ncycle
      if (wwv_time(itww).eq.-1.E+20) wwv_time(itww)=wwv_time(i)

! Read wave amplitude

      if (wwagrd.eq.1.) then
        ierr=nf_fread (wwag(START_2D_ARRAY,i), ncidwave, wwa_id,
     &                                           ww_rec, r2dvar)
# ifdef OPENACC
        call update_device(wwag(:,:,i))
# endif
      else
        ierr=nf_get_var1_FTYPE (ncidwave,wwa_id,ww_rec,wwap(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wind-wave amplitude', ww_rec
        goto 99                                         !--> ERROR
      endif

! Read wave direction

      if (wwdgrd.eq.1.) then
        ierr=nf_fread (wwdg(START_2D_ARRAY,i), ncidwave, wwd_id,
     &                                           ww_rec, r2dvar)
# ifdef OPENACC
        call update_device(wwdg(:,:,i))
# endif
      else
        ierr=nf_get_var1_FTYPE (ncidwave,wwd_id,ww_rec,wwdp(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wind-wave direction', ww_rec
        goto 99                                         !--> ERROR
      endif

! Read wave period

      if (wwpgrd.eq.1.) then
        ierr=nf_fread (wwpg(START_2D_ARRAY,i), ncidwave, wwp_id,
     &                                           ww_rec, r2dvar)
# ifdef OPENACC
        call update_device(wwpg(:,:,i))
# endif
      else
        ierr=nf_get_var1_FTYPE (ncidwave,wwp_id,ww_rec,wwpp(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wind-wave period', ww_rec
        goto 99                                         !--> ERROR
      endif

! Read orbital velocity

# ifdef MUSTANG
      if (wwugrd.eq.1.) then
        ierr=nf_fread (wwug(START_2D_ARRAY,i), ncidwave, wwu_id,
     &                                           ww_rec, r2dvar)
#  ifdef OPENACC
        call update_device(wwug(:,:,i))
#  endif
      else
        ierr=nf_get_var1_FTYPE (ncidwave,wwu_id,ww_rec,wwup(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wave orbital velocity', ww_rec
        goto 99                                         !--> ERROR
      endif
# endif

! Read breaking dissipation

# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
      if (wwebgrd.eq.1.) then
        ierr=nf_fread (wweb(START_2D_ARRAY,i), ncidwave, wweb_id,
     &                                           ww_rec, r2dvar)
#  ifdef OPENACC
        call update_device(wweb(:,:,i))
#  endif
      else
        ierr=nf_get_var1_FTYPE (ncidwave,wweb_id,ww_rec,wwebp(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wave breaking dissipation', ww_rec
        goto 99                                         !--> ERROR
      endif
# endif

! Read frictional dissipation

# if defined MRL_WCI && defined WAVE_OFFLINE_FRICTION
      if (wwedgrd.eq.1.) then
        ierr=nf_fread (wwedd(START_2D_ARRAY,i), ncidwave, wwed_id,
     &                                           ww_rec, r2dvar)
#  ifdef OPENACC
        call update_device(wwedd(:,:,i))
#  endif
     else
        ierr=nf_get_var1_FTYPE (ncidwave,wwed_id,ww_rec,wwedp(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wave frictional dissipation', ww_rec
        goto 99                                         !--> ERROR
      endif
# endif

! Read frictional dissipation

# if defined MRL_WCI && defined WAVE_OFFLINE_ROLLER
      if (wwergrd.eq.1.) then
        ierr=nf_fread (wwer(START_2D_ARRAY,i), ncidwave, wwer_id,
     &                                           ww_rec, r2dvar)
#  ifdef OPENACC
        call update_device(wwer(:,:,i))
#  endif
      else
        ierr=nf_get_var1_FTYPE (ncidwave,wwer_id,ww_rec,wwerp(i))
      endif

      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wave roller dissipation', ww_rec
        goto 99                                         !--> ERROR
      endif
# endif


      itww=i
      MPI_master_only write(stdout,'(6x,A,1x,A,1x,F10.4,1x,I4)')
     &  'GET_WWAVE --','Read wind wave input data for time =', cff
# ifdef USE_CALENDAR
     &  -ww_origin_date_in_sec*sec2day
# endif
# ifdef MPI
     &                                                      , mynode
# endif
      if (ntww.gt.1) goto 1
      if (ntww.eq.1) return
!
!
! Sort out error messages: The following portion of the code is
!===== === ===== ========= not accessed unless something goes wrong.
!
  3   format(/,' GET_WWAVE - unable to find forcing variable: ',a,
     &                          /,15x,'in forcing NetCDF file: ',a)
  4   write(stdout,5) wave_file(1:lstr)
  5   format(/,' GET_WWAVE - unable to open forcing NetCDF file:',
     &                                                        1x,a)
      goto 99

  6   format(/,' GET_WWAVE - error while reading variable: ',a,2x,
     &                                      ' at TIME index = ',i4)
  7   write(stdout,8) ww_rec, ntww, wave_file(1:lstr), tdays,
     &                                wwv_time(itww)*sec2day
# ifdef USE_CALENDAR
     &  -ww_origin_date_in_sec*sec2day
# endif
  8   format(/,' GET_WWAVE - ERROR: requested time record ',I4,
     &       1x,'exceeds the last available', /,14x,'record ',I4,
     &       1x,'in forcing NetCDF file: ',a, /,14x,'TDAYS = ',
     &       g12.4,2x,'last available wwv_time = ',g12.4)

  99  may_day_flag=2

      return
      end
!
!=====================================================================
!
!                             SET_WWAVE
!
!=====================================================================
!
      subroutine set_wwave (tile)
      implicit none
# include "param.h"
      integer tile
# include "compute_tile_bounds.h"
!
      call set_wwave_tile (Istr,Iend,Jstr,Jend)
      return
      end

      subroutine set_wwave_tile (Istr,Iend,Jstr,Jend)
!
!--------------------------------------------------------------------
!  Set-up wave amplitude, direction, period for current tile.
!--------------------------------------------------------------------
!
      implicit none

      integer Istr,Iend,Jstr,Jend, i,j,k, it1,it2
      real    fac,fac1,fac2,fac11,fac12,fac21,fac22,fac31,fac32,
     &                      fac41,fac42,fac51,fac52,fac61,fac62,
     &        wwA, wwD, wwP, wwe, Dstp,cfrq,hrms
# ifdef MUSTANG
     &       , wwU, fac71, fac72
# endif
      real    cff,cff1,cff2, kh,khd,kw, eps, wramp,m2p,isr2
# ifdef BBL_OFFLINE
     &                                            , ab, orb_vel
# endif
      parameter (eps = 1.0d-10, m2p=1.2804973111, isr2=0.70710678)
!
# include "param.h"
# include "grid.h"
# include "forces.h"
# include "scalars.h"
# include "ocean2d.h"
# include "ocean3d.h"
!
# include "compute_extended_bounds.h"
!
! Set coefficients for interpolation. Check that for the next time
! step [when time=time+dt] both weights will still be positive, and
! if not, set synchro_flag to signal that new data should be read
! from an appropriate netCDF input file (master thread only).
! After that, either load time-invariant data, or interpolate in time.
! Complain about error and signal to quit, if interpolation is
! needed, but not possible.
!
!
      it1=3-itww
      it2=itww
      fac=time+0.5*dt
      fac1=wwv_time(it2)-fac
      fac2=fac-wwv_time(it1)
!
! Load time-invariant wind-wave data.
! Time interpolation is not performed in this case.
!
      if (ww_cycle.lt.0.) then
        if (FIRST_TIME_STEP) then

            wwA=wwa_scale*wwap(itww)
            wwD=wwd_scale*wwdp(itww)
            wwP=wwp_scale*wwpp(itww)
# ifdef MUSTANG
            wwU=wwu_scale*wwup(itww)
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
            wwEpb=wweb_scale*wwebp(itww)
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_FRICTION
            wwEpd=wwed_scale*wwedp(itww)
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_ROLLER
            wwEpr=wwer_scale*wwerp(itww)
# endif

!$acc parallel loop if(compute_on_device) default(present)          
            do j=JstrR,JendR
              do i=IstrR,IendR
                Awave(i,j)=(wwagrd)*wwa_scale*wwag(i,j,itww) +
     &                     (1.-wwagrd)*wwA
                Dwave(i,j)=(wwdgrd)*wwd_scale*wwdg(i,j,itww) +
     &                     (1.-wwdgrd)*wwD

                Pwave(i,j)=(wwpgrd)*wwp_scale*wwpg(i,j,itww) +
     &                     (1.-wwpgrd)*wwP
# ifdef MUSTANG
                Uwave(i,j)=(wwugrd)*wwu_scale*wwug(i,j,itww) +
     &                     (1.-wwugrd)*wwU
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
                Ebwave(i,j)=(wwebgrd)*wweb_scale*wweb(i,j,itww) +
     &                     (1.-wwebgrd)*wwEpb
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_FRICTION
                Edwave(i,j)=(wwedgrd)*wwed_scale*wwed(i,j,itww) +
     &                     (1.-wwedgrd)*wwEpd
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_ROLLER
                Erwave(i,j)=(wwergrd)*wwer_scale*wwer(i,j,itww) +
     &                     (1.-wwergrd)*wwEpr
# endif
              enddo
            enddo
        endif
!
! Time-interpolate wave data from gridded or point data.
! Make sure that for the next time step [when time=time+dt]
! time+dt is still between wwv_time(it1) and wwv_time(it2);
! and if not, set synchro_flag to signal that the new forcing data
! should be read from the netCDF input file (master thread only).
!
      elseif (fac1.ge.0. .and. fac2.ge.0.) then
       if (ZEROTH_TILE .and. fac1.lt.dt) synchro_flag=.TRUE.

        fac=wwa_scale/(fac1+fac2)
        fac11=fac*fac1
        fac12=fac*fac2
        wwA=fac11*wwap(it1)+fac12*wwap(it2)

        fac=wwd_scale/(fac1+fac2)
        fac21=fac*fac1
        fac22=fac*fac2
        wwD=fac21*wwdp(it1)+fac22*wwdp(it2)

        fac=wwp_scale/(fac1+fac2)
        fac31=fac*fac1
        fac32=fac*fac2
        wwP=fac31*wwpp(it1)+fac32*wwpp(it2)
# ifdef MUSTANG
        fac=wwu_scale/(fac1+fac2)
        fac71=fac*fac1
        fac72=fac*fac2
        wwU=fac71*wwup(it1)+fac72*wwup(it2)
# endif

# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
        fac=wweb_scale/(fac1+fac2)
        fac41=fac*fac1
        fac42=fac*fac2
        wwEpb=fac41*wwebp(it1)+fac42*wwebp(it2)
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_FRICTION
        fac=wwed_scale/(fac1+fac2)
        fac51=fac*fac1
        fac52=fac*fac2
        wwEpd=fac51*wwedp(it1)+fac52*wwedp(it2)
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_ROLLER
        fac=wwer_scale/(fac1+fac2)
        fac61=fac*fac1
        fac62=fac*fac2
        wwEpr=fac61*wwerp(it1)+fac62*wwerp(it2)
# endif
        
!$acc parallel loop if(compute_on_device) default(present)
        do j=JstrR,JendR
          do i=IstrR,IendR
            Awave(i,j)=wwagrd*(fac11*wwag(i,j,it1)+fac12*wwag(i,j,it2))
     &                +(1.-wwagrd)*wwA
            Dwave(i,j)=wwdgrd*(fac21*wwdg(i,j,it1)+fac22*wwdg(i,j,it2))
     &                +(1.-wwdgrd)*wwD

            Pwave(i,j)=wwpgrd*(fac31*wwpg(i,j,it1)+fac32*wwpg(i,j,it2))
     &                +(1.-wwpgrd)*wwP
# ifdef MUSTANG
            Uwave(i,j)=wwugrd*(fac71*wwug(i,j,it1)+fac72*wwug(i,j,it2))
     &                +(1.-wwugrd)*wwU
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
            wepb(i,j)=wwebgrd*(fac41*wweb(i,j,it1)+fac42*wweb(i,j,it2))
     &                +(1.-wwebgrd)*wwEpb
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_FRICTION
            wepd(i,j)=wwedgrd*(fac51*wwed(i,j,it1)+fac52*wwed(i,j,it2))
     &                +(1.-wwedgrd)*wwEpd
# endif
# if defined MRL_WCI && defined WAVE_OFFLINE_ROLLER
            wepr(i,j)=wwergrd*(fac61*wwer(i,j,it1)+fac62*wwer(i,j,it2))
     &                +(1.-wwergrd)*wwEpr
# endif
          enddo
        enddo
      endif
!
!--------------------------------------------------------------------
!  Convert wave direction Dwave to Cartesian convention
!--------------------------------------------------------------------
!
!  Wave direction conventions:
!  - Cartesian convention (CROCO):
!          the direction to where the vector points,
!          measured counterclockwise from the positive x-axis
!          of the model grid (must account for grid angle)
!  - Meteorological convention (i.e., ECMWF data ...):
!          the direction where the waves come from,
!          measured clockwise from geographic North.
!  - Oceanographic convention:
!          the direction where the waves travel to,
!          measured clockwise from geographic North.
!
# if defined CURVGRID || !defined SHOREFACE
!$acc parallel loop if(compute_on_device) default(present)
      do j=jstrR,jendR
        do i=istrR,iendR
#  if !defined SHOREFACE && !defined MUSTANG
!         Dwave(i,j)=0.5*pi+Dwave(i,j) ! Dwave=0 when waves to N (oceano)
          Dwave(i,j)=1.5*pi-Dwave(i,j) ! Dwave=0 when waves from N (meteo)
#  endif
#  ifdef CURVGRID
#   ifdef MUSTANG
          Dwave(i,j)=Dwave(i,j)-angler(i,j)*rad2deg
#   else
          Dwave(i,j)=Dwave(i,j)-angler(i,j)
#   endif
#  endif
        enddo
      enddo
# endif /* CURVGRID || !SHOREFACE */
!
!--------------------------------------------------------------------
! Simulate nearshore modifications for low resolution data
!
!   Mekong application provides an example that can be used
!   for other muddy areas (still needs work).
!--------------------------------------------------------------------
!
# if !defined MRL_WCI || defined MEKONG
!$acc parallel loop if(compute_on_device) default(present)
      do j=jstrR,jendR
        do i=istrR,iendR
          Dstp=max(h(i,j)+z_w(i,j,N),1.e-3)
#  ifdef MEKONG
!
! First estimate wavenumber
!
          hrms=2*Awave(i,j)
          cfrq=2.0*pi/max(Pwave(i,j),0.1D0)
          khd=Dstp*cfrq*cfrq/g
          kh = sqrt(    khd*khd + khd/(1.0 + khd*(0.6666666666
     &                  +khd*(0.3555555555 + khd*(0.1608465608
     &                  +khd*(0.0632098765 + khd*(0.0217540484
     &                                  +khd*0.0065407983)))))) )
          kw=kh/Dstp
!
! Strong mudbank dissipation formulation: frequency and depth-dependent
! dissipation rate is a function akin to that proposed by Bretschneider
! and Reid (1954) and Grosskopf (1980): longer and higher waves are more
! highly dissipated.
!
          Awave(i,j)=Awave(i,j)*
     &                (.5*(1-tanh(1.-2.*kw**.9*Dstp/hrms**0.2)))**2
#  else
!
! Assume surfzone dissipation due to depth-induced wave breaking
!
          Awave(i,j)=min(0.275*Dstp, Awave(i,j))
#  endif
        enddo
      enddo
# endif /* !MRL_WCI || MEKONG */
!
!--------------------------------------------------------------------
! Compute variables for wave-current interactions (mrl_wci routine)
! --> store in same arrays as in OW_COUPLING case
!--------------------------------------------------------------------
!
# ifdef MRL_WCI
!$acc parallel loop if(compute_on_device) default(present)
      do j=jstrR,jendR
        do i=istrR,iendR
          whrm(i,j)=2*Awave(i,j)                  ! Amp to Hrms conversion
          wfrq(i,j)=2.0*pi/max(Pwave(i,j),0.1D0)  ! Pwave peak period
          wdrx(i,j)=cos(Dwave(i,j))
          wdre(i,j)=sin(Dwave(i,j))
#  ifdef MASKING
          wfrq(i,j)=wfrq(i,j)*rmask(i,j)
          whrm(i,j)=whrm(i,j)*rmask(i,j)
          wdrx(i,j)=wdrx(i,j)*rmask(i,j)
          wdre(i,j)=wdre(i,j)*rmask(i,j)
#   ifdef WAVE_OFFLINE_BREAKING
          wepb(i,j)=wepb(i,j)*rmask(i,j)
#   endif
#   ifdef WAVE_OFFLINE_FRICTION
          wepd(i,j)=wepd(i,j)*rmask(i,j)
#   endif
#   ifdef WAVE_OFFLINE_ROLLER
          wepr(i,j)=wepr(i,j)*rmask(i,j)
#   endif
#  endif
        enddo
      enddo
# endif /* MRL_WCI */
!
!--------------------------------------------------------------------
! Unable to set-up wave fields:
! Complain about the error and signal to quit (ONE THREAD ONLY).
!--------------------------------------------------------------------
!
      if ((ww_cycle.gt.0.) .and. (fac1.lt.0..or.fac2.lt.0.)) then
        if (ZEROTH_TILE) then
          write(stdout,1) 'wwv_time', tdays, wwv_time(it2)*sec2day
# ifdef USE_CALENDAR
     &  -ww_origin_date_in_sec*sec2day
# endif
  1       format(/,' SET_WWAVE_TILE - current model time exceeds',
     &             ' ending value for variable: ',a,/,14x,'TDAYS = ',
     &                                      g12.4,2x,'TEND = ',g12.4)
          may_day_flag=2
        endif
      endif
      return
      end
#else
      subroutine get_wwave_empty
      return
      end
#endif /* BBL && !ANA_WWAVE */
