      subroutine input(filnam,frm,nblock,nskip,npts,rdata)
      character*(*) filnam,frm
      parameter (nrecl=30000)
      dimension rdata(*),recrd(nrecl)
      if(((frm.eq.'U').or.(frm.eq.'u')).and.(nblock.eq.0))goto 100
      if((frm.eq.'F').or.(frm.eq.'f'))goto 200
      open (unit=10,file=filnam,status='old',form='unformatted')
      nrecs=nskip/nblock
      nlefs=nskip-nblock*nrecs
      do 10 i=1,nrecs
   10   read(10,end=11)(recrd(j),j=1,nblock)
      goto 12
   11 npts=0
      return
   12 k=1
      if(nlefs.eq.0)goto 30
      read(10,end=300)(recrd(i),i=1,nblock)
      do 20 i=nlefs+1,min(nblock,npts+nlefs)
        rdata(k)=recrd(i)
   20   k=k+1
   30 continue
      nread=npts-k+1
      nrec=nread/nblock
      nleft=nread-nblock*nrec
      do 40 j=1,nrec
        read(10,end=300)(recrd(i),i=1,nblock)
        do 35 i=1,nblock
          rdata(k)=recrd(i)
   35     k=k+1
   40 continue
      if(nleft.eq.0)goto 300
      read(10,end=300)(recrd(i),i=1,nblock)
      do 50 i=1,nleft
        rdata(k)=recrd(i)
   50   k=k+1
      goto 300
  100 continue
      open (unit=10,file=filnam,status='old',form='unformatted')
      do 110 i=1,nskip
  110   read(10,end=111)
      goto 112
  111 npts=0
      return
  112 continue
      do 120 k=1,npts
  120   read(10,end=300)rdata(k)
      goto 300
  200 continue
      open (unit=10,file=filnam,status='old')
      do 210 i=1,nskip
  210   read(10,*,end=211)
      goto 212
  211 npts=0
      return
  212 continue
      do 220 k=1,npts
  220   read(10,*,end=300)rdata(k)
  300 npts=k-1
      close (unit=10)
      return
      end
      subroutine binput(filnam,imode,nskip,nread,nar,ncomp,nser,
     $           iindex,datum,rdata,ar,cvar)
c imode=1 means no header except possibly ar filter coefficients
c imode=2 means ascii ts format so that data read begins after
c '>INFO_END' is detected and cvar is returned with
c latitude, longitude, and elevation (if present)
      include 'parameters.h'
      parameter (narm=100)
      character*(*) filnam,cvar(3)
      character line*80
      dimension rdata(nptsm,*),datum(*),iindex(*),
     $          ar(narm+1,*)
      open (unit=10,file=filnam,status='old')
      if(nar.lt.0)then
        do 10 n=1,10000
          read(10,'(a)')line
          if(line(1:10).eq.'#AR_FILTER')then
            do 6 j=1,nser
              ar(1,iindex(j))=1.
              do 5 jj=1,-nar/5
    5         read(10,'(1x,5g15.7)')(ar((jj-1)*5+i,
     $                               iindex(j)),i=1,5)
              nl=-nar-5*(-nar/5)
              if(nl.ne.0)then
                read(10,'(1x,5g15.7)')
     $              (ar(5*(-nar/5)+i,iindex(j)),i=1,nl)
              endif
    6       continue
            goto 11
          endif
   10     continue
      endif
   11 if(imode.eq.1)goto 50
      do 15 i=1,3
   15   cvar(i)=' '
      do 20 i=1,10000
        read(10,'(a)',end=90)line
   20   if(line(1:1).eq.'>')goto 22
   22 do 30 i=1,10000
        if(line(2:9).eq.'INFO_END')goto 50
        if(line(2:9).eq.'LATITUDE')then
          i1=index(line,':')
          i2=index(line,'=')
          if(i1.ne.0)then
            call strlen(line(i1+1:80),j1,j2)
            cvar(1)=line(i1+j1-1:i1+j2+1)
          elseif(i2.ne.0)then
            call strlen(line(i2+1:80),j1,j2)
            cvar(1)=line(i2+j1-1:i2+j2+1)
          endif
        elseif(line(2:10).eq.'LONGITUDE')then
          i1=index(line,':')
          i2=index(line,'=')
          if(i1.ne.0)then
            call strlen(line(i1+1:80),j1,j2)
            cvar(2)=line(i1+j1-1:i1+j2+1)
          elseif(i2.ne.0)then
            call strlen(line(i2+1:80),j1,j2)
            cvar(2)=line(i2+j1-1:i2+j2+1)
          endif
        elseif(line(2:10).eq.'ELEVATION')then
          i1=index(line,':')
          i2=index(line,'=')
          if(i1.ne.0)then
            call strlen(line(i1+1:80),j1,j2)
            cvar(3)=line(i1+j1-1:i1+j2+1)
          elseif(i2.ne.0)then
            call strlen(line(i2+1:80),j1,j2)
            cvar(3)=line(i2+j1-1:i2+j2+1)
          endif
        endif
   30   read(10,'(a)')line
      nread=0
      return
   50 do 70 i=1,nskip
   70   read(10,*,end=90)
      do 80 i=1,nread
        read(10,*,end=81)(datum(j),j=1,ncomp)
        do 80 j=1,nser
   80     rdata(i,j)=datum(iindex(j))
      close (unit=10)
      return
   81 nread=i-1
      close (unit=10)
      return
   90 nread=0
      close (unit=10)
      return
      end
      integer function i1mach(i)
      integer i
