# QSUB -lT 80 -lt 60 -q prem -lM 2Mw -eo -s /bin/csh -r getndcs
date
ja
 
cd $TMPDIR
#set QSUB_REQNAME = $0.$$
set ms1 = "/ISRADAR/EHP781102A"
set ms2 = "/ISRADAR/GPI600101A"
set ms3 = "/ISRADAR/IMF631127A"
set ms4 = "/ISRADAR/AEI780101A"
set ms5 = "/ISRADAR/DST570101A"
set ofl = "vishnu.hao:/d/foster/tgcminp/S${QSUB_REQNAME}out"
set ofp = "vishnu.hao:/d/foster/tgcminp/S${QSUB_REQNAME}cgm"
set ofx = "vishnu.hao:/d/foster/tgcminp/S${QSUB_REQNAME}lex"
 
#       Put COS blocked read files in a semi-volitile directory
set fld = "/usr/tmp/$user"
set lf1 = $fld"/"$ms1:t
set lf2 = $fld"/"$ms2:t
set lf3 = $fld"/"$ms3:t
set lf4 = $fld"/"$ms4:t
set lf5 = $fld"/"$ms5:t
if (! -e $fld) mkdir $fld
if (! -e $lf1) msread $lf1 $ms1
if (! -e $lf2) msread $lf2 $ms2
if (! -e $lf3) msread $lf3 $ms3
if (! -e $lf4) msread $lf4 $ms4
if (! -e $lf5) msread $lf5 $ms5
if ($status != 0) exit $status
 
cat << EOF >! in$$
  0 1 0 0 0 1
'$lf1' '$lf2' '$lf3' '$lf4' '$lf5'
1985 0919 0000    0 1985 1006 2400    0 262
EOF
#----------INPUTS ARE:--------------------------------------------------
#   READ IGETHP,IGETGPI,IGETIMF,IGETAE,IGETDST , ICALCT
#   READ NAMEHP, NAMGPI, NAMIMF, NAMAE, NAMDST
#   READ IBEGT, IENDT, LEXBDAY  (or done; i.e., repeat this line for each
#                                time interval to obtain indices)
#  Variable definitions:
#   IGETHP  = 1,0 if do,not get Estimated Hemispheric Power data
#   IGETGPI = 1,0 if do,not get GeoPhysical Indices (SS#,F107,Kp,Ap,..)
#   IGETIMF = 1,0 if do,not get Interplanet Mag Fld data (IMF,VSW,DSW)
#   IGETAE  = 1,0 if do,not get AE hourly data
#   IGETDST = 1,0 if do,not get Dst
#   ICALCT  = 1,0 if do,not calculate crosstail potential and perform
#             lexical write
#   NAMEHP...NAMDST = Name of input files containing indices (apostrophy
#             delimited) must be defined when option is selected; e.g.,
#             if not requesting ae, set IGETAE=0 and NAMAE=' '.
#   IBEGT,IENDT = yyyy mmdd hhmm cs yyyy mmdd hhmm cs defines the UT time
#             interval for this frames plot values
#   LEXBDAY = UT day number in lexical output of the begin day; only
#             relevent when ICALCT=1.
#  Current data availability:
#   HP from  2 Nov 78 to 19 Jul 86
#   GPI from 1 Jan 60 to 31 Jan 92
#   hourly IMF from 27 Nov 63 to 11 May 91
#   hourly AE  from  1 Jan 78 to 31 Dec 86
#   hourly DST from  1 Jan 57 to 31 Dec 90
 
cat << "EOF" >! src$$.f
      PROGRAM GETNDCS
C          This program will retrieve Dave Evan's Hemispheric Power Input
C          values from I.S. database file EHP781102A, the 3 hour values
C          of Kp and ap and daily values of AP, F10.7 flux and sunspot
C          number from I.S. database file GPI600101A, and the hourly
C          Interplanetary Magnetic Field Solar Wind data from NCAR I.S.
C          database IMF631127A.  In addition, it will write out those data for
C          those times that can be used (i.e. when data were recorded) for
C          Hemispheric Power, Crosstail Potential, and the IMF.  The format
C          of the output will allow lexical reads to be made--the first
C          number is the UT date, the second number is the UT hour, the
C          third the UT minute, and the last either POWER, CTPOTEN, or IMF.
C          It will also plot these data versus time.
C
C          This version has been installed on a Cray, and uses the version
C          of rdblk6 which calls rdtape; it does not require the longer
C          work array; see MXLWRK below.
C
C          File linkages:
C          Done as part of the jcl of the batch job:
C                   stdin    = fort.5  - job control = JCTLUN
C                   stdout   = fort.6  - print output and diagnostics = MSGUN
C                   fort.61  = fort.61 - lexical write output file (not
C                                        explicitly opened).
C                   gmeta    = fort.2  - graphics metacode plot output
C          Input data files (NCAR I.S. binary format with Cray block envelope)
C          NOTE:  The fortran unit number file name linkage must be done prior
C          to execution because the OPNBLK entry in the Cray version of
C          RDBLK6 does not attempt to open the named files.
C                   EHP781102A ( = fort.11 ) = IUNEHP
C                   GPI600101A ( = fort.12 ) = IUNGPI
C                   IMF631127A ( = fort.13 ) = IUNIMF
C                   AEI780101A ( = fort.14 ) = IUNAE
C                   DST570101A ( = fort.15 ) = IUNDST
 
C          Additional source files referenced but not included here (all
C          files should be in directory ~/pgms/is/src for this version):
C            rdblk6.f cbfrd.f  timdif.f hours.f  newtim.f cvti2r.f getehp.f
C            gethpd.f getgpd.f getimf.f getae.f  getdst.f gtmtic.f gaxtic.f
C            gnumb.f  pltyvt.f parsb.f
 
C          Define unit number to read job control parameters:
      PARAMETER (JCTLUN=5)
 
C          Declarations for RDBLK6 which are the same for all get__ subs
      PARAMETER (NEOFS=9999 , MSGUN=6 , IOCKSM=0 ,
     +           LMWD=64 , MXVALS=10000 ,
     +           MXWDS=MXVALS*16/LMWD , MXLWRK= 17                   )
CUNIX+           MXWDS=MXVALS*16/LMWD , MXLWRK= 17+20 + 2*512*64/LMWD)
CUNIX      version w/o CBFRD does not need extra work array space
 
C          Declarations for getehp:
      PARAMETER (MXNPWR = 2000)
C                MXNPWR = maximum number of values to be returned by GETEHP
      DIMENSION ITMEHP(4,MXNPWR) ,  HPEHP(MXNPWR) , NDXEHP(MXNPWR) ,
     +            IQLEHP(MXNPWR)
C          Declarations for rdblk6 specifically for getehp
      CHARACTER*80 NAMEHP
      PARAMETER   (IUNEHP=11)
      DIMENSION LREHP(MXVALS) , IUPEHP(MXVALS) , IBKEHP(MXWDS) ,
     +         IWKEHP(MXLWRK)
 
C          Declarations for getgpi:
      PARAMETER (MXNGPI = 40 , MXNG3H = MXNGPI*8)
C                MXNGPI = maximum number of values to be retrieved by GETGPI
      DIMENSION ITM24G(4,MXNGPI) , APDG(MXNGPI) ,  F107G(MXNGPI) ,
     +            F107AG(MXNGPI) , SSNG(MXNGPI) , ISTGPI(MXNGPI) ,
     +          ITM03G(4,MXNG3H) ,   XKPG(MXNG3H)  , APG(MXNG3H) ,
     +         IKP(8) , IAP(8) , I3HT(8) , RF107(MXNGPI+80)
      DATA I3HT/0130 , 0430 , 0730 , 1030 , 1330 , 1630 , 1930 , 2230/
C          Declarations for rdblk6 specifically for getgpi
      CHARACTER*80 NAMGPI
      PARAMETER   (IUNGPI=12)
      DIMENSION LRGPI(MXVALS) , IUPGPI(MXVALS) , IBKGPI(MXWDS) ,
     +         IWKGPI(MXLWRK)
 
C          Declarations for getimf:
      PARAMETER (MXNBY = 750)
C                MXNBY = maximum number of hourly IMF's to be retreived
C                        by GETIMF
      DIMENSION ITMIMF(4,MXNBY) ,
     +              BX(MXNBY)   ,     BY(MXNBY) ,     BZ(MXNBY) ,
     +           BYGSE(MXNBY)   ,  BZGSE(MXNBY) , IQLIMF(MXNBY) ,
     +             SWD(MXNBY)   ,    SWV(MXNBY) , ISTIMF(MXNBY)
C          Declarations for rdblk6 specifically for getimf
      CHARACTER*80 NAMIMF
      PARAMETER   (IUNIMF=13)
      DIMENSION LRIMF(MXVALS) , IUPIMF(MXVALS) , IBKIMF(MXWDS) ,
     +         IWKIMF(MXLWRK)
 
C          Declarations for getae:
      PARAMETER (MXNAE = 750)
      DIMENSION ITMAE(4,MXNAE) , AE(MXNAE), ISTAE(MXNAE)
C          Declarations for rdblk6 specifically for getae
      CHARACTER*80 NAMAE
      PARAMETER   (IUNAE=14)
      DIMENSION LRAE(MXVALS) , IUPAE(MXVALS) , IBKAE(MXWDS) ,
     +         IWKAE(MXLWRK)
 
C          Declarations for getdst:
      PARAMETER (MXNDST = 750)
      DIMENSION ITMDST(4,MXNDST) , DST(MXNDST) , ISTDST(MXNDST)
C          Declarations for rdblk6 specifically for getdst
      CHARACTER*80 NAMDST
      PARAMETER   (IUNDST=15)
      DIMENSION LRDST(MXVALS) , IUPDST(MXVALS) , IBKDST(MXWDS) ,
     +         IWKDST(MXLWRK)
 
C          Miscellaneous declarations
      DIMENSION IBEGT(4) , IENDT(4) , ITMP1(4) , ITMP2(4) ,
     +          IBTHR(4) , IETHR(4) , IBTG(4) , IETG(4)
      DOUBLE PRECISION ICSDIF , IDNHCS , I1HR , I1DA , M40DA
      PARAMETER (MISS=-32767 , RMISS=-32767. , I1HR=360000D0  ,
     +          TMISS=-32766.9 , NOPRT=99999 , I1DA=24D0*I1HR )
 
C          Declarations for parameters found from AE
      DIMENSION CPFAE(MXNAE,2) , NOACP(2) , TIMAE(MXNAE), AENM(MXNAE),
     |          HPFAE(MXNAE)
C          Declarations for parameters found from Hp (P)
      DIMENSION   TIMHP(MXNPWR)  ,  CPFHP(MXNPWR,2), NOHCP(2)
C          Declarations for parameters found from Kp
      DIMENSION    HPFKP(MXNG3H) , CPFKP(MXNG3H,3),  NOKCP(3),
     |             TIMKP(MXNG3H)
C          Declarations for parameters found from IMF (and some others)
      DIMENSION  TMNOG(MXNBY)   ,  SWNOG(MXNBY) ,  SDNOG(MXNBY) ,
     |             POT(MXNBY,3) ,  POTF(MXNBY,3) , NOICP(3), NOICPF(3)
 
C          Get and print desired options and local input file names
      READ (JCTLUN,*) IGETHP,IGETGPI,IGETIMF,IGETAE,IGETDST,ICALCT
