!
      program rdgpi
      implicit none
      character(len=80) :: gpi_ncfile
      integer :: iyd0,iyd1
!
      gpi_ncfile = ' '
      write(gpi_ncfile,"('gpi_1979-2000.nc')")
      call rdfile(gpi_ncfile)
!
!     iyd0 = 1982001
!     iyd1 = 1982030
!     iyd0 = 1982001
!     iyd1 = 1982365
!     iyd0 = 1982075
!     iyd1 = 1982085
!     iyd0 = 1991150
!     iyd1 = 1991170
!     iyd0 = 1997080
!     iyd1 = 1997105
!     iyd0 = 1993080
!     iyd1 = 1993100
      iyd0 = 1993080
      iyd1 = 1993085
      call opngks
      call pltxy(iyd0,iyd1)
      call clsgks
      end
!-----------------------------------------------------------------------
      subroutine pltxy(iyd0,iyd1)
      implicit none
      include "gpi.h"
!
! Args:
      integer,intent(in) :: iyd0,iyd1
      character(len=80) :: toplab
!
! Local:
      integer :: i,ii,iii,ndays,loc0,loc1,icount
      integer,external :: ixfind
      character(len=32) :: xlab
      real :: xcoords(mxgpidays),power(8,mxgpidays),
     |  ctpoten(8,mxgpidays),x3hrinc(mxgpidays*8)
!
      loc0 = ixfind(gpi_iyd,ngpidays,iyd0,icount)
      loc1 = ixfind(gpi_iyd,ngpidays,iyd1,icount)
      if (loc0==0.or.loc1==0) then
        write(6,"(/,'>>> pltxy: iyd0=',i8,' loc0=',i4)") iyd0,loc0
        write(6,"(/,'           iyd1=',i8,' loc1=',i4)") iyd1,loc1
      endif
      ndays = loc1-loc0+1
      write(6,"('pltxy: iyd0=',i8,' iyd1=',i8,' ndays=',i4)")
     |  iyd0,iyd1,ndays
!
      do i=1,ndays
        xcoords(i) = float(i)
      enddo
!
! Daily f107:
      write(toplab,"('F10.7 CM DAILY FLUX (',i4,' DAYS FROM ',i8,
     +  ' TO ',i8,')')") ndays,iyd0,iyd1
      write(xlab,"(i4,' DAYS ',i8,' TO ',i8)") ndays,iyd0,iyd1
      call agsetc('LABEL/NAME.','B')
      call agseti('LINE/NUMBER.',-100)
      call agsetc('LINE/TEXT.',trim(xlab))
      call agsetc ('LABEL/NAME.','T')
      call agseti ('LINE/NUMBER.', 100)
      call agseti ('LINE/MAX.',80)
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','F107D')
      call ezxy(xcoords,gpi_f107d(loc0:loc1),ndays,trim(toplab))
!
! Average f107d:
      write(toplab,"('F10.7 CM AVERAGE FLUX (',
     |  i4,' DAYS FROM ',i8,' TO ',i8,')')") ndays,iyd0,iyd1
      write(xlab,"(i4,' DAYS ',i8,' TO ',i8)") ndays,iyd0,iyd1
      call agsetc('LABEL/NAME.','B')
      call agseti('LINE/NUMBER.',-100)
      call agsetc('LINE/TEXT.',trim(xlab))
      call agsetc ('LABEL/NAME.','T')
      call agseti ('LINE/NUMBER.', 100)
      call agseti ('LINE/MAX.',80)
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','F107A')
      call ezxy(xcoords,gpi_f107a(loc0:loc1),ndays,trim(toplab))
!
      iii = 0
      do i=1,ndays
        do ii=1,8 ! hours 0,3,6,...,24
          iii = iii+1
          x3hrinc(iii) = float(iii)
        enddo
      enddo
!
! Power:
      do i=1,ndays
        do ii=1,8
          power(ii,i) = max(0.,-2.78+9.33*gpi_kp(ii,loc0+i-1))
        enddo
      enddo 
      write(toplab,"('HEMISPHERIC POWER',
     |  '(',i4,' DAYS FROM ',i8,' TO ',i8,')')") ndays,iyd0,iyd1
      write(xlab,"(i4,' DAYS ',i8,' TO ',i8)") ndays,iyd0,iyd1
      call agsetc('LABEL/NAME.','B')
      call agseti('LINE/NUMBER.',-100)
      call agsetc('LINE/TEXT.',trim(xlab))
      call agsetc ('LABEL/NAME.','T')
      call agseti ('LINE/NUMBER.', 100)
      call agseti ('LINE/MAX.',80)
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','HPOWER')
      call ezxy(x3hrinc,power,ndays*8,trim(toplab))