c
c    i1mach( 1) = the standard input unit.
c    i1mach( 2) = the standard output unit.
c    i1mach( 3) = the standard punch unit.
c    i1mach( 4) = the standard error message unit.
c    i1mach( 5) = the number of bits per integer storage unit.
c    i1mach( 6) = the number of characters per character storage unit.
c    integers have form sign ( x(s-1)*a**(s-1) + ... + x(1)*a + x(0) )
c    i1mach( 7) = a, the base.
c    i1mach( 8) = s, the number of base-a digits.
c    i1mach( 9) = a**s - 1, the largest magnitude.
c    floats have form  sign (b**e)*( (x(1)/b) + ... + (x(t)/b**t) )
c               where  emin .le. e .le. emax.
c    i1mach(10) = b, the base.
c  single-precision
c    i1mach(11) = t, the number of base-b digits.
c    i1mach(12) = emin, the smallest exponent e.
c    i1mach(13) = emax, the largest exponent e.
c  double-precision
c    i1mach(14) = t, the number of base-b digits.
c    i1mach(15) = emin, the smallest exponent e.
c    i1mach(16) = emax, the largest exponent e.
c
      integer imach(16), output, sc, small(2)
      save imach, sc
      real rmach
      equivalence (imach(4),output), (rmach,small(1))
      integer i3, j, k, t3e(3)
      data t3e(1) / 9777664 /
      data t3e(2) / 5323660 /
      data t3e(3) / 46980 /
c  this version adapts automatically to most current machines,
c  including auto-double compilers.
c  to compile on older machines, add a c in column 1
c  on the next line
      data sc/0/
c  and remove the c from column 1 in one of the sections below.
c  constants for even older machines can be obtained by
c          mail netlib@research.bell-labs.com
c          send old1mach from blas
c  please send corrections to dmg or ehg@bell-labs.com.
c
c     machine constants for the honeywell dps 8/70 series.
c
c      data imach( 1) /    5 /
c      data imach( 2) /    6 /
c      data imach( 3) /   43 /
c      data imach( 4) /    6 /
c      data imach( 5) /   36 /
c      data imach( 6) /    4 /
c      data imach( 7) /    2 /
c      data imach( 8) /   35 /
c      data imach( 9) / o377777777777 /
c      data imach(10) /    2 /
c      data imach(11) /   27 /
c      data imach(12) / -127 /
c      data imach(13) /  127 /
c      data imach(14) /   63 /
c      data imach(15) / -127 /
c      data imach(16) /  127 /, sc/987/
c
c     machine constants for pdp-11 fortrans supporting
c     32-bit integer arithmetic.
c
c      data imach( 1) /    5 /
c      data imach( 2) /    6 /
c      data imach( 3) /    7 /
c      data imach( 4) /    6 /
c      data imach( 5) /   32 /
c      data imach( 6) /    4 /
c      data imach( 7) /    2 /
c      data imach( 8) /   31 /
c      data imach( 9) / 2147483647 /
c      data imach(10) /    2 /
c      data imach(11) /   24 /
c      data imach(12) / -127 /
c      data imach(13) /  127 /
c      data imach(14) /   56 /
c      data imach(15) / -127 /
c      data imach(16) /  127 /, sc/987/
c
c     machine constants for the univac 1100 series.
c
c     note that the punch unit, i1mach(3), has been set to 7
c     which is appropriate for the univac-for system.
c     if you have the univac-ftn system, set it to 1.
c
c      data imach( 1) /    5 /
c      data imach( 2) /    6 /
c      data imach( 3) /    7 /
c      data imach( 4) /    6 /
c      data imach( 5) /   36 /
c      data imach( 6) /    6 /
c      data imach( 7) /    2 /
c      data imach( 8) /   35 /
c      data imach( 9) / o377777777777 /
c      data imach(10) /    2 /
c      data imach(11) /   27 /
c      data imach(12) / -128 /
c      data imach(13) /  127 /
c      data imach(14) /   60 /
c      data imach(15) /-1024 /
c      data imach(16) / 1023 /, sc/987/
c
      if (sc .ne. 987) then
