module mod_readwaves
c --- Dany's/Jon's wave reading files are in here now [Tim: 2012.10.18]
c --- So far there are
c --- 1a. wavewatch 3  from IFREMER (WAVES_WW3 flag)
c --- 1b. wavewatch 3  from IFREMER - Arctic grid + waves-in-ice (WAVES_WW3_ARCTIC flag)
c --- 2. WAM North Sea from met.no (WAVES_WAMNSEA flag)
c --- 3. WAM era-i from ECMWF (WAVES_ERAI flag)
c --- 4. old WAM (no flag)
c --- swh, mwp, mwd are global variables that are now stored in mod_common_wavesice.F
c --- also need to set waves_timeres (interval between wave readings)
c --- WAVES_SKIP flag can also be used to skip readings

#if defined (WAVES)

      use mod_xc!HYCOM communication interface
      implicit none

      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
     &  fldmask     ! wave mask
#if defined (WAVES_TEST_READWAVES) && defined (WAVES_WW3_ARCTIC)
      !to dump ice concentration used by WW3
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)  :: icec_interp
#endif

      integer, parameter :: 
     &    unit_mwd   = 892
     &   ,unit_mwp   = 893
     &   ,unit_swh   = 894

#if defined (WAVES_SKIP)
      !!can skip some readings to speed up code
      real, parameter :: waves_timeres = 24.0
      !!NB waves_timeres needs to be <=24 and an integral multiple of
      !!the usual readtime of the particular wave model used
      !!(WW3, WAM_NSEA or WAM);
#else
#if defined (WAVES_WW3) || defined (WAVES_WW3_ARCTIC)
      ! set time spacing between calls to readwaves
      real, parameter :: waves_timeres = 3.0 ![h] wavewatch3 gives waves each 3h;
#else
      real, parameter :: waves_timeres = 6.0 ![h] other models give waves each 6h;
#endif
#endif/*WAVES_SKIP*/

#if defined (WAVES_WW3) || defined (WAVES_WW3_ARCTIC)
      ! set time spacing in netcdf file
      real, parameter :: waves_timeres_nc = 3.0 ![h] wavewatch3 gives waves each 3h;
#else
      real, parameter :: waves_timeres_nc = 6.0 ![h] other models give waves each 6h;
#endif

#if defined (WAVES_WW3) || defined(WAVES_WAMNSEA) || defined (WAVES_WW3_ARCTIC)
      !!pivot points only read once then stored as global variables
      integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
     &   ipiv_all,jpiv_all
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
     &   dist_all
#endif

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      contains

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine read_waves(dtime_waves)
      !! Wrapper which diverts caller to appropriate subroutine for the
      !! input wave model;
      !! swh,mwp,mwd interpolated and read into mod_common_wavesice
      !! Also does some tests (WAVES_TEST_READWAVES flag);

      use mod_common_wavesice,only: mwp,mwd
#if defined (WAVES_TEST_READWAVES)
     &   ,swh!,dfloe
     &   ,dump_netcdf
      !use mod_xc
      use mod_hycom_nersc!->ncdraft,rdncrec,rdncrec3d
#endif
      implicit none
      real*8, intent(in)  :: dtime_waves
      real                 :: dir,ddir
      integer              :: i,j
#if defined (WAVES_TEST_READWAVES)
      integer                       :: iyear, ihour, iday
      character(len=4)              :: cyy
      character(len=3)              :: cdd,cdfrac
      character(len=12)             :: dtime_waves_human!!yyyy_ddd.xxx, where dd starts at 0
      character(len=80)             :: ncfil2
      character(len=*), parameter   :: outdir = 'WAVES/test_input/'
      !!
      real, dimension(itdm,jtdm) :: fld
      logical, parameter         :: CHK_XTRA = .true.
      !!
      include 'common_blocks.h'
#endif

      !!IF CHANGING WAVE MODEL MAKE A FLAG, AND MAKE A SPECIAL WAVE
      !! READING SUBROUTINE BELOW;
#if defined (WAVES_WAMNSEA)
      call read_waves_wamnsea(dtime_waves)
#elif defined (WAVES_WW3)
      call read_waves_ww3(dtime_waves)
#elif defined (WAVES_WW3_ARCTIC)
      call read_waves_ww3_arctic(dtime_waves)
#elif defined (WAVES_ERAI)
      call read_waves_erai(dtime_waves)
#else
      call read_waves_ecmwf(dtime_waves)
#endif

      !!Make sure mwd is in [0,360]:
      do j=1,jj
      do i=1,ii
         dir   = mwd(i,j)
         if (mwp(i,j).gt.0.) then
            if (dir.gt.360.) then
               ddir     = dir-360.
               mwd(i,j) = dir-ceiling(ddir/360.)*360.
            elseif (dir.lt.0.) then
               ddir     = -dir
               mwd(i,j) = dir+ceiling(ddir/360.)*360.
            end if
         end if
      end do
      end do

#if defined (WAVES_TEST_READWAVES)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !!dump waves into a netcdf file
      !!set its name:
      call forday(dtime_waves,yrflag, iyear,iday,ihour)

      write(cyy  ,'(i4.4)') iyear
      write(cdd,'(i3.3)') iday-1!!start at iday=0
      write(cdfrac,'(i3.3)') int(1.0e3*(ihour/24.0))
      dtime_waves_human = cyy//'_'//cdd//'.'//cdfrac

#if defined (WAVES_WAMNSEA)
      ncfil2=outdir//'waves_in_wamnsea_'
     &   //dtime_waves_human//'.nc'
#elif (WAVES_WW3)
      ncfil2=outdir//'waves_in_ww3_'
     &   //dtime_waves_human//'.nc'
#elif (WAVES_WW3_ARCTIC)
      ncfil2=outdir//'waves_in_ww3a_'
     &   //dtime_waves_human//'.nc'
#elif (WAVES_ERAI)
      ncfil2=outdir//'waves_in_erai_'
     &   //dtime_waves_human//'.nc'
#else
      ncfil2=outdir//'waves_in_ecmwf_'
     &   //dtime_waves_human//'.nc'
#endif
      if (mnproc==1)write(*,*)"ncfil2: ",ncfil2

      !!set up lon,lat,...
      if (mnproc==1)print*,'TW: dumping input waves into a netcdf file'
      call xcaget(fld,plon,0)
      call ncdraft(trim(ncfil2),fld,'lon','clobber')
      call xcaget(fld,plat,0)
      call ncdraft(trim(ncfil2),fld,'lat','')

      !!Check if wave fields have been correctly read and dump them in a netcdf file
      call dump_netcdf(trim(ncfil2),swh,'swh')
      call dump_netcdf(trim(ncfil2),mwp,'mwp')
      call dump_netcdf(trim(ncfil2),mwd,'mwd')
#if defined (WAVES_WW3_ARCTIC)
      !dump ice conc used by WW3
      call dump_netcdf(trim(ncfil2),icec_interp,'icec')
#endif

      !!Check final wave mask
      if (CHK_XTRA) then
         call dump_netcdf(trim(ncfil2),fldmask,'mask')
      end if
      if (mnproc==1)print*,'TW: finished dumping input waves'
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      end subroutine read_waves
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


