!
      module see_flux_module
      use input_module,only: tempdir,see_ncfile
      implicit none
#ifdef SUN
#include "netcdf3.inc"
#else
#include "netcdf.inc"
#endif
      integer :: 
     |  nstruct,    ! # of structure_elements (unlimited dim)
     |  ndate,      ! # dates
     |  nwave       ! # bins
      integer,allocatable :: date(:)  ! date(ndate)
      real(kind=4),allocatable :: 
     |  wave1(:,:),    ! lower boundary of each bin (wave1(nwave,nstruct))
     |  wave2(:,:),    ! upper boundary of each bin (wave2(nwave,nstruct))
     |  spflux(:,:,:)  ! flux data (nwave,ndate,nstruct)
!
! Current SEE flux, interpolated to current model time:
      real,allocatable :: seeflux(:) ! seeflux(nwave) 
      real,allocatable :: sflux0(:),sflux1(:)
      contains
!-----------------------------------------------------------------------
      subroutine rdseeflux
!
! Read SEE flux nc file. This is called once per run from init.
!
      integer :: i,ncid,istat,ndims,nvars,ngatts,id_unlim,len,itype,
     |  nwave1,nwave2,ndateflux,nwaveflux
      integer ::
     |  idim_date, idv_date,        ! id's for date dimension and var
     |  idim_wave1, idv_wave1,      ! id's for wave1 dimension and var
     |  idim_wave2, idv_wave2,      ! id's for wave2 dimension and var
     |  idim1_spflux, idim2_spflux, ! dim id's for spflux
     |  idv_spflux
      character(len=120) :: diskfile
      character(len=120) :: attname,atttext
!
      write(6,"(/,72('-'),/,'Read SEE flux data:')")
      call mkdiskflnm(see_ncfile,diskfile)
      call getms(see_ncfile,diskfile,tempdir,' ')
!
! Open netcdf dataset:
      istat = nf_open(diskfile,NF_NOWRITE,ncid) 
      if (istat /= NF_NOERR) then
        write(6,"('>>> rdseeflux: error opening see_ncfile ',a)")
     |    trim(diskfile)
        call handle_ncerr(istat,'error opening see_ncfile')
      else
        write(6,"('Opened netCDF file ',a,' istat=',i3,' ncid=',
     |    i3)") trim(diskfile),istat,ncid
      endif
!
! Get number of dims, vars, atts. Print global attributes to stdout.
      istat = nf_inq(ncid,ndims,nvars,ngatts,id_unlim)
      write(6,"(/,'Global attributes from file ',a,':')") 
     |  trim(diskfile)
      do i=1,ngatts
        attname = ' ' ; atttext = ' '
        istat = nf_inq_attname(ncid,NF_GLOBAL,i,attname)
        istat = nf_inq_att(ncid,NF_GLOBAL,attname,itype,len)
        if (itype==NF_CHAR) then
          istat = nf_get_att_text(ncid,NF_GLOBAL,attname,atttext)
          write(6,"('  ',a,': ',a)") trim(attname),trim(atttext)
        endif
      enddo
!
! Get length of unlimited dimension (nstruct):
      istat = nf_inq_dimlen(ncid,id_unlim,nstruct)
!
! Get ndate dimension (dim1_DATE):
      istat = nf_inq_dimid(ncid,'dim1_DATE',idim_date)
      istat = nf_inq_dimlen(ncid,idim_date,ndate)
!
! Get nwave1 dimension (dim1_WAVE1):
      istat = nf_inq_dimid(ncid,'dim1_WAVE1',idim_wave1)
      istat = nf_inq_dimlen(ncid,idim_wave1,nwave1)
!
! Get nwave2 dimension (dim1_WAVE2):
      istat = nf_inq_dimid(ncid,'dim1_WAVE2',idim_wave2)
      istat = nf_inq_dimlen(ncid,idim_wave2,nwave2)
!
! Get dimensions for spflux:
      istat = nf_inq_dimid(ncid,'dim1_SP_FLUX',idim1_spflux)
      istat = nf_inq_dimlen(ncid,idim1_spflux,nwaveflux)
      istat = nf_inq_dimid(ncid,'dim2_SP_FLUX',idim2_spflux)
      istat = nf_inq_dimlen(ncid,idim2_spflux,ndateflux)
