character ans*1,ch1(4)*1 character infile*65,pntfile*12,sta*8,dent*8,ch14*4,ch24*4 equivalence (ch14,ch1(1)) logical exists data n,m,ns,ifred/4*0/,ch14,ch24/2*' '/,flat,flon/2*0.0/ print 100 100 format (' MAPINDEX, Plouff. Program to list in file MAPINDEX.', 1 'PNT standard map names',/,' for gravity data in plouff ', 2 'format. MAPINDEX.PNT presumably will be sorted',/,' by ', 3 'map name in a 1- by 2-degree quadrangle.',/,' DO you want an', 4 ' explanation of how map names are created?') read 501, ans 501 format (a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') print 101 101 format (' Map names follow the style of 1:250,000 quads in Alas', 1 'ka (but 1x2D):',/,' First digit is A, B, C, or D for ', 2 'rows of 15-minute maps from south to north.',/,' Second digi', 3 't is numbers 1 to 8 for columns of maps from east to west.',/, 4 ' For 7-1/2-minute maps, append NE, NW, SE, or SW.',/,' Sta', 5 'tions along N or S edges of interior maps will be listed ', 6 'with the',/,' adjacent map.') print 110 110 format (' You can avoid printing data outside the desired quad', 1 'rangle by first',/,' using program PULLRECT to select data.') pntfile='MAPINDEX.PNT' inquire (file=pntfile,exist=exists) if (exists) then print 102, pntfile 102 format (' File ',a12,' already exists.',/,' Do you want ', 1 'to STOP to save it?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') stop end if print 103 103 format (' TYPE the name of your data file in plouff-format:') read 565, infile 565 format (a65) inquire (file=infile,exist=exists) if (.not. exists) then print 104 104 format (' That file was not found. Try once more.') print 103 read 565, infile inquire (file=infile,exist=exists) if (.not. exists) then print 105 105 format (' **STOP. That file was not found.') stop end if end if print 106 106 format (' There is an option to read pairs of stations from ', 1 'program REDUND.',/,' Are there pairs of redundant stations ', 2 'in this file?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') ifred=1 23 print 107 107 format (' TYPE latitude (integer degrees) of south edge of the', 1 ' 1- by 2-degree sheet:') read (5,*,err=23) lsd if (lsd .lt. 0 .or. lsd .ge. 90) go to 23 24 print 108 108 format (' TYPE longitude (integer degrees) of east edge of ', 1 'the sheet:') read (5,*,err=24) led led=iabs(led) if ((led-2*(led/2)) .ne. 0) then print 109 109 format (' Your response must be an even number.') go to 24 end if open (9,file=infile,form='formatted',status='old',blank='zero') open (7,file=pntfile,form='formatted',status='unknown') if (ifred .eq. 1) then write (7,702) 702 format (2(' MAP STATION LATITUDE LONGITUDE ELEV')) else write (7,703) 703 format (' MAP STATION LATITUDE LONGITUDE ELEV') end if c Start loop to read input. Increment line number. 2 n=n+1 c A line or first of a pair read (9,900,err=88,end=99) dent,ktd,ltm00,knd,lnm00,h 900 format (bz,a8,1x,i2,i4,1x,i3,i4,f6.1) c Find mapname (4 digits in array ch1) and convert to decimal minutes. call mapind (lsd,led,ktd,ltm00,knd,lnm00,ch14,glat,glon) c Inside quad if characters 3 and 4 are not blank. if (ch1(4) .ne. ' ') ns=ns+1 if (ifred .ne. 1) then c No pairs option write (7,701) ch14,dent,ktd,glat,knd,glon,h 701 format (2x,a4,1x,a8,i3,f6.2,i4,f6.2,f8.1) go to 2 end if c Second of nearly coincident pair of stations n=n+1 read (9,900,err=88,end=98) sta,ltd,ltm00,lnd,lnm00,el m=m+1 c Again find map name in case coordinates are slightly different. call mapind (lsd,led,ltd,ltm00,lnd,lnm00,ch24,flat,flon) if (ch14 .eq. ch24) ch24=' ' ih=h+sign(0.5,h) iel=el+sign(0.5,el) c ***Following line takes 2*39=78 columns. 78 write (7,700) ch14,dent,ktd,glat,knd,glon,ih,ch24,sta,ltd,flat, 1 lnd,flon,iel 700 format (2(1x,a4,1x,a8,i3,f6.2,i4,f6.2,i6)) go to 2 88 print 588, n 588 format (' **STOP. Bad format on input line',i5) go to 999 98 n=n-1 print 598, n,m 598 format (' Data ended at line',i5,' before second of pair',i4, 1 ' was read.') go to 999 99 n=n-1 print 111, n 111 format (i8,' input lines were read.') if (ifred .eq. 1) then print 112, m,ns 112 format (i8,' station pairs were read.',/,i8,' pairs of sta', 1 'tions are inside the quadrangle.') write (7,702) else print 113, ns 113 format (i8,' stations are inside the quadrangle.') write (7,703) end if write (7,799) infile 799 format (' Input file=',a65) print 114 114 format (' To sort the output by map name, for example, ', 1 'TYPE:',/,5x,'sort MAPINDEX.PNT -o MAPINDEX.SRT') 999 close (7) close (9) stop end subroutine mapind (lsd,led,ltd,ltm00,lnd,lnm00,ch4,flat,flon) character ch(4)*1,ltdig(4)*1,lndig(8)*1,ch4*4,temp*4 equivalence (temp,ch(1)) c Data variables retain initial values data ltdig/'A','B','C','D'/ data lndig/'1','2','3','4','5','6','7','8'/ temp=' ' c Convert coordinates to 0.01 minute. Relative to S, E map edges. ltmap=6000*(ltd-lsd)+ltm00 if (ltmap .gt. 6000) ch(1)='N' c South edge of map is part of map to south. East, east. if (ltmap .le. 0) ch(1)='S' lnmap=6000*(lnd-led)+lnm00 if (lnmap .gt. 12000) ch(2)='W' if (lnmap .le. 0) ch(2)='E' c Return if outside quadrangle. if (temp .ne. ' ') go to 9 c Latitude reference in units of 15 minutes. ltindex=ltm00/1500 ch(1)=ltdig(ltindex+1) ch(3)='N' if ((ltm00-ltindex*1500) .lt. 751) ch(3)='S' c Subtracted even degrees for longitudes. lnindex=lnmap/1500 ch(2)=lndig(lnindex+1) ch(4)='W' if ((lnmap-lnindex*1500) .lt. 751) ch(4)='E' 9 flat=0.01*float(ltm00) flon=0.01*float(lnm00) ch4=temp return end