*DK GETGPI
      subroutine getgpi
c
c 5/94 (B. Foster):
c
c Get geophysical indices for current time, i.e. f107d, f107a, and Kp(8):
c   These data are in mss file /FOSTER/tgcminp/79001-92031.gpi.dat. 
c   Available dates are 79001 to 92031 (4779 days). This is an unformatted 
c   direct access f77 file, with one record per day containing 
c   iyd,f107d,f107a,rkp(8). (88 bytes/rec)
c Note this routine assumes the daily fluxes are at 12 ut and that the
c   8 3-hourly Kp values for each day are at ut 1.5, 4.5,..., 19.5, 22.5.
c
c If needed, hpower and ctpoten are calculated from Kp.
c (as of 5/94, byimf is not calculated in this routine; it must be provided
c  by the user in lexical reads in INPUT)
c
*CA PARAMS
*CA CONS
*CA STRT
      PARAMETER (NAURPX=59,NLEX=600)
      COMMON/AURDAT/ AURPS(NLEX,2,NAURPX), NPTS(NAURPX),IPR, NAURP, 
     1 PARAMV(NAURPX), JSWOLDA, JDITH, JDIDK
      parameter(spval=-999.) ! special value (must be same as in INPUT)
      parameter(irecl=88) ! iyd,f107d,f107a,rkp(8) * 8 (rec length in bytes)
      common/ingpi/ f107,f107a,ctpoten,hpower,byimf
      dimension fkp_prv(8),fkp_cur(8),fkp_nxt(8)
      character*8 diskfile
      character*40 mssfile
      character*128 errmsg
      data iydmin,iydmax/79001,92031/ ! first and last days of available data
      data mssfile/'/FOSTER/tgcminp/79001-92031.gpi.dat     '/
      data diskfile/'gpi.dat '/, lu/97/, mxrec/4779/
      data if107,if107a,ictpoten,ihpower,ibyimf/0,0,0,0,0/
      data f107_prv,f107_cur,f107_nxt/spval,spval,spval/
      data f107a_prv,f107a_cur,f107a_nxt/spval,spval,spval/
      data iyd_prv,iyd_cur,iyd_nxt/-1,-1,-1/
      save fkp_prv,fkp_cur,fkp_nxt 
c
c Statement function for linear interpolation:
c
      finterp(f0,f1,isec0,isec1,isec) = f0+(f1-f0)*float(isec-isec0)/
     +                                  float(isec1-isec0)
c
c There are 5 input parameters that may be defined by this routine:
c   f107,f107a,ctpoten,hpower, and byimf. If any of these were not
c   defined in INPUT, i.e., not provided by the user, then they 
c   were initialized to spval in INPUT before first iteration, and
c   will be defined from the I.S. database in this routine. Save
c   this info for later iterations: 
c (as of 5/94, byimf is always provided by the user, i.e., ibyimf will = 0)
c
      if (ifrst.eq.0) then
        if (f107.eq.spval) if107 = 1
        if (f107a.eq.spval) if107a = 1
        if (ctpoten.eq.spval) ictpoten = 1
        if (hpower.eq.spval) ihpower = 1
        if (ibyimf.eq.spval) ibyimf = 1
      endif
c
c If all 5 input parameters were defined by user, we have nothing
c   to do here:
c
      if (if107.eq.0.and.if107a.eq.0.and.ictpoten.eq.0.and.
     +    ihpower.eq.0.and.ibyimf.eq.0) return
c
c Determine current calendar day and time:
c
      iyd = (iyear-1900)*1000+iiday
      if (iyd.lt.iydmin.or.iyd.gt.iydmax-1) then
        write(6,"(' ')")
        write(6,"('>>> GETGPI: iyd outside NMC data ',
     +    'range: '/12x'desired iyd=',i5,', NMC data available for ',i5,
     +    ' to ',i5)") iyd,iyd0,iyd1
        stop 'GETGPI'
      endif
      nsec = ifix(amod(float(iter)*c(4),86400.))
      ih = nsec/3600
      nsec = mod(nsec,3600)
      im = nsec/60
      is = mod(nsec,60)
      isec_cur = ih*3600+im*60+is ! current ut in seconds
