C     DTCMAKEDR MAKES DTC INDEX FILE AT SELECTED TIMES FROM DST STORMS
C
C     Dst values input in DSTVAL
C     Make sure Dst values at least 10 days before and after
C       requested DTC value computation time specified below
C
C     PROGRAM DEVELOPED BY Bruce R. Bowman - June 2008
C
C     DTCMAKEDR       Driver generates dTc file
C                    *Outputs dTc file
C       DTCDRV        Returns dTc value for input time
C         DTCAP       Computes dTc from Jacchia equation using ap
C           SOLFLUX   Returns 3-hour ap value for input time
C                    *Reads solar flux file
C         DSTDTC      Returns dTc from Dst storm integration of dTc
C           DSTSTM    Returns Dst storm event times
C             DSTBEG  Selects start of Dst storm
C             DSTMAX  Computes time of storm start
C             DSTMIN  Computes time of storm minimum
C             DSTREC  Computes time of recovery slope change
C             DSTEND  Computes time of storm end
C           DSTVAL    Returns Dst value from input time
C                    *Reads Dst file and loads Dst common
C           DSTERR    Checks for Dst error
C           TMINP     Converts input time in yymmdd... to D1950
C           TMOUTD    Converts D1950 to yyddd
C
C     * Following files used:
C
C       DTCFILE.TXT  output unit  8 in DTCMAKEDR contains dTc values
C       DSTFILE.TXT  input  unit 12 in DSTVAL    contains Dst values
C       DSTVAL.GAP   output unit 13 in DSTVAL    contains Dst gap time
C       SOLRESAP.TXT input  unit 20 in SOLFLUX   contains flux & ap
C
C
C     NOTE: Added two small sections in DSTBEG.for subroutine to correct
C           for small storms per Bruce Bowman May 8 email. WKT [May 9 2017]
C
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
C
      CHARACTER*24 FOUT
C
      DIMENSION IDTC(24)