!
!     write(6,"(/,'rdseeflux: nstruct=',i3,' ndate=',i3,' nwave1=',i3,
!    |  ' nwave2=',i3,' nwaveflux=',i3,' ndateflux=',i3)")
!    |  nstruct,ndate,nwave1,nwave2,nwaveflux,ndateflux
!
! Check dimensions:
      if (ndate /= ndateflux) then
        write(6,"('>>> rdseeflux: ndate /= ndateflux: ndate=',i4,
     |    ' ndateflux=',i4)") ndate,ndateflux
        stop 'rdseeflux'
      endif
      if (nwave1 /= nwave2) then
        write(6,"('>>> rdseeflux: nwave1 /= nwave2: nwave1=',i4,
     |    ' nwave2=',i4)") nwave1,nwave2
        stop 'rdseeflux'
      endif
      nwave = nwave1
      if (nwaveflux /= nwave) then
        write(6,"('>>> rdseeflux: nwaveflux /= nwave: nwaveflux=',i4,
     |    ' nwave=',i4)") nwaveflux,nwave
        stop 'rdseeflux'
      endif
!
! Allocate and read date variable:
      if (allocated(date)) deallocate(date)
      allocate(date(ndate),stat=istat)
      if (istat /= 0) write(6,"('>>> rdseeflux: error allocating ',
     |  ' date var: ndate=',i4)") ndate
      istat = nf_inq_varid(ncid,'DATE',idv_date)
      istat = nf_get_var_int(ncid,idv_date,date)
!     write(6,"('date=',/,(8i8))") date
      write(6,"('DATE min,max=',2i8)") minval(date),maxval(date)
!
! Allocate and read wave1:
      nwave = nwave1
      if (allocated(wave1)) deallocate(wave1)
      allocate(wave1(nwave,nstruct),stat=istat)
      if (istat /= 0) write(6,"('>>> rdseeflux: error allocating ',
     |  ' wave1 var: nwave=',i4,' nstruct=',i4)") nwave,nstruct
      istat = nf_inq_varid(ncid,'WAVE1',idv_wave1)
      istat = nf_get_var_real(ncid,idv_wave1,wave1)
!     write(6,"('wave1=',/,(6e12.4))") wave1
      write(6,"('WAVE1 min,max=',2e12.4)") minval(wave1),maxval(wave1)
!
! Allocate and read wave2:
      if (allocated(wave2)) deallocate(wave2)
      allocate(wave2(nwave,nstruct),stat=istat)
      if (istat /= 0) write(6,"('>>> rdseeflux: error allocating ',
     |  ' wave2 var: nwave=',i4,' nstruct=',i4)") nwave,nstruct
      istat = nf_inq_varid(ncid,'WAVE2',idv_wave2)
      istat = nf_get_var_real(ncid,idv_wave2,wave2)
!     write(6,"('wave2=',/,(6e12.4))") wave2
      write(6,"('WAVE2 min,max=',2e12.4)") minval(wave2),maxval(wave2)
!
! Allocate and read spflux:
      if (allocated(spflux)) deallocate(spflux)
      allocate(spflux(nwave,ndateflux,nstruct),stat=istat)
      if (istat /= 0) write(6,"('>>> rdseeflux: error allocating ',
     |  ' spflux var: nwave=',i4,' ndateflux=',i4,' nstruct=',i3)") 
     |  nwave,ndateflux,nstruct
      istat = nf_inq_varid(ncid,'SP_FLUX',idv_spflux)
      istat = nf_get_var_real(ncid,idv_spflux,spflux)
!     write(6,"('spflux=',/,(6e12.4))") spflux
      write(6,"('SP_FLUX min,max=',2e12.4)") 
     |  minval(spflux),maxval(spflux)
!
! Allocate array for current see flux:
      if (allocated(seeflux)) deallocate(seeflux)
      allocate(seeflux(nwave),stat=istat)
      if (istat /= 0) then
        write(6,"('>>> rdseeflux: error allocating ',
     |    ' seeflux var: nwave=',i4)") nwave
      endif
      if (allocated(sflux0)) deallocate(sflux0)
      allocate(sflux0(nwave),stat=istat)
      if (istat /= 0) then
        write(6,"('>>> rdseeflux: error allocating ',
     |    ' sflux0 var: nwave=',i4)") nwave
      endif
      if (allocated(sflux1)) deallocate(sflux1)
      allocate(sflux1(nwave),stat=istat)
      if (istat /= 0) then
        write(6,"('>>> rdseeflux: error allocating ',
     |    ' sflux1 var: nwave=',i4)") nwave
      endif
!
      write(6,"('Finished reading SEE flux data',/,72('-'))")
      end subroutine rdseeflux
!-------------------------------------------------------------------
      subroutine getseeflux(iyear,iday,iutsec,iprint)