C          IGETHP  = 1,0   if do,not get HP data
C          IGETGPI = 1,0   if do,not get GPI data (SS#,F107,Kp,Ap,..)
C          IGETIMF = 1,0   if do,not get IMF data (IMF,VSW,DSW,CPOTEN)
C          IGETAE  = 1,0   if do,not get AE hourly data
C          IGETDST = 1,0   if do,not get Dst
C          ICALCT  = 1,0   if do,not calculate crostail potential and lexical write
      WRITE (MSGUN,'(''1 PROGRAM GETNDCS - OPTIONS SELECTED:'',/,
     +   ''    IGETHP IGETGPI IGETIMF  IGETAE IGETDST  ICALCT'',/,6I8)')
     +         IGETHP,IGETGPI,IGETIMF, IGETAE,IGETDST, ICALCT
      READ (JCTLUN,*) NAMEHP,NAMGPI,NAMIMF,NAMAE,NAMDST
      IF (IGETHP .EQ. 1) THEN
        CALL PARSB(NAMEHP,IBH,IEH)
        WRITE(MSGUN,'(''    Hemispheric power data are read from '',A)
     +                ') NAMEHP(IBH:IEH)
      ENDIF
      IF (IGETGPI .EQ. 1) THEN
        CALL PARSB(NAMGPI,IBG,IEG)
        WRITE(MSGUN,'(''   GPI (SS#,F107,Kp,Ap,..) are read from '',A)
     +                ') NAMGPI(IBG:IEG)
      ENDIF
      IF (IGETIMF .EQ. 1) THEN
        CALL PARSB(NAMIMF,IBI,IEI)
        WRITE(MSGUN,'(''           Hourly IMF data are read from '',A)
     +                ') NAMIMF(IBI:IEI)
      ENDIF
      IF (IGETAE .EQ. 1) THEN
        CALL PARSB(NAMAE,IBA,IEA)
        WRITE(MSGUN,'(''         Hourly ae indices are read from '',A)
     +                ') NAMAE(IBA:IEA)
      ENDIF
      IF (IGETDST .EQ. 1) THEN
        CALL PARSB(NAMDST,IBD,IED)
        WRITE(MSGUN,'(''           Hourly DST data are read from '',A)
     +                ') NAMDST(IBD:IED)
      ENDIF
 
C          Open GKS for plotting
      CALL OPNGKS
 
C          Loop on each time interval requested
    2 CONTINUE
C          Get and print desired time interval or terminate
      IF (ICALCT .EQ. 1) THEN
        READ (JCTLUN,*,END=1000) (IBEGT(I),I=1,4),(IENDT(J),J=1,4)
     +                          , LEXBDAY
        WRITE (MSGUN,'(
     +        '' Begin time (year moday hrmin centisec) ='',4I5,/,
     +        ''   End time                             ='',4I5,/,
     +        ''              and lexical day (LEXBDAY) ='',I3,
     +                   '' matches begin day'',2I5)')
     +     (IBEGT(I),I=1,4),(IENDT(J),J=1,4),LEXBDAY,(IBEGT(K),K=1,2)
      ELSE
        READ (JCTLUN,*,END=1000) (IBEGT(I),I=1,4),(IENDT(J),J=1,4)
        WRITE (MSGUN,'(
     +        '' Begin time (year moday hrmin centisec) ='',4I5,/,
     +        ''   End time                             ='',4I5)')
     +     (IBEGT(I),I=1,4),(IENDT(J),J=1,4)
      ENDIF
C          Check input time
      NOK = 0
      IF(IBEGT(1) .LT. 1940 .OR. IBEGT(1) .GT. 2100)NOK = 1
      IF(IBEGT(2) .LT. 0101 .OR. IBEGT(2) .GT. 1231)NOK = 1
      IF(IBEGT(3) .LT.    0 .OR. IBEGT(3) .GT. 2400)NOK = 1
      IF(IBEGT(4) .LT.    0 .OR. IBEGT(4) .GT. 5999)NOK = 1
      IF(IENDT(1) .LT. 1940 .OR. IENDT(1) .GT. 2100)NOK = 1
      IF(IENDT(2) .LT. 0101 .OR. IENDT(2) .GT. 1231)NOK = 1
      IF(IENDT(3) .LT.    0 .OR. IENDT(3) .GT. 2400)NOK = 1
      IF(IENDT(4) .LT.    0 .OR. IENDT(4) .GT. 5999)NOK = 1
      IF(NOK .GT. 0)GO TO 9000
C          Determine NDAS the number of days indicies requested
C                and NHRTOT the number of hours difference
      CALL TIMDIF(IBEGT,IENDT,ICSDIF)
      IF (ICSDIF .LT. 0D0) GO TO 9050
      NHRTOT = IDINT( (ICSDIF-1D0) / I1HR ) + 1
      DO 3 I=1,2
      IBTG(I) = IBEGT(I)
    3 IETG(I) = IENDT(I)
      IBTG(3) = 0
      IBTG(4) = 0
      IETG(3) = 2400
      IETG(4) = 0
      CALL TIMDIF(IBTG,IETG,ICSDIF)
      NDAS   = IDINT( (ICSDIF-1D0) / I1DA ) + 1
 
C          Get Ae indicies
      IF(IGETAE .EQ. 0) GO TO 100
C          Open NCAR I.S. format (binary) for read
      CALL OPNBLK(IUNAE,NAMAE(IBA:IEA),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,
     +            MXWDS,IWKAE,MXLWRK)
      WRITE(6,'('' --------- AE INDEX ----------'')')
      IF(NHRTOT .GT. MXNAE)GO TO 9250
      NOKAE = 0
      ITMP1(1) = IBEGT(1)
      ITMP1(2) = IBEGT(2)
      ITMP1(3) = 100*(IBEGT(3)/100) + 30
      ITMP1(4) = 0
      ITMP2(4) = 0
      DO 10 I=1,NHRTOT
      ITMAE(1,I) = ITMP1(1)
      ITMAE(2,I) = ITMP1(2)
      ITMAE(3,I) = ITMP1(3)
      ITMAE(4,I) = 0
      ISTAT = NOPRT
      CALL GETAE(MSGUN,ITMP1,324,AE(I) , ISTAT ,
     +           IWKAE,IBKAE,IUPAE,LRAE)
      ISTAE(I) = ISTAT
      IF (AE(I) .GT. TMISS) NOKAE = NOKAE + 1
C         Increment date
      ITMP2(1) = ITMP1(1)
      ITMP2(2) = ITMP1(2)
      ITMP2(3) = ITMP1(3)
   10 CALL NEWTIM(ITMP2,ITMP1,I1HR,MSGUN)
 
C          Print Ae values except when entire date is missing:
      WRITE(6,"('  Time (mid-point)  Ae(hourly)'
     +   ,/,    '  year  mmdd  hrmn      nT')")
      LY  = 0
      LMD = 0
      DO 20 I=1,NHRTOT
      IF(ISTAE(I) .EQ. 6)THEN
        WRITE(6,9600)(ITMAE(J,I),J=1,3)
        GO TO 20
      ENDIF
      IF(ISTAE(I) .NE. 0 .AND.
     +        ITMAE(1,I) .EQ. LY .AND. ITMAE(2,I) .EQ. LMD)GO TO 20
      LY  = ITMAE(1,I)
      LMD = ITMAE(2,I)
      IF(ISTAE(I) .EQ. 2)WRITE(6,9300)(ITMAE(J,I),J=1,2)
      IF(ISTAE(I) .EQ. 3)WRITE(6,9400)(ITMAE(J,I),J=1,2)
      IF(ISTAE(I) .EQ. 0)WRITE(6,'(3I6,F9.0)')(ITMAE(J,I),J=1,3),AE(I)
   20 CONTINUE
      WRITE(6,'('' '',/,'' '')')
  100 CONTINUE
 
C          Get Hemispheric Power Data
      NVEHP = 0
      IF(IGETHP .EQ. 0) GO TO 200
C          Open NCAR I.S. format (binary) for read
      CALL OPNBLK(IUNEHP,NAMEHP(IBH:IEH),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,
     +            MXWDS,IWKEHP,MXLWRK)
      WRITE(6,'('' ------ ESTIMATED HEMISPHERIC POWER INPUT -------'')')
      CALL GETEHP(MSGUN,IBEGT,IENDT,MXNPWR,NVEHP,
     +            ITMEHP,HPEHP,NDXEHP,IQLEHP,ISTEHP,
     +            IWKEHP,IBKEHP,IUPEHP,LREHP)
      WRITE(6,'(''          Requested Time Interval        # Values''/,
     +          '' -------Begin------- --------End--------   Found '',/,
     +            8I5,I8)')  (IBEGT(I),I=1,4),(IENDT(J),J=1,4),NVEHP
      IF(NVEHP .GT. 0)
     +  WRITE (6,'(''  UT Date/Time  Power(GW) PwrIndex  PwrQualifier'',
     +             /,(3I5,F8.1,I7,I13))') ((ITMEHP(I,J),I=1,3),
     +                          HPEHP(J),NDXEHP(J),IQLEHP(J),J=1,NVEHP)
      WRITE(6,'('' '',/,'' '')')
  200 CONTINUE
 
C          Get GPI (Kp, ap, AP, F10.7, Sunspot Number)
      IF (IGETGPI .EQ. 0) GO TO 300
C          Open NCAR I.S. format (binary) for read
      CALL OPNBLK(IUNGPI,NAMGPI(IBG:IEG),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,
     +            MXWDS,IWKGPI,MXLWRK)
      WRITE(6,'('' ------------- GEOPHYSICAL INDICIES ------------'')')
      IF(NDAS .GT. MXNGPI)GO TO 9100
C          Since 81 day average F10.7 is desired, compute it here using
C          GETGPD to fetch a days worth of values for each call rather
C          than using GETGPI which would force a rewind for each new
C          81 day average.
      NTOTD = NDAS + 80
      IBSAV = 41
      IESAV = NDAS + 40
      M40DA = -I1DA* 40D0
      CALL NEWTIM(IBEGT,ITMP1,M40DA,MSGUN)
      ITMP1(3) = 0
      ITMP1(4) = 0
      ITMP2(3) = 0
      ITMP2(4) = 0
      ISAVD = 0
      NOKAPD = 0
      NOKF10 = 0
      NOKF1A = 0
      NOKSSN = 0
      NOKKP  = 0
      NOKAP3 = 0
      DO 220 I=1,NTOTD
      ISTAT = NOPRT
      CALL GETGPD(MSGUN,ITMP1,IKP,IAP,IAPD,IF107,IFQFR,ISSN,
     +             ISTAT,IWKGPI,IBKGPI,IUPGPI,LRGPI)
C          Always save F10.7 values for 81 day averages
      RF107(I) = RMISS
      IF(IF107 .NE. MISS)RF107(I) = FLOAT(IF107) / 10.
      IF(I .GE. IBSAV .AND. I .LE. IESAV)THEN
C          Save daily average values; hr-min is fixed at 1200,
C          representing the mean time of the daily average values.
        ISAVD = ISAVD + 1
        ISTGPI(ISAVD) = ISTAT
        ITM24G(1,ISAVD) = ITMP1(1)
        ITM24G(2,ISAVD) = ITMP1(2)
        ITM24G(3,ISAVD) = 1200
        ITM24G(4,ISAVD) = 0000
        APDG  (ISAVD) = FLOAT(IAPD)
        F107G (ISAVD) = RF107(I)
        SSNG  (ISAVD) = FLOAT(ISSN)
        IF(IAPD  .NE. MISS) NOKAPD = NOKAPD + 1
        IF(IF107 .NE. MISS) NOKF10 = NOKF10 + 1
        IF(ISSN  .NE. MISS) NOKSSN = NOKSSN + 1
C          Get 3hr average indicies; hr-min is fixed to 0130,0430...2230
C          representing the mean time of day for each index value.
        IB3 = (ISAVD-1)*8 + 1
        IE3 = IB3 + 7
        J = 0
        DO 210 I3H=IB3,IE3
        J = J + 1
        ITM03G(1,I3H) = ITMP1(1)
        ITM03G(2,I3H) = ITMP1(2)
        ITM03G(3,I3H) = I3HT(J)
        ITM03G(4,I3H) = 0000
        APG (I3H) = FLOAT(IAP(J))
        XKPG(I3H) = RMISS
        IF(IKP(J) .NE. MISS)THEN
          XKPG(I3H) = FLOAT(IKP(J)) / 10.
          NOKKP = NOKKP + 1
        ENDIF
        IF(IAP(J) .NE. MISS) NOKAP3 = NOKAP3 + 1
  210   CONTINUE
      ENDIF
C         Increment date
      ITMP2(1) = ITMP1(1)
      ITMP2(2) = ITMP1(2)
  220 CALL NEWTIM(ITMP2,ITMP1,I1DA,MSGUN)
C          Compute F10.7 solar flux 81 day average from saved values
      DO 240 I=1,NDAS
      SUM = 0.
      CNT = 0.
      DO 230 J=I,I+80
      IF(RF107(J) .LT. TMISS)GO TO 230
      SUM = SUM + RF107(J)
      CNT = CNT + 1.
  230 CONTINUE
      F107AG(I) = RMISS
      IF(CNT .GT. 40.)THEN
        F107AG(I) = SUM / CNT
        NOKF1A = NOKF1A + 1
      ENDIF
  240 CONTINUE
C          Establish the number of 3hr GPI assigned
      N3HRG = NDAS*8
C          Print GPI values except when entire date is missing:
      WRITE(6,'(''  Time (mid-point)  Daily   81day   Daily  Daily'',/,
     +          ''  year  mmdd  hrmn   SS#    F10.7   F10.7    AP'')')
      DO 250 I=1,NDAS
      IF(ISTGPI(I) .EQ. 2)WRITE(6,9300)(ITM24G(J,I),J=1,2)
      IF(ISTGPI(I) .EQ. 3)WRITE(6,9400)(ITM24G(J,I),J=1,2)
      IF(ISTGPI(I) .EQ. 6)WRITE(6,9600)(ITM24G(J,I),J=1,2)
      IF(ISTGPI(I) .EQ. 0)WRITE(6,'(3I6,F7.0,2F8.1,F7.0)')
     +        (ITM24G(J,I),J=1,3),SSNG(I),F107AG(I),F107G(I),APDG(I)
  250 CONTINUE
      WRITE(6,'('' '',/,
     +          ''  Time (mid-point)      3hr    3hr'',/,
     +          ''  year  mmdd  hrmn       Kp     ap'')')
      DO 260 I=1,NDAS
      I1 = (I-1)*8 + 1
      IF(ISTGPI(I) .EQ. 2)WRITE(6,9300)(ITM03G(J,I1),J=1,2)
      IF(ISTGPI(I) .EQ. 3)WRITE(6,9400)(ITM03G(J,I1),J=1,2)
      IF(ISTGPI(I) .EQ. 6)WRITE(6,9600)(ITM03G(J,I1),J=1,2)
      IF(ISTGPI(I) .NE. 0)GO TO 260
      WRITE(6,'(3I6,F9.1,F7.0)')(ITM03G(J,I1),J=1,3),XKPG(I1),APG(I1)
      I2 = I1+1
      I8 = I1+7
      WRITE(6,'(12X,I6,F9.1,F7.0)')(ITM03G(3,J),XKPG(J),APG(J),J=I2,I8)
  260 CONTINUE
      WRITE(6,'('' '',/,'' '')')
  300 CONTINUE
 
C          Get IMF and solar wind data
      IF (IGETIMF .EQ. 0) GO TO 400
C          Open NCAR I.S. format (binary) for read
      CALL OPNBLK(IUNIMF,NAMIMF(IBI:IEI),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,
     +            MXWDS,IWKIMF,MXLWRK)
      WRITE(6,'('' ----------------------------------- IMF AND SOLAR'',
     +                '' WIND ------------------------------------'')')
      IF(NHRTOT .GT. MXNBY)GO TO 9200
      NOKSD = 0
      NOKSV = 0
      NOKBX = 0
      NOKBY = 0
      NOKBZ = 0
      SDMIN = 0.
      SDMAX = 0.
      SVMIN = 1.E6
      SVMAX = 0.
      BMIN  = 1.E6
      BMAX  =-1.E6
      ITMP1(1) = IBEGT(1)
      ITMP1(2) = IBEGT(2)
      ITMP1(3) = 100*(IBEGT(3)/100) + 30
      ITMP1(4) = 0
      ITMP2(4) = 0
      DO 310 I=1,NHRTOT
      ITMIMF(1,I) = ITMP1(1)
      ITMIMF(2,I) = ITMP1(2)
      ITMIMF(3,I) = ITMP1(3)
      ITMIMF(4,I) = 0
      ISTAT = NOPRT
      CALL GETIMF(MSGUN,ITMP1,
     +            BX(I),BY(I),BZ(I),BYGSE(I),BZGSE(I),SWDEN,SWSPD,
     +            IQLIMF(I),ISTAT,IWKIMF,IBKIMF,IUPIMF,LRIMF)
      ISTIMF(I)= ISTAT
      SWD(I)   = SWDEN/1000000.
      SWV(I)   = SWSPD/1000.
      IF(SWDEN .LT. TMISS)SWD(I) = RMISS
      IF(SWSPD .LT. TMISS)SWV(I) = RMISS
C          Determine plot stuff:  # ok values and domain
      IF(SWD(I) .GT. TMISS)THEN
        NOKSD = NOKSD + 1
        SDMAX = AMAX1(SWD(I),SDMAX)
      ENDIF
      IF(SWV(I) .GT. TMISS)THEN
        NOKSV = NOKSV + 1
        SVMIN = AMIN1(SWV(I),SVMIN)
        SVMAX = AMAX1(SWV(I),SVMAX)
      ENDIF
      IF(BX(I) .GT. TMISS)THEN
        NOKBX = NOKBX + 1
        IF(BX(I) .LT. BMIN) BMIN = BX(I)
        IF(BX(I) .GT. BMAX) BMAX = BX(I)
      ENDIF
      IF(BY(I) .GT. TMISS)THEN
        NOKBY = NOKBY + 1
        IF(BY(I) .LT. BMIN) BMIN = BY(I)
        IF(BY(I) .GT. BMAX) BMAX = BY(I)
      ENDIF
      IF(BZ(I) .GT. TMISS)THEN
        NOKBZ = NOKBZ + 1
        IF(BZ(I) .LT. BMIN) BMIN = BZ(I)
        IF(BZ(I) .GT. BMAX) BMAX = BZ(I)
      ENDIF
C         Increment date
      ITMP2(1) = ITMP1(1)
      ITMP2(2) = ITMP1(2)
      ITMP2(3) = ITMP1(3)
  310 CALL NEWTIM(ITMP2,ITMP1,I1HR,MSGUN)
C          Print IMF values except when entire date is missing:
      WRITE(6,'(
     +''      Mean     ------- GSM [nT] -------  --- GSE [nT] ---   '',
     +                                ''    Solar Wind'',/,
     +'' ---- Time ----     Bx       By       Bz       By       Bz   '',
     +                          ''   [#/cm3]  [km/s]    Qualifier'')')
      LY  = 0
      LMD = 0
      DO 320 I=1,NHRTOT
      IF(ISTIMF(I) .NE. 0 .AND.
     +        ITMIMF(1,I) .EQ. LY .AND. ITMIMF(2,I) .EQ. LMD)GO TO 320
      LY  = ITMIMF(1,I)
      LMD = ITMIMF(2,I)
      IF(ISTIMF(I) .EQ. 2)WRITE(6,9300)(ITMIMF(J,I),J=1,2)
      IF(ISTIMF(I) .EQ. 3)WRITE(6,9400)(ITMIMF(J,I),J=1,2)
      IF(ISTIMF(I) .EQ. 4)WRITE(6,9500)(ITMIMF(J,I),J=1,2)
      IF(ISTIMF(I) .EQ. 0)WRITE(6,'(3I5,5F9.2,F12.1,F9.1,I7)')
     +        (ITMIMF(J,I),J=1,3) , BX(I),BY(I),BZ(I),
     +        BYGSE(I),BZGSE(I),SWD(I),SWV(I),IQLIMF(I)
  320 CONTINUE
      WRITE(6,'('' '',/,'' '')')
  400 CONTINUE
 
C          Get DST indicies
      IF(IGETDST .EQ. 0)GO TO 500
C          Open NCAR I.S. format (binary) for read
      CALL OPNBLK(IUNDST,NAMDST(IBD:IED),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,
     +            MXWDS,IWKDST,MXLWRK)
      WRITE(6,'('' --------- DST INDEX ---------'')')
      IF(NHRTOT .GT. MXNDST)GO TO 9260
      NOKDST = 0
      ITMP1(1) = IBEGT(1)
      ITMP1(2) = IBEGT(2)
      ITMP1(3) = 100*(IBEGT(3)/100) + 30
      ITMP1(4) = 0
      ITMP2(4) = 0
      DO 410 I=1,NHRTOT
      ITMDST(1,I) = ITMP1(1)
      ITMDST(2,I) = ITMP1(2)
      ITMDST(3,I) = ITMP1(3)
      ITMDST(4,I) = 0
      ISTAT = NOPRT
      CALL GETDST(MSGUN,ITMP1,DST(I),  ISTAT ,
     +                 IWKDST,IBKDST,IUPDST,LRDST)
      ISTDST(I) = ISTAT
      IF(DST(I) .GT. TMISS) NOKDST = NOKDST + 1
C         Increment date
      ITMP2(1) = ITMP1(1)
      ITMP2(2) = ITMP1(2)
      ITMP2(3) = ITMP1(3)
  410 CALL NEWTIM(ITMP2,ITMP1,I1HR,MSGUN)
C          Print DST values except when entire date is missing:
      WRITE(6,'(''  Time (mid-point) DST(hourly)'',/,
     +          ''  year  mmdd  hrmn      nT '')')
      LY  = 0
      LMD = 0
      DO 420 I=1,NHRTOT
      IF(ISTDST(I) .EQ. 6)THEN
        WRITE(6,9600)(ITMDST(J,I),J=1,3)
        GO TO 420
      ENDIF
      IF(ISTDST(I) .NE. 0 .AND.
     +        ITMDST(1,I) .EQ. LY .AND. ITMDST(2,I) .EQ. LMD)GO TO 420
      LY  = ITMDST(1,I)
      LMD = ITMDST(2,I)
      IF(ISTDST(I) .EQ. 2)WRITE(6,9300)(ITMDST(J,I),J=1,2)
      IF(ISTDST(I) .EQ. 3)WRITE(6,9400)(ITMDST(J,I),J=1,2)
      IF(ISTDST(I) .EQ. 0)WRITE(6,'(3I6,F9.0)')
     +                                 (ITMDST(J,I),J=1,3),DST(I)
  420 CONTINUE
      WRITE(6,'('' '',/,'' '')')
  500 CONTINUE
 
C          Done reading indicies
 
C          If requested calculate other parameters from the indicies
      IF (ICALCT .EQ. 0) GO TO 599
      IF (IGETAE .EQ. 1) THEN
       NOKAE = 0
       DO 510 I=1,NHRTOT
       HPFAE(I) = RMISS
       CPFAE(I,1) = RMISS
       CPFAE(I,2) = RMISS
       IF(AE(I) .GT. TMISS) THEN
C  Calculate hemispheric power from AE using Maeda's formula for 1 hr AE
        HPFAE(I) = 4.28 + 0.0624 * AE(I)
C  Calculate CP from P using CP=max(25,-48.05+35.82*ln(P))
        CPFAE(I,1) = AMAX1( -48.05 + 35.82*ALOG(HPFAE(I)), 25. )
C  Calculate CP from AE using Reiff, 1981, JGR, p 7639 (corr coeff 55%)
        CPFAE(I,2) = 41. + 0.11 * AE(I)
C  Refind NOKAE so as to define AENM and TIMAE
        NOKAE = NOKAE + 1
        AENM(NOKAE) = AE(I)
        TIMAE(NOKAE) = MOD(ITMAE(2,I),100)*24 + ITMAE(3,I)/100 +
     |    FLOAT( MOD(ITMAE(3,I),100)) / 60.
       ENDIF
  510  CONTINUE
      ENDIF
 
      IF (IGETHP .EQ. 1 .AND. NVEHP .GT. 0) THEN
       DO 520 I=1,NVEHP
       TIMHP(I) = MOD(ITMEHP(2,I),100)*24 + ITMEHP(3,I)/100 +
     |    FLOAT( MOD(ITMEHP(3,I),100)) / 60.
       CPFHP(I,1) = 16. + 1.3*HPEHP(I)
       CPFHP(I,2) = AMAX1( 25., -48.05+35.82*ALOG(HPEHP(I)) )
  520  CONTINUE
      ENDIF
 
      IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) THEN
        DO 530 I=1,N3HRG
        TIMKP(I) = MOD(ITM03G(2,I),100)*24 + ITM03G(3,I)/100 +
     +                            FLOAT( MOD(ITM03G(3,I),100)) / 60.
C  Median value of P for Kp=0 from Maeda plots is about 5 GW.
        HPFKP(I) = AMAX1(3.0, -2.78+9.33*XKPG(I) )
        CPFKP(I,1) = AMAX1(-4.21+22.47*ALOG(HPFKP(I)), 10.)
        CPFKP(I,2) = AMAX1(-48.05+35.82*ALOG(HPFKP(I)), 25.)
C  From Reiff (71% correlation coeff), Reiff, JGR, 1981, p7639 (refers
C    to correlation coeff but not to formula
        CPFKP(I,3) = 29. + 11. * XKPG(I)
  530   CONTINUE
      ENDIF
 
C          Must do this last, because will fill in IMF gaps with other
C          parameters found previously, i.e.,  CP(AE,P,Kp)
      NUMSW = 0
      DO 535 I=1,NHRTOT
      DO 535 K=1,3
      POT (I,K) = RMISS
  535 POTF(I,K) = RMISS
      IF (IGETIMF .EQ. 1 .AND. NOKBZ .GT. 0) THEN
        NUMSW = 0
        DO 540 I=1,NHRTOT
        DO 538 K=1,3
        POT (I,K) = RMISS
  538   POTF(I,K) = RMISS
C          Save all non-missing solar wind speed and density separately
        IF(SWV(I) .GT. TMISS)THEN
          NUMSW = NUMSW + 1
          TMNOG(NUMSW) = MOD(ITMIMF(2,I),100)*24 + ITMIMF(3,I)/100 +
     +                            FLOAT( MOD(ITMIMF(3,I),100)) / 60.
          SWNOG(NUMSW) = SWV(I)
          SDNOG(NUMSW) = SWD(I)
        ENDIF
  540   CONTINUE
        DO 550 I=1,NHRTOT
        TM = MOD(ITMIMF(2,I),100)*24 + ITMIMF(3,I)/100 +
     +                      FLOAT( MOD(ITMIMF(3,I),100)) / 60.
        SWIND = SWV(I)
        SWDEN = SWD(I)
        IF (SWIND.LT.TMISS .AND. BZ(I).GT.TMISS .AND. NUMSW.GT.0) THEN
C          Interpolate Vsw and Nsw values if within 24 hours
          SWIND = TERP(TMNOG,SWNOG,TM,NUMSW,24.,IFLGWSP)
          SWDEN = TERP(TMNOG,SDNOG,TM,NUMSW,24.,IFLGDSP)
          IF (IFLGWSP .EQ. 1) SWIND = SWV(I)
          IF (IFLGDSP .EQ. 1) SWDEN = SWD(I)
        ENDIF
        FHP = RMISS
        FAE = RMISS
        FKP = RMISS
C  Find interpolated P value if within 1 hour
        IF (NVEHP .GT. 0) FHP = TERP(TIMHP,HPEHP,TM,NVEHP,1.0,IFLGSP)
C  Find interpolated AE value if within 2 hours
        IF (IGETAE .EQ. 1 .AND. NOKAE .GT. 0)
     +                    FAE = TERP(TIMAE,AENM ,TM,NOKAE,2. ,IFLGAE)
C  Find interpolated Kp value if within 3 hours
        IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0)
     +                    FKP = TERP(TIMKP,XKPG ,TM,N3HRG,3. ,IFLGSP)
C  FINDPOT calculates the crosstail potential
        CALL FINDPOT(BX(I),BY(I),BZ(I),SWDEN,SWIND,FHP,FAE,FKP,
     +                                    POT(I,1),POT(I,2),POT(I,3))
  550   CONTINUE
C         Subroutine FILPOT filters the data for the crosstail potential
        CALL FILPOT(POT,POTF,NHRTOT,0,POT)
        CALL FILPOT(POT(1,2),POTF(1,2),NHRTOT,0,POT)
        CALL FILPOT(POT(1,3),POTF(1,3),NHRTOT,1,POT)
      ENDIF
 
C  Print out indicies and calculated values on 6 and lexical reads on 61
C
      WRITE (6,'(''1'')')
      WRITE (61,'(''C  lexical day '',I3,'' = year'',I5,'' month'',I3,
     +  '' day'',I3)') LEXBDAY,IBEGT(1),IBEGT(2)/100,MOD(IBEGT(2),100)
      IMODA1 = MOD(IBEGT(2),100)
      LEXADD = LEXBDAY - IMODA1
 
      IF (IGETDST .EQ. 1 .AND. NOKDST .GT. 0) THEN
       CALL LEXWRT (' ','DST =', 0, DST, ITMDST, NHRTOT, 4, LEXADD)
      ENDIF
 
      IF (IGETAE .EQ. 1 .AND. NOKAE .GT. 0) THEN
       WRITE(6,"('  year  mmdd  hrmn  Ae(hourly)  calc P  CP (P)',3X,
     +  'CP (AE)')")
       DO 560 I=1,NHRTOT
       IF(AE(I) .GT. TMISS)WRITE(6,'(3I6,F9.0,3F9.1)')(ITMAE(J,I),J=1,3)
     |  ,AE(I),HPFAE(I),CPFAE(I,1),CPFAE(I,2)
  560  CONTINUE
       CALL LEXWRT ('C  Hourly AE (nT)', 'AE =', 0, AE, ITMAE, NHRTOT,
     |  4, LEXADD)
       CALL LEXWRT ('C  P from AE (nT), P=4.28+0.0624*AE', 'CALC P =',
     |  1, HPFAE, ITMAE, NHRTOT, 4, LEXADD)
       CALL LEXWRT ('C  CP from P(AE), CP=max(25,-48.05+35.82*ln(P)',
     |  'CALC CP =', 1, CPFAE(1,1), ITMAE, NHRTOT, 4, LEXADD)
       CALL LEXWRT ('C  CP from AE, CP=41+0.11*AE',
     |  'CALC CP =', 1, CPFAE(1,2), ITMAE, NHRTOT, 4, LEXADD)
      ENDIF
 
      IF (IGETHP .EQ. 1 .AND. NVEHP .GT. 0) THEN
       WRITE (6,"(1X/1X,' year mmdd hrmn   P(GW)   CP(P) CP(lnP)')")
       DO 570 I=1,NVEHP
       WRITE (6,"(1X,3I5,3F8.1)") (ITMEHP(J,I),J=1,3),HPEHP(I),
     |  CPFHP(I,1),CPFHP(I,2)
  570  CONTINUE
       CALL LEXWRT ('C  Estimated hemispheric power (GW)', 'POWER =', 1,
     |  HPEHP, ITMEHP, NVEHP, 4, LEXADD)
       CALL LEXWRT ('C  CP from P, CP = 16 + 1.3*P', 'CALC CP =', 1,
     |  CPFHP(1,1), ITMEHP, NVEHP, 4, LEXADD)
       CALL LEXWRT ('C  CP from P, CP = MAX(25.,-48.05+35.82*LN(P))',
     |  'CALC CP =', 1, CPFHP(1,2), ITMEHP, NVEHP, 4, LEXADD)
      ENDIF
 
      IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) THEN
       WRITE (6,"(1X/1X,' year mmdd hrmn      Kp   P(Kp) CP lnPs',
     |  '  CP lnP CP (Kp)')")
       DO 580 I=1,N3HRG
       WRITE (6,"(1X,3I5,5F8.1)") (ITM03G(J,I),J=1,3),XKPG(I),
     |  HPFKP(I),(CPFKP(I,K),K=1,3)
  580  CONTINUE
       CALL LEXWRT (' ','SUNSPOTS =', 0, SSNG, ITM24G, NDAS, 4,LEXADD)
       CALL LEXWRT (' ','F107A =', 1, F107AG, ITM24G, NDAS, 4, LEXADD)
       CALL LEXWRT (' ','F107 =', 1, F107G, ITM24G, NDAS, 4, LEXADD)
       CALL LEXWRT (' ','APD =', 0, APDG, ITM24G, NDAS, 4, LEXADD)
       CALL LEXWRT (' ','KP =', 1, XKPG, ITM03G, N3HRG, 4, LEXADD)
       CALL LEXWRT ('C  P from Kp, P = max(0.1,-2.78+9.33*Kp)',
     |  'CALC P =', 1, HPFKP, ITM03G, N3HRG, 4, LEXADD)
       CALL LEXWRT ('C  CP from P(Kp), CP = max(10.,-4.21+22.47*LN(P))',
     |  'CALC CP =', 1, CPFKP(1,1), ITM03G, N3HRG, 4, LEXADD)
       CALL LEXWRT ('C  CP from P(Kp), CP = max(25.,-48.05+35.82*LN(P))'
     |  ,'CALC CP =', 1, CPFKP(1,2), ITM03G, N3HRG, 4, LEXADD)
       CALL LEXWRT ('C  CP from Kp, CP = 29 + 11 * Kp',
     |  'CALC CP =', 1, CPFKP(1,3), ITM03G, N3HRG, 4, LEXADD)
      ENDIF
 
      IF (IGETIMF .EQ. 1 .AND. NOKBZ .GT. 0) THEN
       WRITE (6,"(1X/1X,' year mmdd hrmn      Bx      By      Bz',5X,
     | 'Vsw     Vden CPuVBQ1 CPus3ng CP vBs3 CP VBQ1 CP s3ng')")
       DO 590 I=1,NHRTOT
       IF (BZ(I) .LT. TMISS) GO TO 590
       WRITE (6,"(1X,3I5,10F8.1)") (ITMIMF(J,I),J=1,3),BX(I),BY(I),
     |  BZ(I),SWV(I),SWD(I),(POT(I,J),J=1,2),(POTF(I,K),K=1,3)
  590  CONTINUE
       CALL LEXWRT (' ','BXIMF =', 1, BX, ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT (' ','BYIMF =', 1, BY, ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT (' ','BZIMF =', 1, BZ, ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT (' ','VSW =', 0, SWV, ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT (' ','DENSW =', 1, SWD, ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT ('C  Filtered CP from vBsin3 IMF with gaps',
     |  'CALC CP =', 1, POTF(1,1), ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT ('C  Filtered CP from vBQ1 IMF with gaps',
     |  'CALC CP =', 1, POTF(1,2), ITMIMF, NHRTOT, 4, LEXADD)
       CALL LEXWRT ('C  Filtered CP from vBsin3 IMF with gaps filled wit
     |h CP(AE,P,Kp)', 'CALC CP =', 1, POTF(1,3), ITMIMF,NHRTOT,4,LEXADD)
      ENDIF
  599 CONTINUE
 
C          Build plot time axis begin (and end) time which is rounded
C          down (and up) to next even hour.
      IBTHR(1) = IBEGT(1)
      IBTHR(2) = IBEGT(2)
      IBTHR(3) = 100*(IBEGT(3)/100)
      IBTHR(4) = 0
      IDNHCS = DBLE(NHRTOT) * I1HR
      CALL NEWTIM(IBTHR,IETHR,IDNHCS,MSGUN)
 
C          Plot Ae
      IF(IGETAE .EQ. 0)GO TO 600
      IF(NOKAE  .EQ. 0)GO TO 600
      AEMIN = 0.
      AEMAX = RMISS
      CALL PLTYVT('Hourly AE Index','[nT]',AE,ITMAE,MXNAE,NHRTOT,
     +             AEMIN,AEMAX,IBTHR,IETHR)
  600 CONTINUE
C          Plot DST
      IF(IGETDST .EQ. 0)GO TO 700
      DSMIN = RMISS
      DSMAX = RMISS
      IF(NOKDST .GT. 0)CALL PLTYVT('Hourly Dst Index','[nT]',
     +                 DST,ITMDST,MXNDST,NHRTOT,DSMIN,DSMAX,IBTHR,IETHR)
  700 CONTINUE
C          Plot IMF Bx By and Bz
      IF(IGETIMF .EQ. 0)GO TO 800
      IF(NOKBX .GT. 0)CALL PLTYVT('Hourly IMF - Bx','[nT]',
     +                BX,ITMIMF,MXNIMF,NHRTOT,BMIN,BMAX,IBTHR,IETHR)
      IF(NOKBY .GT. 0)CALL PLTYVT('Hourly IMF - By','[nT]',
     +                BY,ITMIMF,MXNIMF,NHRTOT,BMIN,BMAX,IBTHR,IETHR)
      IF(NOKBZ .GT. 0)CALL PLTYVT('Hourly IMF - Bz','[nT]',
     +                BZ,ITMIMF,MXNIMF,NHRTOT,BMIN,BMAX,IBTHR,IETHR)
      IF(NOKSV .GT. 0)CALL PLTYVT('Hourly solar wind speed','[km/s]',
     +                SWV,ITMIMF,MXNIMF,NHRTOT,SVMIN,SVMAX,IBTHR,IETHR)
      IF(NOKSD .GT. 0)CALL PLTYVT('Hourly solar wind density','[#/cm3]',
     +                SWD,ITMIMF,MXNIMF,NHRTOT,SDMIN,SDMAX,IBTHR,IETHR)
  800 CONTINUE
 
C          Lenhart tape indicies Kp, ap, Ap, F10.7, Sunspot Number
      IF(IGETGPI .EQ. 0)GO TO 900
C          Build plot time axis begin and end times for GPI parameters
C          which extend backwards to 0 UT and forwards to 24 UT
      IBTG(1) = ITM24G(1,1)
      IBTG(2) = ITM24G(2,1)
      IBTG(3) = 0
      IBTG(4) = 0
      IETG(1) = ITM24G(1,NDAS)
      IETG(2) = ITM24G(2,NDAS)
      IETG(3) = 2400
      IETG(4) = 0
      ADMIN = 0.
      ADMAX = RMISS
      IF(NOKAPD .GT. 0)CALL PLTYVT('Daily Ap Index',' ',
     +                 APDG,ITM24G,MXNGPI,NDAS,ADMIN,ADMAX,IBTG,IETG)
      F1MIN = RMISS
      F1MAX = RMISS
      IF(NOKF10 .GT. 0)CALL PLTYVT('Daily 10.7 cm Solar Flux',
     +                             '[W/m2/Hz*10-23]',
     +                F107G,ITM24G,MXNGPI,NDAS,F1MIN,F1MAX,IBTG,IETG)
      FAMIN = RMISS
      FAMAX = RMISS
      IF(NOKF1A .GT. 0)CALL PLTYVT('81 Day 10.7 cm Solar Flux',
     +                             '[W/m2/Hz*10-23]',
     +               F107AG,ITM24G,MXNGPI,NDAS,FAMIN,FAMAX,IBTG,IETG)
      SSMIN = RMISS
      SSMAX = RMISS
      IF(NOKSSN .GT. 0)CALL PLTYVT('Daily Sunspot Number',' ',
     +                 SSNG,ITM24G,MXNGPI,NDAS,SSMIN,SSMAX,IBTG,IETG)
      XKMIN = 0.
      XKMAX = 9.
      IF(NOKKP  .GT. 0)CALL PLTYVT('3-hour Kp',' ',
     +                 XKPG,ITM03G,MXNG3H,N3HRG,XKMIN,XKMAX,IBTG,IETG)
      A3MIN = 0.
      A3MAX = RMISS
      IF(NOKAP3 .GT. 0)CALL PLTYVT('3-hour ap',' ',
     +                  APG,ITM03G,MXNG3H,N3HRG,A3MIN,A3MAX,IBTG,IETG)
  900 CONTINUE
 
C          Determine maximum power (GW) for consistent plot axis labeling
C          for requested power parameters:  Evans HP, HP from Kp
      GWMIN = 0.
      GWMAX = RMISS
      NOKEHP = NVEHP
      IF (ICALCT .EQ. 1) THEN
        IF (IGETHP .EQ. 1) THEN
          NOKEHP = 0
          DO 910 I=1,NVEHP
          IF(HPEHP(I) .LT. TMISS)GO TO 910
          NOKEHP = NOKEHP + 1
          GWMAX = AMAX1(HPEHP(I),GWMAX)
  910     CONTINUE
        ENDIF
        IF (IGETGPI .EQ. 1) THEN
          NOKCHP = 0
          IF (NOKKP .GT. 0) THEN
            DO 920 I=1,N3HRG
            IF (HPFKP(I) .GT. TMISS) THEN
              NOKCHP = NOKCHP + 1
              GWMAX = AMAX1(HPFKP(I),GWMAX)
            ENDIF
  920       CONTINUE
          ENDIF
        ENDIF
        IF (IGETAE .EQ. 1) THEN
          NOKAHP = 0
          DO 930 I=1,NHRTOT
          IF (HPFAE(I) .GT. TMISS) THEN
            NOKAHP = NOACHP + 1
            GWMAX = AMAX1(HPFAE(I),GWMAX)
          ENDIF
  930     CONTINUE
        ENDIF
      ENDIF
 
C          Plot the "GW" power units parameters
      IF(IGETHP .EQ. 1 .AND. NOKEHP .GT. 0)
     +                 CALL PLTYVT('Estimated Hemispheric Power','[GW]',
     +                HPEHP,ITMEHP,MXNPWR,NVEHP,GWMIN,GWMAX,IBTHR,IETHR)
 
C         Plot the calculated parameters if requested
      IF (ICALCT .EQ. 0) GO TO 980
      IF(IGETGPI .EQ. 1 .AND. NOKCHP .GT. 0)
     +  CALL PLTYVT('H. Power = max(0.1,-2.78+9.33Kp)','[GW]',
     +    HPFKP,ITM03G,MXNG3H,N3HRG,GWMIN,GWMAX,IBTHR,IETHR)
      IF(IGETAE .EQ. 1 .AND. NOKAHP .GT. 0)
     +             CALL PLTYVT('H. Power = 4.28+0.0624AE','[GW]',
     +                HPFAE,ITMAE,MXNAE,NHRTOT,GWMIN,GWMAX,IBTHR,IETHR)
 
C          Determine maximum potential (KV) for consistent plot axis
C          labeling for requested parameters: 3 versions of CALPOT
      VKMIN = 0.
      VKMAX = RMISS
      NOKCP(1) = 0
      NOKCP(2) = 0
      NOKCP(3) = 0
      IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) THEN
        DO 940 K=1,3
        DO 940 I=1,N3HRG
        IF (CPFKP(I,K) .GT. TMISS) THEN
          NOKCP(K) = NOKCP(K) + 1
          VKMAX = AMAX1(CPFKP(I,K),VKMAX)
        ENDIF
  940   CONTINUE
      ENDIF
      NOACP(1) = 0
      NOACP(2) = 0
      IF (IGETAE .EQ. 1) THEN
        DO 950 K=1,2
        DO 950 I=1,NHRTOT
        IF(CPFAE(I,K) .GT. TMISS)THEN
          NOACP(K) = NOACP(K) + 1
          VKMAX = AMAX1(CPFAE(I,K),VKMAX)
        ENDIF
  950   CONTINUE
      ENDIF
      NOHCP(1) = 0
      NOHCP(2) = 0
      IF (IGETHP .EQ. 1) THEN
        DO 960 K=1,2
        DO 960 I=1,NVEHP
        IF(CPFHP(I,K) .GT. TMISS)THEN
          NOHCP(K) = NOHCP(K) + 1
          VKMAX = AMAX1(CPFHP(I,K),VKMAX)
        ENDIF
  960   CONTINUE
      ENDIF
      NOICP(1) = 0
      NOICP(2) = 0
      NOICP(3) = 0
      NOICPF(1) = 0
      NOICPF(2) = 0
      NOICPF(3) = 0
      IF (IGETIMF .EQ. 1) THEN
        DO 970 K=1,3
        DO 970 I=1,NHRTOT
        IF(POT(I,K) .GT. TMISS)THEN
          NOICP(K) = NOICP(K) + 1
          VKMAX = AMAX1(POT(I,K),VKMAX)
        ENDIF
        IF(POTF(I,K) .GT. TMISS)THEN
          NOICPF(K) = NOICPF(K) + 1
          VKMAX = AMAX1(POTF(I,K),VKMAX)
        ENDIF
  970   CONTINUE
      ENDIF
 
C          Plot the "KV" potential units parameters
      IF(IGETGPI .EQ. 1)THEN
        IF(NOKCP(1).GT. 0)CALL PLTYVT('P(Kp)  CP = max{10,-4+22Ln(P)}',
     +  '[KV]',CPFKP,ITM03G,MXNGPI,N3HRG,VKMIN,VKMAX,
     +  IBTHR,IETHR)
        IF(NOKCP(2).GT.0)CALL PLTYVT('P(Kp)  CP = max{25,-48+36Ln(P)}',
     +  '[KV]',CPFKP(1,2),ITM03G,MXNGPI,N3HRG,VKMIN,VKMAX,
     +  IBTHR,IETHR)
        IF(NOKCP(3) .GT. 0)CALL PLTYVT('Crosstail = 29+11Kp',
     +  '[KV]',CPFKP(1,3),ITM03G,MXNGPI,N3HRG,VKMIN,VKMAX,
     +  IBTHR,IETHR)
      ENDIF
      IF(IGETHP .EQ. 1)THEN
        IF(NOHCP(1) .GT. 0)CALL PLTYVT('Crosstail = 16+1.3P',
     +  '[KV]',CPFHP(1,1),ITMEHP,MXNPWR,NVEHP,VKMIN,VKMAX,IBTHR,IETHR)
        IF(NOHCP(2).GT.0)CALL PLTYVT('Crosstail = max{25,-48+36Ln(P)}',
     +  '[KV]',CPFHP(1,2),ITMEHP,MXNPWR,NVEHP,VKMIN,VKMAX,IBTHR,IETHR)
      ENDIF
      IF(IGETAE .EQ. 1)THEN
        IF(NOACP(2) .GT. 0)CALL PLTYVT('Crosstail = 41+.11AE',
     +  '[KV]',CPFAE(1,2),ITMAE,MXNAE,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR)
        IF(NOACP(1).GT.0)CALL PLTYVT('P(AE)  CP = max{25,-48+36Ln(P)}',
     +  '[KV]',CPFAE(1,1),ITMAE,MXNAE,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR)
      ENDIF
      IF(IGETIMF .EQ. 1)THEN
        IF(NOICPF(1) .GT. 0)CALL PLTYVT('Crosstail = vBsin3 filtered',
     +  '[KV]',POTF(1,1),ITMIMF,MXNIMF,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR)
        IF(NOICPF(2) .GT. 0)CALL PLTYVT('Crosstail = vBQ1 filtered',
     +  '[KV]',POTF(1,2),ITMIMF,MXNIMF,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR)
        IF(NOICPF(3) .GT. 0)CALL PLTYVT('Crosstail = vBsin3 filtered gap
     +s filled with CP(AE,P,Kp)',
     +  '[KV]',POTF(1,3),ITMIMF,MXNIMF,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR)
      ENDIF
  980 CONTINUE
 
      GO TO 2
 
 1000 CONTINUE
 
C          Close GKS plotting unit
      CALL CLSGKS
      GO TO 9999
 
C          Error trap diagnostics
 9000 WRITE(MSGUN,'('' Bad time input.  Begin date: '',4I5,/,
     +              ''                    End date: '',4I5)')IBEGT,IENDT
      STOP
 9050 WRITE(MSGUN,'('' Beginning date later than ending date input'',/,
     +              ''         Begin date: '',4I5,/,
     +              ''           End date: '',4I5)')IBEGT,IENDT
      STOP
 9100 WRITE(MSGUN,'('' Dimension exceeded.  Increase MXNGPI to'',I8)')
     +                                                        NDAS
      STOP
 9200 WRITE(MSGUN,'('' Dimension exceeded.  Increase MXNBY to'',I8)')
     +                                                        NHRTOT
      STOP
 9250 WRITE(MSGUN,'('' Dimension exceeded.  Increase MXNAE to'',I8)')
     +                                                        NHRTOT
 9260 WRITE(MSGUN,'('' Dimension exceeded.  Increase MXNDST to'',I8)')
     +                                                        NHRTOT
      STOP
 9300 FORMAT(' Requested date',2I5,' is before first available')
 9400 FORMAT(' Requested date',2I5,' is after last available')
 9500 FORMAT(' No data available for date',2I5)
 9600 FORMAT(' Indicies dataset corrupted?  Data should be',
     +                                     ' continuous yet',/,
     +       '         somehow read past requested time',3I5)
 
C          Normal termination point
 9999 CONTINUE
      END
 
      SUBROUTINE LEXWRT (COMMENT, PNAME, IOF5, ARR, ITM, NAR, IX,LEXADD)
C
C  LEXWRT writes out the lexical read parameters from the GETNDCS job.
C
C   COMMENT = the comment about the parameter to follow
C   PNAME = the parameter name of the lexical read
C   IOF5 = the number of digits seen past the decimal point in the
C          format of F5.IOF5  (is generally 0, 1, or 2).
C   ITM(IX,NAR) = the mid-point time of the measurement (year,
C                 month/day, UThr/min, csec)
C   ARR = the array of the parameter in question
C   NAR = the number of points in the arrays ARR and ITM
C   IX = the first dimension of ITM (usually 4)
C   LEXADD = the number to be added to the day for the purposes of the
C            lexical write.  (If 0, the day of the lexical day is the
C            same as the day of the month.  Could also be arranged so
C            the lexical day starts at 1 or that the lexical day is
C            the day number of the year.)
C
      CHARACTER *(*) COMMENT,PNAME
      CHARACTER *80 ALINE
      DIMENSION ARR(NAR), ITM(IX,NAR)
      DIMENSION NDT(1000), NHR(1000), NMIN(1000), AREV(1000)
      DATA TMISS/-32766.9/
C
      IF (NAR .GT. 1000) THEN
       WRITE (6,"(1X,'LEXWRT:  Need to increase dimension to NAR=',I5)")
     |  NAR
       STOP
      ENDIF
      IF (IOF5 .LT. 0 .OR. IOF5 .GT. 2) THEN
       WRITE (6,"(1X,'LEXWRT:  Need to change IOF5 to 0,1,2.  =',I5)")
     |  IOF5
       STOP
      ENDIF
      NAREV = 0
      DO 100 I=1,NAR
         IF (ARR(I) .LT. TMISS) GO TO 100
           NAREV = NAREV + 1
           AREV(NAREV) = ARR(I)
           NDT(NAREV)=MOD(ITM(2,I),100) + LEXADD
C          NDT represents the UT date
           NHR(NAREV)=ITM(3,I)/100
C          NHR represents the UT hour
           NMIN(NAREV)=MOD(ITM(3,I),100)
C          NMIN represents the UT minute
  100 CONTINUE
      IF (NAREV .LE. 0) RETURN
      NC = LEN(COMMENT)
      WRITE(61,"(A)") COMMENT(1:NC)
      NP = LEN(PNAME)
      WRITE(61,"(1X,A)") PNAME(1:NP)
      NLIN = (NAREV+3) / 4
      NLAST = NAREV - (NLIN-1)*4
      NINTW = NAREV - NLAST
      NW =  NLAST * 16
      IF (IOF5 .EQ. 0) THEN
        IF (NLIN-1 .GT. 0) THEN
          DO 110 NL=1,NLIN-1
          I1 = 1 + (NL-1)*4
          I2 = I1 + 3
  110     WRITE(61,"(1X,4(I3,',',2(I2,','),F5.0,',') )")
     |     (NDT(I),NHR(I),NMIN(I),AREV(I),I=I1,I2)
        ENDIF
        WRITE (ALINE,"(1X,4(I3,',',2(I2,','),F5.0,',') )")
     |     (NDT(I),NHR(I),NMIN(I),AREV(I),I=1+NINTW,NAREV)
      ENDIF
      IF (IOF5 .EQ. 1) THEN
        IF (NLIN-1 .GT. 0) THEN
          DO 120 NL=1,NLIN-1
          I1 = 1 + (NL-1)*4
          I2 = I1 + 3
  120     WRITE(61,"(1X,4(I3,',',2(I2,','),F5.1,',') )")
     |     (NDT(I),NHR(I),NMIN(I),AREV(I),I=I1,I2)
        ENDIF
        WRITE (ALINE,"(1X,4(I3,',',2(I2,','),F5.1,',') )")
     |     (NDT(I),NHR(I),NMIN(I),AREV(I),I=1+NINTW,NAREV)
      ENDIF
      IF (IOF5 .EQ. 2) THEN
        IF (NLIN-1 .GT. 0) THEN
          DO 130 NL=1,NLIN-1
          I1 = 1 + (NL-1)*4
          I2 = I1 + 3
  130     WRITE(61,"(1X,4(I3,',',2(I2,','),F5.2,',') )")
     |     (NDT(I),NHR(I),NMIN(I),AREV(I),I=I1,I2)
        ENDIF
        WRITE (ALINE,"(1X,4(I3,',',2(I2,','),F5.2,',') )")
     |     (NDT(I),NHR(I),NMIN(I),AREV(I),I=1+NINTW,NAREV)
      ENDIF
      WRITE (61,"(A)") ALINE(1:NW)
      RETURN
      END
 
      SUBROUTINE FINDPOT (BX,BY,BZ,SWD,SWV,P,AE,RKP,POT1,POT2,POT3)
C  Find potential from the vBsin3(theta/2) formula from Reiff and
C   Luhmann, 1986
C
      DATA RMISS/-32767./, TMISS/-32766.9/
C
      IF (BX.LT.TMISS .OR. BY.LT.TMISS .OR. BZ.LT.TMISS .OR.
     | SWV.LT.TMISS) THEN
       POT1 = RMISS
       POT2 = RMISS
C  In the absence of IMF data, use formula for CP from AE, P, or Kp
C   (in that order of importance)
       POT3 = RMISS
       IF (RKP .GT. TMISS) POT3 = 29. + 11. * RKP
       IF (P .GT. TMISS) POT3 = AMAX1(25., -48.05 + 35.82 * ALOG(P) )
       IF (AE .GT. TMISS) POT3 = 41. + 0.11 * AE
       RETURN
      ENDIF
      B = SQRT(BX**2 + BY**2 + BZ**2)
      THETA = ACOS(BZ/B)
      BU = AMIN1(B,9.)
      POT1 = 12.7 + 0.041 * SWV * BU * ( SIN(THETA/2.) )**3
      CALL CPVBQ1 (BX,BY,BZ,SWD,SWV,POT2)
      POT3 = 12.7 + 0.041 * SWV * BU * ( SIN(THETA/2.) )**3
      RETURN
      END
 
      SUBROUTINE CPVBQ1 (BX,BY,BZ,SWD,SWV,POT)
C     This subroutine finds the crosstail potential from the Reiff
C     formula for Crooker Compression vBQ1, REIFF and LUHMAN, 1986.
 
      IF ((BX .EQ. -32767.) .OR. (BY .EQ. -32767.) .OR.
     +    (BZ .EQ. -32767.) .OR. (SWD .EQ. -32767.) .OR.
     +    (SWV .EQ. -32767.)) THEN
           POT = -32767.
           RETURN
      ENDIF
 
      B=SQRT(BX**2+BY**2+BZ**2)
      BM=AMIN1(B,7.)
      BT=SQRT(BY**2+BZ**2)
      BE=SQRT(2.*0.002098*SWD*SWV**2+BT**2)
      BS=2.*SQRT(BT*BE)
      BS=AMIN1(BS,BE)
      X=0.
      C8=BS/BE
      C9=BZ/BT
      IF (C8 .GE. C9) X=SWV*BM*C8*(C8-C9)/SQRT(1.+C8**2-2.*C8*C9)
      POT=25.+0.062*X
      RETURN
      END
 
      SUBROUTINE FILPOT(POT,POTF,NUM,IFLGAP,RPOT)
C  This subroutine filters the potential across NUM consecutive hours
C   and stores the filtered potentials in POTF
C  IFLGAP = 0,1  If do,not have filled gaps in IMF POT (If = 1, need to
C                examine RPOT as well as POT)
C  All of the following variables should be dimensioned to MXNBY
C   in the calling routine
      DIMENSION POT(*),POTF(*),RPOT(*)
      DATA TMISS/-32766.9/
 
      POTF(1)=POT(1)
      IF (NUM .EQ. 1) RETURN
      POTF(2)=POT(2)
      IF (POT(2) .GT. TMISS  .AND. POT(1) .GT. TMISS)
     |   POTF(2) = POT(2)*0.57 + POT(1)*0.43
      IF (NUM .EQ. 2) GO TO 101
      DO 100 I=3,NUM
           IF ((POT(I-2) .GT. TMISS) .AND. (POT(I-1) .GT. TMISS)
     +         .AND. (POT(I) .GT. TMISS)) THEN
                    POTF(I)=POT(I)*0.57+POT(I-1)*0.29+POT(I-2)*0.14
           ELSE IF ((POT(I-2) .GT. TMISS) .AND. (POT(I) .GT. TMISS))
     +               THEN
                    POTF(I)=POT(I)*0.57+POT(I-2)*0.43
           ELSE IF ((POT(I-1) .GT. TMISS) .AND. (POT(I) .GT. TMISS))
     +               THEN
                    POTF(I)=POT(I)*0.57+POT(I-1)*0.43
           ELSE
                    POTF(I)=POT(I)
           ENDIF
  100 CONTINUE
  101 IF (IFLGAP .NE. 1) RETURN
C  Have filled in gaps --- do not blend over original data when it
C   comes out of a gap
      IF (RPOT(2) .LT. TMISS) GO TO 150
      POTF(2)=POT(2)
      IF (RPOT(2) .GT. TMISS  .AND. RPOT(1) .GT. TMISS)
     |   POTF(2) = POT(2)*0.57 + POT(1)*0.43
  150 CONTINUE
      IF (NUM .EQ. 2) RETURN
      DO 200 I=3,NUM
        IF (RPOT(I) .LT. TMISS) GO TO 200
C
           IF ((RPOT(I-2) .GT. TMISS) .AND. (RPOT(I-1) .GT. TMISS)
     +         .AND. (RPOT(I) .GT. TMISS)) THEN
                    POTF(I)=POT(I)*0.57+POT(I-1)*0.29+POT(I-2)*0.14
           ELSE IF ((RPOT(I-2) .GT. TMISS) .AND. (RPOT(I) .GT. TMISS))
     +               THEN
                    POTF(I)=POT(I)*0.57+POT(I-2)*0.43
           ELSE IF ((RPOT(I-1) .GT. TMISS) .AND. (RPOT(I) .GT. TMISS))
     +               THEN
                    POTF(I)=POT(I)*0.57+POT(I-1)*0.43
           ELSE
                    POTF(I)=POT(I)
           ENDIF
  200 CONTINUE
      RETURN
      END
 
      FUNCTION TERP (X,Y,XI,N,XMIN,IFLGSP)
C     ****
C     ****     PERFORMS LINEAR INTERPOLATION IN TABLE OF N (X,Y) POINTS
C     ****     X(N) IS ARRAY OF MONOTONICALLY INCREASING ABSCISSAE
C     ****     Y(N) IS ARRAY OF CORRESPONDING ORDINATES
C     ****     XI IS ABSCISSA AT WHICH INTERPOLATION IS REQUIRED
C     ****     N IS NUMBER OF POINTS IN TABLE
C     ****     INTERPOLATED VALUE IS RETURNED IN TERP
C     ****
C     ****
      DIMENSION X(1),Y(1)
C     ****
C     ****     SEARCH FOR INTERVAL IN WHICH XI FALLS
C     ****
C  XMIN =        The minimum spacing in X that XI is required to be
C                within bounds
C  IFLGSP = 0,1  Is flag indicating whether or not there is a large gap
C                (greater than XMIN) near the time XI
C
      IFLGSP = 0
      IF(N.EQ.1)THEN
           TERP=Y(1)
      ELSE
           NN=1
           DO 1 I=1,N-1
           IF(XI.GT.X(I))NN=I
    1      CONTINUE
C          ****
C          ****     PERFORM LINEAR INTERPOLATION
C          ****
           TERP=((X(NN+1)-XI)*Y(NN)+(XI-X(NN))*Y(NN+1))/(X(NN+1)-X(NN))
      IF (XI .LT. X(1)) TERP = Y(1)
      IF (XI .GT. X(N)) TERP = Y(N)
      DELTNP = ABS(X(NN+1)-XI)
      DELTNN = ABS(XI-X(NN))
      IF (DELTNP .GT. XMIN .AND. DELTNN .GT. XMIN) IFLGSP = 1
      ENDIF
      RETURN
      END
 
      SUBROUTINE RDBLK6(IWRK,LREC,IBLK,IUNPK,ISTAT)
C          This subroutine reads NCAR I.S. Data Base Tape formatted data
C          sets.  See "The NCAR I.S. Data Base Tape Format" for a format
C          description.  When called RDBLK6 returns a status flag and the
C          next logical record in the dataset; i.e. sequential access.
C          Alternate entries, OPNBLK, CLSBLK and REWBLK provide a mechanism
C          for establishing connection of a fortran unit to a file name,
C          severing the connection, and rewinding the dataset.  It is
C          necessary to call OPNBLK prior to calling RDBLK6.  Closing and
C          rewinding are at the users discretion.  The status flag provides
C          an indication when there is nothing left to read (or sooner,
C          depending on how the NEOFS termination condition is specified).
C          INPUTS:
C            IWRK  = Work array provided by the user.  See entry OPNBLK for
C                    dimension.
C            LREC  = Although this will contain a logical record on return
C                    the user must provide the array dimensioned to MXVALS.
C                    MXVALS = 10000 is usually enough, but RDBLK6 will
C                    complain if more space is needed.
C            IBLK  = User provided work array, dimensioned to MXWDS.
C                    MXWDS = MXVALS*16/LMWD, where LMWD = computer word size
C                    in bits.  (IBLK contains a tape block).
C            IUNPK = User provided work array, dimensioned to MXVALS.  (IUNPK
C                    contains an unpacked version of IBLK, with the length
C                    prefix and checksum suffix removed).
C          RETURNS:
C            LREC  = Logical record.  See the "NCAR I.S. Database Tape Format"
C                    for a complete description of the record layout.  The
C                    prologue (the first 16 values) will always contain
C                    integers.  If the prologue indicates that the logical
C                    record is type "data", the rest of LREC is also integers.
C                    If, however, LREC contains a "header" or "catalogue"
C                    record, elements starting at position 41 each contain
C                    two ASCII characters in the right most 16 bits.  (Use
C                    subroutine HDRPRT to print the contents of these LRECS).
C            ISTAT = Status of last block read.  This value should be
C                    interpreted differently, depending on how NEOFS was
C                    specified (in the OPNBLK call).
C                  = 0 => Normal.  LREC contains the next logical record
C                         from fortran unit IBUN (IBUN is defined in OPNBLK).
C                  = 1 => EOF (end of file) encountered.
C                         If NEOFS>0, then NEOFS files have been read from
C                         IBUN (and presumably the user wishes to terminate).
C                         If NEOFS=0, then another call to RDBLK6 is in order
C                         as there may be more data to read in IBUN.
C                  = 2 => Parity error.
C                         If NEOFS>0, the user will not see this value as the
C                         RDBLK6 will automatically print a warning and
C                         proceed to read the next block block.
C                         If NEOFS=0, then the user must call RDBLK6 again to
C                         try to read the next block.
C                  = 3 => EOD (end of data) encountered.  Another call to
C                         RDBLK6 will result in an error (unless IBUN is
C                         repositioned at the beginning via entry REWBLK).
C
C          Alternate entrys:
C             OPNBLK(IBUN,NAME,NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,MXWDS,IWRK,LWRK)
C                    OPNBLK must be called prior to calling RDBLK6 in order
C                    to establish connection between IBUN and NAME as well
C                    as some other houskeeping tasks.  All arguments are
C                    inputs:
C                    IBUN  = Fortran unit number to use for reading file NAME.
C                    NAME  = Character string containing the name of the
C                            input dataset.
C                    NEOFS = Number of EOFs to encounter, before RDBLK6 sends
C                            a completion status message.  There are two
C                            ways to use this termination warning:
C                            (1) NEOFS = 0.  Control is returned to the user
C                                for all conditions that ISTAT is not 0 and
C                                LREC is not updated.  This gives the user
C                                the most control; e.g., this is useful when
C                                desiring to examine the last LREC of each
C                                file in NAME or it is useful when desiring to
C                                count files in NAME.
C                            (2) NEOFS > 0.  When an EOF is encountered, the
C                                first LREC of the next file returned.  This
C                                This continues until NEOFS files have been
C                                read; at which point LREC is not updated and
C                                ISTAT=1.  If the EOD is encountered before
C                                NEOFS files are read, ISTAT = 3.  Note that
C                                attempting to read more than is there will
C                                result in an error, so the user should always
C                                check ISTAT after a RDBLK6 call.
C                    IOCKSM= Checksum option; on=1, off=0.  If on, block
C                            checksums will be verified.
C                    MSGUN = Fortran unit to use for any diagnostic messages.
C                    LMWD  = Word size in bits; e.g., 32.
C                    MXVALS= The dimension of arrays, LREC and IUNPK.
C                    MXWDS = The dimension of array IBLK.
C                    IWRK  = Work array provided by the user.
C                    LWRK  = Length of IWRK array provided by the user.  The
C                            dimension of IWRK is a function of the machine
C                            word size:
C                              dimension = 37 + 2 * 512 * 64 / LMWD
C
C             CLSBLK(IWRK)
C                    CLSBLK may be called to break the connection between
C                    IBUN and NAME.  This is handy when desiring to use
C                    IBUN for multiple purposes in the same program.
C
C             REWBLK(IWRK)
C                    REWBLK positions IBUN at the beginning of NAME, thus
C                    providing a way to rewind NAME.
C
C          HISTORY:
C          RDBLK6 is similar to RDBLK5 except that most parameters
C          passed to RDBLK6 are contained in an array called IWRK.
C          This makes RDBLK6 recallable using different Fortran units.
C          The same characteristics for RDBLK5 are true for RDBLK6:
C          There are more call arguments which allow LMWD and all
C          array sizes to set externally.  Rather than using system
C          specific calls to print diagnostics in the log file during
C          batch execution these prints have been moved to unit MSGUN.
C          The alternate entries were added (Jan 90) to provide a more
C          portable interface with UNIX systems in mind.  This version
C          was installed on a 32-bit UNIX machine using my version of rdtape
C          called cbfrd.  Next it was installed back on a COS Cray, and
C          most recently, on UNICOS.  The biggest UNICOS change, was to
C          add logic to set status to EOD encountered based apon reading
C          sequential EOF's.
C
C          The following list of the parameters describes the IWRK array
C          elements used by RDBLK6.  The variable names following the IWRK
C          element numbers are in some cases vestiges of RDBLK5; in those
C          instances a description is provided (as it is not available
C          elsewhere):
C            IWRK( 1) = IBUN
C            IWRK( 2) = NEOFS
C            IWRK( 3) = MSGUN
C            IWRK( 4) = LMWD
C            IWRK( 5) = MXVALS
C            IWRK( 6) = MXWDS
C            IWRK( 7) = FIRST  = Switch (on=1, off=0) to initialize pointers
C                                and counters in RDBLK6.
C            IWRK( 8) = IOCKSM
C            IWRK( 9) = POINTR = Pointer into IUNPK, indicating the begining
C                                of the next lrec.
C            IWRK(10) = NVALS  = The number of assigned values in IUNPK.
C            IWRK(11) = LPREC  = The number of assigned values in IBLK (the
C                                first 16 bit value in IBLK).
C            IWRK(12) = NPREC  = The number of blocks read.
C            IWRK(13) = NWARN  = Number of warning messages issued about
C                                blocks containing zero filled 16-bit values
C                                at the end of the block (before the checksum)
C            IWRK(14) = NLMSGS = The number of blocks which had non-zero
C                                16-bit values after the checksum.
C            IWRK(15) = NPERR  = The number of blocks with parity errors.
C            IWRK(16) = EOFCNT = The number of EOFs encountered.
C            IWRK(17) = NLRECS = The number of logical records read.
C            IWRK(18) to IWRK(x) = Yet another work array that is used by
C                                subroutine  CBFRD (which does the actual
C                                read).  The additional space needed by CBFRD
C                                depends on LMWD and may be calculated as
C                                   dim = 20 + 2 * 512 * 64 / LMWD.
C                                So IWRK must be dimensioned to 17 + dim;
C                                e.g., if LMWD=32, x=2085 or when
C                                         LMWD=16, x=4133.
C
C          Declarations for formal arguments (call variables)
      DIMENSION  LREC(*) , IUNPK(*) , IWRK(*) , IBLK(*)
      CHARACTER*(*) NAME
C
C          Local declarations:
      PARAMETER (MXNPER=10)
      DIMENSION IDUMP(10)
      INTEGER  REMANS,CKSUM,CKSMRD
      LOGICAL NEDMOR
C
C          UNICOS version (curently UNICOS does not return EOD to IOWAIT
C          so this is a patch to keep count of the number of sequential
C          EOF's (no intervening data) for each possible fortran unit:
      SAVE      NSEOF                           ! <- Cray UNICOS version
      DIMENSION NSEOF(100)                      ! <- Cray UNICOS version
C
C          Skip alternate entry points code
      GO TO 5
C
C          Initialize arrays, pointers, etc. for reading:
      ENTRY OPNBLK(IBUN,NAME,NEOFS,IOCKSM,MSGUN,LMWD,MXVALS,MXWDS,
     +             IWRK,LWRK)
      IWRK( 1) = IBUN
      IWRK( 2) = NEOFS
      IWRK( 3) = MSGUN
      IWRK( 4) = LMWD
      IWRK( 5) = MXVALS
      IWRK( 6) = MXWDS
      IWRK( 7) = 1
      IWRK( 8) = IOCKSM
      ICK = MXVALS*16/LMWD
      IF(MXWDS .LT. ICK)GO TO 8100
CUNIX ICK = 37 + 65536 / LMWD
      ICK = 17                                  ! <- Cray UNICOS version
      IF(LWRK  .LT. ICK)GO TO 8200
CUNIX LWK = LWRK - 17
CUNIX CALL CBFOPN(IBUN,NAME,MSGUN,LMWD,1,IWRK(18),LWK)
      OPEN (IBUN,FILE=NAME,FORM='UNFORMATTED')  ! <- Cray UNICOS version
      RETURN
C
C          Close IBUN
      ENTRY  CLSBLK(IWRK)
CUNIX CALL CBFCLS(IWRK(1),IWRK(18))
      CLOSE (IWRK(1))                           ! <- Cray UNICOS version
      RETURN
C
C          Rewind IBUN
      ENTRY REWBLK(IWRK)
      IWRK(7) = 1
CUNIX CALL CBFREW(IWRK(1),IWRK(18))
      REWIND IWRK(1)                            ! <- Cray UNICOS version
      RETURN
C
C          If the FIRST switch is on, initialize and read a block
    5 NEDMOR = .FALSE.
      IF(IWRK(7) .EQ. 1)THEN
        IWRK(7) = 0
        DO 6 I=10,17
    6   IWRK(I) = 0
        NSEOF(IWRK(1)) = 0                      ! <- Cray UNICOS version
        GO TO 100
      ENDIF
C
C          Determine the number of 16-bit values left in IUNPK:
      REMANS = IWRK(10) - IWRK(9) + 1
      IF(REMANS .LT. 0)GO TO 8900
      IF(REMANS .EQ. 0)GO TO  100
      IF(REMANS .LE. 3)THEN
C          If the remaining values are 0's go get another block; else
C          assume they are the beginning of an lrec spanning a block
        DO 7 I=IWRK(9),IWRK(9)+REMANS-1
        IF(IUNPK(I) .NE. 0)GO TO 8
   7    CONTINUE
        IWRK(13) = IWRK(13) + 1
        GO TO 100
   8    CONTINUE
      ENDIF
C          If only a partial LREC remains in IUNPK then fill LREC
C          with the REMANS, determine how many more values are needed,
C          turn on the 'I-need more' flag and go get another block.
      IF(REMANS .LT. IUNPK(IWRK(9)))THEN
        DO 10 I=1,REMANS
   10   LREC(I) = IUNPK(IWRK(9)+I-1)
        NNEDED  = IUNPK(IWRK(9))-REMANS
        NEDMOR  = .TRUE.
        GO TO 100
      ENDIF
C          At least one full LREC exists in IUNPK.  Therefore, assign
C          LREC from IUNPK, move POINTR to the beginning of the next
C          possible LREC and return.  First, check for a reasonable
C          logical record length before using it.
   20 IF(IUNPK(IWRK(9)) .LE. 0)GO TO 9600
      DO 30 I=1,IUNPK(IWRK(9))
   30 LREC(I) = IUNPK(IWRK(9)+I-1)
      IWRK(9)  = IWRK(9) + LREC(1)
      IWRK(17) = IWRK(17) + 1
      RETURN
C
C          Initialize block buffer and pointer
  100 DO 110,I=1,IWRK(6)
  110 IBLK(I) = 0
      DO 120,I=1,IWRK(5)
  120 IUNPK(I) = 0
      IWRK(9) = 1
C          Read a block (physical record) containing a variable number
C          of logical records:
C          UNIX VERSION:  CBFRD replaces calls to RDTAPE and IOWAIT
CUNIX CALL CBFRD(IWRK(1),IBLK,IWRK(6),LBLK,NUB,ISTAT,IWRK(18))
      CALL RDTAPE (IWRK(1),1,0,IBLK,IWRK(6)) ! <- Cray UNICOS version
      CALL IOWAIT (IWRK(1),ISTAT,LBLK)       ! <- Cray UNICOS version
      IF(ISTAT .EQ. 1 .AND.
     +   NSEOF(IWRK(1)) .GT. 0) ISTAT = 3    ! <- Cray UNICOS version
C          RDTAPE reads an entire block from unit IWRK(1) into array IBLK
C          which has been dimensioned to IWRK(6); there is no formatting
C          and variable length blocks are allowed.
C          IOWAIT obtains the status of the previous RDTAPE record read
C          request.
C          The status of unit IBUN is returned as ISTAT, an integer flag,
C          and LBLK, the number of words read).
C            ISTAT = 0 => okay read of LBLK words into IBLK
C                  = 1 => end of file encountered, nothing new in IBLK
C                  = 2 => parity error encountered, nothing new in IBLK
C                  = 3 => end of data encountered, nothing new in IBLK
C
C          Normal and error condition branches from ISTAT values:
      IF(ISTAT .EQ. 1)THEN
        IWRK(16) = IWRK(16) + 1
        IF(NEDMOR)GO TO 9100
        NSEOF(IWRK(1)) = NSEOF(IWRK(1)) + 1     ! <- Cray UNICOS version
        IWRK(9) = IWRK(10) + 1
        IF(IWRK(2)  .LT. 1)RETURN
C>>>> (patch to prt prol of last lrec in each file)
C     CALL LRPRT(IWRK(3),IWRK(3),IWRK(17),'PROL',LREC)
C>>>>
        IF(IWRK(16) .LT. IWRK(2))GO TO 100
        WRITE(IWRK(3),'('' RDBLK6:  Done reading from unit'',I3,
     +         ''; # EOFs read = '',I5)')IWRK(1),IWRK(16)
        RETURN
      ENDIF
      IF(ISTAT .EQ. 2)THEN
        IWRK(9)  = IWRK(10) + 1
        IWRK(12) = IWRK(12) + 1
        NSEOF(IWRK(1)) = 0                      ! <- Cray UNICOS version
        IF(IWRK(2) .LT. 1)RETURN
        GO TO 9700
      ENDIF
      IF(ISTAT .EQ. 3)THEN
        IF(NEDMOR)GO TO 9200
        WRITE(IWRK(3),'('' RDBLK6:  EOD encountered on unit'',I3,
     +            '' after reading'',I5,'' files.'')')IWRK(1),IWRK(16)
        RETURN
      ENDIF
      IF(ISTAT .LT. 0 .OR. ISTAT .GT. 3)STOP
      NSEOF(IWRK(1)) = 0                        ! <- Cray UNICOS version
      IWRK(12) = IWRK(12) + 1
C
C          Good block read; now strip off 16-bit values into array IUNPK.
C          First, get LPREC and CKSMRD (CHECKSUM) and calculate NVALS
C          based upon LPREC and test LBLK VS. LPREC.  Also if CKSMOP
C          is .TRUE. calculate and compare checksums.
C      IWRK(11) LPREC = # 16-bit values read as the first value in the
C                       block.
C      IWRK(10) NVALS = # 16-bit values in the block minus the length
C                       prefix and the checksum suffix.
C               LBLK  = # computer words found in the block.
C         (UNIX VERSION:  NUB = number of unused bits in IBLK(LBLK)
C              CKSMRD = The 16-bit XOR checksum for the block.
C          IUNPK does not contain the prefix, IWRK(11), or suffix, CKSMRD
  500 IF(LBLK .GE. IWRK(6))GO TO 9000
C
C          Subroutine GBYTE moves bits (16 starting at offset 0) from
C          array IBLK into word LPREC.  The 16 bits are right justified
C          in LPREC that it is a normal integer.  Hence, GBYTE
C          unpacks an integer value from an array.
      CALL GBYTE(IBLK,IWRK(11),0,16)
C
C          It is likely that there is a partial word remaining at the
C          end of a block, since the logical records use a word length
C          (16-bit) different than the computer being used.  The amount
C          of extra bits left over is a function of the word size of the
C          computer creating the block and the word size of the computer
C          reading it (the only known word size is the read).  It is good
C          practice to see if end of a block has a truncated LREC or is
C          merely zero filled.  If non-zero values are found, then print
C          a message, otherwise continue processing.  IDIFF is the max
C          number of 16-bit values left over at the end of the block.
C          UNIX VERSION relies on NUB
      IDIFF = (LBLK*IWRK(4) - IWRK(11)*16)/16
CUNIX IDIFF = (LBLK*IWRK(4) - NUB - IWRK(11)*16)/16
C>>>>      Patch for reading padded blocks (fixed length zero fill)
      IMXOK = 3000
C>>>> IMXOK = MAX0 (4,2*IWRK(4)/16)
      IF (IDIFF .LT. 0 .OR. IDIFF .GT. IMXOK) GO TO 9300
      IF (IDIFF .GT. 0) THEN
        NNE0 = 0
        IOFF = 16*IWRK(11)
        DO 510 I=1,IDIFF
        ITEMP = 0
        CALL GBYTE(IBLK,ITEMP,IOFF,16)
        IF (ITEMP .NE. 0) NNE0 = NNE0+1
  510   IOFF = IOFF + 16
        IF (NNE0 .GT. 0) THEN
          IWRK(14) = IWRK(14) + 1
          IF (IWRK(14) .LT. 10)
     +    WRITE(IWRK(3),'('' RDBLK6:  Warning:  Block'',I5,'' had'',I5,
     +        '' 16-bit values after the cksum;'',I3,'' are not 0'')')
     +         IWRK(12),IDIFF,NNE0
        ENDIF
      ENDIF
      IWRK(10) = IWRK(11) - 2
C          Subroutine GBYTES is the plural of GBYTE; it unpacks IWRK(10)
C          values from IBLK into IUNPK, 16 bits at a time starting with
C          an initial offset into IBLK of 16 with an iteration bit skip of 0.
C          This GBYTES call does the unpacking of all lrecs in the prec
      CALL GBYTES(IBLK,IUNPK,16,16,0,IWRK(10))
      IF(IWRK(8) .EQ. 1)THEN
C          This call gets the 16-bit checksum at the end of the IBLK
        CALL GBYTE(IBLK,CKSMRD,16*(IWRK(11)-1),16)
        CKSUM=IWRK(11)
        DO 550 I=1,IWRK(10)
  550   CKSUM=CKSUM .XOR. IUNPK(I)
        IF(CKSUM .NE. CKSMRD)GO TO 9400
      ENDIF
C
C          Convert any negatives (2's complement notation) examples:
C             ABSOLUTE VALUE    MAPS TO
C            ----------------   --------
C                 0 = 0       ->  0
C                 1 = 2**0    ->  1
C             32767 = 2**15-1 ->  32767  max positive value
C             32768 = 2**15   -> -32768  min negative value
C             32769 = 2**15+1 -> -32767
C             65535 = 2**16-1 -> -1      Max absolute val before convert
      DO 600 N=1,IWRK(10)
      IF(IUNPK(N) .GE. 32768)IUNPK(N) = IUNPK(N) - 65536
  600 CONTINUE
C
C          If more values from this IBLK are needed to complete the LREC
C          then assign these values from IUNPK to LREC, turn off the
C          'NEDMOR' switch, and set POINTR to the beginning of the next
C          LREC in IUNPK; if NEDMOR is false then go back to the normal
C          LREC assignment logic (the DO 20 loop).
      IF(NEDMOR)THEN
        DO 700 I=1,NNEDED
  700   LREC(REMANS+I) = IUNPK(I)
        IWRK(9) = NNEDED + 1
        NEDMOR = .FALSE.
        IWRK(17) = IWRK(17) + 1
        RETURN
      ENDIF
      GO TO 20
C
C          Error trap messages:
 8100 WRITE(MSGUN,'('' RDBLK6:  MXVALS and MXWDS do not agree.  MXWDS'',
     +                                  '' is too small.'',/,
     +  ''          MXWDS >= MXVALS*16/LMWD.  Yet MXWDS,MXVALS='',2I8)')
     +                                            MXWDS,MXVALS
 8200 WRITE(MSGUN,'('' RDBLK6:  Insufficient work array size'',
     +              '' specified.  LWRK input ='',I8,/,
     +          ''          but LWRK must be at least ='',I8)')LWRK,ICK
      STOP
 8900 WRITE(IWRK(3),8910) REMANS
 8910 FORMAT(' RDBLK6:  Logic error: REMANS are < 0.  REMANS = ',I10)
      STOP
 9000 WRITE(IWRK(3),9010)IWRK(12),IWRK(16)+1
 9010 FORMAT(' RDBLK6:  Block/file',I6,'/',I4,' might be larger than ',
     +        'presently dimensioned; increase IWRK(5)')
      STOP
 9100 WRITE(IWRK(3),9110)IWRK(12),IWRK(16)+1,(LREC(I),I=1,LREC(1))
 9110 FORMAT(' RDBLK6: Found only a partial lrec at the end of block',
     +                  I5,' of file',I4,' when an EOF was read.',/,
     +       ' Logical records are not allowed to span EOFs.  This',
     +              ' lrec contained:',/,(10I8))
      STOP
 9200 WRITE(IWRK(3),9210)IWRK(12),IWRK(16)+1,(LREC(I),I=1,LREC(1))
 9210 FORMAT(' RDBLK6:  Found only a partial lrec at the end of block',
     +                 I5,' of file',I4,' when the EOD was read.',/,
     +       ' This lrec contained:',/,(10I8))
      STOP
 9300 N16S = LBLK*IWRK(4)/16
      WRITE(IWRK(3),9310)IWRK(12),IWRK(16)+1,IWRK(4),LBLK,IWRK(11)
 9310 FORMAT(' RDBLK6:  Bogus block length read in block',I8,' of file',
     +I4,/,'       Number of computer (',I2,'-bit) words returned =',
     +I8,/,'    Number of 16-bit words specified in the block =',I8,/,
     +  ' The block contents as 16-bit integers:',44X,'Last Index')
      IOFF = 0
      DO 9330 I=1,N16S,10
      DO 9320 J=1,10
 9320 IDUMP(J) = 0
      N2G = MIN0(10,N16S-I-1)
      CALL GBYTES(IBLK,IDUMP,IOFF,16,0,N2G)
      WRITE(IWRK(3),'(10I8,''     >'',I5)') IDUMP , I+N2G-1
 9330 IOFF = IOFF + N2G*16
      STOP
 9400 IB=IWRK(10)-39
      IF(IB .LT. 0)IB=1
      WRITE(IWRK(3),9410)IWRK(12),IWRK(16)+1,CKSMRD,CKSUM,
     +                 (IUNPK(I),I=IB,IWRK(10))
 9410 FORMAT(' RDBLK6:  Check-sum read does not match that computed',
     +                  ' for block',I5,' of file',I4,/,
     +          '                read = ',I8,/,
     +          '            computed = ',I8,/,
     +          ' The end of the block contained:',/,(10I8))
      WRITE(IWRK(3),9420)IWRK(11),(IUNPK(I),I=1,IWRK(10))
 9420 FORMAT(' RDBLK6:  First 16-bit value in the block = ',I8,/,
     +       ' The rest of the block after unpacking:',/,(10I8))
      STOP
 9600 WRITE(IWRK(3),9610)IWRK(12),IWRK(16)+1,REMANS,IUNPK(IWRK(9)),
     +IWRK(9),(IUNPK(I),I=IWRK(9)+1,REMANS)
 9610 FORMAT(' RDBLK6:  Bogus logical record length encountered in',
     +                  ' block ',I5,' of file ',I4,/,
     +          '       Number of values left in block = ',I8,/,
     +          '                          Value found = ',I8,/,
     +          '                    at block position = ',I8,/,
     +          '       Block contents after this value:',/,(10I8))
      STOP
 9700 IWRK(15)=IWRK(15) + 1
      WRITE(IWRK(3),9710)IWRK(12),IWRK(16)+1
 9710 FORMAT(' RDBLK6:  Parity error encountered in block',I5,
     +                                         ' of file',I4)
      IF(IWRK(15) .LT. MXNPER)GO TO 100
      STOP
      END
      SUBROUTINE PARSB(STR,LB,LE)
C          This routine deterimines the location of the first, LB, and
C          last, LE, non-blank characters in the string STR.  If STR is
C          entirely blank, LE = 0 (and LB = 1).
      CHARACTER*(*) STR
 
      NCH = LEN(STR)
      LB = 1
      DO 10 I=1,NCH
      IF (STR(I:I) .NE. ' ') GO TO 20
      LB = LB + 1
   10 CONTINUE
   20 LE = NCH
      DO 30 I=NCH,1,-1
      IF (STR(I:I) .NE. ' ') GO TO 40
      LE = LE - 1
   30 CONTINUE
   40 CONTINUE
      IF (LE .EQ. 0) LB = 1
 
      RETURN
      END
      SUBROUTINE TIMDIF(BEGT,ENDT,IDIF)
C          This computes the time difference in units of centi-secs
C          (i.e. 100th secs) as ENDT - BEGT, where ___T are integer dim 4
C          with elements yyyy, mmdd, hhmm and centiseconds respectively.
 
      INTEGER BEGT(4),ENDT(4),BEGMO,BEGDA,BEGHR,BEGMI,
     +                        ENDMO,ENDDA,ENDHR,ENDMI
      DOUBLE PRECISION IDIF
 
      BEGMO = BEGT(2)/100
      BEGDA = BEGT(2)     - BEGMO*100
      BEGHR = BEGT(3)/100
      BEGMI = BEGT(3)     - BEGHR*100
 
      ENDMO = ENDT(2)/100
      ENDDA = ENDT(2)     - ENDMO*100
      ENDHR = ENDT(3)/100
      ENDMI = ENDT(3)     - ENDHR*100
 
C          Calculate the number centi-seconds difference:
      IDIF = DBLE( (ENDT(4)-BEGT(4)) )
 
C          Include the number minutes difference in units of csecs
      IDIF = IDIF + DBLE( (ENDMI-BEGMI) ) * 6000D0
 
C          Include hours difference
      IDIF = IDIF + DBLE( (ENDHR-BEGHR) ) * 360000D0
 
C          Calculate days and years difference:
      CALL HOURS(BEGT(1),BEGMO,BEGDA,IBEGH)
      CALL HOURS(ENDT(1),ENDMO,ENDDA,IENDH)
C          Subroutine HOURS computes the number of hours since 1920;
C          e.g. 0100 1 Jan 1921 entered would return IBEGH = 1.  HOURS
C          is part of the NCAR library routines.
      IDIF = IDIF + DBLE( (IENDH-IBEGH) ) * 360000D0
 
      RETURN
      END
      SUBROUTINE HOURS (KYR,KMO,KDAY,NHRS)
C
C DIMENSION OF           NONE
C ARGUMENTS
C
C LATEST REVISION        JANUARY 1985
C
C PURPOSE                GIVEN THE YEAR, MONTH AND DAY OF MONTH, THIS
C                        ROUTINE COMPUTES THE NUMBER OF HOURS UP TO
C                        THIS DAY AFTER 1920.
C
C                        BY MAKING TWO CALLS, THIS ROUTINE CAN BE USED
C                        TO FIND THE NUMBER OF HOURS BETWEEN TWO DATES.
C
C                        USE PACKAGE DATE FOR THE REVERSE PROCEDURE.
C
C USAGE                  CALL HOURS (KYR,KMO,KDAY,NHRS)
C
C ARGUMENTS
C
C ON INPUT               KYR
C                          YEAR (4 DIGITS, EG. 1962).
C
C                        KMO
C                          MONTH.
C
C                        KDAY
C                          DAY OF THE MONTH.
C
C ON OUTPUT              NHRS
C                          NUMBER OF HOURS UP TO THIS DAY AFTER 1920.
C
C SPECIAL CONDITIONS     .  WORKS FOR THE PERIOD 1921 TO 2100.
C                        .  FORTRAN STANDARD DATA STATEMENT LOADS THE
C                           ARRAY MO.
C
C I/O                    NONE
C
C PRECISION              SINGLE
C
C REQUIRED LIBRARY       Q8QST4, WHICH IS LOADED BY DEFAULT ON NCAR'S
C ROUTINES               CRAY MACHINES.
C
C LANGUAGE               FORTRAN
C
C HISTORY                WRITTEN BY ROY JENNE AT NCAR IN JANUARY 1966.
C
C PORTABILITY            FORTRAN 66
C ********************************************************************
      DIMENSION       MO(11)
      DATA MO(1),MO(2),MO(3),MO(4),MO(5),MO(6),MO(7),MO(8),MO(9),MO(10),
     *MO(11) /31,59,90,120,151,181,212,243,273,304,334/
      KYRF = KYR-1921
C
C DAYS UP TO THIS YR. MUST HAVE FIXED PT TRUNCATION OF KYRF/4
C
      LDAY = KYRF*365+KYRF/4
      IF (KMO .EQ. 1) GO TO 102
      IF (KMO .EQ. 2) GO TO 101
      IF (((KYR/4)*4-KYR) .NE. 0) GO TO 101
C
C COUNT DAY OF 29 FEB THIS YEAR
C
      LDAY = LDAY+1
  101 LDAY = LDAY+MO(KMO-1)
C
C NOW GET HOURS UP TO THIS DAY OF THIS MONTH
C
  102 NHRS = 24*(LDAY+KDAY-1)
      RETURN
C
C REVISION HISTORY---
C
C JANUARY 1978     DELETED REFERENCES TO THE  *COSY  CARDS AND
C                  ADDED REVISION HISTORY
C-----------------------------------------------------------------------
      END
      SUBROUTINE NEWTIM(IOLDT,NEWT,IDIF,MSGUN)
C          This subroutine creates a new time which is IDIF
C          centiseconds different from the old time.
C          IS date time units format; i.e. IOLDT & NEWT have 4 elements:
C          yyyy, mmdd, hhmm, and centiseconds.  IDIF is the number of
C          centiseconds difference.
C          NOTES:
C          - Only absolute idif's less than 365 days are allowed caurrently.
C          - Times >= 24__ are always converted to 00__;  i.e. the
C            largest hhmm returned will be 2359.
C          - MSGUN is the fortran unit number to print any diagnostic
C            messages.
C          - IDIF is double precision to accomodate short word (32 bit)
C            machines.
C          Formal argument declarations:
      DOUBLE PRECISION IDIF
      DIMENSION IOLDT(4),NEWT(4)
 
C          Local declarations
      DOUBLE PRECISION NDIFDA , NDIFHR , NDIFMI
      DIMENSION LMSUM(12),LYS(12),NLS(12),LENM(12)
      SAVE LYS,NLS,LENM
      DATA NLS / 0,31,59,90,120,151,181,212,243,273,304,334/
      DATA LYS / 0,31,60,91,121,152,182,213,244,274,305,335/
      DATA LENM/31,28,31,30 ,31, 30, 31, 31, 30, 31, 30, 31/
 
      IF (DABS(IDIF) .GT. 3153600000D0) GO TO 9200
      IF ( IDIF .LT. 0D0) THEN
        M6000 = 6000
        M60   = 60
        M24   = 24
        ISGN  = -1
      ELSE
        M6000 = -6000
        M60   = -60
        M24   = -24
        ISGN  = 1
      ENDIF
 
C          Reset arrays for leap year:
      IF (MOD(IOLDT(1),4) .EQ. 0) THEN
        DO 10 I=1,12
   10   LMSUM(I)=LYS(I)
        LENM(2)=29
        LASTDA=366
      ELSE
        DO 20 I=1,12
   20   LMSUM(I)=NLS(I)
        LENM(2)=28
        LASTDA=365
      ENDIF
 
C          Determine # of days, hours, mins and seconds to be added
      IDIFDA = IDINT (IDIF    /   8640000D0)
      NDIFDA = DBLE (IDIFDA)
      IDIFHR = IDINT((IDIF-NDIFDA*8640000D0)   /   360000D0)
      NDIFHR = DBLE (IDIFHR)
      IDIFMI = IDINT((IDIF-NDIFDA*8640000D0-NDIFHR*360000D0)/6000D0)
      NDIFMI = DBLE (IDIFMI)
      IDIFCS = IDINT((IDIF-NDIFDA*8640000D0-NDIFHR*360000D0-
     +                                                NDIFMI*6000D0))
 
C          Convert old time to day # of year & reformat hhmm
      IOMO = IOLDT(2)/100
      IODA = IOLDT(2)-IOMO*100
      IOJD = LMSUM(IOMO)+IODA
      IOHR = IOLDT(3)/100
      IOMI = IOLDT(3)-IOHR*100
 
C          Now build new time starting with centi-seconds:
      NEWCS  = IOLDT(4) + IDIFCS
      ICARRY = 0
      IF (NEWCS .GE. 6000 .OR. NEWCS .LT. 0) THEN
        NEWCS  = NEWCS + M6000
        ICARRY = ISGN
      ENDIF
      NEWT(4) = NEWCS
 
C          Make new minutes & hours
      NEWMI  = IOMI  + IDIFMI + ICARRY
      ICARRY = 0
      IF (NEWMI .GE. 60 .OR. NEWMI .LT. 0) THEN
        NEWMI  = NEWMI + M60
        ICARRY = ISGN
      ENDIF
      NEWHR  = IOHR + IDIFHR + ICARRY
      ICARRY = 0
      IF (NEWHR .GE. 24 .OR. NEWHR .LT. 0) THEN
        NEWHR  = NEWHR+M24
        ICARRY = ISGN
      ENDIF
      NEWT(3) = NEWHR*100 + NEWMI
 
C          Make new day number of year
      NEWJD = IOJD + IDIFDA + ICARRY
      NEWT(1) = IOLDT(1)
      IF (NEWJD .GT. LASTDA .OR. NEWJD .LE. 0) THEN
        NEWT(1) = IOLDT(1) + ISGN
        IF (IDIF .GE. 0D0) NEWJD = NEWJD-LASTDA
C          Reset arrays for leap year:
        IF (MOD(NEWT(1),4) .EQ. 0) THEN
          DO 11 I=1,12
   11     LMSUM(I) = LYS(I)
          LENM(2) = 29
          LASTDA = 366
        ELSE
          DO 21 I=1,12
   21     LMSUM(I) = NLS(I)
          LENM(2) = 28
          LASTDA = 365
        ENDIF
        IF (IDIF .LT. 0D0) NEWJD = NEWJD + LASTDA
      ENDIF
C          Now convert day number of year to mmdd
      DO 30 I=1,12
      ITEST = NEWJD - LMSUM(I)
      IF (ITEST .GT. 0 .AND. ITEST .LE. LENM(I)) GO TO 40
   30 CONTINUE
      STOP
   40 NEWMO = I
      NEWDA = NEWJD - LMSUM(I)
 
      NEWT(2) = NEWMO*100 + NEWDA
 
      RETURN
 
 9200 WRITE(MSGUN,'(''0Subroutine NEWTIM has not been tested for '',
     +                              ''making new times more'',/,
     +              '' than 365 days from the old time.'')')
      STOP
      END
      SUBROUTINE CVTI2R(MSGUN,TIME,ITIME)
C          This subroutine converts ITIME to TIME, where
C            MSGUN = Fortran unit for diagnostic messages
C            TIME  = double precision version of time:  yyddd.fraction
C            ITIME = integer version of time: (1) = year 4-digit
C                                             (2) = month/day
C                                             (3) = hour/minute
C                                             (4) = centisecond
C          MISS = Missing data value
 
C          Formal argument declarations
      DIMENSION ITIME(4)
      DOUBLE PRECISION TIME
 
C          Local declarations
      DOUBLE PRECISION FRACT , RMISS
      PARAMETER (MISS=-32767 , RMISS=-32767.D0)
      DIMENSION JDAZ(12) , LENM(12)
      SAVE JDAZ , LENM
      DATA JDAZ/ 0,31,28,31,30,31,30,31,31,30,31,30/
      DATA LENM/31,28,31,30,31,30,31,31,30,31,30,31/
 
      TIME = RMISS
      IF (ITIME(1) .EQ. MISS) RETURN
 
      YY = FLOAT( ITIME(1)-1900 )
      IF (YY .LT. 0.) GO TO 9100
 
      JDAZ(3) = 28
      LENM(2) = 28
      LENYR   = 365
      IF (MOD(ITIME(1),4) .EQ. 0) THEN
        JDAZ(3) = 29
        LENM(2) = 29
        LENYR   = 366
      ENDIF
 
      NMON = ITIME(2)/100
      NDAY = ITIME(2) - NMON*100
      IF (NMON .LT. 0 .OR. NMON .GT. 12) GO TO 9100
      IF (NDAY .GT. LENM(NMON)) GO TO 9100
      IDDD = 0
      DO 100 I=1,NMON
  100 IDDD = IDDD + JDAZ(I)
      DDD = FLOAT(IDDD + NDAY)
      IF (DDD .LT. 0. .OR. DDD .GT. 366.) GO TO 9100
 
      HR = FLOAT( ITIME(3)/100)
      RM = FLOAT( ITIME(3)) - HR*100.
      SS = FLOAT( ITIME(4)) / 100.
      IF (HR .LT. 0. .OR. HR .GT.  24.) GO TO 9100
      IF (RM .LT. 0. .OR. RM .GT.  59.) GO TO 9100
      IF (SS .LT. 0. .OR. SS .GT.  59.99) GO TO 9100
      FRACT = (DBLE(HR) + DBLE(RM + SS/60.)/60.D0 )/24.D0
 
      IF (IDDD .EQ. LENYR) THEN
        IF (IDINT(FRACT+1.E-10) .EQ. 1) THEN
          YY    = YY + 1
          DDD   = 1.
          FRACT = FRACT - 1.D0
        ENDIF
      ENDIF
 
      TIME = DBLE (YY*1000. + DDD) + FRACT
 
      RETURN
 
C          Error trap messages
 9100 WRITE(MSGUN,'('' CVTI2R: warning goofy date input:'',4I9)')
     +                                     (ITIME(I),I=1,4)
      RETURN
      END
      SUBROUTINE GETEHP(MSGUN,IBEGT,IENDT,MXNPWR,
     +                  NVALS,ITIM,HPWR,NDX,IQL,ISTAT,
     +                  IWRK,IBLK,IUNPK,LREC)
C          This routine returns Dave Evans hemispheric power input data
C          from the vsn, EHP781102A.  NVALS values for the time interval
C          IBEGT to IENDT are returned in arrays ITIM, HPWR, NDX and IQL.
C          INPUTS:
C             MSGUN     = Fortran unit number to print diagnostic msgs
C             IBEGT(1)  = Begin UT year, four digits
C             IBEGT(2)  = Begin UT month and day, e.g. 31 Jan is 0131
C             IBEGT(3)  = Begin UT hour and minute (0000 to 2400)
C             IBEGT(4)  = Begin UT centi-second
C             IENDT(1-4)= End date and time of interval for which
C                         parameters are returned
C             MXNPWR    = Dimension of arrays provided to this routine,
C                         IY, IMD, ITIM , HPWR, NDX, and IQL
C             ISTAT     = 99999 on input, implies only fatal errors
C                         will cause diagnostic messages to be printed,
C                         nonfatal errors will be signified only by the
C                         returned value of ISTAT (thus allowing user
C                         control of most diagnostic print messages).
C                         SEE ALSO RETURN ISTAT values below.
C             IWRK      = Work array for RDBLK6 (dimensioned to 17)
C             IBLK      = Block array for RDBLK6 (see RDBLK6 for dim)
C             IUNPK     = Unpacked version of IBLK for RDBLK6
C             LREC      = Current lrec array for RDBLK6
C          RETURNS:
C             NVALS = Number of values found
C             ITIM  = U.T. Date and time (same 4 integer format as IBEGT)
C             HPWR  = Hemispheric power input (units = GWatts)
C             NDX   = Dave Evans hemispheric power index
C             IQL   = HPWR qualifier; contains the measurement
C                     satellite ID in the tens digit and the
C                     hemisphere of measurement in the units digit:
C                     if the tens digit is 0, TIROS-N made the
C                     measurement, otherwise the digit is the number
C                     of the NOAA satellite; if the units digit is 0,
C                     the northern hemisphere was sampled, if the
C                     the units digit is 1 it is the southern
C                     hemisphere; e.g., 61 => NOAA-6 in Southern Hem.
C             ISTAT = Status flag
C                   = 0 => ok
C                   = 1 => Bad time input
C                   = 2 => Request time before first available
C                   = 3 => Request time after last available
 
C          Formal parameter argument declarations:
      DIMENSION ITIM(4,MXNPWR), HPWR(MXNPWR), NDX(MXNPWR), IQL(MXNPWR),
     +          IBEGT(4), IENDT(4) ,
     +         IWRK(*), IBLK(*) , IUNPK(*) , LREC(*)
 
C          Local declarations
      DOUBLE PRECISION  IDIF                          , I1DA
      PARAMETER (MISS=-32767 , MXNDAZ=40 , MXNVLD=200 , I1DA=8640000D0 ,
     +           NOPRT=99999)
C          MXNDAZ = Maximum number of different dates (yr/mm/dd)
C          MXNVLD = Maximum number of values per day
      DIMENSION IYI(MXNDAZ) , IMDI(MXNDAZ) , IT(4) , JT(4) ,
     +          ITIMD(MXNVLD),HPWRD(MXNVLD),NDXD(MXNVLD),IQLD(MXNVLD)
C          (arrays with D suffix hold a days worth of values)
      NVALS = 0
 
C          Check input time
      NOK = 0
      IF(IBEGT(1) .LT. 1940 .OR. IBEGT(1) .GT. 2100)NOK = 1
      IF(IBEGT(2) .LT. 0101 .OR. IBEGT(2) .GT. 1231)NOK = 1
      IF(IBEGT(3) .LT.    0 .OR. IBEGT(3) .GT. 2400)NOK = 1
      IF(IBEGT(4) .LT.    0 .OR. IBEGT(4) .GT. 5999)NOK = 1
      IF(IENDT(1) .LT. 1940 .OR. IENDT(1) .GT. 2100)NOK = 1
      IF(IENDT(2) .LT. 0101 .OR. IENDT(2) .GT. 1231)NOK = 1
      IF(IENDT(3) .LT.    0 .OR. IENDT(3) .GT. 2400)NOK = 1
      IF(IENDT(4) .LT.    0 .OR. IENDT(4) .GT. 5999)NOK = 1
      IF(NOK .GT. 0)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(IBEGT(I),I=1,4),
     +                                       (IENDT(J),J=1,4)
        ISTAT = 1
        RETURN
      ENDIF
C          Check that input end time is later than begin time:
      CALL TIMDIF(IBEGT,IENDT,IDIF)
      IF(IDIF .LT. 0D0)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(IBEGT(I),I=1,4),
     +                                       (IENDT(J),J=1,4)
        ISTAT = 1
        RETURN
      ENDIF
 
C          Determine how many different dates are requested and build
C          year, IYI, and month/day, IMDI, arrays
      NDATES  = 1
      IYI(1)  = IBEGT(1)
      IMDI(1) = IBEGT(2)
      IT(1)   = IBEGT(1)
      IT(2)   = IBEGT(2)
      IT(3)   = 0
      IT(4)   = 0
      IF(IDIF .LE. I1DA .AND. IBEGT(2) .EQ. IENDT(2))GO TO 100
   10 CALL NEWTIM(IT,JT,I1DA,MSGUN)
      NDATES = NDATES + 1
      IF(NDATES .GT. MXNDAZ)GO TO 9200
      IYI (NDATES) = JT(1)
      IMDI(NDATES) = JT(2)
      IF(JT(1).EQ.IENDT(1) .AND. JT(2).EQ.IENDT(2))GO TO 100
      IT(1)        = JT(1)
      IT(2)        = JT(2)
      GO TO 10
 
C          Loop once for each unique date requested
  100 JST = 0
      DO 1000 I=1,NDATES
C          Get a day's worth of data:
      IF(IWRK(7) .EQ. 1) GO TO 200
      IF(IYI(I) .EQ. LREC(5) .AND. IMDI(I) .EQ. LREC(6))GO TO 300
  200 JSTAT = ISTAT
      CALL GETHPD(MSGUN,IYI(I),IMDI(I),NVALSD,ITIMD,HPWRD,NDXD,IQLD,
     +             JSTAT,IWRK,IBLK,IUNPK,LREC)
      IF(JSTAT .NE. 0)JST = JSTAT
      IF(NVALSD .GT. MXNVLD)GO TO 9300
 
C          Assign the fraction of the days data needed
  300 IF(I .EQ. 1)THEN
C          This is first date requested; do not assign values before
C          requested begin time (hour and minute)
        DO 350 J=1,NVALSD
        IF(ITIMD(J) .GE. IBEGT(3))GO TO 400
  350   CONTINUE
        GO TO 1000
  400   ISKP  = J - 1
        NVALS = NVALSD - ISKP
        IF(NDATES .EQ. 1)THEN
C          Only one date to fill, so must also check end time and omit
C          assignment of any values after the end time:
          DO 450 J=NVALSD,1,-1
          IF(ITIMD(J) .LE. IENDT(3))GO TO 500
  450     CONTINUE
          GO TO 9300
  500     NOMIT = NVALSD - J
          NVALS = NVALS - NOMIT
        ENDIF
        DO 600 J=1,NVALS
        ITIM(1,J) = IYI (I)
        ITIM(2,J) = IMDI(I)
        ITIM(3,J) = ITIMD(ISKP + J)
        ITIM(4,J) = 0
        HPWR(J) = HPWRD(ISKP + J)
        NDX (J) = NDXD (ISKP + J)
  600   IQL (J) = IQLD (ISKP + J)
 
      ELSEIF(I .NE. NDATES)THEN
C          Not first or last date requested, so assign the whole day
        DO 700 J=1,NVALSD
        ITIM(1,J+NVALS) = IYI (I)
        ITIM(2,J+NVALS) = IMDI(I)
        ITIM(3,J+NVALS) = ITIMD(J)
        ITIM(4,J+NVALS) = 0
        HPWR(J+NVALS) = HPWRD(J)
        NDX (J+NVALS) = NDXD (J)
  700   IQL (J+NVALS) = IQLD (J)
        NVALS = NVALS + NVALSD
 
      ELSE
C          This is the last date requested; do not assign values after
C          requested end time (hour and minute)
        DO 750 J=NVALSD,1,-1
        IF(ITIMD(J) .LT. IENDT(3))GO TO 800
  750   CONTINUE
        GO TO 1000
  800   LASTOK = J
        DO 900 J=1,LASTOK
        ITIM(1,J+NVALS) = IYI (I)
        ITIM(2,J+NVALS) = IMDI(I)
        ITIM(3,J+NVALS) = ITIMD(J)
        ITIM(4,J+NVALS) = 0
        HPWR(J+NVALS) = HPWRD(J)
        NDX (J+NVALS) = NDXD (J)
  900   IQL (J+NVALS) = IQLD (J)
        NVALS = NVALS + LASTOK
      ENDIF
 1000 CONTINUE
      ISTAT = JST
      RETURN
 
C          Error trap messages
 9100 FORMAT('0GETEHP:  Bad time input; begin time =',4I5,/,
     +       '                            end time =',4I5)
 9200 WRITE(MSGUN,'(''0GETEHP:  Dimension exceeded, increase MXNDAZ '',
     +              ''to at least '',I7)')NDATES
      STOP
 9300 WRITE(MSGUN,'(''0GETEHP:  Dimension exceeded, increase MXNVLD '',
     +              ''to at least '',I5)')NVALSD
      STOP
      END
      SUBROUTINE GETHPD(MSGUN,IYR,IMMDD,NVALSD,ITIMD,HPWRD,NDXD,IQLD,
     +                  ISTAT,IWRK,IBLK,IUNPK,LREC)
C          This routine returns 1 day of values for the input year, IYR
C          (4 digit) and month-day, IMMDD.  Returned values are:
C                NVALSD = number of values in the following arrays:
C                ITIMD  = U.T. hour and minute of day
C                HPWRD  = Hemispheric power input (units = GWatts)
C                NDXD   = Dave Evans hemispheric power index
C                IQLD   = HPWRD qualifier; contains the measurement
C                         satellite id in the tens digit and the
C                         hemisphere of measurement in the units digit:
C                         if the tens digit is 0, TIROS-N made the
C                         measurement, otherwise the digit is the number
C                         of the NOAA satellite; for hemispheres the
C                         units digit is 0 if northern hemisphere and 1
C                         for the southern hemisphere; e.g., 61 => NOAA-6
C                         in the southern hemisphere.
C                ISTAT  = Status flag
C                       = 99999 ON INPUT, implies only fatal errors
C                         will cause diagnostic messages to be printed,
C                         nonfatal errors will be signified only by the
C                         returned value of ISTAT (thus allowing user
C                         control of most diagnostic print messages).
C                         ISTAT VALUES ON RETURN:
C                       = 0 => ok
C                       = 2 => Request time before first available
C                       = 3 => Request time after last available
C          Notes:
C          - MSGUN is the fortran unit number to print diagnostic
C            messages; see also ISTAT above!
C          - It is assumed that input IYEAR and IMMDD are reasonable
C            since GETEHP does this checking.
C          - The following arrays are used by RDBLK6, which must be
C            dimensioned appropriately by the initial calling routine:
C              IWRK  = Work array (contains fortran unit to read
C                      indicies, etc - see RDBLK6)
C              IBLK  = Packed physical record buffer.
C              IUNPK = Unpacked block buffer.
C              LREC  = Logical record.
C          Formal parameter argument declarations:
C          (the arrays in the first line should be dimmed NVALSD)
      DIMENSION ITIMD(*),HPWRD(*),NDXD(*),IQLD(*),
     +          IWRK(*),IBLK(*),IUNPK(*),LREC(*)
 
C          Local declarations:
      DOUBLE PRECISION IDIF , IOLDIF  , I1DA , JDIF
      PARAMETER (MPAR=6 , NOPRT=99999 , I1DA=8640000D0)
      DIMENSION IT(4) , LBEGT(4) , IOLDB(4)
      LOGICAL ATBEGN
      SAVE    ATBEGN,LBEGT,IOLDB,ISTRD
      DATA    ATBEGN,LBEGT,IOLDB,ISTRD
     +        /.TRUE. , 1940,0101,0,0, 1940,0101,0,0 , 0/
 
C          See if already read past the requested date:
      IT(1) = IYR
      IT(2) = IMMDD
      IT(3) = 0
      IT(4) = 0
      CALL TIMDIF(LBEGT,IT,IDIF)
      IF (IDIF .EQ. 0D0) GO TO 1000
      IF (IDIF .GT. 0D0) THEN
C          Must read more but first see if already at EOD
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)IYR,IMMDD,(LREC(I),I=5,6)
        ISTAT = 3
        RETURN
      ENDIF
C          Already read past the request date
      IF(ATBEGN)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)(IT(I),I=1,2),
     +                                     (LREC(J),J=5,6)
        ISTAT = 2
        RETURN
      ENDIF
C          Maybe no data are available for request date
      CALL TIMDIF(IOLDB,LBEGT,IOLDIF)
C          IOLDIF = time of last lrec read - previous one
      IF (IOLDIF .GT. I1DA) THEN
        NVALSD = 0
        RETURN
      ENDIF
      CALL REWBLK(IWRK)
      IWRK(7) = 1
      ATBEGN  = .TRUE.
 
C          Read a new lrec, compute & test the new difference
  100 NVALSD = 0
      CALL RDBLK6(IWRK,LREC,IBLK,IUNPK,ISTRD)
      IF(ISTRD .EQ. 0 .AND. LREC(2) .NE. 1002)GO TO 100
      IOLDB(1) = LBEGT(1)
      IOLDB(2) = LBEGT(2)
      LBEGT(1) = LREC(5)
      LBEGT(2) = LREC(6)
      CALL TIMDIF(LBEGT,IT,JDIF)
      IF (JDIF .GT. 0D0) THEN
        ATBEGN = .FALSE.
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(IT(I),I=1,2),
     +                                       (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ELSEIF (JDIF .LT. 0D0) THEN
        IF(ATBEGN)THEN
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)(IT(I),I=1,2),
     +                                        (LREC(J),J=5,6)
          ISTAT = 2
        ENDIF
        RETURN
      ENDIF
 
C          Assign new last begin time read and values, then return
 1000 NVALSD   = LREC(16)
      JOFF=16
      DO 1100 I=1,NVALSD
      JOFF=JOFF+MPAR
      ITIMD(I) =       LREC(3+JOFF)
      HPWRD(I) = FLOAT(LREC(4+JOFF)) / 10.
      NDXD (I) =       LREC(5+JOFF)
 1100 IQLD (I) =       LREC(6+JOFF)
      RETURN
 
C          Error trap messages:
 9100 FORMAT('0GETHPD:  Requested date',2I5,' is after last available:
     +',2I5)
 9200 FORMAT('0GETHPD:  Requested date',2I5,' is before first available
     +: ',2I5)
      END
      SUBROUTINE GETGPD(MSGUN,ITIME,IKP,IAP,IAPD,IF107,IFQFR,ISSPT,
     +                  ISTAT,IWRK,IBLK,IUNPK,LREC)
C          This routine returns 1 days (1 logical records) values for
C          the UT date and hour to arrays IKP and IAP (3 hourly)
C          and daily values of IAPD, IF107, IFQFR, ISSPT; where,
C                        IKP(8) =  8 3-hour samples of Kp for the day
C                        IAP(8) =  8 3-hour samples of ap
C                        IAPD   =  dialy mean AP
C                        IF107  =  F10.7 Solar Flux (daily)
C                        IFQFR  =  F10.7 Solar Flux Qualifier (daily)
C                        ISSPT  =  Sun spot number (daily)
C          Notes
C          - Although ITIME is in standard I.S. time format (four
C            integers holding year, monthday, hourminute, and centisec)
C            only the the year, month and day are used here.
C          - It is assumed that input time has reasonable values since
C            GETGPI normally does this checking.
C          - MSGUN is the fortran unit number to print diagnostic
C            messages.
C          - ISTAT  = Status flag
C                   = 99999 ON INPUT, implies only fatal errors
C                     will cause diagnostic messages to be printed,
C                     nonfatal errors will be signified only by the
C                     returned value of ISTAT (thus allowing user
C                     control of most diagnostic print messages).
C                     ISTAT VALUES ON RETURN:
C                   = 0 => ok
C                   = 2 => Request time before first available
C                   = 3 => Request time after last available
C                   = 6 => Indicies dataset corrupted; should be
C                          continuous, yet somehow read past ITIME
C          - The following arrays are used by RDBLK6, which must be
C            dimensioned appropriately by the initial calling routine:
C              IWRK  = Work array (contains fortran unit to read
C                      indicies, etc - see RDBLK6)
C              IBLK  = Packed physical record buffer.
C              IUNPK = Unpacked block buffer.
C              LREC  = Logical record.
 
C          Formal parameter argument declarations:
      DIMENSION ITIME(4), IKP(8), IAP(8), IWRK(*),IBLK(*), IUNPK(*),
     +          LREC(*)
 
C          Local declarations:
      PARAMETER (MPAR=6 , NOPRT=99999 , MISS=-32767)
      DIMENSION IT(4) , LBEGT(4)
      DOUBLE PRECISION IDIF , JDIF
      LOGICAL ATBEGN
      SAVE    ATBEGN,LBEGT,ISTRD
      DATA    ATBEGN,LBEGT,ISTRD/.TRUE. , 1940,0101,0,0 , 0/
 
      DO 10 I=1,8
      IKP(I) = MISS
   10 IAP(I) = MISS
      IAPD   = MISS
      IF107  = MISS
      IFQFR  = MISS
      ISSPT  = MISS
 
C          See if already read past the requested date:
      IT(1) = ITIME(1)
      IT(2) = ITIME(2)
      IT(3) = 0
      IT(4) = 0
      CALL TIMDIF(LBEGT,IT,IDIF)
      IF (IDIF .EQ. 0D0) GO TO 1000
      IF (IDIF .GT. 0D0) THEN
C          Must read more but first see if already at EOD
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(ITIME(I),I=1,2),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ENDIF
      CALL REWBLK(IWRK)
      IWRK(7) = 1
      ATBEGN  = .TRUE.
 
C          Read a new lrec, compute & test the new difference
  100 CALL RDBLK6(IWRK,LREC,IBLK,IUNPK,ISTRD)
      IF(ISTRD .EQ. 0 .AND. LREC(2) .NE. 1002)GO TO 100
      LBEGT(1)=LREC(5)
      LBEGT(2)=LREC(6)
      CALL TIMDIF(LREC(5),IT,JDIF)
      IF (JDIF .GT. 0D0) THEN
C          Must read more but first see if already at EOD
        ATBEGN = .FALSE.
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)ITIME(1),ITIME(2),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ELSEIF (JDIF .LT. 0D0) THEN
C          After request time
        IF(ATBEGN)THEN
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)ITIME(1),ITIME(2),
     +                                        (LREC(J),J=5,6)
          ISTAT = 2
        ELSE
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9500)ITIME(1),ITIME(2)
          ISTAT = 6
        ENDIF
        RETURN
      ENDIF
 
C          Assign new last begin time read and values, then return
C          (JOFF = lprol + 2*jpar = 16 + 8)
 1000 JOFF=24
      DO 1100 I=1,8
      JOFF=JOFF+MPAR
      IKP(I)=LREC(  5  + JOFF)
 1100 IAP(I)=LREC(  6  + JOFF)
      IAPD  =LREC( 21 )
      IF107 =LREC( 22 )
      IFQFR =LREC( 23 )
      ISSPT =LREC( 24 )
      ISTAT = 0
      RETURN
 
C          Error trap messages:
 9100 FORMAT('0GETGPD:  Requested date',2I5,' is after last available:
     +',2I5)
 9200 FORMAT('0GETGPD:  Requested date',2I5,' is before first available
     +: ',2I5)
 9500 FORMAT('0GETGPD:  Indicies dataset corrupted?  Data should be',
     +                                        ' continuous yet',/,
     +       '           somehow read past requested date',2I5)
      END
      SUBROUTINE GETIMF(MSGUN,ITIME,BXGSM,BYGSM,BZGSM,BYGSE,BZGSE,
     +                  SWDEN,SWSPD,IQL, ISTAT, IWRK,IBLK,IUNPK,LREC)
C          This routine returns OMNI tape hourly IMF and Solar Wind data
C          from the vsn, IMF631127A.  If not available for a requested
C          time, value(s) are returned as "-32767." (the standard missing
C          data code) except for the IMF/Solar Wind Qualifier.  Although
C          the standard time format is input as ITIME, only the year,
C          month, day and hour which are used here.  It is assumed that
C          any minute of the hour is a request for the same hourly mean
C          value.  (Centiseconds are ignored):
C                   0000 - 0059    returns same value
C                   0100 - 0159       "      "   "
C                        .            .      .   .
C                   2300 - 2400*      "      "   "
C          *  When ITIME(3)=2400, this is considered to be an "end" time
C             and the last index of the day is returned.  Although one
C             could argue this is inconsistant, it is easy enough to
C             input 0000 on the next date for the next hours indicies.
C          INPUTS:
C            ITIME(1) = Year, four digits Universal time
C            ITIME(2) = Month and day, e.g. 31 Jan is 0131
C            ITIME(3) = Hour and minute of day
C            ITIME(4) = Centisecond (hundredth second of minute)
C            MSGUN    = Fortran unit number to print diagnostic messages
C            ISTAT    = 99999 on input, implies only fatal errors
C                       will cause diagnostic messages to be printed,
C                       nonfatal errors will be signified only by the
C                       returned value of ISTAT (thus allowing user
C                       control of most diagnostic print messages).
C                       SEE ALSO RETURN ISTAT values below.
C            IWRK     = Work array used by RDBLK6 (contains fortran unit
C                       to read indicies, etc - see RDBLK6)
C            IBLK     = Packed physical record buffer - see RDBLK6.
C            IUNPK    = Unpacked block buffer - see RDBLK6.
C            LREC     = Logical record - see RDBLK6.
C          RETURNS:
C            ISTAT = Status flag
C                  = 0 => ok
C                  = 1 => Bad time input
C                  = 2 => Request time before first available
C                  = 3 => Request time after last available
C                  = 4 => Nothing available for requested date
C            BXGSM = X-component IMF(Geocentric Solar Magnetospheric)
C            BYGSM = Y-component IMF GSM Coordinates
C            BXGSM = Z-component IMF GSM
C            BYGSE = Y-component IMF Geocentric Solar Ecliptic
C            BZGSE = Z-component IMF GSE
C                      (Note Bx GSE is not explicitly returned since
C                      it is equivalent to Bx GSM)
C            IQL   = IMF/Solar Wind Qualifier, code 2236, contains
C                    satellite id codes for the measurement of the IMF
C                    and solar wind:  The 1000's and 100's digit
C                    comprise the IMF id and the 10's and 1's digit is
C                    the solar wind id.  Id codes follow:
C                    ---Satellite Name ---               --- ID ---
C                       no measurement                        0
C                       IMP 1 (Expl 18)                      18
C                       IMP 3 (Expl 28)                      28
C                       IMP 4 (Expl 34)                      34
C                       IMP 5 (Expl 41)                      41
C                       IMP 6 (Expl 43)                      43
C                       IMP 7 (Expl 47)                      47
C                       IMP 8 (Expl 50)                      50
C                       AIMP 1 (Expl 33)                     33
C                       AIMP 2 (Expl 35)                     35
C                       HEOS 1 and HEOS 2                     1
C                       VELA 3                                3
C                       OGO 5                                 5
C                       Merged LANL VELA speeds (7/64-3/71)  99
C                       Merged LANL IMP T,N,V (including     98
C                          all IMP 8 LANL plasma)
C                       ISEE 1                               11
C                       ISEE 2                               12
C                       ISEE 3                               13
C                       Prognoz10                            10
      DIMENSION ITIME(4) , IWRK(*) , IBLK(*) , IUNPK(*) , LREC(*)
 
C          Local declarations:
      DOUBLE PRECISION IDIF , IOLDIF , JDIF           , I1DA
      PARAMETER (MPAR=12 , MISS=-32767 , RMISS=-32767., I1DA=8640000D0,
     +           NOPRT=99999)
      DIMENSION IT(4) , LBEGT(4) , IOLDB(4)
      LOGICAL ATBEGN
      SAVE    ATBEGN,LBEGT,IOLDB,ISTRD
      DATA    ATBEGN,LBEGT,IOLDB,ISTRD
     +       /.TRUE. , 1940,0101,0,0 , 1940,0101,0,0 , 0/
 
C          Initialize returned values to missing:
      BXGSM = RMISS
      BYGSM = RMISS
      BZGSM = RMISS
      BYGSE = RMISS
      BZGSE = RMISS
      SWDEN = RMISS
      SWSPD = RMISS
      IQL   = 0
 
C          Check input time
      NOK =0
      IF(ITIME(1) .LT. 1940 .OR. ITIME(1) .GT. 2100)NOK = 1
      IF(ITIME(2) .LT. 0101 .OR. ITIME(2) .GT. 1231)NOK = 1
      IF(ITIME(3) .LT.    0 .OR. ITIME(3) .GT. 2400)NOK = 1
      IF(NOK .GT. 0)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9000)(ITIME(I),I=1,4)
        ISTAT = 1
        RETURN
      ENDIF
 
C          See if already read past the requested date:
      IT(1) = ITIME(1)
      IT(2) = ITIME(2)
      IT(3) = 0
      IT(4) = 0
      CALL TIMDIF(LBEGT,IT,IDIF)
      IF (IDIF .EQ. 0D0) GO TO 1000
      IF (IDIF .GT. 0D0) THEN
C          Must read more but first see if already at EOD
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(ITIME(I),I=1,2),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ENDIF
C          Already read past the request date
      IF(ATBEGN)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)(IT(I),I=1,2),
     +                                     (LREC(J),J=5,6)
        ISTAT = 2
        RETURN
      ENDIF
C          Maybe no data are available for request date
      CALL TIMDIF(IOLDB,LBEGT,IOLDIF)
C          IOLDIF = time of last lrec read - previous one
      IF(IOLDIF .GT. I1DA)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9300)(IT(I),I=1,2)
        ISTAT = 4
        RETURN
      ENDIF
      CALL REWBLK(IWRK)
      IWRK(7) = 1
      ATBEGN  = .TRUE.
 
C          Read a new lrec, compute & test the new difference
  100 CALL RDBLK6(IWRK,LREC,IBLK,IUNPK,ISTRD)
      IF(ISTRD .EQ. 0 .AND. LREC(2) .NE. 1002)GO TO 100
      IOLDB(1) = LBEGT(1)
      IOLDB(2) = LBEGT(2)
      LBEGT(1) = LREC(5)
      LBEGT(2) = LREC(6)
      CALL TIMDIF(LBEGT,IT,JDIF)
      IF (JDIF .GT. 0D0) THEN
C          Must read more
        ATBEGN = .FALSE.
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(IT(I),I=1,2),
     +                                       (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ELSEIF (JDIF .LT. 0D0) THEN
C          Already read past the request date
        IF(ATBEGN)THEN
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)(IT(I),I=1,2),
     +                                       (LREC(J),J=5,6)
          ISTAT = 2
        ELSE
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9300)(IT(I),I=1,2)
          ISTAT = 4
        ENDIF
        RETURN
      ENDIF
 
C          Determine offset into LREC appropriate for assigning values
C          for beginning hour, IHR
 1000 IHR = ITIME(3)/100
      IF(IHR .EQ. 24)IHR = 23
      IOFF = 33 + MPAR*IHR
      IF(LREC(IOFF) .NE. MISS) BXGSM = FLOAT(LREC(IOFF))/100.
      IOFF = IOFF + 1
      IF(LREC(IOFF) .NE. MISS) BYGSM = FLOAT(LREC(IOFF))/100.
      IOFF = IOFF + 1
      IF(LREC(IOFF) .NE. MISS) BZGSM = FLOAT(LREC(IOFF))/100.
      IOFF = IOFF + 1
      IF(LREC(IOFF) .NE. MISS) BYGSE = FLOAT(LREC(IOFF))/100.
      IOFF = IOFF + 1
      IF(LREC(IOFF) .NE. MISS) BZGSE = FLOAT(LREC(IOFF))/100.
      IOFF = IOFF + 1
      IF(LREC(IOFF) .NE. MISS) SWDEN = FLOAT(LREC(IOFF))*100000.
      IOFF = IOFF + 1
      IF(LREC(IOFF) .NE. MISS) SWSPD = FLOAT(LREC(IOFF))*100.
      IOFF = IOFF + 1
                                 IQL =       LREC(IOFF)
      ISTAT = 0
      RETURN
 
C          Error trap messages:
 9000 FORMAT('0GETIMF:  Bad time input; time =',4I5)
 9100 FORMAT('0GETIMF:  Requested date',2I5,' is after last available: '
     +,2I5)
 9200 FORMAT('0GETIMF:  Requested date',2I5,' is before first available:
     + ',2I5)
 9300 FORMAT('0GETIMF:  No data available for date',2I5)
      END
      SUBROUTINE GETAE(MSGUN,ITIME,IPARCOD,VALS, ISTAT ,
     +                 IWRK,IBLK,IUNPK,LREC)
C          GETAE returns 1-minute or hourly Ae, Al, Au, or Ao for the
C          requested time.  The A_ indicies are read from the I.S.
C          data base vsn, AEI780101A.  If not available for a requested
C          time, value(s) are returned as "-32767." (the standard missing
C          data code).  Although the standard time format is input as
C          ITIME, the minute and centisecond are ignored; i.e., any
C          minute of the hour is interpreted as a request for that hours
C          value(s).  If 1-minute samples are requested, 60 values will
C          be returned but a request for an hourly average results in
C          only one returned value:
C                   0000 - 0059    returns same hourly value(s)
C                   0100 - 0159       "      "    "      "
C                        .            .      .    .      .
C                   2300 - 2400*      "      "    "      "
C          *  When ITIME(3)=2400, this is considered to be an "end" time
C             and the last value(s) of the day are returned.  Although
C             one could argue this is inconsistant, it is easy enough to
C             input 0000 on the next date for the next hours indicies.
C
C          INPUTS:
C            MSGUN    = Fortran unit to print diagnostic messages; see
C                       also ISTAT.
C            ITIME(1) = Year, four digits Universal time
C            ITIME(2) = Month and day, e.g. 31 Jan is 0131
C            ITIME(3) = Hour and minute of day
C            ITIME(4) = Centisecond (hundredth second of minute)
C            IPARCOD  = The parameter code for the desired index;
C                       recognized codes include:
C                        320 = Ae Index (1 or 2.5 min sample*)
C                        321 = Al Index (1 or 2.5 min sample*)
C                        322 = Au Index (1 or 2.5 min sample*)
C                        323 = Ao Index (1 or 2.5 min sample*)
C                        324 = Ae Index (hourly mean)
C                        325 = Al Index (hourly mean)
C                        326 = Au Index (hourly mean)
C                        327 = Ao Index (hourly mean)
C                              * Always 1 minute samples after 1977
C            ISTAT    = 99999 on input, implies only fatal errors
C                       will cause diagnostic messages to be printed,
C                       nonfatal errors will be signified only by the
C                       returned value of ISTAT (thus allowing user
C                       control of most diagnostic print messages).
C                       SEE ALSO RETURN ISTAT values below.
C            IWRK     = Work array used by RDBLK6 (contains fortran unit
C                       to read indicies, etc - see RDBLK6)
C            IBLK     = Packed physical record buffer - see RDBLK6.
C            IUNPK    = Unpacked block buffer - see RDBLK6.
C            LREC     = Logical record - see RDBLK6.
C          RETURNS:
C            ISTAT    = Status flag
C                     = 0 => ok
C                     = 1 => Bad time input
C                     = 2 => Request time before first available
C                     = 3 => Request time after last available
C                     = 5 => Bad parameter code input
C                     = 6 => Indicies dataset corrupted; should be
C                            continuos, yet somehow read past ITIME
C            VALS     = The index value(s) for IPARCOD at ITIME; units
C                       are nT (nanoTeslas) or gammas.
C          Formal parameter declarations:
      DIMENSION ITIME(4) , VALS(60) , IWRK(*),IBLK(*),IUNPK(*),LREC(*)
 
C          Local declarations:
      PARAMETER (RMISS=-32767. , NOPRT=99999)
      DIMENSION LBEGT(4) , IT(4)
      DOUBLE PRECISION IDIF , JDIF
      LOGICAL ATBEGN
      SAVE    ATBEGN,LBEGT,ISTRD
      DATA    ATBEGN,LBEGT,ISTRD/.TRUE. , 1940,0101,0000,0000 , 0/
 
C          Initialize returned values array to missing
      DO 10 I=1,60
   10 VALS(I) = RMISS
 
C          Check input time
      NOK = 0
      IF(ITIME(1) .LT. 1940 .OR. ITIME(1) .GT. 2100)NOK = 1
      IF(ITIME(2) .LT. 0101 .OR. ITIME(2) .GT. 1231)NOK = 1
      IF(ITIME(3) .LT.    0 .OR. ITIME(3) .GT. 2400)NOK = 1
      IF(NOK .GT. 0)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9000)(ITIME(I),I=1,4)
        ISTAT = 1
        RETURN
      ENDIF
 
C          See if already read past the requested date and hour:
      IT(1) =  ITIME(1)
      IT(2) =  ITIME(2)
      IT(3) = (ITIME(3)/100) * 100
      IF (IT(3) .EQ. 2400) IT(3) = 2300
      IT(4) = 0
      CALL TIMDIF(LBEGT,IT,IDIF)
      IF (IDIF .EQ. 0D0) GO TO 1000
      IF (IDIF .GT. 0D0) THEN
C          Must read more but first see if already at EOD
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(ITIME(I),I=1,3),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ENDIF
      CALL REWBLK(IWRK)
      IWRK(7) = 1
      ATBEGN  = .TRUE.
 
C          Read a new lrec, compute & test the new difference
  100 CALL RDBLK6(IWRK,LREC,IBLK,IUNPK,ISTRD)
      IF(ISTRD .EQ. 0 .AND. LREC(2) .NE. 1002)GO TO 100
      LBEGT(1) = LREC(5)
      LBEGT(2) = LREC(6)
      LBEGT(3) = LREC(7)
      CALL TIMDIF(LREC(5),IT,JDIF)
      IF (JDIF .GT. 0D0) THEN
C          must read more
        ATBEGN = .FALSE.
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(ITIME(I),I=1,3),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ELSEIF (JDIF .LT. 0D0) THEN
C          after request time
        IF(ATBEGN)THEN
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)(ITIME(I),I=1,3),
     +                                         (LREC(J),J=5,6)
          ISTAT = 2
        ELSE
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9500)(ITIME(I),I=1,3)
          ISTAT = 6
        ENDIF
        RETURN
      ENDIF
 
C          Assign value(s) for requested parameter code
 1000 IOFF = 0
      IF(IPARCOD .EQ. 320)IOFF = 31
      IF(IPARCOD .EQ. 321)IOFF = 32
      IF(IPARCOD .EQ. 322)IOFF = 33
      IF(IPARCOD .EQ. 323)IOFF = 34
      IF(IPARCOD .EQ. 324)IOFF = 14
      IF(IPARCOD .EQ. 325)IOFF = 15
      IF(IPARCOD .EQ. 326)IOFF = 16
      IF(IPARCOD .EQ. 327)IOFF = 17
      IF(IOFF .EQ. 0)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9400)IPARCOD
        ISTAT = 5
        RETURN
      ENDIF
      NVALS = 1
      IF(IPARCOD .LT. 324)NVALS = 60
      DO 1100 I=1,NVALS
      IOFF = IOFF + 8
 1100 VALS(I) = FLOAT( LREC(IOFF) )
      ISTAT = 0
      RETURN
 
C          Error trap messages
 9000 FORMAT('0GETAE:  Bad time input; time =',4I5)
 9100 FORMAT('0GETAE:  Requested time',3I5,' is after last available: ',
     + 2I5)
 9200 FORMAT('0GETAE:  Requested time',3I5,' is before first available:
     +',2I5)
 9400 FORMAT('0GETAE:  Does not recognize parameter code:',I7,/,
     +       '         Known parameter codes include:',/,
     +       '           320 = Ae Index (1 or 2.5 min sample*)  ',/,
     +       '           321 = Al Index (1 or 2.5 min sample*)  ',/,
     +       '           322 = Au Index (1 or 2.5 min sample*)  ',/,
     +       '           323 = Ao Index (1 or 2.5 min sample*)  ',/,
     +       '           324 = Ae Index (hourly mean)           ',/,
     +       '           325 = Al Index (hourly mean)           ',/,
     +       '           326 = Au Index (hourly mean)           ',/,
     +       '           327 = Ao Index (hourly mean)           ',/,
     +       '                 * Always 1 minute samples after 1977 ')
 9500 FORMAT('0GETAE:  Indicies dataset corrupted?  Data should be',
     +                                        ' continuous yet',/,
     +       '         somehow read past requested time',3I5)
      END
      SUBROUTINE GETDST(MSGUN,ITIME,DST, ISTAT ,
     +                 IWRK,IBLK,IUNPK,LREC)
C          GETDST returns the hourly DST index for the UT time input,
C          ITIME.  The indicies are read from the NCAR I.S. database vsn
C          DST570101A.  If not available for a requested time, DST is
C          returned as "-32767." (the standard missing data code).
C          Although the standard time format is input as ITIME, the
C          minute and centisecond are ignored; i.e., any minute of the
C          hour is interpreted as a request for that hour's value:
C                   0000 - 0059    returns same hourly value
C                   0100 - 0159       "      "    "      "
C                        .            .      .    .      .
C                   2300 - 2400*      "      "    "      "
C          *  When ITIME(3)=2400, this is considered to be an "end" time
C             and the last value(s) of the day are returned.  Although
C             one could argue this is inconsistant, it is easy enough to
C             input 0000 on the next date for the next hours indicies.
C
C          INPUTS:
C            MSGUN    = Fortran unit to print diagnostic messages; see
C                       also ISTAT.
C            ITIME(1) = Year, four digits Universal time
C            ITIME(2) = Month and day, e.g. 31 Jan is 0131
C            ITIME(3) = Hour and minute of day
C            ITIME(4) = Centisecond (hundredth second of minute)
C            ISTAT    = 99999 on input, implies only fatal errors
C                       will cause diagnostic messages to be printed,
C                       nonfatal errors will be signified only by the
C                       returned value of ISTAT (thus allowing user
C                       control of most diagnostic print messages).
C                       SEE ALSO RETURN ISTAT values below.
C            IWRK     = Work array used by RDBLK6 (contains fortran unit
C                       to read indicies, etc - see RDBLK6)
C            IBLK     = Packed physical record buffer - see RDBLK6.
C            IUNPK    = Unpacked block buffer - see RDBLK6.
C            LREC     = Logical record - see RDBLK6.
C          RETURNS:
C            ISTAT    = Status flag
C                     = 0 => ok
C                     = 1 => Bad time input
C                     = 2 => Request time before first available
C                     = 3 => Request time after last available
C                     = 6 => Indicies dataset corrupted; should be
C                            continuos, yet somehow read past ITIME
C            DST      = The DST value for time ITIME; units are nT
C                       (nanoTeslas) or gammas.
C          Formal parameter declarations:
      DIMENSION ITIME(4) , IWRK(*),IBLK(*),IUNPK(*),LREC(*)
 
C          Local declarations
      PARAMETER (MISS=-32767 , RMISS=-32767. , NOPRT=99999 ,
     +           IOFF=26,MPAR=5)
C                IOFF=LPROL + 2*JPAR + 2*MPAR
      DIMENSION LBEGT(4) , IT(4)
      DOUBLE PRECISION IDIF , JDIF
      LOGICAL ATBEGN
      SAVE    ATBEGN,LBEGT,ISTRD
      DATA    ATBEGN,LBEGT,ISTRD/.TRUE. , 1940,0101,0000,0000 , 0/
 
C          Initialize returned values array to missing
      DST = RMISS
 
C          Check input time
      NOK = 0
      IF(ITIME(1) .LT. 1940 .OR. ITIME(1) .GT. 2100)NOK = 1
      IF(ITIME(2) .LT. 0101 .OR. ITIME(2) .GT. 1231)NOK = 1
      IF(ITIME(3) .LT.    0 .OR. ITIME(3) .GT. 2400)NOK = 1
      IF(NOK .GT. 0)THEN
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9000)(ITIME(I),I=1,4)
        ISTAT = 1
        RETURN
      ENDIF
 
C          See if already read past the requested date and hour:
      IT(1) = ITIME(1)
      IT(2) = ITIME(2)
      IT(3) = 0
      IT(4) = 0
      CALL TIMDIF(LBEGT,IT,IDIF)
      IF(IDIF .EQ. 0D0)GO TO 1000
      IF(IDIF .GT. 0D0)THEN
C          Must read more but first see if already at EOD
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(ITIME(I),I=1,2),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ENDIF
      CALL REWBLK(IWRK)
      IWRK(7) = 1
      ATBEGN  = .TRUE.
 
C          Read a new lrec, compute & test the new difference
  100 CALL RDBLK6(IWRK,LREC,IBLK,IUNPK,ISTRD)
      IF(ISTRD .EQ. 0 .AND. LREC(2) .NE. 1002)GO TO 100
      LBEGT(1) = LREC(5)
      LBEGT(2) = LREC(6)
      CALL TIMDIF(LREC(5),IT,JDIF)
      IF(JDIF .GT. 0D0)THEN
C          must read more
        ATBEGN = .FALSE.
        IF(ISTRD .NE. 3)GO TO 100
        IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9100)(ITIME(I),I=1,2),
     +                                        (LREC(J),J=5,6)
        ISTAT = 3
        RETURN
      ELSEIF(JDIF .LT. 0D0)THEN
C          after request time
        IF(ATBEGN)THEN
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9200)(ITIME(I),I=1,2),
     +                                         (LREC(J),J=5,6)
          ISTAT = 2
        ELSE
          IF(ISTAT .NE. NOPRT)WRITE(MSGUN,9500)(ITIME(I),I=1,2)
          ISTAT = 6
        ENDIF
        RETURN
      ENDIF
 
C          Compute position of DST value for the requested hour
 1000 IHR = ITIME(3)/100
      IF(IHR .EQ. 24)IHR = 23
      IPOS = IOFF + IHR*MPAR
      DST = FLOAT( LREC(IPOS) )
      ISTAT = 0
      RETURN
 
C          Error trap messages:
 9000 FORMAT('0GETDST:  Bad time input; time =',4I5)
 9100 FORMAT('0GETDST:  Requested time',2I5,' is after last available: '
     +,2I5)
 9200 FORMAT('0GETDST:  Requested time',2I5,' is before first available:
     + ',2I5)
 9500 FORMAT('0GETDST:  Indicies dataset corrupted?  Data should be',
     +                                        ' continuous yet',/,
     +       '          somehow read past requested time',2I5)
      END
      SUBROUTINE GTMTIC (MSGUN,ITIM,MXNMJR,NMJR,PMJR,LBL,NLBC,NMNR)
C          This routine determines tick marks for time axis labelling.
C          It is set up to return the following options:
C              Tmax-Tmin   Labelled ticks every       Unlabelled:
C             -----------  --------------------     ----------------
C                <     0 -> (error stop)
C              0 ->  1 hr        1/4 hr (units=hr)   5 min (MINR_=2)
C                ->  2 hr        1/2 hr              5 min (MINR_=5)
C                -> 12 hr          1 hr             1/4 hr (MINR_=3)
C                -> 24 hr          2 hr              1 hr  (MINR_=1)
C                -> 48 hr          4 hr              1 hr         3
C                -> 72 hr          6 hr              1 hr         5
C                ->  6 days       12 hr              2 hr         5
C                -> 12 days        1 day (units=day) 1/4 day      3
C                -> 24 days        2 day             1/2 day      3
C                -> 48 days        4 day             1   day      3
C                -> 96 days        8 day             2   day      3
C                -> 1.5 yr     1st&15th of month       none       0
C                 > 1.5 yr  -> (error stop)
C          INPUTS:
C            MSGUN  = Fortran unit number to print any error diagnostics
C            ITIM   = Begin and end time (8 element integer array with
C                     4 elements for each, in form yyyy mmdd hhmm csec)
C            MXNMJR = Maximum number of labelled ticks allowed (array
C                     dimension of PMJR and LBL).  For most cases 13
C                     is sufficient, however, intervals longer than 96
C                     days use 2 per month (e.g., for 1 yr provide 24).
C          RETURNS:
C            NMJR   = Number of computed major tick marks.
C            PMJR   = Fractional positions (domain: 0. -> 1.) of computed
C                     tick marks which are to be labelled.
C            LBL    = Labels for major ticks.
C            NLBC   = Maximum number of characters per label (2 or 4)
C            NMNR   = Number of unlabelled ticks to be inserted between
C                     major (labelled) ticks.
C          Declarations for formal arguments:
      DIMENSION   PMJR(MXNMJR) , ITIM(8)
      CHARACTER*4  LBL(MXNMJR)
 
C          Local declarations:
      DOUBLE PRECISION IDT , IDCS , ITOTCS , ITHCS
      PARAMETER (MXYL=24)
      DIMENSION IOLDT(4), NEWT(4) , IYMD(MXYL)
      CHARACTER*4 LYMD(MXYL)
      CHARACTER*3 UNITS
      SAVE IYMD , LYMD
      DATA IYMD /  101 , 115 ,  201 , 215 ,  301 , 315 ,  401 , 415 ,
     +             501 , 515 ,  601 , 615 ,  701 , 715 ,  801 , 815 ,
     +             901 , 915 , 1001 ,1015 , 1101 ,1115 , 1201 ,1215 /
      DATA LYMD /'1Jan',' 15','1Feb',' 15','1Mar',' 15','1Apr',' 15',
     +           '1May',' 15','1Jun',' 15','1Jul',' 15','1Aug',' 15',
     +           '1Sep',' 15','1Oct',' 15','1Nov',' 15','1Dec',' 15'/
 
C          Determine which option is to be used (IDIF and IDT are in
C          units of centiseconds)
      CALL TIMDIF (ITIM(1),ITIM(5),ITOTCS)
      IF (ITOTCS .LT. 0D0) THEN
        GO TO 9100
      ELSEIF (ITOTCS .LE.   360000D0) THEN
C                                         1 hr and  1/4 hr
        IDT  =    90000D0
        NMNR =        2
        NLBC =        4
        UNITS= 'qtr'
      ELSEIF (ITOTCS .LE.   720000D0) THEN
C                                         2 hr and  1/2 hr
        IDT  =   180000D0
        NMNR =        5
        NLBC =        4
        UNITS= 'haf'
      ELSEIF (ITOTCS .LE.  4320000D0) THEN
C                                        12 hr and  1 hr
        IDT  =   360000D0
        NMNR =        3
        NLBC =        2
        UNITS= 'hr '
      ELSEIF (ITOTCS .LE.  8640000D0) THEN
C                                        24 hr and  2 hr
        IDT  =   720000D0
        NMNR =        1
        NLBC =        2
        UNITS= 'hr '
      ELSEIF (ITOTCS .LE. 17280000D0) THEN
C                                        48 hr and  4 hr
        IDT  =  1440000D0
        NMNR =        3
        NLBC =        2
        UNITS= 'hr '
      ELSEIF (ITOTCS .LE. 25920000D0) THEN
C                                        72 hr and  6 hr
        IDT  =  2160000D0
        NMNR =        5
        NLBC =        2
        UNITS= 'hr '
      ELSEIF (ITOTCS .LE. 51840000D0) THEN
C                                         6 days and 12 hrs
        IDT  =  4320000D0
        NMNR =        5
        NLBC =        2
        UNITS= 'hr '
      ELSEIF (ITOTCS .LE.103680000D0) THEN
C                                        12 days and 1 day
        IDT  =  8640000D0
        NMNR =        3
        NLBC =        2
        UNITS= 'day'
      ELSEIF (ITOTCS .LE.207360000D0) THEN
C                                        24 days and  2 day
        IDT  = 17280000D0
        NMNR =        3
        NLBC =        2
        UNITS= 'day'
      ELSEIF (ITOTCS .LE.414720000D0) THEN
C                                        48 days and  1 day
        IDT  = 34560000D0
        NMNR =        3
        NLBC =        2
        UNITS= 'day'
      ELSEIF (ITOTCS .LE.829440000D0) THEN
C                                        96 days and  2 day
        IDT  = 69120000D0
        NMNR =        3
        NLBC =        2
        UNITS= 'day'
      ELSEIF (ITOTCS .LE.4743360000D0) THEN
C                                      1.5yr (549 days) and  no minor ticks
C       IDT  = __________ (not used)
        NMNR =        0
        NLBC =        4
        UNITS= 'yr'
      ELSE
        GO TO 9200
      ENDIF
 
C          Process year differently because of varying time differences
C          between labels:
      IF (UNITS .EQ. 'yr') GO TO 1000
 
C          Round up to find the value and position of the first label
C          To round up, start with a 00:00 UT of the begin date and
C          test that time against the begin time; if it is less
C          than the begin time, add IDT to it and test again.  Continue
C          adding and testing until it is greater.
      IOLDT(1) = ITIM(1)
      IOLDT(2) = ITIM(2)
      IOLDT(3) = 0
      IOLDT(4) = 0
      NEWT (1) = ITIM(1)
      NEWT (2) = ITIM(2)
      NEWT (3) = 0
      NEWT (4) = 0
      CALL TIMDIF (IOLDT,ITIM,IDCS)
      IF (IDCS .LE. 0D0) GO TO 200
  100 CALL NEWTIM (IOLDT,NEWT,IDT,MSGUN)
      CALL TIMDIF (NEWT,ITIM,IDCS)
      IF (IDCS .LE. 0D0) GO TO 200
      DO 150 I=1,4
  150 IOLDT(I) = NEWT(I)
      GO TO 100
 
  200 IF (UNITS .EQ. 'qtr' .OR. UNITS .EQ. 'haf') THEN
        ILB = NEWT(3)
        WRITE (LBL(1),'(I4)') ILB
      ELSEIF (UNITS .EQ. 'hr ') THEN
        ILB = NEWT(3) / 100
        WRITE (LBL(1),'(I2,2X)') ILB
      ELSEIF (UNITS .EQ. 'day') THEN
        ILB = NEWT(2) - (NEWT(2)/100) * 100
        WRITE (LBL(1),'(I2,2X)') ILB
      ENDIF
      CALL TIMDIF (ITIM,NEWT,ITHCS)
      PMJR(1) = SNGL( ITHCS/ITOTCS )
      NMJR    = 1
 
C          Build the rest of the labels
      DO 400 I=2,MXNMJR
      DO 300 J=1,4
  300 IOLDT(J) = NEWT(J)
      CALL NEWTIM (IOLDT,NEWT,IDT,MSGUN)
      CALL TIMDIF (NEWT,ITIM(5),IDCS)
      IF (IDCS .LT. 0D0) GO TO 401
      IF (UNITS .EQ. 'qtr' .OR. UNITS .EQ. 'haf') THEN
        ILB = NEWT(3)
        IF (ILB .EQ. 0) ILB = 2400
        WRITE (LBL(I),'(I4)') ILB
      ELSEIF (UNITS .EQ. 'hr ') THEN
        ILB = NEWT(3) / 100
        IF (ILB .EQ. 0) ILB = 24
        WRITE (LBL(I),'(I2,2X)') ILB
      ELSEIF (UNITS .EQ. 'day') THEN
        ILB = NEWT(2) - (NEWT(2)/100) * 100
        WRITE (LBL(I),'(I2,2X)') ILB
      ENDIF
      CALL TIMDIF (ITIM,NEWT,ITHCS)
      PMJR(I) = SNGL( ITHCS/ITOTCS )
      NMJR = NMJR + 1
  400 CONTINUE
  401 CONTINUE
 
      RETURN
 
C          Build labels for year intervals
 1000 IBMO = ITIM(2)/100
      IBDA = ITIM(2) - IBMO*100
      IF (IBMO .LT. 1 .OR. IBMO .GT. 12) GO TO 9300
 
C          Find index for the first month/day label on or after begin date
      I = IBMO * 2
      IF (IBDA .LT. 15) I = I - 1
      IF (IBDA .EQ. 15 .AND. ITIM(3) .LT. 1200) I = I - 1
      LBL(1) = LYMD(I)
 
C          Compute the fractional position of the first label
      IOLDT(1) = ITIM(1)
      IOLDT(2) = ITIM(2)
      IOLDT(3) = 0
      IOLDT(4) = 0
      NEWT (1) = ITIM(1)
      NEWT (2) = IYMD(I)
      NEWT (3) = 0
      NEWT (4) = 0
      CALL TIMDIF (IOLDT,ITIM,IDCS)
      IF (IDCS .LE. 0D0) IDCS = 0D0
      PMJR(1) = SNGL( IDCS/ITOTCS )
      NMJR    = 1
 
C          Build the rest of the labels
 1100 DO 1110 J=1,4
 1110 IOLDT(J) = NEWT(J)
      I = I + 1
      IF (I .GT. MXYL) THEN
        I = 1
        NEWT(1) = NEWT(1) + 1
      ENDIF
      NEWT(2) = IYMD(I)
      CALL TIMDIF (NEWT,ITIM(5),IDCS)
      IF (IDCS .LT. 0D0) GO TO 1200
      NMJR = NMJR + 1
      IF (NMJR .GT. MXNMJR) GO TO 1200
      CALL TIMDIF (ITIM,NEWT,ITHCS)
      PMJR(NMJR) = SNGL( ITHCS/ITOTCS )
      LBL(NMJR) = LYMD(I)
      GO TO 1100
 
 1200 CONTINUE
 
      RETURN
 
C          Error trap diagnostics:
 9100 WRITE(MSGUN,9150)(ITIM(I),I=1,8)
 9150 FORMAT(' GTMTIC is not smart enough to calculate ',
     +                                        'time tick marks for',
     +       ' a reversed scale; Tbeg must be < Tend but',/,
     +       '                   Tbeg =',4I5,' and',/,
     +       '                   Tend =',4I5)
      STOP
 9200 WRITE(MSGUN,9250)(ITIM(I),I=1,8)
 9250 FORMAT(' GTMTIC is not smart enough to calculate ',
     +                                        'time tick marks when',
     +       ' the range exceeds 1 1/2 yr; Tbeg =',4I5,' and',/,
     +       '                             Tend =',4I5)
      STOP
 9300 WRITE (MSGUN,'('' GTMTIC: bogus dates input '',8I8)') ITIM
      STOP
      END
      SUBROUTINE GAXTIC(MSGUN,YMIN,YMAX,MXNMJR,NMJR,PMJR,LBL,LLBL,NMNR)
C          This routine determines tick marks and labels for axis
C          labeling.  If this works as well as it ought to, credit should
C          go to Dave Kennison who condescended to give me a copy of CPNUMB
C          (renamed here to GNUMB) from CONPACK, Version 3.  Dec 89.
C          INPUTS:
C            MSGUN  = Fortran unit number to print any error diagnostics
C            YMIN   = Smallest value to plot
C            YMAX   = Biggest value to plot.  If YMAX = YMIN, then YMIN
C                     and YMAX will be changed.  If YMAX < YMIN, inverted
C                     scale is assumed; e.g., for a Y-axis, the biggest
C                     tick label value would be on the bottom.
C            MXNMJR = Maximum number of labelled ticks allowed (array
C                     dimension of PMJR and LBL).
C          RETURNS:
C            NMJR   = Number of computed major tick marks.
C            PMJR   = Fractional positions (domain: 0. -> 1.) of computed
C                     tick marks which are to be labelled.
C            LBL    = Labels for major ticks.  USER should anticipate
C                     up to MXNMJR labels of maximum length 14 characters
C                     Error diagnotics will be generated.
C            LLBL   = Number of assigned characters in corresponding LBL.
C            NMNR   = Number of unlabelled ticks to be inserted between
C                     major (labelled) ticks.
C          Declarations for formal arguments:
      DIMENSION    PMJR(MXNMJR) , LLBL(MXNMJR)
      CHARACTER*(*) LBL(MXNMJR)
 
C          Local declarations:
      PARAMETER (NSTEP=3)
      DIMENSION TSTP(NSTEP) , NMINR (NSTEP)
      SAVE      TSTP,NMINR
      DATA      TSTP / 1., 2., 5./
      DATA      NMINR/ 4 , 1 , 4 /
C          Labelling constraints:
C            - Logarithmic scaling is not considered.
C            - Ymin may be less than smallest labelled tick mark.
C            - Ymax may be greater than largest labelled tick mark.
C            - MXNMJR is expected to be reasonable, i.e., more than 15
C              labels becomes crowded; moreover, this scheme will require
C              extending the values of i another decade if MXNMJR exceeds
C              200; but I have an error trap at 50.  (See logic after
C              statement label 200 for definition of i.)
      IF(MXNMJR .GT. 50 .OR. MXNMJR .LT. 2)GO TO 9100
C          The following formula is used to describe label intervals;
C          it says the number of labelled tick intervals is less than
C          or equal to the range:
C
C               (ndy)(step)(10**m) <= R
C
C          where ndy  = number of label intervals; 0 < ndy < MXNMJR.
C                step = 1,2, or 5; a nice label interval.
C                m    = an integer.
      MXNDY = MXNMJR - 1
      RMXNDY = REAL( MXNDY )
      R = YMAX - YMIN
      ISGN = 1
      SGN  = 1.
      IF(R .LT. 0.)THEN
        R    = -R
        ISGN = -1
        SGN  = -1.
      ENDIF
      IF(R .EQ. 0.)THEN
C          Have a constant field to plot; so change Ymin and Ymax to one
C          order of magnitude less/more than current value.
        M = 0
        IF(YMIN .NE. 0.)M = INT( LOG10( ABS(YMIN) ) ) - 1
        DY   = 10.**M
        OYMN = YMIN
        OYMX = YMAX
        YMIN = YMIN - DY
        YMAX = YMAX + DY
        WRITE(MSGUN,9200)OYMN,OYMX,YMIN,YMAX
        R    = YMAX - YMIN
      ENDIF
C          Determine m as the next integer above the minimum possible
C          value.  Note: log R - log(MXNDY*5) =< m =< log R.  So unless
C          MXNDY is huge, m can only be one or two integers.  Find the
C          smallest m here (after statement 100 the next m will be used):
      TSTM = LOG10( R / ( RMXNDY*TSTP(NSTEP) ) )
      M = INT( TSTM )
      IF( TSTM - REAL(M) .GT. 0.) M = M + 1
C          Find step and dy by trial until error.  Since most labels
C          possible is preferred, use ndy = MXNDY and increment step
C          to make new test dy values until MXNDY * dy >= range:
      RDT = R/10.**M
      DO 100 I=1,NSTEP
      IF( RMXNDY*TSTP(I) .GE. RDT )GO TO 200
  100 CONTINUE
C          Must try the next integer for m:
      M = M + 1
      RDT = RDT/10.
      DO 150 I=1,NSTEP
      IF( RMXNDY*TSTP(I) .GE. RDT )GO TO 200
  150 CONTINUE
      STOP
  200 STEP = TSTP(I)
      NMNR = NMINR(I)
      DY   = STEP * 10.**M
 
C          Find smallest label value by finding an integer, i, such that
C          (i-1)(dy) < Ymin & (i)(dy) >= Ymin; i.e., i >= Ymin/dy
      SML = YMIN
      IF(ISGN .EQ. -1)SML = YMAX
      SDY = SML / DY
      IF(SDY .GE. 0.)THEN
        RI  = REAL( INT( SDY + 0.99999999) )
      ELSE
        RI  = REAL( INT( SDY ) )
      ENDIF
      YS = RI*DY
      IF( YS-DY .EQ. SML )YS = SML
 
C          Find exact ndy and biggest label value. This could not be done
C          precisely until Ys was established because there may be up to
C          2dy discrepancy between Ymax-Ymin and Yb-Ys.
      BIG = YMAX
      IF(ISGN .EQ. -1)BIG = YMIN
      CK = BIG - YS
      DO 250 I=MXNDY,1,-1
      IF( REAL(I)*DY .LE. CK )GO TO 300
  250 CONTINUE
  300 NDY  = I
      NMJR = NDY + 1
      YB   = YS + REAL(NDY) * DY
 
C          Determine the format of labels to be built; see comments at
C          beginning of gnumb for parameters description.
      IF(DY .GT. 1.E4 .OR. R .GT. 1.E6 .OR. ABS(YB) .GT. 1.E5 .OR.
     +  (YS .NE. 0. .AND. ABS(YS) .LT. 1.E-5)  )THEN
C          Use exponential form for labels.  Use difference between
C          largest and smallest label value scaled by minimum power of
C          ten in labels to determine how many mantissa digits to display
        YDT  = (YB-YS) / 10.**M
        IF(YDT .LE. 99.9)THEN
          NDGD = 3
          LMSD = M+1
        ELSEIF(YDT .LE. 999.9)THEN
          NDGD = 4
          LMSD = M+2
        ELSE
          NDGD = 3
          LMSD =-10000
        ENDIF
        IEXP =-1
        IOMA = 0
        IODP = 0
        IF(STEP .EQ. 1.)THEN
C          Maybe no mantissa is needed since dy is a power of ten.
C          Must also check to see if a label is a power of ten:
          IF(YS .EQ. 0.)THEN
            IOMA = 1
            IODP = 1
          ELSEIF( AMOD(LOG10(ABS(YS)),1.)  .LT. 0.001 )THEN
            IOMA = 1
            IODP = 1
          ENDIF
        ENDIF
        IOTZ = 0
      ELSE
C          Use non-exponential form; labels will be in range +/- .00001
C          to 100,000.  Determine LMSD, the left most significant digit,
C          based apon the larger of log(Ys) and log(Yb).  (See GNUMB for
C          definition of LMSD).
        YSL = 0.
        IF(YS  .NE. 0.)YSL = LOG10( ABS(YS) )
        IF(YSL .GE. 0.)THEN
          IYSL = INT(YSL)
        ELSE
          IYSL = INT(YSL) - 1
        ENDIF
        YBL = 0.
        IF(YB  .NE. 0.)YBL = LOG10( ABS(YB) )
        IF(YBL .GE. 0.)THEN
          IYBL = INT(YBL)
        ELSE
          IYBL = INT(YBL) - 1
        ENDIF
        LMSD = MAX (IYSL,IYBL)
        NDGD = LMSD - M + 1
        IEXP = 6
        IOMA = 0
        IODP = 0
        IF(AMOD(DY,1.) .LT. 0.01 .AND. AMOD(YS,1.) .LT. 0.01) IODP = 1
        IOTZ = 0
      ENDIF
 
C          Build the labels and their positions
      LMAX = LEN(LBL(1))
      Y1 = YS
      IF(ISGN .EQ. -1)Y1 = YB
      YI = Y1 - SGN*DY
      DO 1000 I=1,NMJR
      YI = YI + SGN*DY
      CALL GNUMB (YI  ,NDGD,LMSD,IEXP, 2  ,'x','10**',' ',1,4,0,
     +                      IOMA,IODP,IOTZ,LBL(I),LLBL(I),NDGS,IEVA)
      IF(LLBL(I) .GT. LMAX )WRITE(MSGUN,9300)LMAX,LLBL(I)
 1000 PMJR(I) = SGN*(YI-YMIN) / R
 
      RETURN
 
C          Error trap diagnostics:
 9100 WRITE(MSGUN,'('' GAXTIC:  MXNMJR is unreasonable!  The maximum'',
     +                           '' number of labelled Y-axis tick'',/,
     +              ''         marks should be a small positive '',
     +                       ''integer, not'',I15)')MXNMJR
      STOP
 9200 FORMAT(' GAXTIC: Warning, changed YMIN and YMAX',/,
     +       '             from:',1P2E11.4,/,
     +       '               to:',  2E11.4)
 9300 FORMAT(' GAXTIC: Warning, truncated an axis label.  Length ',
     +                                 'provided = ',I10,/,
     +       '                                             Lengt',
     +                                 'h needed = ',I10)
      END
      SUBROUTINE GNUMB (VALU,NDGD,LMSD,IEXP,LEXP,CEX1,CEX2,CEX3,LEX1,
     +                  LEX2,LEX3,IOMA,IODP,IOTZ,CBUF,NBUF,NDGS,IEVA)
C
      CHARACTER*(*) CEX1,CEX2,CEX3,CBUF
C          This routine was given to me by Dave Kennison.  It came from
C          CONPACK, Version 3 - NCAR Graphics, where it was called
C          CPNUMB.  I have made changes for implementation here:
C            (1) the common block has been omitted, and
C            (2) the routine, CPINRC, has been incorporated in the
C                if(first)then block at the beginning.
C            (3) some declarations from cpcomn
C          The peculiar indentation is a result of iftran!
C
C This subroutine expresses the value of a real number in a character
C form.  Depending on the values of the arguments, an exponential form
C (for example, "1.36E-2") or a no-exponent form (for example, ".0136")
C may be used.  The arguments are as follows:
C
C VALU is the real number whose value is to be expressed.
C
C NDGD is the desired number of significant digits to be used in the
C character expression of the number.
C
C LMSD is a flag indicating how the leftmost significant digit of VALU
C is to be determined.  VALU may be written in the form
C
C   ... D(3) D(2) D(1) D(0) . D(-1) D(-2) D(-3) D(-4) ...
C
C where, for all integer values of I, D(I) is an integer between 0 and
C 9, inclusive.  There exists an integer ILFT such that D(ILFT) is non-
C zero and, for all I greater than ILFT, D(I) is zero.  The leftmost
C significant digit of VALU is considered to occur in the position
C MAX(ILFT,LMSD).
C
C LMSD may be used to achieve consistency in expressing the values of a
C group of numbers.  For example, suppose that, with NDGD = 3 and LMSD
C = -10000, we get the numbers
C
C   5.00, 10.0, 15.0, ..., 95.0, 100., 105.              (no exponents)
C   5.00E0, 1.00E1, 1.50E1, ..., 9.50E1, 1.00E2, 1.05E2  (exponents)
C
C By resetting LMSD to 2 (which is the position of the leftmost non-zero
C digit in the whole group), we can get instead
C
C   5., 10., 15., ..., 95., 100., 105.                   (no exponents)
C   0.05E2, 0.10E2, 0.15E2, ..., 0.95E2, 1.00E2, 1.05E2  (exponents)
C
C Whether one prefers to see numbers like those in the first set or the
C second set is to some extent a matter of preference.  The second set
C includes fewer extraneous zeroes and allows the values with exponents
C to be compared with each other more easily.  Note that, in the case of
C the exponential form, LMSD may be viewed as specifying the minimum
C exponent value to be used.  Use LMSD = -10000 to indicate that no
C attempt should be made to force consistency.
C
C IEXP specifies how it is to be decided whether to use the exponential
C form or not, as follows:  If IEXP is less than or equal to zero, the
C exponential form is used, no matter what.  If IEXP is greater than
C zero, the no-exponent form is used if the length of the resulting
C string is less than or equal to IEXP; otherwise, the form resulting
C in the shorter string is used.
C
C LEXP is set less than or equal to zero if exponents are to be written
C in their shortest possible form (plus signs are omitted and the fewest
C digits required to express the value of the exponent are used).  LEXP
C is set greater than zero if exponents are to be written in a manner
C more nearly consistent with one another (the exponent is written with
C either a plus sign or a minus sign and the value of LEXP is the
C desired minimum number of digits to be used, leading zeroes being
C supplied to pad the exponent to the desired length).
C
C CEX1 and CEX2 are character strings to be used in the exponential form
C between the mantissa and the exponent.  If IOMA is non-zero, and, as
C a result, a mantissa exactly equal to one is omitted, CEX1 is omitted
C as well.  Blanks are treated as null strings.  Some possibilities are
C 1) CEX1='E' and CEX2=' ' (or vice-versa), which gives a sort of E
C format (in which case IOMA should not be set non-zero), 2) CEX1='x'
C and CEX2='10**', which gives numbers like "1.36453x10**13", and 3)
C CEX1=':L1:4' and CEX2='10:S:', which generates the function codes
C necessary to make the utility PLCHHQ write the number in exponential
C form.
C
C CEX3 is a character string to be used in the exponential form after
C the exponent.  This will usually be a blank, which is treated as a
C null string; an exception is when function codes for PLCHHQ are being
C generated, in which case it is desirable to use ':N:', in order to
C return to normal level.
C
C LEX1, LEX2, and LEX3 are the lengths to be assumed for the character
C strings CEX1, CEX2, and CEX3 in making decisions about the length of
C the exponential form and the no-exponent form.  (Note that these are
C not the actual lengths of the strings CEX1, CEX2, and CEX3.  If, for
C example, CEX1, CEX2 and CEX3 contain the function codes for PLCHHQ
C mentioned above, use LEX1=1, LEX2=2, and LEX3=0.)
C
C IOMA specifies whether or not it is permissible to omit, from the
C exponential form, mantissas of the form "1" or "1." which are not
C necessary to express the value (as, for example, in "1.x10**2").  If
C IOMA is non-zero, such mantissas are omitted; the part of the exponent
C given by CEX1 (probably the "x" above) is also omitted (thus changing
C "1.x10**2" into "10**2").  Such omission takes place even if IODP
C (which see, below) is zero.
C
C IODP specifies whether or not it is allowed to omit a decimal point
C which is unnecessary (as for example, in "23487.").  If IODP is
C non-zero, such decimal points are omitted.
C
C IOTZ specifies whether or not it is allowed to omit trailing zeroes.
C If IOTZ is non-zero, trailing zeroes are omitted.
C
C CBUF is a character buffer in which the character string is returned.
C If this buffer is not long enough to hold all the characters, no error
C results; the extra characters are simply lost.  This is potentially
C useful, since the object of the call may be simply to obtain the
C number of significant digits and the exponent value.
C
C NBUF is an output parameter; it says how many characters have been
C put into the character buffer CBUF.
C
C NDGS is an output parameter; it contains the number of significant
C digits which were used to express the value of VALU.
C
C IEVA is another output parameter; it is the power to which 10 must be
C raised to obtain a scale factor which will reduce VALU to the range
C from .1 to 1.  That is, the expression "VALU/10.**IEVA" is guaranteed
C (subject to round-off problems) to be greater than or equal to .1 and
C less than 1.  Another way of interpreting IEVA is that it specifies
C the position preceding the leftmost significant digit of VALU (where
C the one's position is numbered 0, the ten's position 1, the hundred's
C position 2, the tenth's position -1, etc.  Thus, the significant
C digits occur in positions IEVA-1 (the leftmost) through IEVA-NDGS
C (the rightmost).
C
C
C Declare all of the CONPACK common blocks.
C
C   .CALL CPCOMN,/$SAVE-COMMON$/1/
C          Declarations from cpcomn used here:
      CHARACTER*500 CTMB
      CHARACTER*8 FRMT
      SAVE NLZF ,  FRMT
      DATA NLZF / 0 /
C
C Declare a variable to hold single characters for testing purposes.
C
      CHARACTER*1 SCHR
 
C          The following is a patch insertion of routine cpinrc, --------
C          which allows me to avoid messing with common blocks:  --------
      LOGICAL FIRST
      SAVE    FIRST,NSDR
      DATA    FIRST / .TRUE. /
      IF(FIRST)THEN
        FIRST = .FALSE.
C CPINRC sets constants that are required by CONPACK and that cannot be
C defined in a DATA statement because determining their values requires
C that code be executed.
C
C Find out how many significant digits a real can represent and use it
C to compute machine constants "epsilon" and "1+epsilon" and to set up
C the format to be used by CPNUMB.
        NSDR=0
    1   NSDR=NSDR+1
        TMP1=10.**(-NSDR)
        TMP2=  1.+TMP1
        TMP3=TMP2+TMP1
        IF (.NOT.(TMP2.EQ.1..OR.TMP3.EQ.TMP2.OR.NSDR.GE.100))GO TO 1
 
        EPSI=10.**(1-NSDR)
        FRMT(1:2)='(E'
        IF (NSDR+8.LE.9) THEN
          FRMT(3:3)=CHAR(ICHAR('0')+NSDR+8)
          ITMP=4
        ELSE
          FRMT(3:3)=CHAR(ICHAR('0')+(NSDR+8)/10)
          FRMT(4:4)=CHAR(ICHAR('0')+MOD(NSDR+8,10))
          ITMP=5
        END IF
        FRMT(ITMP:ITMP)='.'
        IF (NSDR.LE.9) THEN
          FRMT(ITMP+1:ITMP+1)=CHAR(ICHAR('0')+NSDR)
          ITMP=ITMP+2
        ELSE
          FRMT(ITMP+1:ITMP+1)=CHAR(ICHAR('0')+(NSDR)/10)
          FRMT(ITMP+2:ITMP+2)=CHAR(ICHAR('0')+MOD(NSDR,10))
          ITMP=ITMP+3
        END IF
        FRMT(ITMP:ITMP)=')'
      ENDIF
C          End of patch insert of cpinrc -----------------------------
 
C
C Find the real lengths of the three parts of the exponent-creating
C string.
C
      LCX1=LEN(CEX1)
      IF (CEX1.EQ.' ') LCX1=0
      LCX2=LEN(CEX2)
      IF (CEX2.EQ.' ') LCX2=0
      LCX3=LEN(CEX3)
      IF (CEX3.EQ.' ') LCX3=0
C
C Find the length of the character buffer and initialize it to blanks.
C
      LBUF=LEN(CBUF)
      CBUF=' '
C
C Use the local I/O routines to generate an E-format representation of
C the number.
C
      WRITE (CTMB(1:NSDR+8),FRMT) VALU
C
C We're about to scan the E-format representation.  Initialize NBUF,
C which is the number of characters put into CBUF, NDGS, which is the
C number of significant digits found in CTMB, IDPT, which is the number
C of the significant digit after which the decimal point was found,
C IEXF, which is a flag indicating whether or not the exponent has been
C found yet, and IRND, which is a rounding flag.
C
      NBUF=0
      NDGS=0
      IDPT=0
      IEXF=0
      IRND=0
C
C Scan the E-format representation.
C
      DO    10 I=1,NSDR+8
C
C If a minus sign is found, and it's not part of the exponent, put it
C into the user's character buffer.  If it is a part of the exponent,
C set the exponent sign.  On the Cray, large exponents will cause the
C 'E' to be omitted, in which case the sign introduces the exponent.
C
        IF (CTMB(I:I).EQ.'-') THEN
          IF (IEXF.EQ.0) THEN
            IF (NDGS.EQ.0) THEN
              NBUF=NBUF+1
              IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)='-'
            ELSE
              IEXF=1
              IESI=-1
              IEVA=0
            END IF
          ELSE
            IESI=-1
          END IF
C
C If a plus sign is found, it can usually just be skipped.  On the Cray,
C large exponents will cause the 'E' to be omitted, in which case the
C sign introduces the exponent.
C
        ELSE IF (CTMB(I:I).EQ.'+') THEN
          IF (IEXF.EQ.0.AND.NDGS.NE.0) THEN
            IEXF=1
            IESI=1
            IEVA=0
          END IF
C
C If a digit is found, and it's not a part of the exponent, copy it to
C the beginning of the temporary buffer; save at most NDGD such digits.
C If a digit is found, and it's part of the exponent, update the value
C of the exponent.
C
        ELSE IF (ICHAR(CTMB(I:I)).GE.ICHAR('0').AND.ICHAR(CTMB(I:I)).LE.
     +IC
     +  HAR('9')) THEN
          IF (IEXF.EQ.0) THEN
            IF (NDGS.EQ.0) THEN
              IF (CTMB(I:I).NE.'0') THEN
                NDGS=1
                SCHR=CTMB(I:I)
                CTMB(1:1)=SCHR
                NZRS=0
                IF (SCHR.EQ.'9') THEN
                  NNNS=1
                ELSE
                  NNNS=0
                END IF
              ELSE
                IDPT=IDPT-1
              END IF
            ELSE IF (NDGS.LT.NDGD) THEN
              NDGS=NDGS+1
              SCHR=CTMB(I:I)
              CTMB(NDGS:NDGS)=SCHR
              IF (SCHR.EQ.'0') THEN
                NZRS=NZRS+1
                NNNS=0
              ELSE
                NZRS=0
                IF (SCHR.EQ.'9') THEN
                  NNNS=NNNS+1
                ELSE
                  NNNS=0
                END IF
              END IF
            ELSE IF (IRND.EQ.0) THEN
              IRND=1+(ICHAR(CTMB(I:I))-ICHAR('0'))/5
            END IF
          ELSE
            IEVA=10*IEVA+ICHAR(CTMB(I:I))-ICHAR('0')
          END IF
C
C If a decimal point is found, record the index of the digit which it
C followed.
C
        ELSE IF (CTMB(I:I).EQ.'.') THEN
          IDPT=NDGS
C
C If an "E" or an "e" is found, reset the flags to start processing of
C the exponent.
C
        ELSE IF (CTMB(I:I).EQ.'E'.OR.CTMB(I:I).EQ.'e') THEN
          IEXF=1
          IESI=1
          IEVA=0
        END IF
C
   10 CONTINUE
C
C If no significant digits were found, or if no exponent was found,
C assume that the number was exactly zero and return a character string
C reflecting that (unless the use of consistent exponents is forced,
C which requires special action).
C
      IF (NDGS.EQ.0.OR.IEXF.EQ.0) THEN
        IF (IEXP.GT.0.OR.LMSD.EQ.-10000) THEN
          CBUF='0'
          NBUF=1
          NDGS=1
          IEVA=0
          RETURN
        ELSE
          NBUF=0
          GO TO    30
        END IF
      END IF
C
C Round the number, take care of trailing zeroes and nines, and compute
C the final number of significant digits.
C
      IF (IRND.LT.2) THEN
        IF (NZRS.NE.0) NDGS=NDGS-NZRS
      ELSE
        IF (NNNS.NE.0) NDGS=NDGS-NNNS
        IF (NDGS.EQ.0) THEN
          IDPT=IDPT+1
          CTMB(1:1)='1'
          NDGS=1
        ELSE
          SCHR=CHAR(ICHAR(CTMB(NDGS:NDGS))+1)
          CTMB(NDGS:NDGS)=SCHR
        END IF
      END IF
C
C Compute the final value of the exponent which would be required if
C the decimal point preceded the first significant digit in CTMB.
C
      IEVA=IESI*IEVA+IDPT
C
C If the leftmost significant digit is to the right of the one the user
C wants, supply some leading zeroes and adjust the parameters giving the
C number of digits in CTMB and the exponent value.  We must provide for
C the possibility that this will reduce the number to zero.
C
      IF (IEVA-1.LT.LMSD) THEN
        NLZS=LMSD-(IEVA-1)
        IF (NLZS.LT.NDGD) THEN
          NDGT=MIN(NDGS+NLZS,NDGD)
          DO    40 I=NDGT,NLZS+1,-1
            SCHR=CTMB(I-NLZS:I-NLZS)
            CTMB(I:I)=SCHR
   40     CONTINUE
          DO    50 I=1,NLZS
            CTMB(I:I)='0'
   50     CONTINUE
          NDGS=NDGT
          IEVA=LMSD+1
        ELSE
          GO TO    30
        END IF
      ELSE
        NLZS=0
      END IF
C
C Control arrives at this block to generate a multi-digit zero.
C
      GO TO    70
   30 CONTINUE
        CTMB(1:1)='0'
        NDGS=1
        NLZS=0
        IEVA=LMSD+1
   70 CONTINUE
C
C Decide how many digits to output.  This depends on whether the user
C wants to omit trailing zeroes or not.
C
      IF (IOTZ.EQ.0) THEN
        NDTO=NDGD
      ELSE
        NDTO=NDGS
      END IF
C
C Compute the lengths of the character strings required for the form
C without an exponent (LWOE) and for the form with an exponent (LWIE).
C In certain cases, the values given are dummies, intended to force the
C use of one form or the other.  Note that leading zeroes are included
C in computing LWOE, even though they may be omitted from the output,
C in order to achieve consistency of sets of labels.
C
      IF (IEXP.GT.0) THEN
        LWOE=NBUF+MAX(NDTO,IEVA)-MIN(IEVA,0)
        IF (IEVA.LE.NLZS.AND.NLZF.NE.0) LWOE=LWOE+1
        IF (IEVA.GE.NDTO.AND.IODP.EQ.0) LWOE=LWOE+1
        IF (LWOE.LE.IEXP) THEN
          LWOE=0
          LWIE=0
        ELSE
          LWIE=NBUF+NDTO+2+LEX1+LEX2+LEX3
          IF (NDTO.EQ.1) THEN
            IF (IOMA.NE.0.AND.CTMB(1:1).EQ.'1') THEN
              LWIE=LWIE-2-LEX1
            ELSE IF (IODP.NE.0) THEN
              LWIE=LWIE-1
            END IF
          END IF
          IF (IEVA-1.LT.0.OR.LEXP.GT.0) LWIE=LWIE+1
          IF (ABS(IEVA-1).GT.9.OR.LEXP.GE.2) LWIE=LWIE+1
          IF (ABS(IEVA-1).GT.99.OR.LEXP.GE.3) LWIE=LWIE+1
          IF (ABS(IEVA-1).GT.999.OR.LEXP.GE.4) LWIE=LWIE+1
        END IF
      ELSE
        LWOE=1
        LWIE=0
      END IF
C
C Depending on the lengths, generate a string without an exponent ...
C
      IF (LWOE.LE.LWIE) THEN
C
        DO    80 I=MIN(IEVA+1,NLZS+1),MAX(NDTO,IEVA)
          IF (I.EQ.IEVA+1) THEN
            IF (I.LE.NLZS+1.AND.NLZF.NE.0) THEN
              NBUF=NBUF+1
              IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)='0'
            END IF
            NBUF=NBUF+1
            IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)='.'
          END IF
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) THEN
            IF (I.GE.1.AND.I.LE.NDGS) THEN
              CBUF(NBUF:NBUF)=CTMB(I:I)
            ELSE
              CBUF(NBUF:NBUF)='0'
            END IF
          END IF
   80   CONTINUE