*        *** check for autodouble ***
         small(2) = 0
         rmach = 1e13
         if (small(2) .ne. 0) then
*           *** autodoubled ***
            if (      (small(1) .eq. 1117925532
     *           .and. small(2) .eq. -448790528)
     *       .or.     (small(2) .eq. 1117925532
     *           .and. small(1) .eq. -448790528)) then
*               *** ieee ***
               imach(10) = 2
               imach(14) = 53
               imach(15) = -1021
               imach(16) = 1024
            else if ( small(1) .eq. -2065213935
     *          .and. small(2) .eq. 10752) then
*               *** vax with d_floating ***
               imach(10) = 2
               imach(14) = 56
               imach(15) = -127
               imach(16) = 127
            else if ( small(1) .eq. 1267827943
     *          .and. small(2) .eq. 704643072) then
*               *** ibm mainframe ***
               imach(10) = 16
               imach(14) = 14
               imach(15) = -64
               imach(16) = 63
            else
               write(*,9010)
               stop 777
               end if
            imach(11) = imach(14)
            imach(12) = imach(15)
            imach(13) = imach(16)
         else
            rmach = 1234567.
            if (small(1) .eq. 1234613304) then
*               *** ieee ***
               imach(10) = 2
               imach(11) = 24
               imach(12) = -125
               imach(13) = 128
               imach(14) = 53
               imach(15) = -1021
               imach(16) = 1024
               sc = 987
            else if (small(1) .eq. -1271379306) then
*               *** vax ***
               imach(10) = 2
               imach(11) = 24
               imach(12) = -127
               imach(13) = 127
               imach(14) = 56
               imach(15) = -127
               imach(16) = 127
               sc = 987
            else if (small(1) .eq. 1175639687) then
*               *** ibm mainframe ***
               imach(10) = 16
               imach(11) = 6
               imach(12) = -64
               imach(13) = 63
               imach(14) = 14
               imach(15) = -64
               imach(16) = 63
               sc = 987
            else if (small(1) .eq. 1251390520) then
*              *** convex c-1 ***
               imach(10) = 2
               imach(11) = 24
               imach(12) = -128
               imach(13) = 127
               imach(14) = 53
               imach(15) = -1024
               imach(16) = 1023
            else
               do 10 i3 = 1, 3
                  j = small(1) / 10000000
                  k = small(1) - 10000000*j
                  if (k .ne. t3e(i3)) go to 20
                  small(1) = j
 10               continue
*              *** cray t3e ***
               imach( 1) = 5
               imach( 2) = 6
               imach( 3) = 0
               imach( 4) = 0
               imach( 5) = 64
               imach( 6) = 8
               imach( 7) = 2
               imach( 8) = 63
               call i1mcr1(imach(9), k, 32767, 16777215, 16777215)
               imach(10) = 2
               imach(11) = 53
               imach(12) = -1021
               imach(13) = 1024
               imach(14) = 53
               imach(15) = -1021
               imach(16) = 1024
               go to 35
 20            call i1mcr1(j, k, 16405, 9876536, 0)
               if (small(1) .ne. j) then
                  write(*,9020)
                  stop 777
                  end if