!
! Ctpoten:
      do i=1,ndays
        do ii=1,8
          ctpoten(ii,i) = 29.+11.*gpi_kp(ii,loc0+i-1)
        enddo
      enddo 
      write(toplab,"('CROSS-CAP POTENTIAL (VOLTS)',
     |  '(',i4,' DAYS FROM ',i8,' TO ',i8,')')") ndays,iyd0,iyd1
      write(xlab,"(i4,' DAYS ',i8,' TO ',i8)") ndays,iyd0,iyd1
      call agsetc('LABEL/NAME.','B')
      call agseti('LINE/NUMBER.',-100)
      call agsetc('LINE/TEXT.',trim(xlab))
      call agsetc ('LABEL/NAME.','T')
      call agseti ('LINE/NUMBER.', 100)
      call agseti ('LINE/MAX.',80)
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','CTPOTEN')
      call ezxy(x3hrinc,ctpoten,ndays*8,trim(toplab))

      end subroutine pltxy
!-----------------------------------------------------------------------
      subroutine rdfile(gpi_ncfile)
!
! Obtain and read gpi_ncfile netcdf data file containing f107d, f107a, 
! and kp. These data are obtained in ascii format from NGDC. See comments 
! below.
!
      implicit none
      include "netcdf.inc"
      include "gpi.h"
!
! Args:
      character(len=*),intent(in) :: gpi_ncfile
!
! Local:
      integer :: istat,ncid,nkp,i
      integer :: id_ndays,id_nkp,idv_iyd,idv_f107d,idv_f107a,idv_kp
!
! Open netcdf file:
      call nc_open(ncid,trim(gpi_ncfile),'OLD','READ')
      if (ncid==0) then
        write(6,"(/,'>>> rdgpi: error opening netcdf gpi file ',
     |    a)") trim(gpi_ncfile)
        stop 'rdgpi'
      endif
!
! Get ndays dimension:
      istat = nf_inq_dimid(ncid,'ndays',id_ndays)
      istat = nf_inq_dimlen(ncid,id_ndays,ngpidays)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting ndays dimension')
      if (ngpidays > mxgpidays) then
        write(6,"(/,'>>> rdgpi: need to increase mxgpidays: ',
     |    'mxgpidays=',i8,' ngpidays=',i8)") mxgpidays,ngpidays
        stop 'rdgpi'
      endif