#if defined(WAVES_WAMNSEA)
      subroutine read_pivots_wamnsea()
      use mod_za  ! HYCOM I/O interface

      implicit none

c --- Data catalogue name
      character(len=*), parameter :: piv_path = './WAVES/pivots/'
      character(len=80)           :: pivfil
      character(len=*), parameter :: vname = 'wamnsea'

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=248, nlat=400
c
      integer  :: i,j,k,m,ix,jy,nopa
      character preambl(5)*79
c
      include 'common_blocks.h'
c

c --- ******  WRITE  ********    
      pivfil   = piv_path//'pivots_'//vname

      !iunit for forfun.F (rdmonthck)
      nopa  = 926
      call zaiopf(trim(pivfil)//'.a','old',nopa)
      if(mnproc.eq.1) then  ! .b file from 1st tile only
         print*,"Read in pivot points for WAM NSEA files"
         write(*,*)"ipiv files: ",trim(pivfil)//'.[ab]'

         !open .b file & read preable before calling rdmonthck
         open (unit=uoff+nopa,file=trim(pivfil)//'.b',
     &     status='old', action='read')
         read (uoff+nopa,'(a79)') preambl
      endif !1st tile
      call preambl_print(preambl)
      if(mnproc==1)print*,'1st rec: '//pivfil//'.a'
      call rdmonthck(util1, nopa,1)
      if(mnproc==1)print*,'2nd rec: '//pivfil//'.a'
      call rdmonthck(util2, nopa,2)
      if(mnproc==1)print*,'3rd rec: '//pivfil//'.a'
      call rdmonthck(util3, nopa,3)
      if(mnproc==1)print*,'ipiv: finished rdmonthck'

      !close .b file after calling rdmonthck
      if (mnproc==1)close(unit=uoff+nopa)
      if(mnproc==1)print*,'ipiv: closed .b file'
      call zaiocl(nopa)
      if(mnproc==1)print*,'ipiv: closed .a file'
      call xcsync(flush_lp)

      do j= 1,jj
      do i= 1,ii
         ipiv_all(i,j) = util1(i,j)
         jpiv_all(i,j) = util2(i,j)
         dist_all(i,j) = util3(i,j)
      enddo!i
      enddo!j

      end subroutine read_pivots_wamnsea

      subroutine read_waves_wamnsea(dtime_waves) 
      ! rewritten from read_waves_ecmwf(dtime_waves)
      use mod_hycom_nersc
      use netcdf
      use mod_za  ! HYCOM I/O interface
      use mod_common_wavesice, only:swh,mwp,mwd
  
      use mod_year_info, only:
     &   year_info, year_day,refyear,datetojulian,juliantodate
      
      implicit none
      real*8 ,  intent(in) :: dtime_waves
c
c --- Data catalogue name
      character(len=*), parameter :: wave_path0='./WAVES/WAMNSEA/'
      character(len=80) :: wave_path ! ='./WAVES/WAMNSEA/YYYY/analysis
c
c --- Variable names - myocean.met.no WAM north sea data 
      character(len=*), parameter :: vmwd='wave_direction'
      character(len=*), parameter :: vmwp='peak_wave_period'
      character(len=*), parameter :: vswh='significant_wave_height'

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=248, nlat=400
c
      integer  :: i,j,k,m,ix,jy,irec,nrecs,maxrec,nlon2,nlat2
      real, dimension(nlon,nlat) :: tmp_data
      integer :: ipiv,jpiv
!     
      real                 :: missing
      integer              :: dimid, varid, ncid
      
      character(len=23)    :: ncname
      character(len=80)    :: ncfil1, ncfil2
      character(len=4)     :: cyy
      character(len=2)     :: cmm, cdm
      character(len=8)     :: cdate0, cdate1
      character(len=10)    :: ctt
      integer, parameter   :: Nfc=7   ! no of records in forecast
      real                 :: time_vec(Nfc),rday1,rday0
      real*8               :: dtime_diff
      !!
      integer           :: jday0,jday1,ryear,rmon,rday
      type(year_info)   :: tt0,tt1

      logical :: file_exists
!
      include 'common_blocks.h'
!
      !WAMNSEA ref time is 1970-01-01 0:00:00 (units = seconds)
      ryear    = 1970
      rmon     = 1
      rday     = 1

      !!get time/date info for current day
      call year_day(dtime_waves,refyear,tt1,yrflag,.false.)          !refyear not used with refflag=.false.
      jday1 = datetojulian(tt1%iyy,tt1%imm,tt1%idm+1,ryear,rmon,rday)!days from ref date
      rday1 = jday1+tt1%ihh/24.
      if(mnproc==1)print*,'TW: wamnsea',jday1,rday1,tt1%ihh
     &                     ,tt1%ihh/waves_timeres_nc

      ncname      = 'wam_nsea.an.'//tt1%cdate//'.nc'
      wave_path   = wave_path0//'/'//tt1%cyy//'/analysis/'
      ncfil1      = trim(wave_path)//ncname
      nrecs       = 24./waves_timeres_nc
      inquire(file=ncfil1,exist=file_exists)

      if (file_exists) then
         !!analysis for current day exists, so continue
         irec  = nint(tt1%ihh/waves_timeres_nc+1)

         !!check correct number of time rec's
         call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
         call ncerr(nf90_inq_dimid(ncid,"time",dimid))
         call ncerr(nf90_inquire_dimension(ncid,dimid,len=maxrec))
         call ncerr(nf90_close    (ncid))
         if (mnproc==1) then
            print*,ncfil1
            print*,maxrec,nrecs
         end if
         if (maxrec/=nrecs) then
            if (mnproc==1) then
               print*,
     &         'readwaves_wamnsea: Dimension mismatch in WAM NSEA data'
               print*,'maxrec :',maxrec,nrecs
            end if
         end if

      else
         !!analysis for current day doesn't exist, so try forecast file
         ncfil2   = wave_path0//'/wam_nsea.fc.latest.nc'
         if (mnproc==1) then
            print*,'analysis file '//ncfil1//' not present'
            print*,'looking for forecast file '//ncfil2
         end if

         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         inquire(file=ncfil2,exist=file_exists)
         if (.not.file_exists) then
            if(mnproc==1)print*,'forecast file '//ncfil2//' not present'
         end if
         call xcstop('(read_waves_wamnsea)')
         stop '(read_waves_wamnsea)'

!        if(mnproc==1)print*,'wamnsea: opening '//ncfil2
         call ncerr(NF90_open(trim(ncfil2),NF90_NOCLOBBER,ncid))
         call ncerr(nf90_inq_dimid(ncid,"time",dimid))
         call ncerr(nf90_inquire_dimension(ncid,dimid,len=maxrec))
         call ncerr(nf90_close    (ncid))
         if (mnproc==1) then
            !!test
            print*,ncfil2
            print*,maxrec,Nfc
         end if
         if (maxrec/=Nfc) then
            if (mnproc==1) then
               print*,
     &         'readwaves_wamnsea: Dimension mismatch in WAM NSEA data'
               print*,'maxrec :',maxrec,Nfc
            end if
         end if

         !!get time values
!        if(mnproc==1)print*,'wamnsea: opening '//ncfil2
         call ncerr(NF90_open(trim(ncfil2),NF90_NOCLOBBER,ncid))
         call ncerr(nf90_inq_varid(ncid,'time',varid))
         call ncerr(nf90_get_var  (ncid,varid,time_vec))
         call ncerr(nf90_close    (ncid))

         rday0 = time_vec(1)/(3600.*24.)              !days since 1970-1-1 of TOMORROW
         call year_day(rday0,ryear,tt0,yrflag,.true.) !1 Jan in "ryear" is ref date when refflag=.true.
         if (mnproc==1) then
            print*,'date of forecast start: ',tt0%cdate
            print*,'current date: '//tt1%cdate
            print*,'current time: '//tt1%ctime
         end if

         !!check if file is usable
         if (rday1.lt.rday0) then
            !!current date is before start of forecast file
            if (mnproc==1) then
               print*,'current date is before start of forecast file'
               print*,'but '//ncfil1//' is not present'
            end if
            call xcstop('(readwaves_wamnsea)')
            stop '(readwaves_wamnsea)'

         else
            !!current date is after start of forecast file
            dtime_diff  = (rday1-rday0)*24.!time measured in days - convert to hours
            irec        = nint(dtime_diff/waves_timeres_nc)+1
            if (mnproc==1) then
               !!test
               print*,'wamnsea days from 1970-1-1 (current day,netcdf)'
     &            ,rday1,rday0
               print*,'time diff (h)',dtime_diff
               print*,'irec',irec
            end if

            if (irec.gt.maxrec) then
               !!too far into future
               if(mnproc==1)then
                  print*,'irec too large:',irec,maxrec
               end if
               call xcstop('(readwaves_wamnsea)')
               stop '(readwaves_wamnsea)'

            else
               !!can use forecast file
               ncfil1   = ncfil2
            end if
         end if!check forecast file
      end if!check analysis file exists

! --- ******  WRITE  ********    
      if (mnproc==1) then
         write(*,*)" "
         write(*,*)"*********************************"
         write(*,*)"***** READ WAVES WAMNSEA ********"
         write(*,*)"*********************************"
         write(*,*)"ncfil,irec:",trim(ncfil1),irec
         write(*,*)"Current date: "//tt1%cdate
         write(*,*)"Current time: "//tt1%ctime
         write(*,*)"*********************************"
         write(*,*)" "
      end if

      !!check grid dimensions in .nc file
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_dimid(ncid,"rlon",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlon2))
      call ncerr(nf90_inq_dimid(ncid,"rlat",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlat2))
      call ncerr(nf90_close    (ncid))
      
      if (nlon/=nlon2.or.nlat/=nlat2) then
         if (mnproc==1) then
            write(lp,'(a)') 'DD :: Dimension mismatch in WAM NSEA data'
            write(lp,'(a,2i5)') 'nlon :',nlon,nlon2
            write(lp,'(a,2i5)') 'nlat :',nlat,nlat2
         end if
      end if
      

