!======================================================================
! 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 ANA_BRY_WKB && defined WKB_WWAVE

      subroutine get_bry_wkb        ! Read side boundary forcing
      implicit none                 ! fields from boundary file

# define wac_west 'cdir_west'
# include "param.h"
# include "scalars.h"
# include "ncscrum.h"
# include "netcdf.inc"
# include "boundary.h"
      real cff
      integer  lstr,lvar,lenstr, ierr, ierr_all
      integer  nf_read_bry_EW, nf_read_bry_NS
!
! Initialization: Check, whether boundary forcing filefor WKB model is
!================ already opened, and if not, open it.  Find and save
! netCDF IDs for relevant variables.  Determine whether there is cycling
! to reuse the input data and find cycling period "brywkb_cycle", set
! initial cycling index "brywkb_ncycle" and record index "brywkb_rec",
! time index "itbrywkb" and both time record bounds to large negative
! artificial values, so that it will trigger the logic in reading
! part below.
! Note that the frequency may be different than for other bry files
!
      ierr=nf_noerr
      lstr=lenstr(brywkb_file)
      if (iic.eq.0 ) then
        if (brywkb_id.eq.-1) then
          ierr=nf_open (brywkb_file(1:lstr), nf_nowrite, brywkb_id)
          if (ierr.ne.nf_noerr) write(stdout,'(/1x,4A/)') 'ERROR ',
     &               'in get_bry_wkb: can not open netCDF file ''',
     &                                   brywkb_file(1:lstr),   '''.'
        endif
        if (ierr.eq.nf_noerr) then
          ierr_all=0
          ierr=nf_inq_varid (brywkb_id, 'brywkb_time',  brywkb_time_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'brywkb_time'
          ierr_all=ierr_all+ierr

# ifdef WKB_OBC_WEST
          ierr=nf_inq_varid (brywkb_id, 'wac_west', wacbry_west_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wac_west'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wkx_west', wkxbry_west_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wkx_west'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wke_west', wkebry_west_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wke_west'
          ierr_all=ierr_all+ierr
# endif
# ifdef WKB_OBC_EAST
          ierr=nf_inq_varid (brywkb_id, 'wac_east', wacbry_east_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wac_east'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wkx_east', wkxbry_east_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wkx_east'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wke_east', wkebry_east_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wke_east'
          ierr_all=ierr_all+ierr
# endif
# ifdef WKB_OBC_SOUTH
          ierr=nf_inq_varid (brywkb_id, 'wac_south', wacbry_south_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wac_south'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wkx_south', wkxbry_south_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wkx_south'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wke_south', wkebry_south_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wke_south'
          ierr_all=ierr_all+ierr
# endif
# ifdef WKB_OBC_NORTH
          ierr=nf_inq_varid (brywkb_id, 'wac_north', wacbry_north_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wac_north'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wkx_north', wkxbry_north_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wkx_north'
          ierr_all=ierr_all+ierr
          ierr=nf_inq_varid (brywkb_id, 'wke_north', wkebry_north_id)
          if (ierr.ne.nf_noerr) write(stdout,1) 'wke_north'
          ierr_all=ierr_all+ierr
# endif
          ierr=ierr_all
          if (ierr.eq.nf_noerr) then
            call set_cycle (brywkb_id, brywkb_time_id, ntbrywkb,
     &           brywkb_cycle, brywkb_ncycle, brywkb_rec)
            itbrywkb=1
            brywkb_time(1)=-1.E+20
            brywkb_time(2)=-1.E+20
          else
            write(stdout,'(8x,4A)') 'ERROR(s) occur while examining',
     &          ' content of netCDF file ''', brywkb_file(1:lstr), '''.'
          endif
        endif
      endif
  1   format(' ERROR in get_bry_wkb: cannot find variable ''',A,'''')
!
! Read data from the file:  Check if model time is bounded by past
!===== ==== ==== === =====  and future data times: if not, increment
! record and cycling indices, flip time index and read a new portion
! of data. Repeat until model time falls between the two data times.
!
      do while (brywkb_time(itbrywkb).lt.time+.5*dt .and. ierr.eq.nf_noerr)
        call advance_cycle (brywkb_cycle,ntbrywkb,brywkb_ncycle,brywkb_rec)
        if (ierr.eq.nf_noerr) then
          ierr=nf_get_var1_FTYPE (brywkb_id, brywkb_time_id, brywkb_rec, cff)
          if (ierr.eq.nf_noerr) then
# ifdef USE_CALENDAR
            call tool_origindate(brywkb_id,brywkb_time_id,
     &                         brywkb_origin_date_in_sec)
            cff=cff+brywkb_origin_date_in_sec*sec2day
# endif
            itbrywkb=min(3-itbrywkb,ntbrywkb)
            brywkb_time(itbrywkb)=cff*day2sec + brywkb_cycle*brywkb_ncycle
            ierr_all=0
# ifdef WKB_OBC_WEST
            ierr=nf_read_bry_EW (wacbry_west_dt(START_1D_ARRAYETA,itbrywkb
     &          ), brywkb_id,wacbry_west_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wac_west'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_EW (wkxbry_west_dt(START_1D_ARRAYETA,itbrywkb
     &          ), brywkb_id,wkxbry_west_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wkx_west'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_EW (wkebry_west_dt(START_1D_ARRAYETA,itbrywkb
     &          ), brywkb_id,wkebry_west_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wke_west'
            ierr_all=ierr_all+ierr
# endif
# ifdef WKB_OBC_EAST
            ierr=nf_read_bry_EW (wacbry_east_dt(START_1D_ARRAYETA,itbrywkb
     &          ), brywkb_id,wacbry_east_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wac_east'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_EW (wkxbry_east_dt(START_1D_ARRAYETA,itbrywkb
     &          ), brywkb_id,wkxbry_east_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wkx_east'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_EW (wkebry_east_dt(START_1D_ARRAYETA,itbrywkb
     &          ), brywkb_id,wkebry_east_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wke_east'
            ierr_all=ierr_all+ierr
# endif
# ifdef WKB_OBC_SOUTH
            ierr=nf_read_bry_NS (wacbry_south_dt(START_1D_ARRAYXI,itbrywkb
     &          ), brywkb_id,wacbry_south_id, brywkb_rec,  r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wac_south'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_NS (wkxbry_south_dt(START_1D_ARRAYXI,itbrywkb
     &          ), brywkb_id,wkxbry_south_id, brywkb_rec,  r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wkx_south'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_NS (wkebry_south_dt(START_1D_ARRAYXI,itbrywkb
     &          ), brywkb_id,wkebry_south_id, brywkb_rec,  r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wke_south'
            ierr_all=ierr_all+ierr
# endif
# ifdef WKB_OBC_NORTH
            ierr=nf_read_bry_NS (wacbry_north_dt(START_1D_ARRAYXI,itbrywkb
     &          ), brywkb_id,wacbry_north_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wac_north'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_NS (wkxbry_north_dt(START_1D_ARRAYXI,itbrywkb
     &          ), brywkb_id,wkxbry_north_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wkx_north'
            ierr_all=ierr_all+ierr
            ierr=nf_read_bry_NS (wkebry_north_dt(START_1D_ARRAYXI,itbrywkb
     &          ), brywkb_id,wkebry_north_id, brywkb_rec, r2dvar)
            if (ierr.ne.nf_noerr) write(stdout,2) 'wke_north'
            ierr_all=ierr_all+ierr
# endif
            ierr=ierr_all
            if (ierr.eq.0) then
              write(stdout,'(6x,A,9x,A,1x,G12.4,1x,I4)')
     &       'GET_BRY_WKB - Read wave boundary data',
     &                                  'for time =', cff
# ifdef USE_CALENDAR
     &  -brywkb_origin_date_in_sec*sec2day
# endif
# ifdef MPI
     &                                              , mynode
# endif
            else
              write(stdout,'(1x,2A,I4/8x,3A)')
     &           'ERROR(s) occur while trying to read record ',
     &            brywkb_rec,  'in file ''',brywkb_file(1:lstr),'''.'
            endif
          else
            write(stdout,2) 'brywkb_time'
          endif
        else
           write(stdout,'(/1x,A,I4,1x,A,I4/7x,4A/7x,2(A,G12.4)/)')
     &    'ERROR in get_bry_all: requested time record ',  brywkb_rec,
     &    'exeeds the last record', ntbrywkb,  'available in netCDF ',
     &    'file ''',   brywkb_file(1:lstr),  '''',  'tdays = ',tdays,
     &    '  but the last available  brywkb_time =',
     &                                   brywkb_time(itbrywkb)*sec2day
# ifdef USE_CALENDAR
     &  -brywkb_origin_date_in_sec*sec2day
# endif
        endif
      enddo
  2   format(' ERROR in get_bry: cannot read variable ''',A,'''')
      if (ierr.ne.nf_noerr) may_day_flag=2
      return
      end

      subroutine set_bry_wkb (tile)
      implicit none
      integer tile
# include "param.h"
# include "compute_tile_bounds.h"
      call set_bry_wkb_tile (Istr,Iend,Jstr,Jend)
      return
      end


      subroutine set_bry_wkb_tile (Istr,Iend,Jstr,Jend)
!
!--------------------------------------------------------------------
!  Set-up all boundary forcing fields
!--------------------------------------------------------------------
!
      implicit none
      integer Istr,Iend,Jstr,Jend, i,j, it1,it2
      real cff, cff1,cff2
# include "param.h"
# include "scalars.h"
# include "boundary.h"
!
# include "compute_extended_bounds.h"

      it1=3-itbrywkb
      it2=itbrywkb
      cff=time+0.5*dt
      cff1=brywkb_time(it2)-cff
      cff2=cff-brywkb_time(it1)
      if (ZEROTH_TILE .and. cff1.lt.dt) synchro_flag=.true.

      if (cff1.ge.0. .and. cff2.ge.0.) then
        cff=1./(cff1+cff2)                        ! interpolate
        cff1=cff1*cff                             ! boundary values
        cff2=cff2*cff                             ! in time
# ifdef WKB_OBC_WEST
        if (WESTERN_EDGE) then
          do j=JstrR,JendR
            wacbry_west(j)=cff1*wacbry_west_dt(j,it1)
     &                    +cff2*wacbry_west_dt(j,it2)
            wkxbry_west(j)=cff1*wkxbry_west_dt(j,it1)
     &                    +cff2*wkxbry_west_dt(j,it2)
            wkebry_west(j)=cff1*wkebry_west_dt(j,it1)
     &                    +cff2*wkebry_west_dt(j,it2)
          enddo
        endif
# endif
# ifdef WKB_OBC_EAST
        if (EASTERN_EDGE) then
          do j=JstrR,JendR
            wacbry_east(j)=cff1*wacbry_east_dt(j,it1)
     &                    +cff2*wacbry_east_dt(j,it2)
            wkxbry_east(j)=cff1*wkxbry_east_dt(j,it1)
     &                    +cff2*wkxbry_east_dt(j,it2)
            wkebry_east(j)=cff1*wkebry_east_dt(j,it1)
     &                    +cff2*wkebry_east_dt(j,it2)
          enddo
        endif
# endif
# ifdef WKB_OBC_SOUTH
        if (SOUTHERN_EDGE) then
          do i=IstrR,IendR
            wacbry_south(i)=cff1*wacbry_south_dt(i,it1)
     &                     +cff2*wacbry_south_dt(i,it2)
            wkxbry_south(i)=cff1*wkxbry_south_dt(i,it1)
     &                     +cff2*wkxbry_south_dt(i,it2)
            wkebry_south(i)=cff1*wkebry_south_dt(i,it1)
     &                     +cff2*wkebry_south_dt(i,it2)
          enddo
        endif
# endif
# ifdef WKB_OBC_NORTH
        if (NORTHERN_EDGE) then
          do i=IstrR,IendR
            wacbry_north(i)=cff1*wacbry_north_dt(i,it1)
     &                     +cff2*wacbry_north_dt(i,it2)
            wkxbry_north(i)=cff1*wkxbry_north_dt(i,it1)
     &                     +cff2*wkxbry_north_dt(i,it2)
            wkebry_north(i)=cff1*wkebry_north_dt(i,it1)
     &                     +cff2*wkebry_north_dt(i,it2)
          enddo
        endif
# endif
      elseif (ZEROTH_TILE) then
        write(stdout,'(/2(1x,A)/3(1x,A,F16.10)/)')
     &   'SET_BRY_WKB_TILE - current model time is out of bounds of',
     &   '''brywkb_time''.',     'BRY_WKB_TSTART=', brywkb_time(it1)*sec2day
# ifdef USE_CALENDAR
     &  -brywkb_origin_date_in_sec*sec2day
# endif
     &   ,'TDAYS=',  tdays,    'BRY_WKB_TEND=', brywkb_time(it2)*sec2day
# ifdef USE_CALENDAR
     &  -brywkb_origin_date_in_sec*sec2day
# endif
        may_day_flag=2
      endif
      return
      end
#else
      subroutine get_bry_wkb_empty
      end
#endif

