c
      subroutine getgpi
c
c 5/94 (B. Foster):
c
c Get geophysical indices for current time, i.e. f107, f107a, and Kp(8):
c   These data are in mss file /FOSTER/tgcminp/79001-94080.gpi.dat.
c   Available dates are 79001 to 94080 (5559 days). This is an unformatted
c   direct access f77 file, with one record per day containing
c   iyd,f107,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
      include "params.h"
      include "cons.h"
      include "strt.h"
      include "sechis.h"
      PARAMETER (NAURPX=55,NLEX=800)
      COMMON/AURDAT/ AURPS(NLEX,2,NAURPX), NPTS(NAURPX),IPR, NAURP,
     1 PARAMV(NAURPX), JSWOLDA, JDITH, JDIDK
! spval is in sechis.h
!     parameter(spval=-999.) ! special value (must be same as in INPUT)
      parameter(irecl=88) ! iyd,f107,f107a,rkp(8) * 8 (rec length in bytes)
      common/ingpi/ f107,f107a,ctpoten,hpower,byimf,fkp
      dimension fkp_p2d(8),fkp_prv(8),fkp_cur(8),fkp_nxt(8),fkp_n2d(8)
      character*8 diskfile
      character*40 mssfile
      character*128 errmsg
      data iydmin,iydmax/79001,94080/ ! first and last days of available data
      data mssfile/'/FOSTER/tgcminp/79001-94080.gpi.dat     '/
      data diskfile/'gpi.dat '/, lu/97/, mxrec/5559/
      data if107,if107a,ictpoten,ihpower,ibyimf/0,0,0,0,0/
      data f107_p2d,f107_prv,f107_cur,f107_nxt,f107_n2d
     +  /spval,spval,spval,spval,spval/
      data f107a_p2d,f107a_prv,f107a_cur,f107a_nxt,f107a_n2d
     +  /spval,spval,spval,spval,spval/
      data iyd_p2d,iyd_prv,iyd_cur,iyd_nxt,iyd_n2d/-1,-1,-1,-1,-1/
      data isecf107/43200/ ! assume ut 12 for daily fluxes
      save fkp_p2d,fkp_prv,fkp_cur,fkp_nxt,fkp_n2d
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 (iyear and iiday are in a model common /STRT/)
c
      iyd_cur = (iyear-1900)*1000+iiday
      if (iyd_cur.lt.iydmin.or.iyd_cur.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_cur,iydmin,iydmax
        stop 'GETGPI'
      endif
      nsec = ifix(amod(float(iter)*c(4),86400.))
      call isec2hms(nsec,ih,im,is)
      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,"(/72('-'))")
        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:len_trim(mssfile))
        write(6,"('Interpolation starts at yyddd ',i5,' hr:min:sec ',
     +    i2,':',i2,':',i2,' (iter=',i8,')')") iyd_cur,ih,im,is,iter
        write(6,"(72('-')/)")
c
c Get data file (assume msread and open only on first iter):
c
        call getms(mssfile,diskfile,tmpdir,' ')
        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
        call rdgpi(lu,iyd_cur,iydrd,f107_cur,f107a_cur,fkp_cur)
        iyd_prv = inciyd(iyd_cur,-1)
        call rdgpi(lu,iyd_prv,iydrd,f107_prv,f107a_prv,fkp_prv)
        iyd_nxt = inciyd(iyd_cur,1)
        call rdgpi(lu,iyd_nxt,iydrd,f107_nxt,f107a_nxt,fkp_nxt)
c
c If getting fluxes, will need 2nd day previous and 2nd day next:
c
        if (if107.gt.0.or.if107a.gt.0) then
          iyd_p2d = inciyd(iyd_cur,-2)
          call rdgpi(lu,iyd_p2d,iydrd,f107_p2d,f107a_p2d,fkp_p2d)
          iyd_n2d = inciyd(iyd_cur,2)
          call rdgpi(lu,iyd_n2d,iydrd,f107_n2d,f107a_n2d,fkp_n2d)
        endif
      endif ! ifrst == 0
c
c Read new day if necessary:
c
      idayit = iter*ifix(c(4))/86400
      if (idayit*86400.eq.iter*ifix(c(4))) then
        if (if107.gt.0.or.if107a.gt.0) then
          iyd_p2d = iyd_prv
          f107_p2d = f107_prv
          f107a_p2d = f107a_prv
          fkp_p2d = fkp_prv
        endif
        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(:)
        if (if107.gt.0.or.if107a.gt.0) then
          iyd_nxt = iyd_n2d
          f107_nxt = f107_n2d
          f107a_nxt = f107a_n2d
          fkp_nxt = fkp_n2d
          iyd_n2d = inciyd(iyd_n2d,1)
          call rdgpi(lu,iyd_n2d,iydrd,f107_n2d,f107a_n2d,fkp_n2d)
        else
          iyd_nxt = inciyd(iyd_nxt,1)
          call rdgpi(lu,iyd_nxt,iydrd,f107_nxt,f107a_nxt,fkp_nxt)
        endif
      endif ! get new next day
c
c Interpolate fluxes:
c If ut = 12, then use value for current day
c If ut < 12, then use previous two days and current day
c If ut > 12, then use next two days and current day
c     subroutine timeterp(d1,d2,d3,d4,n,ipos,isec_dat,isec_req,fout,
c    +  iprnt,ier)
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)
     +      call timeterp(f107_p2d,f107_prv,f107_cur,f107_nxt,1,
     +        2,isecf107,isec_cur,f107,0,ier)
          if (if107a.gt.0)
     +      call timeterp(f107a_p2d,f107a_prv,f107a_cur,f107a_nxt,1,
     +        2,isecf107,isec_cur,f107a,0,ier)
        else ! interp between current and next days
          if (if107.gt.0)
     +      call timeterp(f107_prv,f107_cur,f107_nxt,f107_n2d,1,
     +        2,isecf107,isec_cur,f107,0,ier)
          if (if107a.gt.0)
     +      call timeterp(f107a_prv,f107a_cur,f107a_nxt,f107a_n2d,1,
     +        2,isecf107,isec_cur,f107a,0,ier)
        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*2
            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)
c
c 6/3/94: Ray wants to set minimum power to 3 gw:
c
        if (ihpower.gt.0) then
          hpower = amax1(3.,-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
      if (ifrst.eq.0) then
        write(6,"(' ')")
        write(6,"('  iter    iyd   hr:min:sec    f107   f107a',
     +    '      Kp ctpoten  hpower')")
        write(6,"(72('-'))")
      endif
      write(6,"(i6,i7,i5,':',i2,':',i2,2x,5f8.3)")
     +  iter,iyd_cur,ih,im,is,f107,f107a,fkp,ctpoten,hpower
      return
c
c i/o error trap:
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'
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c