C
        IF (IEVA.GE.NDTO.AND.IODP.EQ.0) THEN
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)='.'
        END IF
C
C ... or a string with an exponent.
C
      ELSE
C
        IF (NDTO.NE.1.OR.CTMB(1:1).NE.'1'.OR.IOMA.EQ.0) THEN
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)=CTMB(1:1)
        END IF
C
        IF (NDTO.NE.1.OR.((CTMB(1:1).NE.'1'.OR.IOMA.EQ.0).AND.IODP.EQ.0)
     +)
     +  THEN
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)='.'
        END IF
C
        DO    90 I=2,NDTO
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) THEN
            IF (I.LE.NDGS) THEN
              CBUF(NBUF:NBUF)=CTMB(I:I)
            ELSE
              CBUF(NBUF:NBUF)='0'
            END IF
          END IF
   90   CONTINUE
C
        IF (LCX1.NE.0.AND.(NDTO.NE.1.OR.CTMB(1:1).NE.'1'.OR.IOMA.EQ.0))
     +TH
     +  EN
          NBUF=NBUF+LCX1
          IF (NBUF.LE.LBUF) CBUF(NBUF-LCX1+1:NBUF)=CEX1
        END IF
C
        IF (LCX2.NE.0) THEN
          NBUF=NBUF+LCX2
          IF (NBUF.LE.LBUF) CBUF(NBUF-LCX2+1:NBUF)=CEX2
        END IF