!
! Get kp dimension (nkp should be 8 for 3-hourly kp):
      istat = nf_inq_dimid(ncid,'nkp',id_nkp)
      istat = nf_inq_dimlen(ncid,id_nkp,nkp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting nkp dimension')
!
! Get gpi_iyd:
      istat = nf_inq_varid(ncid,'year_day',idv_iyd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting f107d var year_day')
      istat = nf_get_var_int(ncid,idv_iyd,gpi_iyd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting variable year_day')
      write(6,"('rdgpi: ndays=',i8,' (year_days ',i8,' to ',i8,')')") 
     |  ngpidays,gpi_iyd(1),gpi_iyd(ngpidays)
!     write(6,"('rdgpi: iyd=',/(8i9.7))") 
!    |  (gpi_iyd(i),i=1,ngpidays)
!
! Get f107d:
      istat = nf_inq_varid(ncid,'f107d',idv_f107d)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting f107d var id')
      istat = nf_get_var_real(ncid,idv_f107d,gpi_f107d)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting variable f107d')
!     write(6,"('rdgpi: f107d=',/(8f8.2))") 
!    |  (gpi_f107d(i),i=1,ngpidays)
!
! Get f107a:
      istat = nf_inq_varid(ncid,'f107a',idv_f107a)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting f107a var id')
      istat = nf_get_var_real(ncid,idv_f107a,gpi_f107a)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting variable f107a')
!     write(6,"('rdgpi: f107a=',/(8f8.2))") 
!    |  (gpi_f107a(i),i=1,ngpidays)
!
! Get kp:
      istat = nf_inq_varid(ncid,'kp',idv_kp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting kp var id')
      istat = nf_get_var_real(ncid,idv_kp,gpi_kp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting variable kp')
!     write(6,"('rdgpi: kp=',/(8f8.2))") 
!    |  (gpi_kp(:,i),i=1,ngpidays)
!
! Get beginning and ending integer dates:
! (these should be gpi_iyd(3) and gpi_iyd(ngpidays-2) respectively,
!  to enable interpolation using previous 2 days or last 2 days
!  at the extremes of the dates available).
!
      istat = nf_get_att_int(ncid,NF_GLOBAL,"yearday_beg",iyd_beg)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting global attribute yearday_beg')
      istat = nf_get_att_int(ncid,NF_GLOBAL,"yearday_end",iyd_end)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdgpi: Error getting global attribute yearday_end')
      write(6,"('rdgpi: iyd_beg=',i8,' iyd_end=',i8)") iyd_beg,iyd_end
!
! Close the file:
      call nc_close(ncid)
      write(6,"('rdgpi: completed read from gpi data file ',a)")
     |  trim(gpi_ncfile)
      end subroutine rdfile
!-----------------------------------------------------------------------
      subroutine nc_open(ncid,flnm,status,rdwr)
!     call nc_open(ncid,gpi_ncfile,'OLD','READ')
      implicit none
      include "netcdf.inc"
!
! Open or create a netcdf dataset flnm, returning netcdf file
!   id handle in ncid.
! If status=='NEW' or 'REPLACE', a new netcdf file is created,
!   otherwise (status=='OLD'), open existing dataset).
! If an existing file is opened and rdwr=='READ', the dataset is
!   opened read-only (NF_NOWRITE), otherwise (rdwr=='WRITE' or 
!   'APPEND') the dataset is opened read/write (NF_WRITE).
! If the open fails, print error msg and return ncid == 0.
!
! Args:
      integer,intent(out) :: ncid    ! netcdf file id
      character(len=*),intent(in) :: flnm,status,rdwr
!
! Local:
      integer :: istat
      character(len=120) :: char120
!
      ncid = 0
      if (status=='NEW'.or.status=='REPLACE') then
!
! Create new dataset (clobber old one if it pre-exists):
!
        istat = nf_create(flnm,NF_CLOBBER,ncid)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_create for netcdf ',
     |      'file ',a,' status=',a)") trim(flnm),status
          call handle_ncerr(istat,char120)
          ncid = 0
        else
          write(6,"('Created netcdf file ',a,' ncid=',i8)") 
     |      trim(flnm),ncid
        endif
!
! Open pre-existing dataset:
! NF_WRITE opens with read-write, where write means append, change, etc
! NF_NOWRITE is read-only.
!
      elseif (status=='OLD') then
        if (rdwr=='WRITE'.or.rdwr=='APPEND') then
          istat = nf_open(flnm,NF_WRITE,ncid)
        elseif (rdwr=='READ') then
          istat = nf_open(flnm,NF_NOWRITE,ncid)
        else
          write(6,"('>>> nc_open: unrecognized rdwr=',a)") trim(rdwr)
          write(6,"('  Will open NF_WRITE (read/write).')")
          istat = nf_open(flnm,NF_WRITE,ncid)
        endif
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_open for netcdf ',
     |      'file ',a,' status=',a)") trim(flnm),status
          call handle_ncerr(istat,char120)
          ncid = 0
        else
          write(6,"('Opened existing netcdf file ',a,' ncid=',i8)") 
     |      trim(flnm),ncid
        endif
!
! Unknown status request:
      else
        write(6,"(/,'>>> nc_open: unrecognized status = ',a)")
     |    status
        ncid = 0
      endif
      end subroutine nc_open
!-------------------------------------------------------------------
      subroutine nc_close(ncid)
      implicit none
      include "netcdf.inc"
!
! Close a netcdf dataset file:
!
! Arg:
      integer,intent(in) :: ncid
!
! Local:
      integer :: istat
!
      istat = nf_close(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error return from nf_close')
      end subroutine nc_close
!-------------------------------------------------------------------
      subroutine handle_ncerr(istat,msg)
!
! Handle a netcdf lib error:
!
      integer,intent(in) :: istat
      character(len=*),intent(in) :: msg
      include "netcdf.inc"
!
      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
!-------------------------------------------------------------------
      integer function ixfind(iarray,idim,itarget,icount)
!
! Search iarray(idim) for itarget, returning first index in iarray
! where iarray(idim)==target. Also return number of elements of
! iarray that == itarget in icount.
!
 
!
! Args:
      integer,intent(in) :: idim,itarget
      integer,intent(in) :: iarray(idim)
      integer,intent(out) :: icount
!
! Local:
      integer :: i
!
      ixfind = 0
      icount = 0
      if (.not.any(iarray==itarget)) return
      icount = count(iarray==itarget)
      do i=1,idim
        if (iarray(i)==itarget) then
          ixfind = i
          exit
        endif
      enddo
      end function ixfind
!-----------------------------------------------------------------------