*              *** cray 1, xmp, 2, and 3 ***
               imach(1) = 5
               imach(2) = 6
               imach(3) = 102
               imach(4) = 6
               imach(5) = 46
               imach(6) = 8
               imach(7) = 2
               imach(8) = 45
               call i1mcr1(imach(9), k, 0, 4194303, 16777215)
               imach(10) = 2
               imach(11) = 47
               imach(12) = -8188
               imach(13) = 8189
               imach(14) = 94
               imach(15) = -8141
               imach(16) = 8189
               go to 35
               end if
            end if
         imach( 1) = 5
         imach( 2) = 6
         imach( 3) = 7
         imach( 4) = 6
         imach( 5) = 32
         imach( 6) = 4
         imach( 7) = 2
         imach( 8) = 31
         imach( 9) = 2147483647
 35      sc = 987
         end if
 9010 format(/' adjust autodoubled i1mach by uncommenting data'/
     * ' statements appropriate for your machine and setting'/
     * ' imach(i) = imach(i+3) for i = 11, 12, and 13.')
 9020 format(/' adjust i1mach by uncommenting data statements'/
     * ' appropriate for your machine.')
      if (i .lt. 1  .or.  i .gt. 16) go to 40
      i1mach = imach(i)
      return
 40   write(*,*) 'i1mach(i): i =',i,' is out of bounds.'
      stop
      end
      subroutine i1mcr1(a, a1, b, c, d)
**** special computation for old cray machines ****
      integer a, a1, b, c, d
      a1 = 16777216*b + c
      a = 16777216*a1 + d
      end
      real function r1mach(i)
      integer i
c
c  single-precision machine constants
c  r1mach(1) = b**(emin-1), the smallest positive magnitude.
c  r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude.
c  r1mach(3) = b**(-t), the smallest relative spacing.
c  r1mach(4) = b**(1-t), the largest relative spacing.
c  r1mach(5) = log10(b)
c
      integer small(2)
      integer large(2)
      integer right(2)
      integer diver(2)
      integer log10(2)
c     needs to be (2) for autodouble, harris slash 6, ...
      integer sc
      save small, large, right, diver, log10, sc
      real rmach(5)
      equivalence (rmach(1),small(1))
      equivalence (rmach(2),large(1))
      equivalence (rmach(3),right(1))
      equivalence (rmach(4),diver(1))
      equivalence (rmach(5),log10(1))
      integer j, k, l, t3e(3)
      data t3e(1) / 9777664 /
      data t3e(2) / 5323660 /
      data t3e(3) / 46980 /
c  this version adapts automatically to most current machines,
c  including auto-double compilers.
c  to compile on older machines, add a c in column 1
c  on the next line
      data sc/0/
c  and remove the c from column 1 in one of the sections below.
c  constants for even older machines can be obtained by
c          mail netlib@research.bell-labs.com
c          send old1mach from blas
c  please send corrections to dmg or ehg@bell-labs.com.
c
c     machine constants for the honeywell dps 8/70 series.
c      data rmach(1) / o402400000000 /
c      data rmach(2) / o376777777777 /
c      data rmach(3) / o714400000000 /
c      data rmach(4) / o716400000000 /
c      data rmach(5) / o776464202324 /, sc/987/
c
c     machine constants for pdp-11 fortrans supporting
c     32-bit integers (expressed in integer and octal).
c      data small(1) /    8388608 /
c      data large(1) / 2147483647 /
c      data right(1) /  880803840 /
c      data diver(1) /  889192448 /
c      data log10(1) / 1067065499 /, sc/987/
c      data rmach(1) / o00040000000 /
c      data rmach(2) / o17777777777 /
c      data rmach(3) / o06440000000 /
c      data rmach(4) / o06500000000 /
c      data rmach(5) / o07746420233 /, sc/987/
c
c     machine constants for the univac 1100 series.
c      data rmach(1) / o000400000000 /
c      data rmach(2) / o377777777777 /
c      data rmach(3) / o146400000000 /
c      data rmach(4) / o147400000000 /
c      data rmach(5) / o177464202324 /, sc/987/
c
      if (sc .ne. 987) then
*        *** check for autodouble ***
         small(2) = 0
         rmach(1) = 1e13
         if (small(2) .ne. 0) then
*           *** autodoubled ***
            if (      small(1) .eq. 1117925532
     *          .and. small(2) .eq. -448790528) then
*              *** ieee big endian ***
               small(1) = 1048576
               small(2) = 0
               large(1) = 2146435071
               large(2) = -1
               right(1) = 1017118720
               right(2) = 0
               diver(1) = 1018167296
               diver(2) = 0
               log10(1) = 1070810131
               log10(2) = 1352628735
            else if ( small(2) .eq. 1117925532
     *          .and. small(1) .eq. -448790528) then
