c
      program rdgpi
c
      include 'rdgpi.h'
      dimension rkp(8),iasf(13),xdayinc(nrec),x3hrinc(nrec*8)
      integer itlexda(4,nrec),itlexhr(4,nrec*8)
      character*56 dskfile,xlab
      character*80 toplab 
      data lugpi/20/, ludat/21/, irecl/88/, iasf/13*1/, iframe/0/
!     data iyd0/79001/
      data iyd0/1979001/
c
      call getinp
c
c Set up gks if making plots:
c
      if (ipltline.gt.0) then
        call opngks
        call gsclip(0)
        call gsasf(iasf)
        call gsfais(1)
      endif
c
c Open file to receive data statments:
c
      if (iwrlex.gt.0) then
        open(ludat,file='rdgpi.lex',err=901)
      endif
c
      call getmss(mssfile)
      call tail(mssfile,dskfile)
      ldsk = lenstr(dskfile)
      open(lugpi,file=dskfile(1:ldsk),form='UNFORMATTED',
     +  access='DIRECT',recl=irecl,status='OLD',err=900)
      write(6,"('Opened file ',a,' lugpi=',i2)") dskfile(1:ldsk),lugpi
c
      irec0 = igetrec(iyd0,iyds(1))
      irec1 = igetrec(iyd0,iyds(2))
      ndays = irec1-irec0+1
      iyd = iyds(1)
      mday1 = iyds(1)-(iyds(1)/1000*1000)
c
c Rec loop:
c Note iydrd is read from the file and checked against the desired
c   iyd to confirm the determined record.
c
      iydcur = iyds(1)
      nr = 0
      write(6,"('rdgpi: reading records ',i5,' to ',i5,' (ndays=',
     +  i4,' from ',i8,' to ',i8,')')") irec0,irec1,ndays,iyds
      write(6,"('  mday1=',i4)") mday1
      do irec=irec0,irec1
        nr = nr+1
        read(lugpi,rec=irec) iydrd,f107d(nr),f107a(nr),rkp
        if (iydrd.ne.iydcur) then ! make sure we got desired iyd
          write(6,"('>>> rdgpi: iydrd != iydcur: iydrd=',i8,' iydcur=',
     +      i8,' irec=',i5,' nr=',i5)") iydrd,iydcur,irec,nr
          stop 'iydrd'
        endif
        call iyd2date(iydrd,imo,ida,iyr)
        iyd(nr) = iydrd
        do i=1,8
          itlexhr(1,(nr-1)*8+i) = iyr			! yyyy (e.g., 1972)
          itlexhr(2,(nr-1)*8+i) = imo*100+ida		! moda (e.g., 0230)
          itlexhr(3,(nr-1)*8+i) = ((i-1)*3+1)*100+30	! midpt uthr+min
          itlexhr(4,(nr-1)*8+i) = 0			! csecs
        enddo
        itlexda(1,nr) = iyr
        itlexda(2,nr) = imo*100+ida
        itlexda(3,nr) = 1200
        itlexda(4,nr) = 0

!       write(6,"('irec=',i6,' nr=',i3,' iydrd=',i5,' ida=',i5,
!    +    ' itlexda(2,nr)=',i6)") irec,nr,iydrd,ida,itlexda(2,nr)
!
!       write(6,"('irec=',i6,' nr=',i3,' iydrd=',i8,' imo=',i5,
!    |    ' ida=',i2,' iyr=',i4,' f107d=',f8.3,' f107a=',f8.3)") 
!    |    irec,nr,iydrd,imo,ida,iyr,f107d(nr),f107a(nr)

c
c Calculate hpower and ctpoten from Kp at hrs 0,3,6,...,24:
c
c       fkp(1,nr) = rkp(1)
c       do i=2,8
c         fkp(i,nr) = .5*(rkp(i-1)+rkp(i)) 
c       enddo
c       fkp(9,nr) = rkp(8)
c       do i=1,9
        do i=1,8
          fkp(i,nr) = rkp(i)
          hpower(i,nr) = amax1(0., -2.78+9.33*fkp(i,nr))
          ctpoten(i,nr) = 29.+11.*fkp(i,nr)