!
! Args:
      integer,intent(in) :: iyear,iday,iutsec,iprint
!
! Local:
      integer :: i,iyd,ixdate,isec0,isec1
      integer,parameter :: isec12 = 12*60*60   ! secs at ut 12
      integer,parameter :: istruct = 1
!
!     write(6,"('getseeflux: iyear=',i5,' iday=',i4,' iutsec=',
!    |  i8)") iyear,iday,iutsec
!
! Check that requested date is available:
      iyd = iyear*1000+iday
      if (iyd < date(1) .or. iyd > date(ndate)) then
        write(6,"(/,'>>> getseeflux: requested year-day is not ',
     |    'available from see_ncfile data file:')")
        write(6,"(4x,'Requested year-day=',i8,' date(1)=',i8,
     |    ' date(ndate)=',i8,/)") iyd,date(1),date(ndate)
        write(6,"(72('-'),/)")
        stop 'SEE'
      endif
!
! Get index to requested date:
      ixdate = 0
      dateloop: do i=2,ndate
        if (iyd <= date(i) .and. iyd >= date(i-1)) then
          ixdate = i
          exit dateloop 
        endif
      enddo dateloop
      if (ixdate == 0) then
        write(6,"('>>> getseeflux: could not bracket iyd ',i5,
     |    ' date=',/,(8i8))") iyd,date
        stop 'SEE'
      endif
!     write(6,"('getseeflux: iyd=',i8,' ixdate=',i8,' date(ixdate)=',
!    |  i8)") iyd,ixdate,date(ixdate)
!
! If current ut is after 12 ut, use next iyd, otherwise use prev iyd:
      if (iutsec==isec12) then      ! no interp necessary
        seeflux(:) = spflux(:,ixdate,istruct)
      elseif (iutsec < isec12) then ! ut < 12 -> use prev iyd
        if (ixdate==1) then
          write(6,"('>>> getseeflux: ixdate=',i2,' -> need ',
     |      'previous day data for interpolation.')") ixdate
          stop 'SEE'
        endif
        isec0 = 0
        isec1 = isec12*(date(ixdate)-date(ixdate-1))
        sflux0(:) = spflux(:,ixdate-1,istruct) ! note spflux is kind=4
        sflux1(:) = spflux(:,ixdate  ,istruct)

!       write(6,"('date(ixdate-1)=',i8,' isec0=',i8,' sflux0=',/,
!    |    (6e12.4))") date(ixdate-1),isec0,sflux0
!       write(6,"('date(ixdate)  =',i8,' isec1=',i8,' sflux1=',/,
!    |    (6e12.4))") date(ixdate),isec1,sflux1

      else                            ! ut > 12 -> use next iyd
        if (ixdate==ndate) then
          write(6,"('>>> getseeflux: ixdate=',i4,' ndate=',i4,
     |      ' -> need next day data for interpolation.')") ixdate,
     |      ndate
          stop 'SEE'
        endif
        isec0 = isec12
        isec1 = isec12+isec12*(date(ixdate+1)-date(ixdate))
        sflux0(:) = spflux(:,ixdate  ,istruct) ! note spflux is kind=4
        sflux1(:) = spflux(:,ixdate+1,istruct)

!       write(6,"('date(ixdate)  =',i8,' isec0=',i8,' sflux0=',/,
!    |    (6e12.4))") date(ixdate),isec0,sflux0
!       write(6,"('date(ixdate+1)=',i8,' isec1=',i8,' sflux1=',/,
!    |    (6e12.4))") date(ixdate+1),isec1,sflux1
        
      endif
      call finterpa(sflux0,sflux1,seeflux,nwave,isec0,isec1,iutsec)
!     write(6,"('getseeflux: iyd=',i8,' iutsec=',i8,' seeflux=',/,
!    |  (6e12.4))") iyd,iutsec,seeflux

      end subroutine getseeflux
!-------------------------------------------------------------------
      subroutine handle_ncerr(istat,msg)
!
! Handle a netcdf lib error:
!
! (Cannot use-associate this from nchist because of circular dependency
!  problems)
!
      integer,intent(in) :: istat
      character(len=*),intent(in) :: msg
!
      write(6,"(/72('-'))")
      write(6,"('>>> Error from netcdf library:')")
      write(6,"(a)") trim(msg)
      write(6,"('istat=',i5)") istat
      write(6,"(a)") nf_strerror(istat)
      write(6,"(72('-')/)")
      return
      end subroutine handle_ncerr
!-----------------------------------------------------------------------
      end module see_flux_module