C
C     SET UP OUTPUT FILE
      KOUNIT = 8
      FOUT = 'DTCFILE.TXT'
      OPEN (KOUNIT,FILE=FOUT,ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
C
C     SET START TO END TIME SPAN
C     DST VALUES MUST INCLUDE AT LEAST 10 DAYS PRIOR TO AND AFTER TIMES:
      IYRBEG  =  97
      DYBEG   = 001.0D0
      IYREND  =  07
      DYEND   = 366.999999D0
C
C     CONVERT TIME SPAN TO DAYS SINCE 1950 DAY 0
      IF (IYRBEG.LT.50) IYRBEG = IYRBEG + 100
      IYY = ((IYRBEG-1)/4-12)
      IYY = (IYRBEG-50)*365 + IYY
      D1950B = IYY*1.D0 + DYBEG
      IF (IYREND.LT.50) IYREND = IYREND + 100
      IYY = ((IYREND-1)/4-12)
      IYY = (IYREND-50)*365 + IYY
      D1950E = IYY*1.D0 + DYEND
C
C     COMPUTE DTC VALUE EVERY HOUR
      DTSTEP = 1.0D0/24.D0
C     INITIALIZE VALUES
      DO 5 IND=1,24
    5 IDTC(IND) = 0
C
C     STEP FROM BEG TO END BY DTSTEP
      TSTEP = D1950B - DTSTEP
   10 TSTEP = TSTEP + DTSTEP
      IF (TSTEP.GT.D1950E) GO TO 90
C
C     COMPUTE DTC VALUE AT TSTEP FROM DST ALGORITHMS
      CALL DTCDRV (TSTEP,DTCVAL)
C
C     OUTPUT DTC DATA
C     FORMAT EXAMPLE:
C   YYDDD DTC(DEG)
CTC 97125 010 015 025 035 060 066 078 105 120 135 150 145 141 132 ...
C
      CALL TMOUTD(TSTEP,NYR,DAY)
      NDAY = DAY
      TFDAY = DAY - NDAY
      IHR = TFDAY*24. + 0.0001
      IND = IHR + 1
      IDTC(IND) = DTCVAL + 0.5
      IF (IND.EQ.24) THEN
        WRITE (KOUNIT,105) NYR,NDAY,(IDTC(I),I=1,24)
  105   FORMAT('DTC',I3,I3,24I4)
        DO 15 IND=1,24
   15   IDTC(IND) = 0
      ENDIF
C
      GO TO 10
C
   90 STOP
      END
C
C***********************************************************************
C
      SUBROUTINE DTCDRV (D1950,DTCVAL)
C
C     COMPUTE DELTA TC AT INPUT TIME AND OUTPUT DTCVAL
C
C     IF STORM FOUND (DST <= -75) THEN USE DST TO COMPUTE DTC
C     DURING DST STORM USE DST THROUGH ENTIRE STORM
C     IF NOT STORM THEN COMPUTE DTC FROM AP
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: D1950 - DAYS SINCE 1950 JAN 0
C
C     OUTPUT: DTCVAL - DTC VALUE (DEG)
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     COMPUTE DTC FROM DST IF STORM PRESENT
      CALL DSTDTC (D1950,DTCVAL,NSTORMS)
C
      IF (NSTORMS.EQ.0) THEN
C       Use 3-HR ap value from Jacchia-70
        CALL DTCAP (D1950,DTCVAL)
      ENDIF
C
      RETURN
      END SUBROUTINE DTCDRV
C
C***********************************************************************
C
      SUBROUTINE DTCAP (D1950,DTCVAL)
C
C     COMPUTE DELTA TC FROM ap AT INPUT TIME AND OUTPUT DTCVAL
C
C     THIS IS FOR LOW 3-HOUR AP ACTIVITY ONLY
C     IF AP > 50 THEN SET AP TO 50
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: D1950 - DAYS SINCE 1950 JAN 0
C
C     OUTPUT: DTCVAL - DTC VALUE (DEG)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     USE 6.7 HR LAG FOR ap INFLUENCE
C
      TIME = D1950 - 0.279D0
      CALL SOLFLUX(TIME,X1,X2,X3,AP3HR)
      IF (AP3HR.GT.50.) AP3HR = 50.
C
C     Use 3-HR ap value from Jacchia-70
      EXPAP = DEXP(-0.08D0 * AP3HR)
      DTCVAL = AP3HR + 100.D0 * (1.D0 - EXPAP)
C
      RETURN
      END SUBROUTINE DTCAP
C
C***********************************************************************
C
      SUBROUTINE DSTDTC (TINP,DTCVAL,NSTORMS)
C
C     COMPUTE DELTA TC AT INPUT TIME TINP
C
C     DETERMINE STORM INTERVAL THAT INCLUDES TINP
C       INTEGRATE DTC FROM STORM START (MAX) TO STORM END
C
C     IF NO STORM THEN DTC = 0 FOR TIME AND NSTORMS = 0
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: TINP - DAYS SINCE 1950 JAN 0
C
C     OUTPUT: DTCVAL - DTC VALUE (DEG)
C             NSTORMS - NUMBER OF STORMS
C                       0  = NO STORM OCCURRING AT TINP
C                       1  = SINGLE STORM OCCURRING AT TINP
C                       >1 = MULTIPLE STORM OCCURRING
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     COMMON DSTDATA HAS DST VALUES LOADED FROM PROGRAM DSTVAL
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
C
C     COMMON DTCTIME HAS DTC VALUES LOADED FROM THIS PROGRAM DSTDTC
      COMMON /DTCTIME/ BEGTIM,ENDTIM,DTC(2400),NUMSTMS
C
C     MAXIMUM STORMS FOR MULTIPLE STORM IS 3
      DIMENSION STMTMAX(3),STMTMIN(3),STMTREC(3),STMTEND(3)
      CHARACTER*3 ATYPE
C
C     SET TIME ON 1HR DST GRID POINT PRIOR TO TINP
C     ASSUMES DST EVERY 1 HOUR
      ITIME = TINP
      IHR = (TINP-ITIME)*24. + 0.0001
      TIME = ITIME + IHR/24.
C
C     INITIALIZE DST VALUES - INPUT ALL DST IF NOT ALREADY READ
      IF (IDSTRD.NE.1) CALL DSTVAL (TIME,IDSTVAL)
C
C     DETERMINE TIME FOR STORM FOR INPUT TIME IN STORM INTERVAL
C       STORM START=TMAX, MIN=TMIN, RECOVERY PT=TREC, AND END=TEND
C     IF NO STORM THEN NSTORMS=0 AND FIRST STORM MAX TIME SET TO
C        TIME + 5 DAYS
C     FOR MULTIPLE STORMS WHERE DT INTEGRATION NEEDS TO START AT FIRST
C        STORM AND CONTINUE THROUGH TO LAST STORM, USE NSTORMS
C        THREE STORMS MAX FOR THIS SCENARIO
C
      CALL DSTSTM (TIME,STMTMAX,STMTMIN,STMTREC,STMTEND,NSTORMS)
C
      IF (NSTORMS.EQ.0) THEN
        INDDTC = 2
        DTC(1) = 0.D0
        DTC(2) = DTC(1)
        BEGTIM = TIME
        ENDTIM = STMTMAX(1)
        IF (ENDTIM.GT.BEGTIM+5.D0) ENDTIM = BEGTIM + 5.D0
        DTCVAL  = 0.D0
        TSTEP = TIME
   10   TSTEP = TSTEP + DTDST
        INDDTC = INDDTC + 1
        DTC(INDDTC) = 0.D0
        IF (TSTEP.LT.ENDTIM) GO TO 10
        RETURN
      ENDIF
C
C     INTEGRATE DTC FROM BEGINNING TO END OF STORM(S)
C
C     INITIALIZE DTC AT START TIME
C     ADD TIME DELAY OF 1-2 HRS BASED ON STORM MAGNITUDE
      TMIN = STMTMIN(1)
      INDX = 1 + (TMIN - DS1950 + 0.0001D0)*24
      IDSTMIN = IDST(INDX)
      DELAY = 0.D0
      IF (IDSTMIN.LE.-250.AND.IDSTMIN.GT.-350) DELAY = 1.D0/24.D0
      IF (IDSTMIN.GT.-250) DELAY = 2.D0/24.D0
      TMAX = STMTMAX(1) + DELAY
      IF (TINP.LT.TMAX) THEN
        DTCVAL = 0.D0
        NSTORMS = 0
        RETURN
      ENDIF
C     START STORM AT CURRENT AP DENSITY
      CALL DTCAP(TMAX,DTCVAL)
      INDDTC = 2
      DTC(1) = DTCVAL
      DTC(2) = DTC(1)
C
      DO 50 I=1,NSTORMS
C       SET VALUES FOR STORM
        TMAX = STMTMAX(I)
        TMIN = STMTMIN(I)
        TREC = STMTREC(I)
        TEND = STMTEND(I)
C
        INDX = 1 + (TMIN - DS1950 + 0.0001D0)*24
        IDSTMIN = IDST(INDX)
        DSTMINV = IDSTMIN
C
C       INTEGRATE MAIN PHASE ( START=MAX TO MIN )
        TSTEP = TMAX
   15   TSTEP = TSTEP + DTDST
        INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
C       SET MAIN PHASE DERIVATIVE BASED ON STORM MAGNITUDE
        SLPMAIN =  -1.5050D-05*(DSTMINV)**2 - 1.0604D-02*DSTMINV - 3.20
        IF (IDSTMIN.LT.-450) SLPMAIN = -1.40
        DERIV = SLPMAIN
C       CHECK FOR MAX VALUE TOO LARGE FROM SUDDEN COMMENCEMENT WAVE
        IDST1 = IDST(INDX)
        IDST0 = IDST(INDX-1)
        IF (IDST1.GT.0) IDST1 = 0
        IF (IDST0.GT.0) IDST0 = 0
C
        DELDST = IDST1 - IDST0
        IF (DELDST.GE.0.D0) THEN
          SFAC = 0.30
          DTCVAL = DTCVAL - SFAC*DERIV*DELDST
        ELSE
          DTCVAL = 0.846D0*DTCVAL
     *             + DERIV*(IDST1 - 0.870D0*IDST0)
        ENDIF
        INDDTC = INDDTC + 1
        DTC(INDDTC) = DTCVAL
C
        IF (TSTEP.LT.TMIN) GO TO 15
C
C       INTEGRATE RECOVERY PHASE ( MIN TO RECOVERY INFLECTION PT )
   20   TSTEP = TSTEP + DTDST
        INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
        DERIV =  0.13D0
        DTCVAL = 1.000D0*DTCVAL
     *           + DERIV*(IDST(INDX) - 0.000D0*IDST(INDX-1))
        INDDTC = INDDTC + 1
        DTC(INDDTC) = DTCVAL
C
        IF (TSTEP.LT.TREC) GO TO 20
C
C       INTEGRATE POST RECOVERY PHASE ( RECOVERY INFLECTION PT TO END )
   30   TSTEP = TSTEP + DTDST
        INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
        DELDST = IDST(INDX) - IDST(INDX-1)
        SLPEND  =  -2.5D0
        DERIV = SLPEND
        IF (DELDST.LT.0.D0) DERIV = SLPMAIN
        DTCVAL = DTCVAL + DERIV*DELDST
        IF (DTCVAL.LT.0.D0) DTCVAL = 0.D0
        INDDTC = INDDTC + 1
        DTC(INDDTC) = DTCVAL
C
        IF (TSTEP.LE.TEND-0.0000001D0) GO TO 30
C
   50 CONTINUE
C
C     ADD TIME DELAY OF 1-2 HRS BASED ON STORM MAGNITUDE
      DELAY = 0.D0
      IF (IDSTMIN.LE.-250.AND.IDSTMIN.GT.-350) DELAY = 1.D0/24.D0
      IF (IDSTMIN.GT.-250) DELAY = 2.D0/24.D0
      BEGTIM = STMTMAX(1) + DELAY
      ENDTIM = TEND + DELAY
      IF (ENDTIM.GT.BEGTIM+5.D0) ENDTIM = BEGTIM + 5.D0
C
C     LOAD DTC AND DST VALUE FOR TIME OF INPUT
C
      INDDST = 1 + (TIME - DS1950 + 0.0000001D0)*24.
      IF (INDDST.LT.2.OR.INDDST.GT.IDSTNM) CALL DSTERR
C     INDDTC=2 STARTS AT BEGTIM, INDDTC=1 SAME VALUE AS AT BEGTIM
      INDDTC = 2. + (TIME - BEGTIM + 0.0000001D0)*24.
      IF (INDDTC.LT.2) INDDTC = 2
      DTCVAL1  = DTC(INDDTC)
      DTCVAL2  = DTC(INDDTC+1)
C
C     INTERPOLATE VALUES
C
      FX = 2. + (TINP - BEGTIM + 0.0000001D0)*24.
      FAC = (FX - INDDTC)/1.D0
      DTCVAL  = DTCVAL1  + FAC*(DTCVAL2-DTCVAL1)
C
      RETURN
      END SUBROUTINE DSTDTC
C
C***********************************************************************
C
      SUBROUTINE DSTSTM (TIME,STMTMAX,STMTMIN,STMTREC,STMTEND,NSTORMS)
C
C     DETERMINE TIMES FOR STORM FOR INPUT TIME IN STORM INTERVAL
C       STORM START=TMAX, MIN=TMIN, RECOVERY PT=TREC, AND END=TEND
C     STORM DEFINED AS DST <=-75 AND (DSTmax - DSTmin) >= 50
C
C     STORM MAIN PHASE IS BETWEEN MAX AND MIN
C     STORM RECOVERY PHASE IS BETWEEN MIN AND END
C
C     STORM BEGINNING IS MAX OF STORM DEFINNING START
C     STORM MINIMUM   IS MIN DST OF STORM
C     STORM RECOVERY  IS REFLECTION POINT WHERE SLOPE CHANGES
C     STORM END       IS BASED ON EQUATION OF DEL TIME VS DSTmin
C
C     SUB DSTBEG DETERMINES BEGINNING (MAX) POINT OF STORM BY
C       ITERATING ON DETERMINING MIN, THEN MAX OF STORM, CHECK
C       IF TIME POINT IN NO STORM AREA WHERE DST PTS AROUND TIME >-40
C     SUB DSTMAX STEPS FROM MIN TIME BACK TO SPECIFIED BEG TIME
C       CHECKING FOR MAX INBETWEEN
C     SUB DSTMIN STEPS FROM MAX TIME FORWARD TO SPECIFIED END TIME
C       CHECKING FOR MIN INBETWEEN
C     SUB DSTEND COMPUTES END TIME FROM MIN BASED ON MAGNITUDE OF MIN
C       CHECKING FOR NEXT STORM POSSIBLY OCCURING BEFORE COMPUTED END
C     SUB DSTREC STEPS FORWARD FROM MIN TO END CHECKING FOR SLOPE
C       CHANGING ENOUGH TO SIGNIFY RECOVERY INFLECTION POINT
C
C     IF NO STORM THEN NSTORMS=0 AND FIRST STORM MAX TIME SET TO
C        TIME + 5 DAYS
C     FOR MULTIPLE STORMS WHERE DT INTEGRATION NEEDS TO START AT FIRST
C        STORM AND CONTINUE THROUGH TO LAST STORM, USE NSTORMS > 0
C        3 STORMS MAX FOR THIS SCENARIO
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT:  TIME - DAYS SINCE 1950 JAN 0
C
C     OUTPUT: STMTMAX - ARRAY OF STORM MAX (START) TIMES
C             STMTMIN - ARRAY OF STORM MIN TIMES
C             STMTREC - ARRAY OF STORM RECOVERY SLOPE CHANGE TIMES
C             STMTEND - ARRAY OF STORM END TIMES
C             NSTORMS - NUMBER OF STORMS (MAX=3)
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
      DIMENSION STMTMAX(3),STMTMIN(3),STMTREC(3),STMTEND(3)
C
C     INITIALIZE
      DO 5 I=1,3
      STMTMAX(1) = 0
      STMTMIN(1) = 0
      STMTREC(1) = 0
   5  STMTEND(1) = 0
C
C     DETERMINE STARTING POINT OF STORM
      CALL DSTBEG (TIME,TSTART,NSTORM)
C
      IF (TSTART.GT.TIME.OR.NSTORM.EQ.0) THEN
C       NO STORM FOUND PRIOR TO CURRENT INPUT TIME
        STMTMAX(1) = TSTART
        NSTORMS = 0
        RETURN
      ENDIF
C
C     DETERMINE MAIN PHASE MIN, THEN MAX TIME
      TEND = TSTART + 2.
      CALL DSTMIN (TSTART,TEND,TMIN,IMIN)
      CALL DSTMAX (TMIN,TSTART,TMAX,IMAX)
C
      IF (TMAX.GT.TSTART) THEN
C       NEW MAX FOUND AFTER TSTART TIME
C       NO STORM FOUND PRIOR TO CURRENT INPUT TIME
        STMTMAX(1) = TSTART
        NSTORMS = 0
        RETURN
      ENDIF
C
C     DETERMINE RECOVERY PT AND STORM END
C     SET STORM END TIME BASED ON DST MIN VALUE
C
      CALL DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
      CALL DSTREC (TMIN,TEND,TREC)
      IF (TIME.GT.TEND) NEWSTM = 1
C
C     SET STORM TIMES FOR 1ST STORM
      STMTMAX(1) = TMAX
      STMTMIN(1) = TMIN
      STMTREC(1) = TREC
      STMTEND(1) = TEND
      NSTORMS = 1
C
      IF (NEWSTM.EQ.1) THEN
C       2ND STORM FOUND
C       DETERMINE MAIN PHASE MIN USING MAX TIME = 1ST STORM END TIME
        TSTART = TEND
        TEND = TSTART + 10.
        CALL DSTMIN (TSTART,TEND,TMIN,IMIN)
        CALL DSTMAX (TMIN,TSTART,TMAX,IMAX)
C
C       DETERMINE RECOVERY PT AND STORM END
C       SET STORM END TIME BASED ON DST MIN VALUE
C
        CALL DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
        CALL DSTREC (TMIN,TEND,TREC)
C
C       SET STORM TIMES FOR 2ND STORM
        STMTMAX(2) = TMAX
        STMTMIN(2) = TMIN
        STMTREC(2) = TREC
        STMTEND(2) = TEND
        NSTORMS = 2
C
        IF (NEWSTM.EQ.1) THEN
C         3RD STORM FOUND
C         DETERMINE MAIN PHASE MIN USING MAX TIME = 2ND STORM END TIME
          TSTART = TEND
          TEND = TSTART + 10.
          CALL DSTMIN (TSTART,TEND,TMIN,IMIN)
          CALL DSTMAX (TMIN,TSTART,TMAX,IMAX)
C
C         DETERMINE RECOVERY PT AND STORM END
C         SET STORM END TIME BASED ON DST MIN VALUE
C
          CALL DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
          CALL DSTREC (TMIN,TEND,TREC)
C
C         SET STORM TIMES FOR 3RD STORM
          STMTMAX(3) = TMAX
          STMTMIN(3) = TMIN
          STMTREC(3) = TREC
          STMTEND(3) = TEND
          NSTORMS = 3
        ENDIF
      ENDIF
C
      RETURN
      END SUBROUTINE DSTSTM
C
C***********************************************************************
C
      SUBROUTINE DSTBEG (TIME,TSTART,NSTORM)
C
C     FIND STARTING PT OF STORM BY HAVING 6 CONSECUTIVE PTS > -40
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C     NSTORM = 0 FOR NO STORM, 1 FOR STORM AT CURRENT TIME
C
C     INPUT: TIME - DAYS SINCE 1950 JAN 0
C
C     OUTPUT: TSTART - STORM START TIME
C             NSTORM - 0=NO STORM, 1=STORM
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
C
C     SET MIN DEL MAG TO BE STORM
      IMAG = 50
      NSTORM = 1
C
C     DETERMINE DST AT CURRENT TIME
      INDX = 1 + (TIME - DS1950 + 0.0001D0)*24
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
      IDSTVAL = IDST(INDX)
C
      IF (IDSTVAL.GE.-40) THEN
C       CHECK FOR 6 PRIOR PTS IN NO TO SMALL STORM REGION
        IPTS = 1
        TSTEP = TIME
   10   TSTEP = TSTEP - DTDST
        INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
        IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
        IDSTVAL = IDST(INDX)
        IF (IDSTVAL.GE.-40) IPTS = IPTS + 1
C       RESET NUM PTS IN NO STORM REGION IF STEPPED OUTSIDE REGION
        IF (IDSTVAL.LT.-50) IPTS = 0
C       STORM FOUND JUST BEFORE INPUT TIME
        IF (IDSTVAL.LE.-75) GO TO 20
        IF (IPTS.LT.6) GO TO 10
        TSAVE = TSTEP
C
C       ALL PTS IN NO TO SMALL STORM REGION
C       STEP FORWARD UNTIL STORM START FOUND
C       FIND FIRST WHERE DST < -75 OR 5 DAYS AFTER TIME
        TSTEP = TIME
   15   TSTEP = TSTEP + DTDST
        INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
C
        IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
        IDSTVAL = IDST(INDX)
        IF (IDSTVAL.GT.-75.AND.TSTEP.LT.TIME+7.D0) GO TO 15
C
        IF (TSTEP.GE.TIME+7.D0) THEN
C         NO STORM FOUND
          NSTORM = 0
C         USE 5 DAYS INSTEAD OF 7 SINCE AT END OF 7 DAYS COULD BE
C           ALREADY BEGINNING ON NEXT STORM
          TSTART = TIME+5.D0
          GO TO 80
        ENDIF
C
C       STORM FOUND (DST < -75) , DETERMINE MAX TIME FOR STORM START
C       CHECK FOR VALID STORM MAGNITUDE
        TBEG = TSAVE - 1.
        CALL DSTMAX (TSTEP,TBEG,TMAX,IMAX)
        TEND = TMAX + 4.
        CALL DSTMIN (TMAX,TEND,TMIN,IMIN)
C
C       DETERMINE MAX AGAIN FROM REAL MIN FOUND
        TBEG = TMIN - 1.5
        CALL DSTMAX (TMIN,TBEG,TMAX,IMAX)
        TMAXSAV = TMAX
C
C       CHECK IF THIS MAX IS PRIOR TO END OF PREVIOUS STORM
C       STEP BACKWARD UNTIL PREVIOUS MAX, MIN, AND THEN END COMPUTED
        TBEG = TMAX - 4.
        CALL DSTMAX (TMAX,TBEG,TMAXP,IMAXP)
        TENDP = TMAX
        CALL DSTMIN (TMAXP,TENDP,TMNP,IMINP)
C
        IDELDST = IMAXP - IMINP
        IF (IMINP.GT.-75.OR.IDELDST.LE.50) THEN
C         MIN TOO SMALL FOR STORM
          GO TO 24
        ENDIF
C
        CALL DSTEND (TMNP,IMINP,IMAXP,TENDP,NEWSTM)
C
        IF (TENDP.GT.TMAX) THEN
C         ADJUST CURRENT MAX POINT TO END OF PREVIOUS STORM
          TMAX = TENDP
          TMAXSAV = TMAXP
          INDX = 1 + (TMAX - DS1950 + 0.0001D0)*24
          IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
          IMAX = IDST(INDX)
C         RECOMPUTE MIN AND END POINTS
          TEND = TMAX + 3.
          CALL DSTMIN (TMAX,TEND,TMIN,IMIN)
          CALL DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
        ENDIF
C
        IDELDST = IMAX - IMIN
        IF (IDELDST.LT.IMAG.OR.IMIN.GE.-50) THEN
C         NO NEXT STORM FOUND
          NSTORM = 0
          TSTART = TMIN
          GO TO 80
        ENDIF
        NSTORM = 1
        TSTART = TMAXSAV
        GO TO 80
      ENDIF
C
C     SMALL TO LARGE STORM IN PROGRESS
C     STEP BACK UNTIL MAX POINT = STORM START FOUND
C     CHECK IF STORM BIG ENOUGH ( DDST > 75 ) TO BE STORM
C
   20 TBEG = TIME
   22 TBEG = TBEG - 5.D0
C
      CALL DSTMAX (TIME,TBEG,TMAX,IMAX)
      TEND = TMAX + 5.D0
      CALL DSTMIN (TMAX,TEND,TMIN,IMIN)
      IF (IMIN.GT.-75) THEN
C       MIN TOO SMALL FOR STORM
        NSTORM = 0
        TSTART = TMIN
        GO TO 80
      ENDIF
C
CCCC
CCCC  CHANGE ADDED MAY 2017
      TSRT = TMIN - 3.0D0
CCCC
CCCC  TSRT = TMIN - 1.5D0
      CALL DSTMAX (TMIN,TSRT,TMAX,IMAX)
      TMAXSAV = TMAX
C     CHECK FOR PREVIOUS STORM NOT RECOVERED YET, LOOK FOR THAT MAX
CCCC
CCCC  CHANGE ADDED MAY 2017
      IF ((IMAX.LT.-50).AND.(TIME-TBEG.LT.20.0)) GO TO 22
CCCC
CCCC  IF (IMAX.LT.-50) GO TO 22
      IDELDST = IMAX - IMIN
      IF (IDELDST.LT.50) THEN
C       MAX-MIN DIFF TOO SMALL FOR STORM
        NSTORM = 0
        TSTART = TMIN
        GO TO 80
      ENDIF
C
      CALL DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
C
C     CHECK IF THIS MAX IS PRIOR TO END OF PREVIOUS STORM
C     STEP BACKWARD UNTIL PREVIOUS MAX, MIN, AND THEN END COMPUTED
      TBEG = TMAX - 4.
      CALL DSTMAX (TMAX,TBEG,TMAXP,IMAXP)
      TENDP = TMAX
      CALL DSTMIN (TMAXP,TENDP,TMNP,IMINP)
C
      IDELDST = IMAXP - IMINP
      IF (IMINP.GT.-75.OR.IDELDST.LE.50) THEN
C       MIN TOO SMALL FOR STORM
        GO TO 24
      ENDIF
C
      CALL DSTEND (TMNP,IMINP,IMAXP,TENDP,NEWSTM)
C
      IF (TENDP.GE.TMAX) THEN
C       ADJUST CURRENT MAX POINT TO END OF PREVIOUS STORM
        TMAX = TENDP
        TMAXSAV = TMAXP
        INDX = 1 + (TMAX - DS1950 + 0.0001D0)*24
        IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
        IMAX = IDST(INDX)
C       RECOMPUTE MIN AND END POINTS
        TEND = TMAX + 3.
        CALL DSTMIN (TMAX,TEND,TMIN,IMIN)
        CALL DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
      ENDIF
C
C     CHECK TO DETERMINE IF STORM OVER BUT DST NOT YET RECOVERED
   24 IF (TEND.LT.TIME) THEN
C       STORM OVER, STEP FORWARD UNTIL NEXT STORM START FOUND
C       FIND WHERE FIRST DST < -75 OR  < -100 LARGE STORM
C         OR 5 DAYS AFTER TIME
C       FOR LARGE STORM DST DOES NOT RECOVER QUICKLY
        IDSTLIM = -75
C       FOR LARGE STORM
        IF (IDELDST.GT.200) IDSTLIM = -100
        TSTEP = TIME
   25   TSTEP = TSTEP + DTDST
        INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
        IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
        IDSTVAL = IDST(INDX)
        IF (IDSTVAL.GE.IDSTLIM.AND.TSTEP.LT.TEND+7.D0) GO TO 25
C
        IF (TSTEP.GE.TEND+7.D0) THEN
C         NO STORM FOUND
          NSTORM = 0
          TSTART = TEND+5.D0
          GO TO 80
        ENDIF
C
C       NEXT STORM FOUND ( DST < -75)  OR DST < -100 LARGE STORM )
C       DETERMINE MAX AND MIN TIMES FOR STORM
        CALL DSTMAX (TSTEP,TEND,TMAX,IMAX)
        TENDX = TMAX + 4.
        CALL DSTMIN (TMAX,TENDX,TMIN,IMIN)