!         write(6,"('  i=',i1,' hour ',i3,' to ',i3,' fkp=',f7.3,
!    |      ' hpower=',f7.3,' ctpoten=',f7.3)") i,(i-1)*3,(i-1)*3+3,
!    |      fkp(i,nr),hpower(i,nr),ctpoten(i,nr)
        enddo
c
!       call inciyd(iydcur,1) ! update next desired iyd
        iydcur = inciyd(iydcur,1)
      enddo
      write(6,"('Completed reads:')")
      write(6,"('iyd=',/(12i6))") (iyd(i),i=1,ndays)
      write(6,"('f107d=',/(8f9.3))") (f107d(i),i=1,ndays)
      write(6,"('f107a=',/(8f9.3))") (f107a(i),i=1,ndays)
      write(6,"('Kp=',/(14f5.2))") ((fkp(i,ii),i=1,8),ii=1,ndays)
      write(6,"('hpower=',/(14f7.3))") 
     |  ((hpower(i,ii),i=1,8),ii=1,ndays)
      write(6,"('ctpoten=',/(14f7.3))") 
     |  ((ctpoten(i,ii),i=1,8),ii=1,ndays)
c
c Set up for line plots:
c xdayinc = days of the year, increasing across year boundaries
c x3hrinc = decimal days every 3-hours, increasing across year boundaries
c
      if (ipltline.gt.0.or.iwrlex.gt.0) then
        xdayinc(1) = float(iyd(1)-iyd(1)/1000*1000)
        x3hrinc(1) = float(iyd(1)-iyd(1)/1000*1000)
        do i=1,ndays
          if (i.gt.1) xdayinc(i) = xdayinc(i-1)+1. ! monot increasing day
          do ii=1,8 ! hours 0,3,6,...,24
            x3hrinc((i-1)*8+ii) = xdayinc(i)+float((ii-1)*3)/24.
          enddo
        enddo
        iyr0 = iyds(1)/1000
        iyr1 = iyds(2)/1000
        if (iyr0.eq.iyr1) then ! no year boundaries
          write(xlab,"('DAY OF ',i4)") iyr0
        else ! at least one year boundary
          write(xlab,"('DAY OF YEARS ',i4,' TO ',i4)")
     +      iyr0,iyr1
        endif
        if (ipltline.gt.0) then
          call agsetc('LABEL/NAME.','B')
          call agseti('LINE/NUMBER.',-100)
          call agsetc('LINE/TEXT.',xlab(1:lenstr(xlab)))
          call agsetc ('LABEL/NAME.','T')
          call agseti ('LINE/NUMBER.', 100)
          call agseti ('LINE/MAX.',80)
          write(6,"(' ')")
          write(6,"('Line plots for the period ',i8,' to ',i8,
     +      ' (',i4,' days)')") iyds,ndays
        endif
        if (ndays.le.1) goto 100
c
c f107d:
c
        ip = 1
        call clearstr(toplab)
        write(toplab,"('F10.7 CM DAILY FLUX (',i4,' DAYS FROM ',i8,
     +    ' TO ',i8,')')") ndays,iyds    
        if (ipltline.gt.0.and.ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','F107D')
          call ezxy(xdayinc,f107d,ndays,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': f107d')") iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'F107 = ',1,
     +      f107d,itlexda,ndays,4,mday1)
          write(6,"('Wrote lex read statements of f107d')")
        endif