*              *** ieee little endian ***
               small(2) = 1048576
               small(1) = 0
               large(2) = 2146435071
               large(1) = -1
               right(2) = 1017118720
               right(1) = 0
               diver(2) = 1018167296
               diver(1) = 0
               log10(2) = 1070810131
               log10(1) = 1352628735
            else if ( small(1) .eq. -2065213935
     *          .and. small(2) .eq. 10752) then
*              *** vax with d_floating ***
               small(1) = 128
               small(2) = 0
               large(1) = -32769
               large(2) = -1
               right(1) = 9344
               right(2) = 0
               diver(1) = 9472
               diver(2) = 0
               log10(1) = 546979738
               log10(2) = -805796613
            else if ( small(1) .eq. 1267827943
     *          .and. small(2) .eq. 704643072) then
*              *** ibm mainframe ***
               small(1) = 1048576
               small(2) = 0
               large(1) = 2147483647
               large(2) = -1
               right(1) = 856686592
               right(2) = 0
               diver(1) = 873463808
               diver(2) = 0
               log10(1) = 1091781651
               log10(2) = 1352628735
            else
               write(*,9010)
               stop 777
               end if
         else
            rmach(1) = 1234567.
            if (small(1) .eq. 1234613304) then
*              *** ieee ***
               small(1) = 8388608
               large(1) = 2139095039
               right(1) = 864026624
               diver(1) = 872415232
               log10(1) = 1050288283
            else if (small(1) .eq. -1271379306) then
*              *** vax ***
               small(1) = 128
               large(1) = -32769
               right(1) = 13440
               diver(1) = 13568
               log10(1) = 547045274
            else if (small(1) .eq. 1175639687) then
*              *** ibm mainframe ***
               small(1) = 1048576
               large(1) = 2147483647
               right(1) = 990904320
               diver(1) = 1007681536
               log10(1) = 1091781651
            else if (small(1) .eq. 1251390520) then
*              *** convex c-1 ***
               small(1) = 8388608
               large(1) = 2147483647
               right(1) = 880803840
               diver(1) = 889192448
               log10(1) = 1067065499
            else
               do 10 l = 1, 3
                  j = small(1) / 10000000
                  k = small(1) - 10000000*j
                  if (k .ne. t3e(l)) go to 20
                  small(1) = j
 10               continue
*              *** cray t3e ***
               call i1mcra(small(1), k, 16, 0, 0)
               call i1mcra(large(1), k, 32751, 16777215, 16777215)
               call i1mcra(right(1), k, 15520, 0, 0)
               call i1mcra(diver(1), k, 15536, 0, 0)
               call i1mcra(log10(1), k, 16339, 4461392, 10451455)
               go to 30
 20            call i1mcra(j, k, 16405, 9876536, 0)
               if (small(1) .ne. j) then
                  write(*,9020)
                  stop 777
                  end if
*              *** cray 1, xmp, 2, and 3 ***
               call i1mcra(small(1), k, 8195, 8388608, 1)
               call i1mcra(large(1), k, 24574, 16777215, 16777214)
               call i1mcra(right(1), k, 16338, 8388608, 0)
               call i1mcra(diver(1), k, 16339, 8388608, 0)
               call i1mcra(log10(1), k, 16383, 10100890, 8715216)
               end if
            end if
 30      sc = 987
         end if
*     sanity check
      if (rmach(4) .ge. 1.0) stop 776
      if (i .lt. 1 .or. i .gt. 5) then
         write(*,*) 'r1mach(i): i =',i,' is out of bounds.'
         stop
         end if
      r1mach = rmach(i)
      return
 9010 format(/' adjust autodoubled r1mach by getting data'/
     *' appropriate for your machine from d1mach.')
 9020 format(/' adjust r1mach by uncommenting data statements'/
     *' appropriate for your machine.')
      end

      subroutine i1mcra(a, a1, b, c, d)
**** special computation for cray machines ****
      integer a, a1, b, c, d
      a1 = 16777216*b + c
      a = 16777216*a1 + d
      end

      double precision function d1mach(i)
      integer i