C
        TBEG = TMIN - 1.5
        IF (TBEG.LT.TEND) TBEG = TEND
        CALL DSTMAX (TMIN,TBEG,TMAX,IMAX)
C
        IF (TMAX.EQ.TEND) THEN
C         START OF THIS NEXT STORM A CONTINUATION OF PREVIOUS STORM
          TSTART = TMAXSAV
          NSTORM = 1
          GO TO 80
        ENDIF
C
        IDELDST = IMAX - IMIN
        IF (IDELDST.LT.IMAG) THEN
C         NO NEXT STORM FOUND
          NSTORM = 0
          TSTART = TMIN
          GO TO 80
        ENDIF
        NSTORM = 1
        TSTART = TMAX
        GO TO 80
      ENDIF
C
C     INPUT TIME IN MIDDLE OF STORM SINCE TEND > TIME
      IDELDST = IMAX - IMIN
      IF (IDELDST.LT.IMAG) THEN
C       NO VALID STORM
        NSTORM = 0
        TSTART = TEND
      ELSE
        NSTORM = 1
        TSTART = TMAXSAV
      ENDIF
C
   80 IF (TIME.LT.TSTART) NSTORM = 0
C
      RETURN
      END SUBROUTINE DSTBEG
C
C***********************************************************************
C
      SUBROUTINE DSTMAX (TINIT,TBEG,TMAX,IMAX)