c --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)"Start reading vswh data, irec=",irec

c --- WAM NSEA waves (every 6 hours)
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_varid(ncid,vswh,varid))
      call ncerr(nf90_get_var  (ncid,varid,tmp_data,start=(/1,1,irec/)))
      call ncerr(nf90_get_att  (ncid,varid,'_FillValue',missing))
      call ncerr(nf90_close    (ncid))

      ! Build the data mask from missing values
      fldmask  = 0.0
      do j=1,jj
      do i=1,ii
         ipiv  = ipiv_all(i,j)
         jpiv  = jpiv_all(i,j)
         if (ipiv.gt.0) then
            !!inside the data grid
            if (.not.(tmp_data(ipiv,jpiv)==missing)) then
               !!data isn't missing
               fldmask(i,j)   = 1.0
            end if
         end if
      enddo
      enddo

      !!get significant wave height
      call rdncrec(ncfil1,vswh,tmp_data,nlon,nlat,irec)

      swh   = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_wamnsea
         !! (this only needs to be done once, in mod_restart.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            swh(i,j) = tmp_data(ipiv,jpiv)
         endif
      enddo !i
      enddo !j


! --- Read in mean wave direction mwd
      call rdncrec(ncfil1,vmwd,tmp_data,nlon,nlat,irec)

      mwd   = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_wamnsea
         !! (this only needs to be done once, in mod_waves_init.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            mwd(i,j) = tmp_data(ipiv,jpiv)+180.0
            !! convert from waves-to direction to waves-from direction;
            !! convention is: 0 = from north, 90 = from east
         endif
      enddo !i
      enddo !j

! --- Read in peak wave period mwp
      call rdncrec(ncfil1,vmwp,tmp_data,nlon,nlat,irec)!applies scale factor etc

      mwp   = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_wamnsea
         !! (this only needs to be done once, in mod_restart.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            mwp(i,j) = tmp_data(ipiv,jpiv)
         endif
      enddo !i
      enddo !j

      end subroutine read_waves_wamnsea

#elif defined(WAVES_WW3)
      subroutine read_pivots_ww3!(dtime_waves)
      use mod_za  ! HYCOM I/O interface

      implicit none

c --- Data catalogue name
      !!wave_path=/work/shared/nersc/msc/WW3
      character(len=*), parameter :: piv_path = './WAVES/pivots/'

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=720, nlat=317
c
      integer  :: i,j,k,m,ix,jy!,irec,maxrec,nlon2,nlat2

! --- Initialize pivot points relating model grid and ww3 grid 
!     integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
!    & ipiv,jpiv
!     real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
!    & ipivreal, jpivreal, dist


      character preambl(5)*79
c
      include 'common_blocks.h'
c

c --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)
     & "Read in ipiv,jpiv from Model to WW3 Global files"

      !wave_path=./WAMNSEA/
      ! use rdmonthch directly
      call zaiopf(piv_path//'pivots_ww3.a', 'old'
     & , 926)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+926,
     & file=piv_path//'pivots_ww3.b',
     &     status='old', action='read')
        read (uoff+926,'(a79)') preambl
        endif !1st tile
      call preambl_print(preambl)
      call rdmonthck(util1, 926,1)
      call rdmonthck(util2, 926,2)
      call rdmonthck(util3, 926,3)
!      call xctilr( util1, 1,1,    1,   1, halo_ps)
! --- call rdmonth(util3, 990)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
      write(*,*)"ipiv file: ",piv_path,'pivots_ww3.a'
      close (unit=uoff+926)
        endif !1st tile
      call zaiocl(926)
      call xcsync(flush_lp)

       do j= 1,jj
         do i= 1,ii
            ipiv_all(i,j) = util1(i,j)
            jpiv_all(i,j) = util2(i,j)
            dist_all(i,j) = util3(i,j)
         enddo !i
      enddo !j

      end subroutine read_pivots_ww3

      subroutine read_waves_ww3(dtime_waves)
      use mod_hycom_nersc
      use netcdf
      use mod_za  ! HYCOM I/O interface
      use mod_year_info, only: year_info, year_day,refyear
      use mod_common_wavesice, only:swh,mwp,mwd

      implicit none
      real*8 ,  intent(in) :: dtime_waves
c      integer,  intent(in) :: slot
      type(year_info) :: tt

c
c --- Data catalogue name
      !!wave_path=/work/shared/nersc/msc/WW3
      character(len=*), parameter :: wave_path = './WAVES/WW3/'
      character(len=17)           :: wave_path_year
c
c --- Variable names - ERA40 1.125 by 1.125 
      character(len=*), parameter :: vmwd='dir'
      character(len=*), parameter :: vmwp='t'
      character(len=*), parameter :: vswh='hs' 

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=720, nlat=317
      integer, parameter :: recs_per_day_ww3 = 8!3h resolution
c
      integer  :: i,j,k,m,ix,jy,irec,maxrec,nlon2,nlat2
      integer  :: mon_num,nrec_ww3,nrec2
      real, dimension(nlon,nlat) :: tmp_data!!TW: to save memory, reduced number of big arrays

! --- Initialize pivot points relating model grid and ww3 grid 
      integer  :: ipiv,jpiv!TW: NB no longer an array

      real               :: missing
      integer            :: dimid, varid, ncid
      integer            :: iyear, ihour, iday, icdm

      character(len=80)  :: ncfil1, ncfil2
      character(len=4)   :: cyy, cyyp1
      character(len=10)  :: ctt
      character(len=2)   :: cmm, cdm 

c
      include 'common_blocks.h'
c
      call forday(dtime_waves,yrflag, iyear,iday,ihour)
!      irec  =  (iday-1)*4 + ihour/6 +1

c --- Set file names
      write(cyy  ,'(i4.4)') iyear
      write(cyyp1,'(i4.4)') iyear+1

! --- Set file names
      call year_day(dtime_waves,refyear,tt,yrflag)

      cmm=''
      cmm(1:2)=tt%cmm(1:2)
      

! --- Calculate the record in the montly file, 8 rec/h (3h)      
      !irec = (tt%idm-1)*8 + ihour/3 +1
      irec = (tt%idm)*8 + ihour/3 +1

! --- Days per month
      mon_num  = tt%imm
      nrec_ww3 = recs_per_day_ww3*tt%totdim(mon_num)

! --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)"*********************************"
      if (mnproc==1)write(*,*)"***** READ WAVES WW3 Global ********"
      if (mnproc==1)write(*,*)"*********************************"
      if (mnproc==1)write(*,*)"iy,id,ih,irec:",iyear,iday,ihour,irec
      if (mnproc==1)write(*,*)"cyy,cmm,idm+1: ",tt%cyy,tt%cmm,tt%idm+1