c
c Print to stdout if first iteration of run:
c
      if (ifrst.eq.0) then
        write(6,"(' ')")
        write(6,"('--------------------------------------------------',
     +    '-----------------------')")
        write(6,"('GETGPI:')")
        write(6,"(' ')")
        write(6,"('Will obtain the following inputs from I.S. (GPI)',
     +    ' database:')")
        if (if107.gt.0) write(6,"('F107',$)")
        if (if107a.gt.0) write(6,"('  F107A',$)")
        if (ihpower.gt.0) write(6,"('  POWER',$)")
        if (ictpoten.gt.0) write(6,"('  CTPOTEN',$)")
        if (ibyimf.gt.0) write(6,"('  BYIMF',$)") ! will not happen yet
        write(6,"(' ')")
        write(6,"('Data is on mss unformatted direct access file ',
     +    a)") mssfile(1:lenstr(mssfile))
        write(6,"('Interpolation starts at yyddd ',i5,' hr:min:sec ',
     +    i2,':',i2,':',i2,' (iter=',i8,')')") iyd,itime,iter
        write(6,"('--------------------------------------------------',
     +    '-----------------------')")
        write(6,"(' ')")
c
c Get data file (assume msread and open only on first iter):
c
        lmssfile = lenstr(mssfile)
        call msread(iier,diskfile,mssfile(1:lmssfile),' ',' ')
        if (iier.ne.3.and.iier.ne.0) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: error from msread: mssfile=',a)")
     +      mssfile(1:lmssfile)
          call mserror(errmsg)
          write(6,"('    mss error: ',a)") errmsg
          write(6,"(' ')")
          stop 'getgpi'
        endif
        open(lu,file=diskfile,access='DIRECT',status='OLD',
     +    form='UNFORMATTED',recl=irecl,err=900,iostat=istat)
c
c Read data for current, previous, and next days:
c
        irec = iyddif(iydmin,iyd_cur)+1
        if (irec.le.0.or.irec.gt.mxrec) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: bad irec for iyd_cur=',i5,' irec=',
     +      i5)") iyd_cur,irec
          write(6,"(' ')")
          stop 'getgpi'
        endif
        read(lu,rec=irec,err=901,iostat=istat) 
     +    iyd,f107d_cur,f107a_cur,fkp_cur
        if (iyd_cur.ne.iyd) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: iyd_cur.ne.iyd: iyd_cur=',i5,
     +      ' iyd=',i5,' irec=',i5)") iyd_cur,iyd,irec
          write(6,"(' ')")
          stop 'getgpi' 
        endif
c
        iyd_prv = inciyd(iyd_cur,-1)
        irec = iyddif(iydmin,iyd_prv)+1
        if (irec.le.0.or.irec.gt.mxrec) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: bad irec for iyd_prv=',i5,' irec=',
     +      i5)") iyd_prv,irec
          write(6,"(' ')")
          stop 'getgpi'
        endif
        read(lu,rec=irec,err=901,iostat=istat) 
     +    iyd,f107d_prv,f107a_prv,fkp_prv
        if (iyd_prv.ne.iyd) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: iyd_prv.ne.iyd: iyd_prv=',i5,
     +      ' iyd=',i5,' irec=',i5)") iyd_prv,iyd,irec
          write(6,"(' ')")
          stop 'getgpi' 
        endif
c
        iyd_nxt = inciyd(iyd_cur,1)
        irec = iyddif(iydmin,iyd_nxt)+1
        if (irec.le.0.or.irec.gt.mxrec) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: bad irec for iyd_nxt=',i5,' irec=',
     +      i5)") iyd_nxt,irec
          write(6,"(' ')")
          stop 'getgpi'
        endif
        read(lu,rec=irec,err=901,iostat=istat) 
     +    iyd,f107d_nxt,f107a_nxt,fkp_nxt
        if (iyd_nxt.ne.iyd) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: iyd_nxt.ne.iyd: iyd_nxt=',i5,
     +      ' iyd=',i5,' irec=',i5)") iyd_nxt,iyd,irec
          write(6,"(' ')")
          stop 'getgpi' 
        endif
      endif ! ifrst == 0
c
c Read new day if necessary:
c
      if (itimeq(iyd_cur,isec_cur,iyd_nxt,0).gt.0) then
        iyd_prv = inciyd(iyd_prv,1)
        f107_prv = f107_cur
        f107a_prv = f107a_cur
        fkp_prv(:) = fkp_cur(:)
        f107_cur = f107_nxt
        f107a_cur = f107a_nxt
        fkp_cur(:) = fkp_nxt(:)
        iyd_nxt = inciyd(iyd_nxt,1)
        irec = iyddif(iydmin,iyd_nxt)+1
        if (irec.le.0.or.irec.gt.mxrec) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: bad irec for iyd_nxt=',i5,' irec=',
     +      i5)") iyd_nxt,irec
          write(6,"(' ')")
          stop 'getgpi'
        endif
        read(lu,rec=irec,err=901,iostat=istat) 
     +    iyd,f107d_nxt,f107a_nxt,fkp_nxt
        if (iyd_nxt.ne.iyd) then
          write(6,"(' ')")
          write(6,"('>>> GETGPI: iyd_nxt.ne.iyd: iyd_nxt=',i5,
     +      ' iyd=',i5,' irec=',i5)") iyd_nxt,iyd,irec
          write(6,"(' ')")
          stop 'getgpi' 
        endif
      endif ! get new next day