C
        ITMP=IEVA-1
C
        IF (ITMP.LT.0.OR.LEXP.GT.0) THEN
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) THEN
            IF (ITMP.LT.0) THEN
              CBUF(NBUF:NBUF)='-'
            ELSE
              CBUF(NBUF:NBUF)='+'
            END IF
          END IF
        END IF
C
        ITMP=MIN(ABS(ITMP),9999)
C
        IF (ITMP.GT.999) THEN
          NTTL=4
          IDIV=1000
        ELSE IF (ITMP.GT.99) THEN
          NTTL=3
          IDIV=100
        ELSE IF (ITMP.GT.9) THEN
          NTTL=2
          IDIV=10
        ELSE
          NTTL=1
          IDIV=1
        END IF
C
        IF (LEXP.GT.0) THEN
          DO   100 I=1,LEXP-NTTL
            NBUF=NBUF+1
            IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)='0'
  100     CONTINUE
        END IF
C
        DO   110 I=1,NTTL
          NBUF=NBUF+1
          IF (NBUF.LE.LBUF) CBUF(NBUF:NBUF)=CHAR(ICHAR('0')+ITMP/IDIV)
          ITMP=MOD(ITMP,IDIV)
          IDIV=IDIV/10
  110   CONTINUE
C
        IF (LCX3.NE.0) THEN
          NBUF=NBUF+LCX3
          IF (NBUF.LE.LBUF) CBUF(NBUF-LCX3+1:NBUF)=CEX3
        END IF
