! $Id: get_wwave.F 1458 2014-02-03 15:01:25Z gcambon $
!
!======================================================================
! CROCO is a branch of ROMS developped at IRD and INRIA, 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 BBL || defined MRL_WCI )

                              ! Read point or gridded wind wave
      subroutine get_wwave    ! height, direction, period at approp.
                              ! time from forcing NetCDF file.
      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.
!

      if (may_day_flag.ne.0) return      !-->  EXIT
      if (itww.eq.0 .or. iic.eq.0) then
        lstr=lenstr(frcname)
        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 (ncidfrc.eq.-1) then
          ierr=nf_open(frcname(1:lstr), nf_nowrite, ncidfrc)
          if (ierr .ne. nf_noerr) goto 4                !--> ERROR 
        endif
      

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

! Check if wave amplitude present, and if so, if field or scalar series
      
        lvar=lenstr(vname(1,indxWWA))
        ierr=nf_inq_varid (ncidfrc,vname(1,indxWWA)(1:lvar),wwa_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, 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),frcname(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 (ncidfrc,vname(1,indxWWD)(1:lvar),wwd_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, 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),frcname(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 (ncidfrc,vname(1,indxWWP)(1:lvar),wwp_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, 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),frcname(1:lstr)
          goto 99                                         !--> ERROR
        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 (ncidfrc,vname(1,indxWEB)(1:lvar),wweb_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, 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),frcname(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 (ncidfrc,vname(1,indxWED)(1:lvar),wwed_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, 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),frcname(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 (ncidfrc,vname(1,indxWER)(1:lvar),wwer_id)
        if (ierr .eq. nf_noerr) then
          ierr=nf_inq_varndims (ncidfrc, 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),frcname(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 (ncidfrc,   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.
     
      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(ncidfrc, ww_tid, ww_rec, cff)
      if (ierr .ne. nf_noerr) then
        write(stdout,6) 'wwv_time', ww_rec 
        goto 99                                         !--> ERROR
      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), ncidfrc, wwa_id,
     &                                           ww_rec, r2dvar)
      else
        ierr=nf_get_var1_FTYPE (ncidfrc,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), ncidfrc, wwd_id,
     &                                           ww_rec, r2dvar)
      else
        ierr=nf_get_var1_FTYPE (ncidfrc,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), ncidfrc, wwp_id,
     &                                           ww_rec, r2dvar)
      else
        ierr=nf_get_var1_FTYPE (ncidfrc,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 breaking dissipation
# if defined MRL_WCI && defined WAVE_OFFLINE_BREAKING
      if (wwebgrd.eq.1.) then
        ierr=nf_fread (wweb(START_2D_ARRAY,i), ncidfrc, wweb_id,
     &                                           ww_rec, r2dvar)
      else
        ierr=nf_get_var1_FTYPE (ncidfrc,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), ncidfrc, wwed_id,
     &                                           ww_rec, r2dvar)
      else
        ierr=nf_get_var1_FTYPE (ncidfrc,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), ncidfrc, wwer_id,
     &                                           ww_rec, r2dvar)
      else
        ierr=nf_get_var1_FTYPE (ncidfrc,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
      write(stdout,'(6x,A,1x,A,1x,g12.4,1x,I4)') 'GET_WWAVE --',
     &              'Read wind wave input data for time =', cff
# 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) frcname(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, frcname(1:lstr), tdays,
     &                                wwv_time(itww)*sec2day
  8   format(/,' GET_WWAVE - ERROR: requested time record ',I4,
     &       1x,'exeeds 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
      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)
# 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
          
            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
# 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)

# 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
        
        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
# 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
      do j=jstrR,jendR
        do i=istrR,iendR
#  ifndef SHOREFACE
!         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
          Dwave(i,j)=Dwave(i,j)-angler(i,j)
#  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 SHOREFACE || defined MEKONG
      do j=jstrR,jendR
        do i=istrR,iendR
#  ifdef MEKONG
! First estimate wavenumber
!
          Dstp=max(h(i,j)+z_w(i,j,N),1.e-3)
          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 
#  endif /* MEKONG */

#  ifndef SHOREFACE
! Simulate surf-zone dissipation due to depth-induced wave breaking
!
          Awave(i,j)=min(0.2*Dstp, Awave(i,j))
#  endif /* SHOREFACE */
        enddo
      enddo
# endif /* !SHOREFACE || MEKONG */
!
!--------------------------------------------------------------------
! Compute variables for wave-current interactions (mrl_wci routine)
! --> store in same arrays as in OW_COUPLING case
!--------------------------------------------------------------------
!
# ifdef MRL_WCI
      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  
  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 */