c
c f107a:
c
        ip = 2
        call clearstr(toplab)
        write(toplab,"('F10.7 CM 81-DAY MEAN FLUX (',i4,' DAYS FROM ',
     +    i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','F107A')
          call ezxy(xdayinc,f107a,ndays,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': f107a')") iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'F107A = ',1,
     +      f107a,itlexda,ndays,4,mday1)
          write(6,"('Wrote lex read statements of f107a')")
        endif
 100    continue
c
c 3-hr Kp:
c
        ip = 3
        call clearstr(toplab)
        write(toplab,"('Kp (every 3 hours) (',i4,' DAYS FROM ',
     +    i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','3-hour Kp')
          call ezxy(x3hrinc,fkp,ndays*8,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': 3-hour Kp')") iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'KP = ',1,
     +      fkp,itlexhr,ndays*8,4,mday1)
          write(6,"('Wrote lex read statements of 3-hr Kp')")
        endif
c
c Daily mean Kp:
c
        if (ndays.le.1) goto 105
        ip = 4
        call clearstr(toplab)
        do i=1,ndays
          fkp_daily(i) = 0.
          do ii=1,8
            fkp_daily(i) = fkp_daily(i)+fkp(ii,i)
          enddo
          fkp_daily(i) = fkp_daily(i)/8.
        enddo
        write(toplab,"('Kp (daily mean) (',i4,' DAYS FROM ',
     +    i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','Kp (daily mean)')
          call ezxy(xdayinc,fkp_daily,ndays,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': daily mean Kp')") iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'KP = ',1,
     +      fkp_daily,itlexda,ndays,4,mday1)
          write(6,"('Wrote lex read statements of daily Kp')")
        endif
 105    continue
c
c 3-hr ctpoten:
c       ctpoten(:,nr) = 29.+11.*fkp(:,nr)
c
        ip = 5
        call clearstr(toplab)
        write(toplab,"('CTPOTEN ((29+11*Kp), every 3 hours) (',
     +    i4,' DAYS FROM ',i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','3-hour ctpoten')
          call ezxy(x3hrinc,ctpoten,ndays*8,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': 3-hour ctpoten (29+11*Kp))')") 
     +      iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'CTPOTEN = ',1,
     +      ctpoten,itlexhr,ndays*8,4,mday1)
          write(6,"('Wrote lex read statements of 3-hr ctpoten')")
        endif
c
c Daily mean Cp:
c
        if (ndays.le.1) goto 110
        ip = 6
        call clearstr(toplab)
        do i=1,ndays
          cp_daily(i) = 0.
          do ii=1,8
            cp_daily(i) = cp_daily(i)+ctpoten(ii,i)
          enddo
          cp_daily(i) = cp_daily(i)/8.
        enddo
        write(toplab,"('CTPOTEN ((29+11*Kp), daily mean) (',
     +    i4,' DAYS FROM ',i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','CTPOTEN (daily mean)')
          call ezxy(xdayinc,cp_daily,ndays,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': daily mean ctpoten (29+11*Kp))')") 
     +      iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'CTPOTEN = ',1,
     +      cp_daily,itlexda,ndays,4,mday1)
          write(6,"('Wrote lex read statements of daily Kp')")
        endif
 110    continue
c
c Line plot 3-hr hpower:
c       hpower(:,nr) = amax1(0., -2.78+9.33*fkp(:,nr))
c
        ip = 7
        call clearstr(toplab)
        write(toplab,"('HEM POWER ((-2.78+9.33*Kp), every 3 hours) (',
     +    i4,' DAYS FROM ',i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','3-hour HPOWER')
          call ezxy(x3hrinc,hpower,ndays*8,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': 3-hour hpower (-2.78+9.33*Kp))')") 
     +      iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'POWER = ',2,
     +      hpower,itlexhr,ndays*8,4,mday1)
          write(6,"('Wrote lex read statements of 3-hr hpower')")
        endif