c
c  double-precision machine constants
c  d1mach( 1) = b**(emin-1), the smallest positive magnitude.
c  d1mach( 2) = b**emax*(1 - b**(-t)), the largest magnitude.
c  d1mach( 3) = b**(-t), the smallest relative spacing.
c  d1mach( 4) = b**(1-t), the largest relative spacing.
c  d1mach( 5) = log10(b)
c
      integer small(2)
      integer large(2)
      integer right(2)
      integer diver(2)
      integer log10(2)
      integer sc, cray1(38), j
      common /d9mach/ cray1
      save small, large, right, diver, log10, sc
      double precision dmach(5)
      equivalence (dmach(1),small(1))
      equivalence (dmach(2),large(1))
      equivalence (dmach(3),right(1))
      equivalence (dmach(4),diver(1))
      equivalence (dmach(5),log10(1))
c  this version adapts automatically to most current machines.
c  r1mach can handle auto-double compiling, but this version of
c  d1mach does not, because we do not have quad constants for
c  many machines yet.
c  to compile on older machines, add a c in column 1
c  on the next line
      data sc/0/
c  and remove the c from column 1 in one of the sections below.
c  constants for even older machines can be obtained by
c          mail netlib@research.bell-labs.com
c          send old1mach from blas
c  please send corrections to dmg or ehg@bell-labs.com.
c
c     machine constants for the honeywell dps 8/70 series.
c      data small(1),small(2) / o402400000000, o000000000000 /
c      data large(1),large(2) / o376777777777, o777777777777 /
c      data right(1),right(2) / o604400000000, o000000000000 /
c      data diver(1),diver(2) / o606400000000, o000000000000 /
c      data log10(1),log10(2) / o776464202324, o117571775714 /, sc/987/
c
c     machine constants for pdp-11 fortrans supporting
c     32-bit integers.
c      data small(1),small(2) /    8388608,           0 /
c      data large(1),large(2) / 2147483647,          -1 /
c      data right(1),right(2) /  612368384,           0 /
c      data diver(1),diver(2) /  620756992,           0 /
c      data log10(1),log10(2) / 1067065498, -2063872008 /, sc/987/
c
c     machine constants for the univac 1100 series.
c      data small(1),small(2) / o000040000000, o000000000000 /
c      data large(1),large(2) / o377777777777, o777777777777 /
c      data right(1),right(2) / o170540000000, o000000000000 /
c      data diver(1),diver(2) / o170640000000, o000000000000 /
c      data log10(1),log10(2) / o177746420232, o411757177572 /, sc/987/
c
c     on first call, if no data uncommented, test machine types.
      if (sc .ne. 987) then
         dmach(1) = 1.d13
         if (      small(1) .eq. 1117925532
     *       .and. small(2) .eq. -448790528) then
*           *** ieee big endian ***
            small(1) = 1048576
            small(2) = 0
            large(1) = 2146435071
            large(2) = -1
            right(1) = 1017118720
            right(2) = 0
            diver(1) = 1018167296
            diver(2) = 0
            log10(1) = 1070810131
            log10(2) = 1352628735
         else if ( small(2) .eq. 1117925532
     *       .and. small(1) .eq. -448790528) then
*           *** ieee little endian ***
            small(2) = 1048576
            small(1) = 0
            large(2) = 2146435071
            large(1) = -1
            right(2) = 1017118720
            right(1) = 0
            diver(2) = 1018167296
            diver(1) = 0
            log10(2) = 1070810131
            log10(1) = 1352628735
         else if ( small(1) .eq. -2065213935
     *       .and. small(2) .eq. 10752) then
*               *** vax with d_floating ***
            small(1) = 128
            small(2) = 0
            large(1) = -32769
            large(2) = -1
            right(1) = 9344
            right(2) = 0
            diver(1) = 9472
            diver(2) = 0
            log10(1) = 546979738
            log10(2) = -805796613
         else if ( small(1) .eq. 1267827943
     *       .and. small(2) .eq. 704643072) then
*               *** ibm mainframe ***
            small(1) = 1048576
            small(2) = 0
            large(1) = 2147483647
            large(2) = -1
            right(1) = 856686592
            right(2) = 0
            diver(1) = 873463808
            diver(2) = 0
            log10(1) = 1091781651
            log10(2) = 1352628735
         else if ( small(1) .eq. 1120022684
     *       .and. small(2) .eq. -448790528) then