c
c Interpolate fluxes:
c
      if (if107.gt.0.or.if107a.gt.0) then
        if (isec_cur.eq.43200) then ! ut=12, no interp necessary:
          if (if107.gt.0) f107 = f107_cur
          if (if107a.gt.0) f107a = f107a_cur
        elseif (isec_cur.lt.43200) then ! interp between prev and curr days
          if (if107.gt.0) 
     +      f107 = finterp(f107_prv,f107_cur,43200,129600,isec_cur)
          if (if107a.gt.0) 
     +      f107a = finterp(f107a_prv,f107a_cur,43200,129600,isec_cur)
        else ! interp between current and next days
          if (if107.gt.0) 
     +      f107 = finterp(f107_cur,f107_nxt,43200,129600,isec_cur)
          if (if107a.gt.0) 
     +      f107a = finterp(f107a_cur,f107a_nxt,43200,129600,isec_cur)
        endif
        if (if107.gt.0) then
          c(61) = f107
          npts(54) = 0
        endif
        if (if107a.gt.0) then
          c(62) = f107a
          npts(55) = 0
        endif
      endif ! doing f107 or f107a
c
c If need hem power or cross-tail potential, interpolate Kp to current
c time and use it for hpower and ctpoten calculations:
c
      if (ihpower.gt.0.or.ictpoten.gt.0.or.ibyimf.gt.0) then
        isecmin = ifix(1.5*3600.)
        isecmax = ifix(22.5*3600.)
        if (isec_cur.ge.isecmin.and.isec_cur.le.isecmax) then ! use curr day
          do i=1,7
            isec0 = ifix((float(i-1)*3.+1.5)*3600.)
            isec1 = isec0 + isecmin
            if (isec_cur.ge.isec0.and.isec_cur.le.isec1) then
              fkp0 = fkp_cur(i)
              fkp1 = fkp_cur(i+1)
              goto 100
            endif
          enddo
 100      continue
        elseif (isec_cur.lt.isecmin) then ! use previous day
          isec0 = -isecmin
          isec1 = isecmin
          fkp0 = fkp_prv(8)
          fkp1 = fkp_cur(1)
        else ! use next day
          isec0 = isecmax
          isec1 = isecmax+ifix(3.*3600.)
          fkp0 = fkp_cur(8)
          fkp1 = fkp_nxt(1)
        endif
        fkp = finterp(fkp0,fkp1,isec0,isec1,isec_cur)
        if (ihpower.gt.0) then
          hpower = amax1(0.,-2.78+9.33*fkp)
          paramv(1) = hpower
          npts(1) = 0
        endif
        if (ictpoten.gt.0) then
          ctpoten = 29.+11.*fkp
          paramv(2) = ctpoten
          npts(2) = 0
        endif
c
c If later calculate byimf, then it goes here:
c
c       if (ibyimf.gt.0) then
c         [put calculation of byimf here]
c         paramv(3) = byimf
c         npts(3) = 0
c       endif
c
      endif ! doing hpower, ctpoten, or byimf
c
c Write to stdout for testing:
c
      write(6,"('--------------------------------------------------',
     +  '-----------------------')")
      write(6,"('getgpi: iter=',i6,' iyd=',i5,' time=',i2,':',i2,':',i2,
     +  /' f107d=',f9.4,' f107a=',f9.4,' fkp=',f7.4,/,' ctpoten=',f8.4,
     +  ' hpower=',f8.4,' byimf=',f8.4)")
     +  iter,iyd,itime,f107,f107a,fkp,ctpoten,hpower,byimf
      return
c
c i/o error traps:
c
 900  continue
      write(6,"(' ')")
      write(6,"('>>> GETGPI: error opening file ',a,': ',
     +  'istat=',i5,' lu=',i2,' irecl=',i4)") diskfile,istat,lu,irecl
      write(6,"(' ')")
      stop 'getgpi'
 901  continue
      write(6,"(' ')")
      write(6,"('>>> GETGPI: error reading record ',i5,
     +  ' from file ',a,': istat=',i4,' lu=',i2)")
     +  irec,diskfile,istat,lu
      write(6,"(' ')")
      stop 'getgpi'
      end