c
c Daily mean Hp:
c
        if (ndays.le.1) goto 115
        ip = 8
        call clearstr(toplab)
        do i=1,ndays
          hp_daily(i) = 0.
          do ii=1,8
            hp_daily(i) = hp_daily(i)+hpower(ii,i)
          enddo
          hp_daily(i) = hp_daily(i)/9.
        enddo
        write(toplab,"('HEM POWER ((-2.78+9.33*Kp), daily mean) (',
     +    i4,' DAYS FROM ',i5,' TO ',i5,')')") ndays,iyds    
        if (ifplt(ip).gt.0) then
          call agsetc('LABEL/NAME.','L')
          call agseti('LINE/NUMBER.',100)
          call agsetc('LINE/TEXT.','HEM POWER (daily mean)')
          call ezxy(xdayinc,hp_daily,ndays,toplab(1:lenstr(toplab)))
          iframe = iframe+1
          write(6,"('Frame ',i2,': daily mean hpower ',
     +      '(-2.78+9.33*Kp))')") iframe 
        endif
        if (iflex(ip).gt.0) then
          write(toplab,"('C ',a)") toplab(1:lenstr(toplab))
          call lexwrt(ludat,toplab(1:lenstr(toplab)),'POWER = ',2,
     +      hp_daily,itlexda,ndays,4,mday1)
          write(6,"('Wrote lex read statements of daily hpower')")
        endif
 115    continue
        if (ipltline.gt.0) call clsgks
        if (iwrlex.gt.0) close(ludat)
      endif
c
c Send files (ipltline is > 0 only if sendcgm has length)
c            (iwrlex   is > 0 only if sendlex has length)
c
      if (ipltline.gt.0) call rcpfile(0,'gmeta',sendcgm)
      if (iwrlex.gt.0) call rcpfile(0,'rdgpi.lex',sendlex)
      stop 'done'
 900  continue
      write(6,"('>>> Error opening file ',a)") dskfile(1:ldsk)
      stop 'open'
 901  continue
      write(6,"('>>> Error opening file rdgpi.lex to receive ',
     +  'data statements')") 
      stop 'open'
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function igetrec(iyd0,iyd)
c
c Given that a file contains one record per day, beginning at
c iyd0, return the record corresponding to iyd
c (iyd0 and iyd are yyddd (year-day) format)
c
      igetrec = 0
      ier = ichiyd(iyd0)
      if (ier.ne.0) then
        write(6,"('>>> igetrec: bad iyd0=',i6)") iyd0
        return
      endif
      ier = ichiyd(iyd)
      if (ier.ne.0) then
        write(6,"('>>> igetrec: bad iyd=',i6)") iyd
        return
      endif
      if (iyd.lt.iyd0) then
        write(6,"('>>> igetrec: bad iyd=',i6,
     +   ' (iyd must be >= iyd0 and iyd0=',i5)')") iyd,iyd0
         return
      endif
      igetrec = iyddif(iyd0,iyd)+1
c     write(6,"('igetrec returning with rec=',i5,' (iyd=',i5,')')")
c    +  igetrec,iyd
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine agchnl(iaxs,vils,chrm,mcim,ncim,ipxm,
     +  chre,mcie,ncie)
      character*(*) chrm,chre
      include 'rdgpi.h'
      data ncalls/0/
c
      if (iaxs.eq.3) then
        ncalls = ncalls+1
        if (ndays.le.7) return
        v = vils
        if (vils.gt.365.) then
          ny = ifix(vils)/365
          iyr = iyds(1)/1000+ny+1900
          if (mod(iyr,4).ne.0) then ! not a leap year
            v = amod(vils,365.)
          else
            if (vils.gt.366.) v = amod(vils,366.)
          endif
c         write(6,"('vils > 365: vils=',f5.1,' v=',f5.1)") vils,v
        endif
        write(chrm,"(i3)") int(v)
        ncim = 5
c       write(6,"('ncalls=',i3,' vils=',f5.0,' chrm=',a)") 
c    +    ncalls,vils,chrm
      endif
      return
      end