c --- Check dimensions against hardcoded ones
c --- Get longitude  & latitude for data  - use land mask file
      wave_path_year=wave_path//cyy//'/'
      ncfil1=wave_path_year//'ww3.'//cyy//cmm//'_hs.nc'

! --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)"name of IN netcdf file: ",ncfil1

      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_dimid(ncid,"longitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlon2))
      call ncerr(nf90_inq_dimid(ncid,"latitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlat2))
      call ncerr(nf90_inq_dimid(ncid,"time",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nrec2))
      call ncerr(nf90_close(ncid))

      if (mnproc==1) write(*,*)'nlon2,nlat2,nrec',nlon2,nlat2,nrec2 
      if (nlon/=nlon2.or.nlat/=nlat2) then
         if (mnproc==1) then
            write(lp,'(a)') 'DD :: Dimension mismatch in WW3Global data'
            write(lp,'(a,2i5)') 'nlon :',nlon,nlon2
            write(lp,'(a,2i5)') 'nlat :',nlat,nlat2
         end if
      end if

      if (nrec_ww3/=nrec2) then
         if (mnproc==1) then
            write(lp,'(a)') 'TW: Corrupted netcdf file:',ncfil1
            write(lp,'(a,2i5)') 'nrec :',nrec_ww3,nrec2
            write(lp,'(a,2i5)') 'irec :',irec
         end if
      end if

c --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)"Start reading Hs data, irec=",irec

c --- WW3Global waves (every 3 hours)
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))!!NEEDED?? (USE SAME FILE FOR lon/lat EARLIER)
      call ncerr(nf90_inq_varid(ncid,vswh,varid))
      call ncerr(nf90_get_var  (ncid,varid,tmp_data,start=(/1,1,irec/)))
      call ncerr(nf90_get_att  (ncid,varid,'_FillValue',missing))
      call ncerr(nf90_close    (ncid))
      
      if (mnproc==1) print*,'TW3 (mrw)'
      ! Build the data mask from missing values
      ! NB do this before rdncrec, since
      !  that applies scale_factor and add_offset
      fldmask  = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_ww3
         !! (this only needs to be done once, in mod_waves_init.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (ipiv.gt.0) then
            !!inside the data grid
            if (.not.(tmp_data(ipiv,jpiv)==missing)) then
               !!data isn't missing
               fldmask(i,j)   = 1.0
            end if
         endif
      enddo !i
      enddo !j

      ! Get swh
      call rdncrec(ncfil1,vswh,tmp_data,nlon,nlat,irec)
      swh   = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_ww3
         !! (this only needs to be done once, in mod_waves_init.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            swh(i,j) = tmp_data(ipiv,jpiv)
         endif
      enddo !i
      enddo !j

! Read in mean wave period data from new file
      !!wave_path=/work/shared/nersc/msc/WW3
      ncfil1=wave_path_year//'ww3.'//cyy//cmm//'_t.nc'
      if (mnproc==1)write(*,*)"Start reading Tp data, irec=",irec
      if (mnproc==1) print*,'ncfil: ',ncfil1

      !call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      !call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      !call ncerr(nf90_inq_varid(ncid,vmwp,varid))
      !call ncerr(nf90_get_var  (ncid,varid,tmp_data,start=(/1,1,irec/)))
      !call ncerr(nf90_get_att  (ncid,varid,'_FillValue',missing))
      !call ncerr(nf90_close    (ncid))

      if (mnproc==1) print*,'TW4a (mnw)'
      call rdncrec(ncfil1,vmwp,tmp_data,nlon,nlat,irec)
      mwp   = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_ww3
         !! (this only needs to be done once, in mod_restart.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            mwp(i,j) = tmp_data(ipiv,jpiv)
         endif
      enddo !i
      enddo !j


! Read in mean wave direction data from new file
      !!wave_path=/work/shared/nersc/msc/WW3
      ncfil1=wave_path_year//'ww3.'//cyy//cmm//'_dir.nc'
      if (mnproc==1)write(*,*)"Start reading mwd data, irec=",irec
      if (mnproc==1) print*,'ncfil: ',ncfil1

      !call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      !call ncerr(nf90_inq_varid(ncid,vmwd,varid))
      !call ncerr(nf90_get_var  (ncid,varid,tmp_data,start=(/1,1,irec/)))
      !call ncerr(nf90_close(ncid))
      
      if (mnproc==1) print*,'TW5a (mnw)'
      call rdncrec(ncfil1,vmwd,tmp_data,nlon,nlat,irec)

      if (mnproc==1) print*,'TW4'
      mwd   = -360.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_ww3
         !! (this only needs to be done once, in mod_restart.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            !! this is already waves-from direction so don't need to add 180
            !! convention is: 0 = from north, 90 = from east
            mwd(i,j) = tmp_data(ipiv,jpiv)
         endif

         if ((i==itest).and.(j==jtest)) then
            print*,'TW: i/j/lon/lat:', i+i0,j+j0,plon(i,j),plat(i,j)
            print*,'TW: swh',swh(i,j)
            print*,'TW: mwp',mwp(i,j)
            print*,'TW: mwd',mwd(i,j)
            print*,'TW: mask',fldmask(i,j)
         end if

      enddo !i
      enddo !j
       
      end subroutine read_waves_ww3

#elif defined(WAVES_WW3_ARCTIC)
      subroutine read_pivots_ww3_arctic
      use mod_za  ! HYCOM I/O interface

      implicit none

c --- Data catalogue name
      !!wave_path=/work/shared/nersc/msc/WW3
      character(len=*), parameter :: piv_path = './WAVES/pivots/'
      character(len=*), parameter :: piv_name = 'pivots_ww3_arctic'

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=450, nlat=400
c
      integer  :: i,j,k,m,ix,jy!,irec,maxrec,nlon2,nlat2

! --- Initialize pivot points relating model grid and ww3 grid 
!     integer, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
!    & ipiv,jpiv
!     real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
!    & ipivreal, jpivreal, dist


      character preambl(5)*79
c
      include 'common_blocks.h'
c

c --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)
     & "Read in ipiv,jpiv from Model to WW3_ARCTIC files"

      !wave_path=./WAMNSEA/
      ! use rdmonthch directly
      call zaiopf(piv_path//piv_name//'.a', 'old'
     & , 926)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+926,
     & file=piv_path//piv_name//'.b',
     &     status='old', action='read')
        read (uoff+926,'(a79)') preambl
        endif !1st tile
      call preambl_print(preambl)
      call rdmonthck(util1, 926,1)
      call rdmonthck(util2, 926,2)
      call rdmonthck(util3, 926,3)
!      call xctilr( util1, 1,1,    1,   1, halo_ps)
! --- call rdmonth(util3, 990)
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
      write(*,*)"ipiv file: ",piv_path,piv_name//'.a'
      close (unit=uoff+926)
        endif !1st tile
      call zaiocl(926)
      call xcsync(flush_lp)

       do j= 1,jj
         do i= 1,ii
            ipiv_all(i,j) = util1(i,j)
            jpiv_all(i,j) = util2(i,j)
            dist_all(i,j) = util3(i,j)
         enddo !i
      enddo !j

!     Create ipiv,jpiv as real varables for output
!     ipivreal = ipiv_all
!     jpivreal = jpiv_all

      end subroutine read_pivots_ww3_arctic

      subroutine read_waves_ww3_arctic(dtime_waves)
      use mod_hycom_nersc
      use netcdf
      use mod_za  ! HYCOM I/O interface
      use mod_common_wavesice, only:swh,mwp,mwd
      use mod_year_info, only:
     &   year_info, year_day,refyear,datetojulian,juliantodate

      implicit none
      real*8 ,  intent(in) :: dtime_waves
c      integer,  intent(in) :: slot
      type(year_info) :: tt

c
c --- Data catalogue name
      !!wave_path=/work/shared/nersc/msc/WW3
      character(len=*), parameter   :: wave_path = './WAVES/WW3_ARCTIC/'
      character(len=50)             :: wave_path_year
      character(len=32)             :: ncname
c
c --- Variable names - ERA40 1.125 by 1.125 
      character(len=*), parameter :: vmwd='dir'
      character(len=*), parameter :: vmwp='fp' !NB peak freq not period
      character(len=*), parameter :: vswh='hs' 

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=450, nlat=400
      integer, parameter :: recs_per_day_ww3 = 8!3h resolution
c
      integer  :: i,j,k,m,ix,jy,irec,nrecs,maxrec,nlon2,nlat2
      integer  :: mon_num,nrec_ww3,nrec2
      real, dimension(nlon,nlat) :: tmp_data,ice_conc
      real,parameter :: icec_min = 0.01 !waves set to 0 if ice_conc>icec_min

! --- Initialize pivot points relating model grid and ww3 grid 
      integer  :: ipiv,jpiv!TW: NB no longer an array

      real               :: missing
      integer            :: dimid, varid, ncid
      integer            :: iyear, ihour, iday, icdm

      character(len=80)  :: ncfil1, ncfil2
      character(len=4)   :: cyy, cyyp1
      character(len=10)  :: ctt
      character(len=2)   :: cmm, cdm 

      integer, parameter   :: Nfc = 49   ! no of records in forecast
      real                 :: time_vec(Nfc),rday1,rday0,fp_ij
      real*8               :: dtime_diff
      integer              :: jday0,jday1,ryear,rmon,rday
      type(year_info)      :: tt0,tt1

      logical :: file_exists

c
      include 'common_blocks.h'
c
      !WW3_ARCTIC ref time is 1990-01-01 0:00:00 (units = days)
      ryear    = 1990
      rmon     = 1
      rday     = 1

      !!get time/date info for current day
      call year_day(dtime_waves,refyear,tt1,yrflag,.false.)          !refyear not used with refflag=.false.
      jday1 = datetojulian(tt1%iyy,tt1%imm,tt1%idm+1,ryear,rmon,rday)!days from ref date
      rday1 = jday1+tt1%ihh/24.
      if(mnproc==1)print*,'TW: ww3_arctic',jday1,rday1,tt1%ihh
     &                     ,tt1%ihh/waves_timeres_nc


      ncname         = 'SWARP_WW3_ARCTIC-12K_'//tt1%cdate//'.nc'
      wave_path_year = wave_path//'/'//tt1%cyy
      ncfil1         = trim(wave_path_year)//'/analysis_m1/'//ncname
      nrecs          = 24./waves_timeres_nc
      inquire(file=ncfil1,exist=file_exists)

      if (file_exists) then
         !!analysis for current day exists, so continue
         if(mnproc==1)print*,'Using WW3 Arctic ncfile: '//ncfil1

         irec  = nint(tt1%ihh/waves_timeres_nc+1)

         !!check correct number of time rec's
         call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
         call ncerr(nf90_inq_dimid(ncid,"time",dimid))
         call ncerr(nf90_inquire_dimension(ncid,dimid,len=maxrec))
         call ncerr(nf90_close    (ncid))

         if (maxrec/=nrecs) then
            if (mnproc==1) then
               print*, 'readwaves_ww3_arctic: '
     &                     //'Dimension mismatch in WW3 ARCTIC data'
               print*,'Obtained/expected no of recs:',maxrec,nrecs
            end if
            !call xcstop('(read_waves_ww3_arctic)')
            !stop '(read_waves_ww3_arctic)'
         end if

      else
         !!analysis for current day doesn't exist, so try forecast file
         ncfil2   = wave_path//'/SWARP_WW3_ARCTIC-12K.fc.latest.nc'
         if (mnproc==1) then
            print*,'analysis file '//ncfil1//' not present'
            print*,'looking for forecast file '//ncfil2
         end if

         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         inquire(file=ncfil2,exist=file_exists)
         if (.not.file_exists) then
            if(mnproc==1)print*,'forecast file '//ncfil2//' not present'
            call xcstop('(read_waves_ww3_arctic)')
            stop '(read_waves_ww3_arctic)'
         end if

!        if(mnproc==1)print*,'wamnsea: opening '//ncfil2
         call ncerr(NF90_open(trim(ncfil2),NF90_NOCLOBBER,ncid))
         call ncerr(nf90_inq_dimid(ncid,"time",dimid))
         call ncerr(nf90_inquire_dimension(ncid,dimid,len=maxrec))
         call ncerr(nf90_close    (ncid))
         if (mnproc==1) then
            !!test
            print*,ncfil2
            print*,maxrec,Nfc
         end if
         if (maxrec/=Nfc) then
            if (mnproc==1) then
               print*,
     &            'WARNING (readwaves_ww3_arctic): '
     &               //'Dimension mismatch in WW3 Arctic data'
               print*,'maxrec :',maxrec,Nfc
            end if
            !call xcstop('(read_waves_ww3_arctic)')
            !stop '(read_waves_ww3_arctic)'
         end if

         !!get time values
!        if(mnproc==1)print*,'wamnsea: opening '//ncfil2
         call ncerr(NF90_open(trim(ncfil2),NF90_NOCLOBBER,ncid))
         call ncerr(nf90_inq_varid(ncid,'time',varid))
         call ncerr(nf90_get_var  (ncid,varid,time_vec(1:maxrec)))
         call ncerr(nf90_close    (ncid))

         rday0 = time_vec(1)                          !days since 1990-1-1 of TODAY
         call year_day(rday0,ryear,tt0,yrflag,.true.) !1 Jan in "ryear" is ref date when refflag=.true.
         if (mnproc==1) then
            print*,'date of forecast start: ',tt0%cdate
            print*,'current date: '//tt1%cdate
            print*,'current time: '//tt1%ctime
         end if

         !!check if file is usable
         if (rday1.lt.rday0) then
            !!current date is before start of forecast file
            if (mnproc==1) then
               print*,'current date is before start of forecast file'
               print*,'but '//ncfil1//' is not present'
            end if
            call xcstop('(readwaves_ww3_arctic)')
            stop '(readwaves_ww3_arctic)'

         else
            !!current date is after start of forecast file
            dtime_diff  = (rday1-rday0)*24.!time measured in days - convert to hours
            irec        = nint(dtime_diff/waves_timeres_nc)+1
            if (mnproc==1) then
               !!test
               print*,'wamnsea days from 1990-01-01'
     &            //' (current day,netcdf)',rday1,rday0
               print*,'time diff (h)',dtime_diff
               print*,'irec',irec
            end if

            if (irec.gt.maxrec) then
               !!too far into future
               if(mnproc==1)then
                  print*,'irec too large:',irec,maxrec
               end if
               call xcstop('(readwaves_ww3_arctic)')
               stop '(readwaves_ww3_arctic)'

            else
               !!can use forecast file
               ncfil1   = ncfil2
            end if
         end if!check forecast file
      end if!check analysis file exists

! --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)"name of IN netcdf file: ",ncfil1

      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_dimid(ncid,"longitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlon2))
      call ncerr(nf90_inq_dimid(ncid,"latitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlat2))
      call ncerr(nf90_close(ncid))

      if (mnproc==1) write(*,*)'nlon2,nlat2',nlon2,nlat2 
      if (nlon/=nlon2.or.nlat/=nlat2) then
         if (mnproc==1) then
            write(lp,'(a)') 'Dimension mismatch in WW3 ARCTIC data'
            write(lp,'(a,2i5)') 'nlon :',nlon,nlon2
            write(lp,'(a,2i5)') 'nlat :',nlat,nlat2
         end if
      end if

c --- ******  WRITE  ********    
      if (mnproc==1)write(*,*)"Start reading Hs data, irec=",irec

c --- WW3 ARCTIC waves (every 3 hours)
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_varid(ncid,vswh,varid))
      call ncerr(nf90_get_var  (ncid,varid,tmp_data,start=(/1,1,irec/)))
      call ncerr(nf90_get_att  (ncid,varid,'_FillValue',missing))
      call ncerr(nf90_close    (ncid))
      
      ! Build the data mask from missing values
      ! NB do this before rdncrec, since
      !  that applies scale_factor and add_offset
      fldmask  = 0.0
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_ww3_arctic
         !! (this only needs to be done once, in mod_waves_init.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (ipiv.gt.0) then
            !!inside the data grid
            if (.not.(tmp_data(ipiv,jpiv)==missing)) then
               !!data isn't missing
               fldmask(i,j)   = 1.0
            end if
         endif
      enddo !i
      enddo !j

      ! Get ice conc
      ! - attenuation applied in ice, but this could be outside
      !   our ice mask
      call rdncrec(ncfil1,'ice',ice_conc,nlon,nlat,irec)
#if defined (WAVES_TEST_READWAVES)
      !for testing: interp ice conc so it can be dumped to netcdf
      call interp_nearest(icec_interp,ice_conc,nlon,nlat)
#endif

      ! Get swh (outside ice)
      call rdncrec(ncfil1,vswh,tmp_data,nlon,nlat,irec)
      where (ice_conc.gt.icec_min) tmp_data = 0.
      call interp_nearest(swh,tmp_data,nlon,nlat)

      ! Read in mean wave period data (outside ice)
      if (mnproc==1)write(*,*)"Start reading Tp data, irec=",irec
      if (mnproc==1) print*,'ncfil: ',ncfil1
      call rdncrec(ncfil1,vmwp,tmp_data,nlon,nlat,irec)
      where (ice_conc.gt.icec_min) tmp_data = 0.
      call interp_nearest(mwp,tmp_data,nlon,nlat)
      where (mwp.gt.0) mwp = 1/mwp !NB convert from freqency to period


      ! Read in mean wave direction data (outside ice)
      ! - this is already waves-from direction so don't need to add 180
      ! - convention is: 0 = from north, 90 = from east
      if (mnproc==1)write(*,*)"Start reading mwd data, irec=",irec
      if (mnproc==1) print*,'ncfil: ',ncfil1
      call rdncrec(ncfil1,vmwd,tmp_data,nlon,nlat,irec)
      where (ice_conc.gt.icec_min) tmp_data = 0.
      call interp_nearest(mwd,tmp_data,nlon,nlat)

      do j= 1,jj
      do i= 1,ii
         if ((i==itest).and.(j==jtest)) then
            print*,'TW: i/j/lon/lat:', i+i0,j+j0,plon(i,j),plat(i,j)
            print*,'TW: swh',swh(i,j)
            print*,'TW: mwp',mwp(i,j)
            print*,'TW: mwd',mwd(i,j)
            print*,'TW: mask',fldmask(i,j)
         end if
      enddo !i
      enddo !j
      end subroutine read_waves_ww3_arctic
       
#elif defined (WAVES_ERAI)
c --- ERA-I 1deg x 1deg global data from 
c --- http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/
c --- 6-hourly time res

      subroutine read_waves_erai(dtime_waves)
      use mod_hycom_nersc
      use netcdf
      use mod_common_wavesice, only:swh,mwp,mwd

      implicit none
      real*8 ,             intent(in)  :: dtime_waves
c
c --- Data catalogue name
      character(len=*), parameter :: wave_path='./WAVES/ERA-I_1deg/'
c
c --- Variable names - ERA40 0.125 by 0.125 
      character(len=*), parameter :: vmwd='mwd'
      character(len=*), parameter :: vmwp='mwp'
      character(len=*), parameter :: vswh='swh'

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=360, nlat=181
c
      integer  :: i,j,k,m,ix,jy,irec,maxrec,nlon2,nlat2
      real, dimension(nlon,nlat) :: tmpmwd, tmpmwp, tmpswh, wavemask

      real               :: eclon(nlon), eclat(nlat)
      real               :: dlon,dlat,flon,flat,llon,llat,rtmp
      real               :: missing,fij
      logical            :: crit
      integer            :: dimid, varid, ncid
      integer            :: iyear, ihour, iday
      character(len=80)  :: ncfil1
      character(len=4)   :: cyy, cyyp1
      character(len=10)  :: ctt
      integer, parameter :: itype=1 ! Sets interpolation type 0:bilinear, 1:bicubic
c
      include 'common_blocks.h'
c
      call forday(dtime_waves,yrflag, iyear,iday,ihour)
      irec  =  (iday-1)*4 + ihour/6 +1!iday starts at 1, not zero - HYCOM day is iday-1

c --- Set file names
      write(cyy  ,'(i4.4)') iyear
      write(cyyp1,'(i4.4)') iyear+1
      !!
      ncfil1=wave_path//'MWP_'//cyy//'.nc'   ! ERA-I
      if (mnproc==1) then
         print*,' '
         print*,'Reading ERA-I 12km waves from:',ncfil1
         print*,'dtime:',dtime_waves
         print*,'Year, HYCOM day, hour:',iyear,iday-1,ihour
         print*,'Record number:',irec
         print*,' '
      end if

c --- Check dimensions against hardcoded ones
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_dimid(ncid,"longitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlon2))
      call ncerr(nf90_inq_dimid(ncid,"latitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlat2))
      !call ncerr(nf90_inq_dimid(ncid,"time",dimid))
      !call ncerr(nf90_inquire_dimension(ncid,dimid,len=maxrec))
      
      if (nlon/=nlon2.or.nlat/=nlat2) then
         if (mnproc==1) then
            write(lp,'(a)') 'DD :: Dimension mismatch in ERAI data'
            write(lp,'(a,2i5)') 'nlon :',nlon,nlon2
            write(lp,'(a,2i5)') 'nlat :',nlat,nlat2
         end if
      end if

c --- Get longitude  & latitude for data  - use land mask file
      call ncerr(nf90_inq_varid(ncid,"longitude",varid))
      call ncerr(nf90_get_var  (ncid,varid,eclon))
      call ncerr(nf90_inq_varid(ncid,"latitude",varid))
      call ncerr(nf90_get_var  (ncid,varid,eclat))
      call ncerr(nf90_close    (ncid))
      flon=eclon(1); flat=eclat(1);
      dlon=eclon(2)-eclon(1); dlat=eclat(2)-eclat(1)
      if (mnproc==1) write(*,*) flon, flat, dlon, dlat

c --- ECMWF waves (every 6 hours)
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_varid(ncid,vmwp,varid))
      call ncerr(nf90_get_var  (ncid,varid,tmpmwp,start=(/1,1,irec/)))
      call ncerr(nf90_get_att  (ncid,varid,'missing_value',missing))
      call ncerr(nf90_close    (ncid))
      
      ! Build the data mask from missing values
      m = 1
      do i=1,nlon
      do j=1,nlat
         if ( tmpmwp(i,j)==missing ) then
            !wavemask(i,j)=0.
            wavemask(i,j)=0./0. !set to NaN
         else
            wavemask(i,j)=1.
         end if
      enddo
      enddo

      call interpug(wavemask,nlon,nlat,flon,flat,dlon,dlat,
     &              fldmask(1-nbdy,1-nbdy),plon,plat,0)

      do j=1-nbdy,jj+nbdy
      do i=1-nbdy,ii+nbdy
         !if (fldmask(i,j).lt.0.99998) then
         !   fldmask(i,j)=0.
         !else
         !   fldmask(i,j)=1.
         !end if
         fij   = fldmask(i,j)
         crit  = ((.not.(fij.lt.0.)).and.(.not.(fij.gt.0.)))!test for NaN
         if (crit) then
            fldmask(i,j)   = 0.
         else
            fldmask(i,j)   = 1.
         end if
      enddo
      enddo

      !mean wave period
      call rdncrec(ncfil1,vmwp,tmpmwp,nlon,nlat,irec)
      call interpug(tmpmwp,nlon,nlat,flon,flat,dlon,dlat,
     &             mwp(1-nbdy,1-nbdy),plon,plat,itype)


      !mean wave direction
      ncfil1=wave_path//'MWD_'//cyy//'.nc'   ! ERA-I
      if(mnproc==1) print*,ncfil1

      call rdncrec(ncfil1,vmwd,tmpmwd,nlon,nlat,irec)
      call interpug(tmpmwd,nlon,nlat,flon,flat,dlon,dlat,
     &             mwd(1-nbdy,1-nbdy),plon,plat,itype)
      !!convention is this is waves-from (0: from north, 90: from east)
      !!TODO: check does this data follow convention?


      !significant wave height
      ncfil1=wave_path//'SWH_'//cyy//'.nc'   ! ERA-I
      if(mnproc==1) print*,ncfil1

      call rdncrec(ncfil1,vswh,tmpswh,nlon,nlat,irec)
      call interpug(tmpswh,nlon,nlat,flon,flat,dlon,dlat,
     &             swh(1-nbdy,1-nbdy),plon,plat,itype)

      mwp   = fldmask*mwp
      mwd   = fldmask*mwd
      swh   = fldmask*swh

      if (mnproc==1) then
         print*,' '
         print*,'Finished reading ERA-I 12km waves'
         print*,' '
      end if
       
      end subroutine read_waves_erai

#else/* use old WAM waves */
      subroutine read_waves_ecmwf(dtime_waves)
      use mod_hycom_nersc
      use netcdf
      use mod_common_wavesice, only:swh,mwp,mwd

      implicit none
      real*8 ,             intent(in)  :: dtime_waves
c
c --- Data catalogue name
      character(len=*), parameter :: wave_path='./WAVES/WAM/'
c
c --- Variable names - ERA40 1.125 by 1.125 
      character(len=*), parameter :: vmwd='mwd'
      character(len=*), parameter :: vmwp='mwp'
      character(len=*), parameter :: vswh='swh' 

c --- Hardcoded grid sizes 
      integer, parameter :: nlon=240, nlat=121
c
      integer  :: i,j,k,m,ix,jy,irec,maxrec,nlon2,nlat2
      real, dimension(nlon,nlat) :: tmpmwd, tmpmwp, tmpswh, wavemask

      real               :: eclon(nlon), eclat(nlat)
      real               :: dlon,dlat,flon,flat,llon,llat,rtmp
      real               :: missing,fij
      logical            :: crit
      integer            :: dimid, varid, ncid
      integer            :: iyear, ihour, iday
      character(len=80)  :: ncfil1
      character(len=4)   :: cyy, cyyp1
      character(len=10)  :: ctt
      integer, parameter :: itype=1 ! Sets interpolation type 0:bilinear, 1:bicubic
c
      include 'common_blocks.h'
c
      call forday(dtime_waves,yrflag, iyear,iday,ihour)
      irec  =  (iday-1)*4 + ihour/6 +1!iday starts at 1, not zero - HYCOM day is iday-1

c --- Set file names
      write(cyy  ,'(i4.4)') iyear
      write(cyyp1,'(i4.4)') iyear+1
      !!
      ncfil1=wave_path//'waves_'//cyy//'_global.nc'   ! ERA40 !
      if (mnproc==1) then
         print*,' '
         print*,'Reading WAM waves from:',ncfil1
         print*,'dtime:',dtime_waves
         print*,'Year, HYCOM day, hour:',iyear,iday-1,ihour
         print*,'Record number:',irec
         print*,' '
      end if

c --- Check dimensions against hardcoded ones
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_dimid(ncid,"longitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlon2))
      call ncerr(nf90_inq_dimid(ncid,"latitude",dimid))
      call ncerr(nf90_inquire_dimension(ncid,dimid,len=nlat2))
      !call ncerr(nf90_inq_dimid(ncid,"time",dimid))
      !call ncerr(nf90_inquire_dimension(ncid,dimid,len=maxrec))
      
      if (nlon/=nlon2.or.nlat/=nlat2) then
         if (mnproc==1) then
            write(lp,'(a)') 'DD :: Dimension mismatch in ERAI data'
            write(lp,'(a,2i5)') 'nlon :',nlon,nlon2
            write(lp,'(a,2i5)') 'nlat :',nlat,nlat2
         end if
      end if

c --- Get longitude  & latitude for data  - use land mask file
      call ncerr(nf90_inq_varid(ncid,"longitude",varid))
      call ncerr(nf90_get_var  (ncid,varid,eclon))
      call ncerr(nf90_inq_varid(ncid,"latitude",varid))
      call ncerr(nf90_get_var  (ncid,varid,eclat))
      call ncerr(nf90_close    (ncid))
      flon=eclon(1); flat=eclat(1);
      dlon=eclon(2)-eclon(1); dlat=eclat(2)-eclat(1)
      if (mnproc==1) write(*,*) flon, flat, dlon, dlat

c --- ECMWF waves (every 6 hours)
      call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      call ncerr(nf90_inq_varid(ncid,vmwp,varid))
      call ncerr(nf90_get_var  (ncid,varid,tmpmwp,start=(/1,1,irec/)))
      call ncerr(nf90_get_att  (ncid,varid,'missing_value',missing))
      call ncerr(nf90_close    (ncid))
      
      ! Build the data mask from missing values
      m = 1
      do i=1,nlon
      do j=1,nlat
         if ( tmpmwp(i,j)==missing ) then
            !wavemask(i,j)=0.
            wavemask(i,j)=0./0. !set to NaN
         else
            wavemask(i,j)=1.
         end if
      enddo
      enddo

      call interpug(wavemask,nlon,nlat,flon,flat,dlon,dlat,
     &              fldmask(1-nbdy,1-nbdy),plon,plat,0)

      do j=1-nbdy,jj+nbdy
      do i=1-nbdy,ii+nbdy
         !if (fldmask(i,j).lt.0.99998) then
         !   fldmask(i,j)=0.
         !else
         !   fldmask(i,j)=1.
         !end if
         fij   = fldmask(i,j)
         crit  = ((.not.(fij.lt.0.)).and.(.not.(fij.gt.0.)))!test for NaN
         if (crit) then
            fldmask(i,j)   = 0.
         else
            fldmask(i,j)   = 1.
         end if
      enddo
      enddo

      call rdncrec(ncfil1,vmwp,tmpmwp,nlon,nlat,irec)
      call interpug(tmpmwp,nlon,nlat,flon,flat,dlon,dlat,
     &             mwp(1-nbdy,1-nbdy),plon,plat,itype)


      !call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      !call ncerr(nf90_inq_varid(ncid,vmwd,varid))
      !call ncerr(nf90_get_var  (ncid,varid,tmpmwd,start=(/1,1,irec/)))
      !call ncerr(nf90_close(ncid))
      call rdncrec(ncfil1,vmwd,tmpmwd,nlon,nlat,irec)
      call interpug(tmpmwd,nlon,nlat,flon,flat,dlon,dlat,
     &             mwd(1-nbdy,1-nbdy),plon,plat,itype)
      !!convention is this is waves-from (0: from north, 90: from east)
      !!TODO: check does this data follow convention?


      !call ncerr(NF90_open(trim(ncfil1),NF90_NOCLOBBER,ncid))
      !call ncerr(nf90_inq_varid(ncid,vswh,varid))
      !call ncerr(nf90_get_var(ncid,varid,tmpswh,start=(/1,1,irec/)))
      !call ncerr(nf90_close(ncid))
      call rdncrec(ncfil1,vswh,tmpswh,nlon,nlat,irec)
      call interpug(tmpswh,nlon,nlat,flon,flat,dlon,dlat,
     &             swh(1-nbdy,1-nbdy),plon,plat,itype)

      mwp   = fldmask*mwp
      mwd   = fldmask*mwd
      swh   = fldmask*swh

      if (mnproc==1) then
         print*,' '
         print*,'Finished reading WAM waves from:',ncfil1
         print*,' '
      end if
       
      end subroutine read_waves_ecmwf
#endif/*flag check to define wave model*/

#if defined (WAVES_WW3) || defined(WAVES_WAMNSEA) || defined (WAVES_WW3_ARCTIC)
      subroutine interp_nearest(data_out,data_in,nlon,nlat)
      !data_out on HYCOM grid; data_in on own grid (size=[nlon,nlat])
      implicit none

      integer, intent(in)  :: nlon,nlat
      real, dimension(nlon,nlat),intent(in) ::
     &   data_in 
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),intent(out)  ::
     &   data_out
      integer  :: i,j,ipiv,jpiv

      data_out = 0.
      do j= 1,jj
      do i= 1,ii
         !! pivots have come from read_pivots_ww3
         !! (this only needs to be done once, in mod_restart.F)
         ipiv = ipiv_all(i,j)
         jpiv = jpiv_all(i,j)
         if (fldmask(i,j).gt.0) then
            !! this is already waves-from direction so don't need to add 180
            !! convention is: 0 = from north, 90 = from east
            data_out(i,j) = data_in(ipiv,jpiv)
         endif
      enddo !i
      enddo !j
      end subroutine interp_nearest
#endif/*flag check for if pivots are needed*/

#endif/* WAVES */
      end module mod_readwaves