*           *** convex c-1 ***
            small(1) = 1048576
            small(2) = 0
            large(1) = 2147483647
            large(2) = -1
            right(1) = 1019215872
            right(2) = 0
            diver(1) = 1020264448
            diver(2) = 0
            log10(1) = 1072907283
            log10(2) = 1352628735
         else if ( small(1) .eq. 815547074
     *       .and. small(2) .eq. 58688) then
*           *** vax g-floating ***
            small(1) = 16
            small(2) = 0
            large(1) = -32769
            large(2) = -1
            right(1) = 15552
            right(2) = 0
            diver(1) = 15568
            diver(2) = 0
            log10(1) = 1142112243
            log10(2) = 2046775455
         else
            dmach(2) = 1.d27 + 1
            dmach(3) = 1.d27
            large(2) = large(2) - right(2)
            if (large(2) .eq. 64 .and. small(2) .eq. 0) then
               cray1(1) = 67291416
               do 10 j = 1, 20
                  cray1(j+1) = cray1(j) + cray1(j)
 10               continue
               cray1(22) = cray1(21) + 321322
               do 20 j = 22, 37
                  cray1(j+1) = cray1(j) + cray1(j)
 20               continue
               if (cray1(38) .eq. small(1)) then
*                  *** cray ***
                  call i1mcry(small(1), j, 8285, 8388608, 0)
                  small(2) = 0
                  call i1mcry(large(1), j, 24574, 16777215, 16777215)
                  call i1mcry(large(2), j, 0, 16777215, 16777214)
                  call i1mcry(right(1), j, 16291, 8388608, 0)
                  right(2) = 0
                  call i1mcry(diver(1), j, 16292, 8388608, 0)
                  diver(2) = 0
                  call i1mcry(log10(1), j, 16383, 10100890, 8715215)
                  call i1mcry(log10(2), j, 0, 16226447, 9001388)
               else
                  write(*,9000)
                  stop 779
                  end if
            else
               write(*,9000)
               stop 779
               end if
            end if
         sc = 987
         end if
*    sanity check
      if (dmach(4) .ge. 1.0d0) stop 778
      if (i .lt. 1 .or. i .gt. 5) then
         write(*,*) 'd1mach(i): i =',i,' is out of bounds.'
         stop
         end if
      d1mach = dmach(i)
      return
 9000 format(/' adjust d1mach by uncommenting data statements'/
     *' appropriate for your machine.')
      end
      subroutine i1mcry(a, a1, b, c, d)
**** special computation for old cray machines ****
      integer a, a1, b, c, d
      a1 = 16777216*b + c
      a = 16777216*a1 + d
      end
      function phase(z)
      parameter(pi=3.1415926535,deg=180./pi)
      complex z
      z1=real(z)
      z2=aimag(z)
      if((z1.ne.0.).and.(z2.ne.0.))then
        phase=deg*atan2(z2,z1)
      else
        phase=0.
      endif
      return
      end
      logical function compare(x1,x2,n)
c  compares two floating point numbers and returns .true.
c  if their mantissas are the same to n significant figures
      compare=.false.
      j=int(log10(x1))
      if(j.ne.int(log10(x2)))return
      ix1=int(10.**(-j+n)*x1)
      ix2=int(10.**(-j+n)*x2)
      if(ix1.ne.ix2)return
      compare=.true.
      return
      end
      subroutine jfinit(ilev,nout,ninp,nref,nr3,nr2,nrr,
     $                  tbw,deltat,nfft,nsctinc,nsctmax,nf1,nfinc,
     $                  nfsect,mfft,uin,ainlin,ainuin,c2threshb,
     $                  c2threshe,nz,c2threshe1,perlo,perhi,nprej,
     $                  frej,npcs,nar,imode,jmode)
      logical ilev
      dimension frej(nprej)
      if(.not.ilev)then
        write(18,10)'#BIRRP Version 5 basic mode output'
      else
        write(18,10)'#BIRRP Version 5 advanced mode output'
      endif
      write(18,11)'#outputs=',nout,' inputs=',ninp,' references=',nref
      if(ilev)then
        write(18,11)'# nr3=',nr3,' nr2=',nr2,' nrr=',nrr
      endif
      write(18,12)'#tbw=',tbw,' deltat=',deltat
      write(18,11)'#nfft=',nfft,' nsctinc=',nsctinc,
     $           ' nsctmax=',nsctmax
      write(18,11)'#nf1=',nf1,' nfinc=',nfinc,' nfsect=',nfsect
      if(ilev)then
        write(18,12)'#mfft=',mfft
        write(18,12)'#uin=',uin,' ainuin=',ainuin
        write(18,12)'#c2threshb=',c2threshb,' c2threshe=',c2threshe
        write(18,12)'#perlo=',perlo,' perhi=',perhi
        write(18,11)'#nprej=',nprej
        write(18,10)'#periods'
        do 1 i=1,nprej
   1      write(18,13)'#',deltat/frej(i)
      else
        write(18,12)'#uin=',uin,' ainlin=',ainlin,' ainuin=',ainuin
        write(18,12)'#c2threshe=',c2threshe
      endif
      write(18,11)'#nz=',nz
      write(18,12)'#c2threshe1=',c2threshe1
      write(18,11)'# npcs=',npcs,' nar=',nar
      write(18,11)'#imode=',imode,' jmode=',jmode
   10 format(a)
   11 format(3(a,i8))
   12 format(3(a,g15.7))
   13 format(1x,g15.7)
      return
      end
      integer function tdiff( time1, time2 )
