character*4 src1,src2,src(75)
      character dat*9,clock*8,cntfile*65,pairfile*65,temfile*12,ans*1
      logical exists
      double precision avel,avog,dn,del(75,75),dog(75,75),dog2(75,75),
     1 del2(75,75)
      dimension num(75,75)
      print 100
  100 format (' Plouff 1-89',/,' Program to read file of pairs of ',
     1'redundant stations from program REDUND and',/'  to print the ',
     2'number of pairs for each combination of 4-digit source codes.',
     3 /,' The original file usually was obtained from a DMA file ',
     4 'and was sorted',/,'  by station name.  The list of ',
     5 'source codes (or first 4 digits of station ',/,'  names) is ',
     6 'obtained from the output of program COUNTDMA with the file ',
     7 'of pairs',/,'  as input.')
      temfile='COMBPAIR.PNT'
      inquire (file=temfile,exist=exists)
      if (exists) then
        print 200
  200   format (' A print file COMBPAIR.PNT already exists.',/,
     1   ' Do you want to STOP to save it?')
        read 501, ans
  501   format (a1)
        if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') stop
      end if
      print 101
  101 format (' TYPE the name of the .PNT file from the CNTSRC or ',
     1 'COUNTDMA programs:')
      read 565, cntfile
  565 format (a65)
      inquire (file=cntfile,exist=exists)
      if (.not. exists) then
        print 201
  201   format (' **STOP. That file does not exist.')
        stop
      end if
      print 102
  102 format (' TYPE the name of the output pairs file from the ',
     1 'REDUND program:')
      read 565, pairfile
      inquire (file=pairfile,exist=exists)
      if (.not. exists) then
        print 202
  202   format ('That file does not exist.  Try once more.')
        print 102
        read 565, pairfile
        inquire (file=pairfile,exist=exists)
        if (.not. exists) then
          print 201
          stop
        end if
      end if
      open (16,file='COMBPAIR.PNT',form='formatted',status='unknown')
      call date(dat)
      call time(clock)
      write (16,160) dat,clock
  160 format (1x,a9,1x,a8)
      write (16,100)
      write (16,161) cntfile
  161 format (' COUNTDMA program file of source codes:',/,3x,a65)
      open (7,file=cntfile,form='formatted',status='old',blank='zero')
      write (16,162) pairfile
  162 format (' REDUND program file of pairs:',/,3x,a65)
      open (9,file=pairfile,form='formatted',status='old',blank='zero')
      n=0
    1 read (7,700,end=2) src1
  700 format (1x,a4)
      n=n+1
      if (n .gt. 75) then
        print 103
        write (16,103)
  103   format (' ***STOP. More than 75 source codes is prohibited.')
        go to 999
      end if
      src(n)=src1
      go to 1
    2 write (16,603) n
      print 603, n
  603 format (' There is a total of',i3,' source codes.')
      if (n .lt. 2) then
        print 104
        write (16,104)
  104   format (' ***STOP. There are too few source codes to continue.')
        go to 999
      end if
      do 5 i=1,n
      do 4 j=1,n
      num(i,j)=0
      del(i,j)=0.0d0
      del2(i,j)=0.0d0
      dog(i,j)=0.0d0
    4 dog2(i,j)=0.0d0
    5 continue
      m=0
    6 read (9,900,end=99) src1,lel1,log1
  900 format (bz,a4,19x,i6,i7)
      read (9,900,end=88) src2,lel2,log2
      m=m+1
      do 7 i=1,n
      if (src1 .eq. src(i)) then
        ii=i
        go to 8
      end if
    7 continue
      go to 77
    8 do 9 j=1,n
      if (src2 .eq. src(j)) then
        jj=j
        go to 10
      end if
    9 continue
      go to 77
   10 num(ii,jj)=num(ii,jj)+1
      el=0.1*(lel2-lel1)
      del(ii,jj)=del(ii,jj)+el
      og=0.01*(log2-log1)
      dog(ii,jj)=dog(ii,jj)+og
      del2(ii,jj)=del2(ii,jj)+el*el
      dog2(ii,jj)=dog2(ii,jj)+og*og
      go to 6
   77 line2=2*m
      line1=line2-1
      print 177, src1,src2,line1,line2
      write (16,177) src1,src2,line1,line2
  177 format (' ***STOP. Source code ',a4,' or ',a4,' in the station ',
     1 'pair on lines',i5,' and',i5,/,'  was not found. Preliminary ',
     2 'results will not be printed.')
      go to 999
   88 print 188, m
      write (16,188) m
  188 format (' ***WARNING. The pairs file ended prematurely after',
     1 i5,' pairs plus one station.',/,' Therefore, the following ',
     2 'results are incomplete or doubtful.')
   99 print 199, m
      write (16,199) m
  199 format (' A total of',i5,' pairs were read.')
      write (16,609)
  609 format (/,' A list of paired source codes follows.  The list ',
     1 'includes the',/,'  number of occurrences, the average ',
     2 'difference of elevations in feet, the',/,'  the standard devi',
     3 'ation of the elevation difference, the average difference',/,
     4 '  of observed gravity in milligals, and the standard ',
     5 'deviation in milligals.',/,' CODE1 CODE2 NUMBER EL2-EL1 ',
     6 'STD DEV  OG2-OG1 STD DEV')
      k=0
      do 12 i=1,n
      do 11 j=1,n
      number=num(i,j)
      if (number .eq. 0) go to 11
      k=k+1
      dn=number
      avel=del(i,j)/dn
      avog=dog(i,j)/dn
      if (number .eq. 1) then
        stdg=0.0
        stde=0.0
      else
        stde=dsqrt(dabs(del2(i,j)-dn*avel*avel)/(dn-1.0d0))
        stdg=dsqrt(dabs(dog2(i,j)-dn*avog*avog)/(dn-1.0d0))
      end if
      if (avel .lt. -99999.9) avel=-99999.9
      if (avel .gt. 999999.9) avel= 999999.7
      write (16,610) src(i),src(j),number,avel,stde,avog,stdg
  610 format (2(2x,a4),i7,2f8.1,f9.2,f8.2)
   11 continue
   12 continue
      print 209, k
      write (16,209) k
  209 format (' There are',i4,' different pairs of source codes.',/,
     1 ' Some pairs may be repeated in reverse order if the ',
     2 'original file was not',/,'  initially sorted by station name.')
  999 print 109
  109 format (' Printout of results are in a file named COMBPAIR.PNT')
      close(16)
      close(7)
      close(9)
      stop
      end