character*4 srcma(80),srcmb(80),srca(9001),srcb(9001) character dat*9,clock*8,tempfile*12 character*65 filea,fileb double precision avel,avog,dn,del(80,80),dog(80,80),dog2(80,80), 1 del2(80,80) dimension num(80,80),nsrca(9001),nsrcb(9001),iela(9001), 1 ielb(9001),ioga(9001),iogb(9001),lata(9001),lona(9001), 2 latb(9001),lonb(9001),ifil(2) logical exists data nred,nredund,na,nb,npair/5*0/,ma,mb/2*1/,ifil/1,2/ print 100 100 format (' COMPARE2. Plouff,3-98. Program to primarily summar', 1 'ize differences of',/,' observed gravity at nearly coinci', 2 'dent stations in two plouff-format files.',/,' Differences ', 3 'are grouped into combinations matched with the first four',/, 4 ' digits of station names. Files preferably should be sorted', 5 ' by station',/,' name and run through programs REDUND (no ', 6 'internally coincident stations)',/,' and COUNTDMA (no more ', 7 'than 80 unique 4-digit source-code names in files).',/, 8 ' Files must have no more than 9000 data points.') print 101 101 format (' The output of this program is in print file COMPARE2.', 1 'PNT') tempfile='COMPARE2.PNT' inquire (file=tempfile,exist=exists) if (exists) then print 102 102 format (' **STOP. File COMPARE2.PNT already exists. Rename ', 1 'it.') stop end if print 103 103 format (' TYPE the name of one file:') read 565, filea 565 format (a65) inquire (file=filea,exist=exists) if (.not. exists) then print 104 104 format (' **STOP. That input file does not exist.') stop end if print 105 105 format (' TYPE the name of the other file:') read 565, fileb inquire (file=fileb,exist=exists) if (.not. exists) then print 104 stop end if print 200 200 format (' Data points from the two sets need not exactly coinc', 1 'ide to be tentatively',/,' assumed to be at the same loca', 2 'tion. The tolerance is expressed in',/,' geographic min', 3 'utes (both longitude and latitude). This program should',/, 4 ' first be run with zero tolerance, to test for copied data.') 9 print 201 201 format (' TYPE the tolerance in minutes:') read (5,*,err=9) d d=abs(d) id=100.0*d+0.5 open (16,file='COMPARE2.PNT',form='formatted',status='new') call date(dat) call time(clock) write (16,160) dat,clock 160 format (1x,a9,1x,a8) write (16,100) write (16,200) write (16,202) d 202 format (' You selected a tolerance of',f6.2,' minutes.') open (7,file=filea,form='formatted',status='old') lines=13 1 na=na+1 if (na .gt. 9001) then print 203, filea write (16,203) filea 203 format (' **STOP. File has more than 9000 stations.',/,5x,a65) lines=lines+2 go to 9999 end if read (7,700,end=10,err=11) srca(na),latd,latm,lond,lonm, 1 iela(na),ioga(na) 700 format (bz,a4,4x,i3,3i4,i6,i7) c Location expressed in 0.01-minute units lata(na)=6000*iabs(latd)+latm lona(na)=6000*iabs(lond)+lonm if (na .eq. 1) then srcma(1)=srca(na) nsrca(1)=1 else do 2 i=1,ma if (srcma(i) .ne. srca(na)) go to 2 nsrca(na)=i go to 1 2 continue c Gets here if no match to previous source code was found ma=ma+1 if (ma .gt. 80) then print 106, filea write (16,106) filea 106 format (' **STOP. There are more than 80 unique 4-digit ', 1 'source codes in file:',/,5x,a65) lines=lines+2 go to 9999 end if nsrca(na)=ma srcma(ma)=srca(na) end if go to 1 10 na=na-1 print 107, ifil(1),filea,na,ma write (16,107) ifil(1),filea,na,ma 107 format (' File',i2,': ',a65,/,' has',i5,' data points and the ', 1 'following',i3,' unique source codes:') print 204, (srcma(k),k=1,ma) write (16,204) (srcma(k),k=1,ma) 204 format (13(1x,a4,']')) lc=ma/13 if (13*lc .ne. ma) lc=lc+1 lines=lines+2+lc open (9,file=fileb,form='formatted',status='old') go to 3 11 print 108, na,filea write (16,108) na,filea 108 format (' **STOP. Format error on line',i5,' of file:',/,5x,a65) lines=lines+2 go to 9999 3 nb=nb+1 if (nb .gt. 9001) then print 203, fileb write (16,203) fileb lines=lines+2 go to 9999 end if read (9,700,end=13,err=14) srcb(nb),latd,latm,lond,lonm, 1 ielb(nb),iogb(nb) latb(nb)=6000*iabs(latd)+latm lonb(nb)=6000*iabs(lond)+lonm if (nb .eq. 1) then nsrcb(1)=1 srcmb(1)=srcb(nb) else do 4 i=1,mb if (srcmb(i) .ne. srcb(nb)) go to 4 nsrcb(nb)=i go to 3 4 continue c Gets here if no match to previous source code was found mb=mb+1 if (mb .gt. 80) then print 106, fileb write (16,106) fileb lines=lines+2 go to 999 end if nsrcb(nb)=mb srcmb(mb)=srcb(nb) end if go to 3 13 nb=nb-1 print 107, ifil(2),fileb,nb,mb write (16,107) ifil(2),fileb,nb,mb print 204, (srcmb(k),k=1,mb) write (16,204) (srcmb(k),k=1,mb) lc=mb/13 if (13*lc .ne. mb) lc=lc+1 lines=lines+2+lc go to 15 14 print 108, nb,fileb write (16,108) nb,fileb lines=lines+2 go to 999 c Initialize 4-digit source array combinations and assoc values 15 do 17 i=1,ma do 16 j=1,mb num(i,j)=0 del(i,j)=0.0d0 del2(i,j)=0.0d0 dog(i,j)=0.0d0 16 dog2(i,j)=0.0d0 17 continue c Loop all data points for both arrays do 30 i=1,na ii=nsrca(i) latm=lata(i) lonm=lona(i) iflag=0 do 29 j=1,nb if (iabs(latb(j)-latm) .gt. id) go to 29 if (iabs(lonb(j)-lonm) .gt. id) go to 29 c Count total duplicates and flag filea data points with any duplicates nred=nred+1 iflag=1 jj=nsrcb(j) num(ii,jj)=num(ii,jj)+1 el=0.1*(ielb(j)-iela(i)) del(ii,jj)=del(ii,jj)+el og=0.01*(iogb(j)-ioga(i)) dog(ii,jj)=dog(ii,jj)+og del2(ii,jj)=del2(ii,jj)+el*el dog2(ii,jj)=dog2(ii,jj)+og*og 29 continue c Keep track of number of stations in filea that have any b duplicates if (iflag .ne. 0) nredund=nredund+1 30 continue 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') lines=lines+6 do 32 i=1,ma do 31 j=1,mb number=num(i,j) if (number .eq. 0) go to 31 npair=npair+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.9 write (16,610) srcma(i),srcmb(j),number,avel,stde,avog,stdg 610 format (2(2x,a4),i7,2f8.1,f9.2,f8.2) lines=lines+1 31 continue 32 continue nred=nred-nredund print 611, npair,nredund,nred write (16,611) npair,nredund,nred 611 format (' There are',i5,' matched pairs of source codes.',i6, 1 ' data points in file 1',/,' have points at the same loca', 2 'tion plus',i4,' extra matches in file 2.') lines=lines+2 999 close(9) 9999 print 109, lines 109 format (' Printout is in a file named COMPARE2.PNT, which has', 1 i5,' lines.') close(16) close(7) stop end