c
c     T D I F F     derives the time difference in seconds
c                   between two times given in format yyyy-mm-dd hh:mm:ss
c
c    Note that this is an INTEGER*4 function to cope with large
c    time differences (up to 68 years)
c
      character*19 time1, time2
      integer year1, month1, day1, hour1, min1, sec1
      integer year2, month2, day2, hour2, min2, sec2
      integer year, ddiff, hdiff, mdiff, sdiff, dayyr
      call date( 0, time1, year1, month1, day1, hour1, min1, sec1 )
      call date( 0, time2, year2, month2, day2, hour2, min2, sec2 )
      ddiff = dayyr(year2,month2,day2 ) + (year2-year1)*365 -
     &        dayyr(year1,month1,day1)
      if(((mod(year1,4).eq.0.and.mod(year1,100).ne.0)
     &    .or.(mod(year1,400).eq.0)).and.(day1.le.59))ddiff=ddiff+1
      hdiff = hour2 - hour1
      mdiff = min2  - min1
      sdiff = sec2 - sec1
      do 1 year = year1+1, year2-1
        if( (mod(year,4).eq.0  .and. mod(year,100).ne.0)
     &                .or. (mod(year,400).eq.0) ) then
          ddiff = ddiff + 1
        endif
    1 continue
      tdiff = ddiff*24*60*60  +  hdiff*60*60  + mdiff*60 + sdiff
      return
      end
      integer function dayyr(year,month,day)
      integer year, month, day, leap, i
      integer daytab(2,12)
      data daytab/31,31, 28,29, 31,31, 30,30, 31,31, 30,30,
     &              31,31, 31,31, 30,30, 31,31, 30,30, 31,31/
 
      if((mod(year,4).eq.0 .and.mod(year,100).ne.0)
     &              .or.(mod(year,400).eq.0)) then
        leap = 2
      else
        leap = 1
      endif
      dayyr = day
      do 1 i = 1, month-1
        dayyr = dayyr + daytab(leap,i)
   1  continue
      return
      end
      subroutine date( conv, cdate, year, month, day, hour, min, sec )
      character cdate*19
      integer conv, year, month, day, hour, min, sec
c
c     converts from date string format to integer format
c
c     conv = 1:  convert from 6 integers to character*19 string
c          = 0:  convert from character*19 string to 6 integers
      if( conv.eq.1) then
        cdate = '0000-00-00 00:00:00'
        write( cdate( 1: 4), '(i4.4)' ) year
        write( cdate( 6: 7), '(i2.2)' ) month
        write( cdate( 9:10), '(i2.2)' ) day
        write( cdate(12:13), '(i2.2)' ) hour
        write( cdate(15:16), '(i2.2)' ) min
        write( cdate(18:19), '(i2.2)' ) sec
      else
        read( cdate( 1: 4), '(i4.4)' ) year
        read( cdate( 6: 7), '(i2.2)' ) month
        read( cdate( 9:10), '(i2.2)' ) day
        read( cdate(12:13), '(i2.2)' ) hour
        read( cdate(15:16), '(i2.2)' ) min
        read( cdate(18:19), '(i2.2)' ) sec
      endif
      return
      end
      