C
C     DSTMAX STEPS BACKWARD FROM TINIT UNTIL EITHER 6 DST PTS > -40 OR
C       TBEG REACHED.  MAX TIME = TMAX SAVED AND OUTPUT
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: TINIT - TIME OF INTEREST IN DAYS SINCE 1950 JAN 0
C            TBEG  - TIME FOR EARLIEST POSSIBLE STORM START
C
C     OUTPUT: TMAX - START TIME OF STORM
C             IMAX - DST VALUE AT TMAX
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
C
      IMAX = -999
      TMAX = TINIT
      TMIN = TINIT - 1.
      IMIN = 999
      IPTS = 0
      TSTEP = TINIT + DTDST
   10 TSTEP = TSTEP - DTDST
      INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
      IDSTVAL = IDST(INDX)
      IF (TMAX.LT.TMIN.AND.IMIN.LT.-75.AND.IMAX.GT.-50
     *    .AND.IDSTVAL.LT.-75) THEN
C       MIN, THEN MAX REACHED, NOW STEPPING TOWARD PREVIOUS STORM MIN
C       STOP SINCE THIS STORMS MAX REACHED
        GO TO 90
      ENDIF
      IF (IDSTVAL.GT.IMAX) THEN
        IMAX = IDSTVAL
        TMAX = TSTEP
      ENDIF
      IF (IDSTVAL.LT.IMIN) THEN
        IMIN = IDSTVAL
        TMIN = TSTEP
        IPTS = 0
      ENDIF