C
      END IF
C
C Done.
C
      RETURN
C
      END
      SUBROUTINE PLTYVT(HDR,YLBL,Y,IT,MXNY,NY,YMIN,YMAX,ITMIN,ITMAX)
C          Plot Y versus T, time. Uses GKS version of NCAR Graphics,
C          and, therefore, writes plot metacode to FT02, also called
C          "GMETA".  OPNGKS and CLSGKS must be done before the first
C          and after the last call to this routine.
C          INPUTS:
C            HDR      = Header label; will be centered above plot on a
C                       single line (no longer than 80 characters).
C            YLBL     = Y-axis label; will be centered on Y-axis with
C                       text runnning from bottom to top on a single line
C            Y        = Y array, of values to be plotted.  Missing points
C                       are excluded; missing values = -32767.
C            IT       = X, or time, array.  Uses standard I.S. 4 integers
C                       notation:  yr, moda, hrmn, centisec.
C            MXNY     = True dimension of Y and second dimension of T.
C            NY       = Number of assigned values in Y and T; must be
C                       less than or equal to MXNY.
C            YMIN     = Minimum Y value to be plotted.  It does not need
C                       to be a "nice" value since rounding will be done
C                       to accomodate nice labels.  Ymin may not equal
C                       the actual minimum, thus, allowing multiple
C                       frames to be produced with the same Y-axis labels
C                       labels.  If YMIN is set to standard missing code,
C                       -32767., YMIN will be selected from the Y array.
C            YMAX     = Maximum Y value to be plotted (see YMIN).  If
C                       YMAX is set to standard missing code, YMAX will
C                       selected from the Y array.
C            ITMIN    = Time (X) axis min value (in standard I.S. four
C                       integer notation.  May be earlier than IT array
C                       first time
C            ITMAX    = Time (X) axis max value (in standard I.S. four
C                       integer notation.  May be later than IT array
C                       last time.
 
C          Formal argument declarations:
      CHARACTER*(*) HDR,YLBL
      DIMENSION  Y(MXNY) , IT(4,MXNY) , ITMIN(4) , ITMAX(4)
 
C          Local declarations:
      PARAMETER (MISS=-32767 , MSGUN=6)
 
      PARAMETER(YTOP=1.00 , YTPLT=0.85 , YBPLT=0.15 , YBTIC=0.10 ,
     +          YBOT=0.00 , YTTIC=0.90 ,
     +          XLFT=0.00 , XLTIC=0.05 , XLPLT=0.15 , XRPLT=0.88 ,
     +          XRTIC=1.00 , XRHT=1.00 )
C          Parameters delineating the various fractions of the frame to
C          be used:
C            YTOP  = Top of frame
C            YTPLT = Top of plot rectangle in frame
C            YBPLT = Bot of plot rectangle in frame
C            YBTIC = Bot of tick mark area for X axis labeling (area
C                    below plot rectangle - for labelled time ticks)
C            YBOT  = Bot of frame.
C            YTTIC = Top of tick mark area for X axis labeling (area
C                    above plot rectangle - for unlabelled time ticks)
C            XLFT  = Left edge of frame
C            XLTIC = Left edge tick mark area to left of plot area.  This
C                    side includes numeric discriptions of "major" ticks
C            XLPLT = Left edge of plot area.
C            XRPLT = Right edge of plot area.
C            XRTIC = Right edge of tick mark area to right of plot area
C                    (no numeric values)
C            XRHT  = Right edge of frame.
 
      PARAMETER(IBIGC=18 , BIGTL=17. , GAP=5. , WSML = 12. ,
     +          ISMLC=12 , SMLTL=10. ,          HSML = 24. ,RES=1024.)
C            IBIGC  = Large character size (1024 resolution)
C            ISMLC  = Small character size
C            BIGTL  = Length of labelled tick marks
C            SMLTL  = Length of unlabelled tick marks
C            GAP    = Distance from tick mark to edge of label
C            WSML   = Width of small character
C            HSML   = Height of small charcter
C            RES    = Resolution of entire plot frame
 
C          Define arrays for holding the position of labelled tick marks
C          and their labels
      PARAMETER (MXNMJRY=16 , MXNMJRT=13)
      DIMENSION   PMJRY(MXNMJRY) , PMJRT(MXNMJRT) , LLBLY(MXNMJRY)
      CHARACTER*14 LBLY(MXNMJRY)
      CHARACTER*4                   LBLT(MXNMJRT)
 
C          Definitions for building the time (X) axis labels and
C          day crossover plot lines
      DOUBLE PRECISION I1DA , I1MIN , IDIF , DBEGT,DENDT , DTJ,DTK,DTDC
      PARAMETER (MXNDCH=300 , I1DA=8640000D0 , I1MIN=6000D0)
      DIMENSION IRNDT(8) , ITMP(4) , ITP2(4) , TDC(MXNDCH)
      CHARACTER *41 HDRT
      CHARACTER *3  MON(12)
      DATA MON /'Jan' , 'Feb' , 'Mar' , 'Apr' , 'May' , 'Jun' , 'Jul',
     +          'Aug' , 'Sep' , 'Oct' , 'Nov' , 'Dec' /
 
C          Define the two RGB color triples to be used:
C          Color for index 0 is black
C          Color for index 1 is white
      PARAMETER (NCLR = 2)
      DIMENSION NDXC(NCLR) , RGB(3,NCLR)
      DATA NDXC/ 0 , 1 /
      DATA RGB / 0.000 , 0.000  , 0.000 ,
     +           1.000 , 1.000  , 1.000 /
 
C          See if Ymin must be determined
      IF (INT(YMIN) .EQ. MISS) THEN
        YMIN = Y(1)
        DO 1 I=2,NY
        IF(Y(I) .LT. YMIN)YMIN = Y(I)
    1   CONTINUE
      ENDIF
 
C          See if Ymax must be determined
      IF (INT(YMAX) .EQ. MISS) THEN
        YMAX = Y(1)
        DO 5 I=2,NY
        IF(Y(I) .GT. YMAX)YMAX = Y(I)
    5   CONTINUE
      ENDIF
 
C          Shut off clipping
      CALL GSCLIP(0)
 
C          Assign each color and its index
      DO 10 J=1,NCLR
   10 CALL GSCR(1,NDXC(J),RGB(1,J),RGB(2,J),RGB(3,J))
 
C          Build and write header label above the plot area
      CALL SET(XLPLT,XRPLT,YTTIC,YTOP,0.,1.,0.,1., 1)
C          Write main label first
      NC = LEN(HDR)
      L  = MIN(80,NC)
      DO 100 I=L,1,-1
      IF(HDR(I:I) .NE. ' ')GO TO 101
  100 NC = NC - 1
  101 IF (NC .GT. 0) CALL WTSTR(0.5,0.60,HDR(1:NC),IBIGC,0,0)
 
C          Write Y-axis label:
      CALL SET(XLFT,XLTIC,YBPLT,YTPLT,0.,1.,0.,1., 1)
      NC = LEN(YLBL)
      L  = MIN(80,NC)
      DO 150 I=L,1,-1
      IF(YLBL(I:I) .NE. ' ')GO TO 151
  150 NC = NC - 1
  151 IF (NC .GT. 0) CALL WTSTR(0.5,0.5,YLBL(1:NC),ISMLC,90,0)
 
C          Y-axis ticks will bound YMIN and YMAX such that nice
C          labels are possible.
      CALL GAXTIC(MSGUN,YMIN,YMAX,MXNMJRY,NMJRY,PMJRY,LBLY,LLBLY,NMNRY)
 
C          Write Y-axis tick marks to left of plot area
      CALL SET(XLTIC,XLPLT,YBPLT,YTPLT,0.,1.,0.,1., 1)
C          Establish fractional positions of:
      SCLX = 1. / ( RES * (XLPLT-XLTIC) )
C          PRED = right edge of all Y-axis ticks
C          PLTL = labelled tick left edge
C          PUTL = unlabelled tick left edge
C          PLBR = label right edge
      PRED = 1.
      PLTL = PRED - BIGTL * SCLX
      PUTL = PRED - SMLTL * SCLX
      PLBR = PLTL - GAP   * SCLX
 
C          Write any unlabelled ticks before first label on the Y-axis
      DMJRY = PMJRY(2) - PMJRY(1)
      IF (NMJRY .EQ. 1) DMJRY = PMJRY(1)
      DMNRY = DMJRY / REAL(NMNRY + 1)
      NBOT  = INT( PMJRY(1)/DMNRY )
      PMNRY0= PMJRY(1) - REAL(NBOT+1) * DMNRY
      PMNRY = PMNRY0
      DO 200 I=1,NBOT
      PMNRY = PMNRY + DMNRY
  200 CALL LINE (PUTL,PMNRY,PRED,PMNRY)
 
C          Write labelled tick marks and labels for the Y-axis
      DO 400 I=1,NMJRY
      CALL LINE (PLTL,PMJRY(I),PRED,PMJRY(I))
      CALL WTSTR(PLBR,PMJRY(I),LBLY(I)(1:LLBLY(I)),ISMLC,0,1)
C          Write unlabelled tick marks which are between labelled ticks:
      IF (I .LT. NMJRY) THEN
        PMNRY = PMJRY(I)
        DO 300 J=1,NMNRY
        PMNRY = PMNRY + DMNRY
  300   CALL LINE (PUTL,PMNRY,PRED,PMNRY)
      ENDIF
  400 CONTINUE
 
C          Write any unlabelled ticks after the last label on the Y-axis:
      NTOP = INT( (1. - PMJRY(NMJRY)) / DMNRY )
      IF (PMJRY(NMJRY)+REAL(NTOP+1)*DMNRY .LE. 1.00001) NTOP = NTOP+1
      PMNRY = PMJRY(NMJRY)
      DO 500 I=1,NTOP
      PMNRY = PMNRY + DMNRY
      CALL LINE (PUTL,PMNRY,PRED,PMNRY)
  500 CONTINUE
 
C          Write Y-axis tick marks to right of plot area - omit labels
      CALL SET(XRPLT,XRTIC,YBPLT,YTPLT,0.,1.,0.,1., 1)
C          Establish fractional positions of:
      SCLX = 1. / ( RES * (XRTIC-XRPLT) )
C          PLED = left edge of right side Y tick marks
C          PLTR = labelled tick right edge
C          PUTR = unlabelled tick right edge
      PLED = 0.
      PLTR = BIGTL * SCLX
      PUTR = SMLTL * SCLX
 
C          Write any unlabelled ticks before first label
      PMNRY = PMNRY0
      DO 510 I=1,NBOT
      PMNRY = PMNRY + DMNRY
  510 CALL LINE (PLED,PMNRY,PUTR,PMNRY)
 
C          Write labelled tick marks and labels for the Y-axis - omit lb
      DO 530 I=1,NMJRY
      CALL LINE (PLED,PMJRY(I),PLTR,PMJRY(I))
C          Write unlabelled tick marks which are between labelled ticks:
      IF (I .LT. NMJRY) THEN
        PMNRY = PMJRY(I)
        DO 520 J=1,NMNRY
        PMNRY = PMNRY + DMNRY
  520   CALL LINE (PLED,PMNRY,PUTR,PMNRY)
      ENDIF
  530 CONTINUE
 
C          Write any unlabelled ticks after the last label on the Y-axis:
      PMNRY = PMJRY(NMJRY)
      DO 540 I=1,NTOP
      PMNRY = PMNRY + DMNRY
      CALL LINE (PLED,PMNRY,PUTR,PMNRY)
  540 CONTINUE
 
C          Build and write X-axis label:
      CALL SET(XLPLT,XRPLT,YBOT,YBTIC,0.,1.,0.,1., 1)
C          Truncate UT begin time down to next even minute
      DO 700 I=1,3
  700 IRNDT(I) = ITMIN(I)
      IRNDT(4) = 0
C          Truncate UT end time up to next even minute
      DO 800 I=5,8
  800 IRNDT(I) = ITMAX(I-4)
      IF (IRNDT(8) .GT. 0) THEN
        CALL NEWTIM(ITMAX,IRNDT(5),I1MIN,MSGUN)
        IRNDT(8) = 0000
      ENDIF
C          If begin time is 24:00, express it as 00:00 (of next day)
      IF (IRNDT(3) .EQ. 2400) THEN
        DO 850 I=1,4
  850   ITMP(I) = IRNDT(I)
        CALL NEWTIM(ITMP,IRNDT(1),I1MIN,MSGUN)
        IRNDT(3) = 2400
      ENDIF
C          If end time is 00:00, express it as 24:00 (of previous day)
      IF (IRNDT(7) .EQ. 0) THEN
        DO 900 I=1,4
  900   ITMP(I) = IRNDT(I+4)
        CALL NEWTIM(ITMP,IRNDT(5),-I1MIN,MSGUN)
        IRNDT(7) = 2400
      ENDIF
      MONB = IRNDT(2)/100
      IDAB = IRNDT(2) - MONB*100
      IHRB = IRNDT(3)/100
      MINB = IRNDT(3) - IHRB*100
      MONE = IRNDT(6)/100
      IDAE = IRNDT(6) - MONE*100
      IHRE = IRNDT(7)/100
      MINE = IRNDT(7) - IHRE*100
      WRITE(HDRT,1100)IRNDT(1),MON(MONB),IDAB,IHRB,MINB,
     +                IRNDT(5),MON(MONE),IDAE,IHRE,MINE
 1100 FORMAT(I4,1X,A3,2I3,':',I2,' to ',I4,1X,A3,2I3,':',I2,' UT')
      IF (HDRT(13:13) .EQ. ' ') HDRT(13:13) = '0'
      IF (HDRT(16:16) .EQ. ' ') HDRT(16:16) = '0'
      IF (HDRT(34:34) .EQ. ' ') HDRT(34:34) = '0'
      IF (HDRT(37:37) .EQ. ' ') HDRT(37:37) = '0'
      CALL WTSTR(0.5,0.5,HDRT(1:41),ISMLC,0,0)
 
C          Determine appropriate X-axis tick marks
      CALL GTMTIC(MSGUN,IRNDT,MXNMJRT,NMJRT,PMJRT,LBLT,NLBLTC,NMNRT)
 
C          Write X-axis tick marks below plot area:
      CALL SET(XLPLT,XRPLT,YBTIC,YBPLT,0.,1.,0.,1., 1)
C          Establish fractional positions of:
      SCLY = 1. / ( RES * (YBPLT-YBTIC) )
      SCLX = 1. / ( RES * (XRPLT-XLPLT) )
C          PTED = top edge of all ticks
C          PLTB = lbld tick bottom edge
C          PUTB = unlbld tick bottom edge
C          PLCY = label center (vertically)
      PTED = 1.
      PLTB = PTED - BIGTL * SCLY
      PUTB = PTED - SMLTL * SCLY
      PLCY = PLTB - (GAP+0.5*HSML)*SCLY
 
C          Write any unlabelled ticks before first label on the X-axis
      IF (NMJRT .LE. 1) GO TO 1201
      DMJRT = PMJRT(2) - PMJRT(1)
      DMNRT = DMJRT / FLOAT(NMNRT + 1)
      NLFT  = IFIX( PMJRT(1)/DMNRT )
      PMNRT0= PMJRT(1) - FLOAT(NLFT+1) * DMNRT
      PMNRT = PMNRT0
      DO 1200 I=1,NLFT
      PMNRT = PMNRT + DMNRT
 1200 CALL LINE (PMNRT,PTED,PMNRT,PUTB)
 1201 CONTINUE
 
C          Write labelled tick marks and labels for the X-axis
      DO 1400 I=1,NMJRT
      CALL LINE (PMJRT(I),PTED,PMJRT(I),PLTB)
C          If lbl has only one character must shift to left
      OFSET = 0.
      IF (LBLT(I)(1:1) .EQ.  ' ') OFSET = 0.5*WSML * SCLX
      IF (LBLT(I)(1:2) .EQ. '  ') OFSET =     WSML * SCLX
      PLCX = PMJRT(I) - OFSET
      CALL WTSTR(PLCX,PLCY,LBLT(I)(1:NLBLTC),ISMLC,0,0)
C          Write unlabelled tick marks which are between labelled ticks:
      IF (I .LT. NMJRT) THEN
        PMNRT = PMJRT(I)
        DO 1300 J=1,NMNRT
        PMNRT = PMNRT + DMNRT
 1300   CALL LINE (PMNRT,PTED,PMNRT,PUTB)
      ENDIF
 1400 CONTINUE
 
C          Write any unlabelled ticks after the last label on the X-axis:
      IF (NMJRT .LE. 1) GO TO 1501
      NRHT  = INT( (1. - PMJRT(NMJRT)) / DMNRT )
      IF (PMJRT(NMJRT)+REAL(NRHT+1)*DMNRT .LE. 1.00001) NRHT = NRHT+1
      PMNRT = PMJRT(NMJRT)
      DO 1500 I=1,NRHT
      PMNRT = PMNRT + DMNRT
      CALL LINE (PMNRT,PTED,PMNRT,PUTB)
 1500 CONTINUE
 1501 CONTINUE
 
C          Write time X-axis ticks on top of plot area
C          (omitting labeling the ticks on the top of the plot area)
      CALL SET(XLPLT,XRPLT,YTPLT,YTTIC,0.,1.,0.,1., 1)
C          Establish fractional positions of:
      SCLY = 1. / ( RES * (YTTIC-YTPLT) )
      SCLX = 1. / ( RES * (XRPLT-XLPLT) )
C          PBED = bot edge of all ticks
C          PLTT = lbld tick top edge
C          PUTT = unlbld tick top edge
C          PLCY = label center (vertically)
      PBED = 0.
      PLTT = BIGTL * SCLY
      PUTT = SMLTL * SCLY
      PLCY = PLTT + (GAP+0.5*HSML)*SCLY
 
C          Write any unlabelled ticks before first label on the X-axis
C          (omitting labeling the ticks on the top of the plot area)
      IF (NMJRT .LE. 1) GO TO 1601
      PMNRT = PMNRT0
      DO 1600 I=1,NLFT
      PMNRT = PMNRT + DMNRT
 1600 CALL LINE (PMNRT,PBED,PMNRT,PUTT)
 1601 CONTINUE
 
C          Write labelled tick marks and labels for the X-axis
C          (omitting labeling the ticks on the top of the plot area)
      DO 1800 I=1,NMJRT
      CALL LINE (PMJRT(I),PBED,PMJRT(I),PLTT)
C          Write unlabelled tick marks which are between labelled ticks:
C          (omit labeling the ticks on the top of the plot area)
      IF (I .LT. NMJRT) THEN
        PMNRT = PMJRT(I)
        DO 1700 J=1,NMNRT
        PMNRT = PMNRT + DMNRT
 1700   CALL LINE (PMNRT,PBED,PMNRT,PUTT)
      ENDIF
 1800 CONTINUE
 
C          Write any unlabelled ticks after the last label on the X-axis:
C          (omitting labeling the ticks on the top of the plot area)
      IF (NMJRT .LE. 1) GO TO 1901
      PMNRT = PMJRT(NMJRT)
      DO 1900 I=1,NRHT
      PMNRT = PMNRT + DMNRT
      CALL LINE (PMNRT,PBED,PMNRT,PUTT)
 1900 CONTINUE
 1901 CONTINUE
 
C          Plot Y vs X.  While doing so, it is necessary to convert
C          time units from 4 integers to a single real
      YDIF = 0.
      CALL CVTI2R(MSGUN,DBEGT,IRNDT(1))
      CALL CVTI2R(MSGUN,DENDT,IRNDT(5))
      BEGT = SNGL (DBEGT)
      ENDT = SNGL (DENDT)
      IB1000 = INT (BEGT) / 1000
      IXOVR  = INT (ENDT) / 1000 - IB1000
      IF (IXOVR .GT. 0) THEN
C          A year has been crossed over.  Must change ENDT for set call
C          Check that never more than one year crossed over
        IF(IXOVR .GT. 1)STOP
        YDIF = 635.
        IF(MOD(IRNDT(1),4) .EQ. 0)YDIF = 634.
        ENDT = ENDT - YDIF
      ENDIF
      CALL SET(XLPLT,XRPLT,YBPLT,YTPLT,BEGT,ENDT,YMIN,YMAX,1)
C          Draw line(s) connecting adjacent non-missing values; if a
C          good value is isolated, just draw a plus
      J = 0
 2100 J = J+1
      I = J-1
      K = J+1
      IF (J .GT. NY) GO TO 2200
      IF (INT(Y(J)) .EQ. MISS) GO TO 2100
      IF (NY .EQ. 1) GO TO 2180
      IF (NY .GT. 2) GO TO 2150
      IF (J .EQ.  2) GO TO 2180
      IF (INT(Y(2)) .EQ. MISS) GO TO 2180
 2150 CONTINUE
      IF (J .EQ.  1 .AND. INT(Y(K)) .EQ. MISS) GO TO 2180
      IF (J .EQ. NY .AND. INT(Y(I)) .EQ. MISS) GO TO 2180
      IF (I .GT.  0 .AND. INT(Y(I)) .EQ. MISS .AND.
     +    K .LE. NY .AND. INT(Y(K)) .EQ. MISS) GO TO 2180
      IF (K .LE. NY .AND. INT(Y(K)) .EQ. MISS) GO TO 2100
      CALL CVTI2R(MSGUN,DTJ,IT(1,J))
      CALL CVTI2R(MSGUN,DTK,IT(1,K))
      TJ = SNGL (DTJ)
      TK = SNGL (DTK)
      IF (INT(TJ)/1000 .GT. IB1000) TJ = TJ - YDIF
      IF (INT(TK)/1000 .GT. IB1000) TK = TK - YDIF
      CALL LINE( TJ,Y(J) , TK,Y(K) )
      GO TO 2190
 2180 CALL CVTI2R(MSGUN,DTJ,IT(1,J))
      TJ = SNGL (DTJ)
      IF (INT(TJ)/1000 .GT. IB1000) TJ = TJ - YDIF
      CALL POINTS( TJ,Y(J),1,-2,0)
 2190 IF (K .LT. NY) GO TO 2100
 2200 CONTINUE
 
C          Draw vertical lines delineating multiple days on the plot if
C          there are more than 1 day changes.  Build times of the day
C          changes and count how many there are, before deciding to draw.
C          Note that day changes at the left and right side of the plot
C          are excluded from this count.
      NDC = 0
      ITMP(1) = IRNDT(1)
      ITMP(2) = IRNDT(2)
      ITMP(3) = 2400
      ITMP(4) = 0
 2210 NDC = NDC + 1
      IF (NDC .GT. MXNDCH) GO TO 9100
      CALL CVTI2R(MSGUN,DTDC,ITMP)
      TDC(NDC) = SNGL (DTDC)
      IF(INT(TDC(NDC))/1000 .GT. IB1000)TDC(NDC) = TDC(NDC) - YDIF
      DO 2220 I=1,4
 2220 ITP2(I) = ITMP(I)
      CALL NEWTIM(ITP2,ITMP,I1DA,MSGUN)
      CALL TIMDIF(ITMP,IRNDT(5),IDIF)
      IF (IDIF .GT. 0D0) GO TO 2210
      IF (NDC .GT. 1) THEN
        DO 2230 I=1,NDC
 2230   CALL LINE( TDC(I),YMIN , TDC(I),YMAX )
      ENDIF
 
C          If the range of Y values includes 0., draw the zero line
      IF (YMIN .LT. 0. .AND. YMAX .GT. 0.) CALL LINE(BEGT,0. , ENDT,0.)
 
C          Draw a perimeter around the plot area
      CALL SET(XLPLT,XRPLT,YBPLT,YTPLT,0.,1.,0.,1.,1)
      CALL LINE(0.,0.,0.,1.)
      CALL LINE(0.,1.,1.,1.)
      CALL LINE(1.,1.,1.,0.)
      CALL LINE(1.,0.,0.,0.)
 
C          Advance frame
      CALL FRAME
      RETURN
 
C          Error trap diagnostics
 9100 WRITE(MSGUN,'('' Dimension exceeded.  Increase MXNDCH!'')')
      STOP
      END
"EOF"
 
cft77 -e mx -l cmpl src$$.f || rcp cmpl ${ofl}cmpl
if ($status != 0) exit
 
#       Link to standard libraries including ncargraphics
segldr -o exe$$ src$$.o -L /lib,/usr/lib,/usr/local/lib \
          -l mss,ncarm,ncaro,net,ncarg,ncarg_gks,ncarg_c,ncarg_loc,X11
 
exe$$ < in$$  > out$$
 
if (! -z out$$) then
  echo "            Returning listing as $ofl:t"
  rcp out$$ $ofl
endif
 
if (-e gmeta) then
  if (! -z gmeta) then
    echo "            Returning plot as $ofp:t"
    rcp gmeta $ofp
  endif
endif
 
if (-e fort.61) then
  if (! -z fort.61) then
    echo "            Returning lexical output as $ofx:t"
    rcp fort.61 $ofx
  endif
endif
 
rm *$$*
ja -s
date