C
C     RESET PTS IN NO STORM OR SMALL STORM IF IT DROPS TOO LOW
      IF (IDSTVAL.LT.-60) IPTS = 0
C     CHECK TO SEE IF LAST MAX STILL INCREASING
      IF (IDSTVAL.GE.-40.AND.IDSTVAL.LE.IMAX) IPTS = IPTS + 1
      IF (IPTS.LT.6.AND.TSTEP.GT.TBEG) GO TO 10
C
   90 IF (TMAX.GT.TMIN) THEN
C       FAILED TO FIND MAX BEFORE MIN, RESTART STEPPING AT TMIN
        IPTS = 0
        IMAX = -999
        TSTEP = TMIN - DTDST
        GO TO 10
      ENDIF
C
      RETURN
      END SUBROUTINE DSTMAX
C
C***********************************************************************
C
      SUBROUTINE DSTMIN (TINIT,TEND,TMIN,IMIN)
C
C     DSTMIN STEPS FORWARD FROM TINIT UNTIL EITHER 2 DST PTS > -40
C       OR MIN REACHED FOLLOWED BY MAX ( MAX-MIN > 75 AND MAX > -75 )
C         AND DOES NOT REACH 2 PTS > -40 BEFORE DESCEND TO ANOTHER MIN
C         OR TEND REACHED
C       MIN TIME = TMIN SAVED AND OUTPUT
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: TINIT - TIME OF INTEREST IN DAYS SINCE 1950 JAN 0
C            TEND  - TIME FOR LATEST POSSIBLE STORM MIN TO OCCUR
C
C     OUTPUT: TMIN - TIME OF STORM MIN
C             IMIN - DST VALUE AT TMIN
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
C
      IMIN = +999
      TMIN = TINIT
      IPTS = 0
      TSTEP = TINIT
   10 TSTEP = TSTEP + DTDST
      INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
      IDSTVAL = IDST(INDX)
      IF (IDSTVAL.LT.IMIN) THEN
        IMIN = IDSTVAL
        TMIN = TSTEP
        IMAX = IMIN
      ENDIF
C
      IF (IDSTVAL.GT.IMAX) THEN
        IMAX = IDSTVAL
        TMAX = TSTEP
      ENDIF
C
C     CHECK FOR MIN AND MAX DIFF MORE THAN 125 TO SHOW MIN REACHED
      IF (IMAX.GT.IMIN+125) GO TO 90
C
C     CHECK FOR IMAX > IMIN + 75 AND IMAX > 75 TO SIGNIFY DONE
      IF (IMAX-IMIN.GT.75.AND.IMAX.GT.-75) THEN
C       PREVIOUS MIN IS STORM EVEN THOUGH MAX NOT REACHED RECOVERY
        GO TO 90
      ENDIF
C
C     CHECK FOR END OF STORM NOT BEGINNING SINCE 6PTS AROUND MAX
C       COULD STILL BE > -50
      IF (IDSTVAL.LT.-40) IPTS = 0
      IF (IMIN.LT.-75) IPTS = IPTS + 1
      IF (IPTS.LT.2.AND.TSTEP+0.00001D0.LT.TEND) GO TO 10
C
   90 RETURN
      END SUBROUTINE DSTMIN
C
C***********************************************************************
C
      SUBROUTINE DSTREC (TMIN,TEND,TREC)
C
C     DSTREC STEPS FORWARD FROM TMIN UNTIL
C       6 DST PTS > -40 OR TEND REACHED
C     IF SLOPE < SLOPE LIMIT THEN SET RECOVERY INFLECTION PT = TREC
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: TMIN - TIME OF STORM MIN
C            TEND - TIME OF STORM END
C
C     OUTPUT: TREC - TIME OF STORM RECOVERY SLOPE CHANGE PT
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
C
C     SLOPE LIMIT IS 100 DST/DAY
      SLPLIM = 100.
      TREC = 0.D0
      IPTS = 0
      IREC = 0
      TSTEP = TMIN  + DTDST
   10 TSTEP = TSTEP + DTDST
      INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
      IDSTVAL = IDST(INDX)
C
C     CHECK FOR SLOPE LESS THAN SLOPE LIMIT FOR REC PT
C     COMPUTE SLOPE BASED ON PTS EITHER SIDE OF CURRENT PT
      SLOPE = (IDST(INDX+1)-IDST(INDX-1))/(2.*DTDST)
      IF (SLOPE.LT.SLPLIM) THEN
        IREC = IREC + 1
        IF (IREC.GE.3) THEN
          TREC = TSTEP - DTDST
          GO TO 90
        ENDIF
      ENDIF
C
      IF (IDSTVAL.GE.-40) IPTS = IPTS + 1
      IF (IPTS.LT.6.AND.TSTEP+0.00001D0.LT.TEND) GO TO 10
      TREC = TEND
C
   90 RETURN
      END SUBROUTINE DSTREC
C
C***********************************************************************
C
      SUBROUTINE DSTEND (TMIN,IMIN,IMAX,TEND,NEWSTM)
C
C     DETERMINE END TIME OF STORM
C     INITIALLY SET STORM END TIME BASED ON DST MIN VALUE
C       OR 6 PTS < -40
C     THEN STEP FROM MIN TO COMPUTED END
C     IF /DST/ DROPS > 75 FROM MAX DURING STEP REACHED THEN NEW STORM
C       SO NEWSTM = 1
C
C     ALL TIMES IN DAYS SINCE 1950 JAN 0
C
C     INPUT: TMIN - TIME OF STORM MIN
C            IMIN - DST MIN VALUE AT TMIN
C            IMAX - DST MAX VALUE AT STORM START
C
C     OUTPUT: TEND - TIME OF STORM END
C             NEWSTM - 0=NO FOLLOWING STORM
C                      1=NEW STORM STARTS BEFORE REAL STORM END
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
C
C     TMINX = TMIN
C     CHECK FOR NEW MIN TIME AT END OF FLAT MIN OCCURRING
      TSTEP = TMIN
   10 TSTEP = TSTEP + DTDST
      INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24.
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
      IDSTVAL = IDST(INDX)
      IF (IDSTVAL.GT.IMIN+15) THEN
C       SET NEW MIN AT PREVIOUS STEP AS END OF FLAT MIN PTS
        TMINX = TSTEP - DTDST
        GO TO 15
      ENDIF
      GO TO 10
C
   15 DELDST = IMAX - IMIN
      TEND = TMINX + 0.0075D0*DELDST
C     SET TIME ON 1HR DST GRID POINT PRIOR TO TEND
      ITIME = TEND
      IHR = (TEND-ITIME)*24. + 0.0000001D0
      TEND = ITIME + IHR/24.
C
C     CHECK FOR STORM ENDING BEFORE THIS OR NEW STORM BEFORE THIS
      NEWSTM = 0
      IPTS = 0
      IMAXST = IMIN
      TSTOP = TEND
      TSTEP = TMINX
   20 TSTEP = TSTEP + DTDST
      INDX = 1 + (TSTEP - DS1950 + 0.0001D0)*24
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
      IDSTVAL = IDST(INDX)
      IF (IDSTVAL.GT.IMAXST) THEN
        IMAXST = IDSTVAL
        TMAXST = TSTEP
      ENDIF
      IF (IDSTVAL-IMAXST.LT.-75) THEN
C       NEW STORM FOUND
        TEND = TMAXST
        NEWSTM = 1
        GO TO 90
      ENDIF
      IF (IDSTVAL.GT.-75) IPTS = IPTS + 1
      IF (IPTS.EQ.6.OR.TSTEP+0.00001D0.GE.TSTOP) THEN
        TEND = TSTEP
        GO TO 90
      ENDIF
      GO TO 20
C
   90 RETURN
      END SUBROUTINE DSTEND
C
C***********************************************************************
C
      SUBROUTINE DSTVAL (T1950,IDSTVAL)
C
C     INPUT T1950 AND READ FILE FOR VALUE FOR DST
C     DST ON 1 HR INCREMENT, SO OUTPUT DST CLOSEST POINT PRIOR TO T1950
C
C     PROGRAM INITIALLY LOADS COMMON OF ALL DST VALUES FROM FILE
C       SUBSEQUENT CALLS OBTAINS VALUES FROM COMMON
C
C     PROGRAM CHECKS FOR MISSING DST VALUES IN FIRST CALL
C       IF MISSING VALUES THEN PROGRAM WRITES TO GAP FILE AND ABORTS
C
C     PROGRAM CHECKS FOR INPUT TIME OUTSIDE OF DST FILE RANGE
C       IF OUTSIDE FILE RANGE PROGRAM ABORTS
C
C
C     INPUT: T1950 - TIME OF INTEREST IN DAYS SINCE 1950 JAN 0
C
C     OUTPUT: IDSTVAL - DST VALUE AT T1950
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      CHARACTER*50 FDST,FGAP
C
      COMMON /DSTDATA/ DS1950,DTDST,IDST(500000),IDSTRD,IDSTNM
      DIMENSION IDSTH(24)
C
C     SET UP I/O FILES
C
      KIUNIT = 12
      KGUNIT = 13
C
      FDST = 'DSTFILE.TXT'
      FGAP = 'DSTVAL.GAP'
C
      OPEN (KIUNIT,FILE=FDST,ACCESS='SEQUENTIAL',STATUS='OLD')
C
      IF (IDSTRD.NE.1) THEN
C
C       INITIALIZE DST FILE AND DATA
        OPEN (KIUNIT,FILE=FDST,ACCESS='SEQUENTIAL',STATUS='OLD')
C
C       INPUT EXAMPLE:
CST9701*01  X219 000-005 000 003 006 009 007 001-003 000 005 002 001 001-002-003-004-006-006-005-003-006-003-009-006-001
C
        READ (KIUNIT,101,END=30) IY,IM,IDAY,KY,(IDSTH(I),I=1,24),IDSTMON
        IHR  = 0
        IMIN = 0
        SEC  = 0.D0
        CALL TMINP(IY,IM,IDAY,IHR,IMIN,SEC,DS1950)
C
        BACKSPACE (KIUNIT)
        IDSTRD = 1
        IDSTNM = 0
        TDATA = 0.D0
        DTDIF = 1.01D0/24.D0
        DTDST = 1.00D0/24.D0
        KNUM = 0
C
   10   READ (KIUNIT,101,END=30) IY,IM,IDAY,KY,(IDSTH(I),I=1,24),IDSTMON
  101   FORMAT(3X,I2,I2,1X,I2,2X,1X,1X,I2,4X,24I4,I4)
C
        DO 20 I=1,24
          KNUM = KNUM + 1
          IHR  = I - 1
          IMIN = 0
          SEC  = 0.D0
          CALL TMINP(IY,IM,IDAY,IHR,IMIN,SEC,D1950)
          DT = D1950 - TDATA
C
C         CHECK FOR DATA GAPS OR MISSING DATA
          IF (KNUM.GT.1.AND.(IDSTH(I).GT.900.OR.DT.GT.DTDIF)) THEN
            OPEN (KGUNIT,FILE=FGAP,ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
            CALL TMOUTD(D1950,NYR,DAY)
            WRITE (KGUNIT,104) D1950,NYR,DAY,IDSTH(I),DS1950
  104       FORMAT(F15.5,I5,F12.5,I10,F15.5)
            STOP
          ENDIF
C         LOAD DST DATA INTO COMMON
          IDST(KNUM)  = IDSTH(I)
          TDATA = D1950
   20   CONTINUE
        GO TO 10
C
   30   IDSTNM = KNUM
      ENDIF
C
      INDX = 1 + (T1950 - DS1950 + 0.0001D0)*24
C
      IF (INDX.LT.1.OR.INDX.GT.IDSTNM) CALL DSTERR
C
      IDSTVAL = IDST(INDX)
C
      RETURN
      END SUBROUTINE DSTVAL
C
C***********************************************************************
C
      SUBROUTINE DSTERR
C
C     ERROR DUE TO REQUESTED DST OUTSIDE DST DATA SPAN
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      WRITE(*,101)
  101 FORMAT('  ERROR IN DSTDRV - TIME OUTSIDE OF DST VALUES ')
      STOP
      END SUBROUTINE DSTERR
C
C***********************************************************************
C
      SUBROUTINE SOLFLUX (TIME,F10,F10B,AP,AP3HR)
C
C     READ SOLAR FLUX FROM FILE SOLRES.DAT FOR INPUT TIME
C     PROGRAM ABORTS IF MISSING DAYS IN FILE
C
C     INPUT: TIME - DAYS SINCE 1950 JAN 0
C
C     OUTPUT: F10
C             F10B
C             AP - DAILY VALUE
C             AP3HR
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /FLUXDATA/ FTIME,FT0,IF10,IF10B,IAP,IAPAVE,IREAD
      DIMENSION IAP(8)
      CHARACTER*50 FLUXUNIT
C
C     SET UP INPUT SOLAR DATA FILE
      KFUNIT = 20
      FLUXUNIT = 'SOLRESAP.TXT'
C
    5 IF (IREAD.NE.1) THEN
C       INITIALIZE FLUX FILE AND DATA
        OPEN (KFUNIT,FILE=FLUXUNIT,ACCESS='SEQUENTIAL',STATUS='OLD')
        READ (KFUNIT,101,ERR=80) IDAY,IF10,IF10B,(IAP(I),I=1,8),IY
  101   FORMAT(I3,2I4,8I4,26X,I2)
        IF (IY.LT.50) IY = IY + 100
        IYY = (IY-1)/4 - 12
        D1950 = (IY-50)*365 + IYY + IDAY
        FTIME = D1950
        FT0 = D1950
        IREAD = 1
      ENDIF
C
 10   IF ((TIME-FTIME).LT.1.0D0.AND.TIME.GE.FTIME) THEN
C       USE FLUX VALUES FROM COMMON
        F10 = IF10
        F10B = IF10B
        IX = (TIME-FTIME)*8. + 1
        AP3HR = IAP(IX)
C       COMPUTE AVE DAILY AP
        SUM = 0.
        ISUM = 0
        DO 15 I=1,8
        IF (IAP(I).LT.299) THEN
          ISUM = ISUM + 1
          SUM = SUM + IAP(I)
        ENDIF
 15     CONTINUE
        AP = 1.
        IF (ISUM.GE.1) AP = SUM/ISUM
        IAPAVE = AP
        GO TO 90
      ENDIF
C
      IF (TIME.GT.FTIME) THEN
C       READ NEXT RECORD
        READ (KFUNIT,101,ERR=80) IDAY,IF10,IF10B,(IAP(I),I=1,8),IY
        IF (IY.LT.50) IY = IY + 100
        IYY = (IY-1)/4 - 12
        D1950 = (IY-50)*365 + IYY + IDAY
        IF (D1950-FTIME.NE.1.D0) GO TO 60
        FTIME = D1950
        GO TO 10
      ENDIF
C
      IF (FTIME-TIME.GT.TIME-FT0) THEN
C       START AT BEGINNING OF FILE
        REWIND (KFUNIT)
        IREAD = 0
        GO TO 5
      ENDIF
C
C     REQUESTED TIME LESS THAN CURRENT FILE TIME
      BACKSPACE (KFUNIT)
      BACKSPACE (KFUNIT)
      READ (KFUNIT,101,ERR=80) IDAY,IF10,IF10B,(IAP(I),I=1,8),IY
      IF (IY.LT.50) IY = IY + 100
      IYY = (IY-1)/4 - 12
      D1950 = (IY-50)*365 + IYY + IDAY
      FTIME = D1950
      GO TO 10
C
   60 WRITE(8,105) IY,IDAY
  105 FORMAT(' DATA GAP IN SOLAR DATA',5X,I2,I4)

   80 STOP
C
   90 RETURN
      END SUBROUTINE SOLFLUX
C
C***********************************************************************
C
      SUBROUTINE TMINP(IYR,IMON,IDAY,IHR,IMIN,SEC,D1950)
C
C     CONVERTS YEAR, MONTH, AND DAY OF MONTH TO DAYS SINCE 1950
C
C     Program: This converts a date in year/month/day//hr/min/sec
C       format into Julian date and D1950
C     inputs: iyr=cal year (2 digit between 1950 and 2049)
C             imon=month of year
C             iday=day (day of month)
C             ihr=hour of the day
C             imin=minute of the hour
C             sec=seconds
C     outputs: D1950
C
C     This is based on Algorithm2: Julian date pg 67 in 'Fundamentals of
C     Astrodaynamics and Application by David Vallado.
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     Correction of 2 digit year, only works between 1950 and 2049
C
      if (iyr .gt. 50) then
        ifullyr= iyr + 1900
      else
        ifullyr= iyr + 2000
      endif
C
      k1 = (imon + 9)/12
      k2 = (7*(ifullyr +k1))/4
      k3 = 275*imon/9
      a4 = ((sec/60.d0 + imin)/60.d0 + ihr)/24.d0

C     Calculation of the Julian date
      DJ = 367*ifullyr - k2 + k3 + iday + 1721013.5d0 + a4
      D1950 = DJ - 2433281.5
C
      RETURN
      END SUBROUTINE TMINP
C
C******************************************************************
C
      SUBROUTINE TMOUTD(D1950,IYR,DAY)
C
C     COMPUTE DAY AND YEAR FROM TIME D1950 (DAYS SINCE 1950)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      IYDAY = D1950
      FRAC = D1950 - IYDAY
      IYDAY = IYDAY + 364
      ITEMP = IYDAY/1461
      IYDAY = IYDAY - ITEMP*1461
      IYR = 1949 + 4*ITEMP
      ITEMP = IYDAY/365
      IF (ITEMP.GE.3) ITEMP = 3
      IYR = IYR + ITEMP
      IYDAY = IYDAY - 365*ITEMP + 1
      IYR = IYR - 1900
      DAY = IYDAY + FRAC
      IF (IYR.GE.100) IYR = IYR - 100
C
      RETURN
      END SUBROUTINE TMOUTD
C