c DMA91 has coordinates in degrees and minutes and is sorted by c latitude and longitude blocks (was pure latitude sort before). c Locations are arranged are in a sequence in one-degree rows from c south to north. For a given latitude row, there are cells of one c degree longitude (negative numbers for degree part) arranged from c west to east. For a given one-degree cell, latitude minutes c carried to 0.01 (without decimal point) increase northward. c Program "testdma91" showed that longitudes range from -174.5 to c -66.1 degrees. Editing showed that latitudes range from 18.9 (only c at -155 longitude) to 71.5 (only -157 to -154 degree bands) degrees. character infile*65,outfile*65,rest*53,source*4,sequence*4,ans*1 logical exists print 100 100 format(' Plouff, 1-94. PULLDMA91.',/,' Program to extract U.S. ', 1 'DMA gravity data from file /bgp/data/ngdc/dma91.dat',/,' The ', 2 'file has data from 18 to 72 degrees latitude and -174 to -66 ', 3 'degrees',/,' longitude. See Donald Plouff for a map of ', 4 'coverage. Your output will be',/,' in abbreviated plouff ', 5 'format through columns including observed gravity.',/, 7 ' Does the above file name specify the correct path?') read 501, ans 501 format(a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') then infile='/bgp/data/ngdc/dma91.dat' else print 101 101 format (' TYPE the pathname of the DMA91 gravity file:') read 565, infile 565 format(a65) end if inquire (file=infile,exist=exists) if (.not. exists) then print 102 102 format (' Pathname does not work. Try once more.') print 101 read 565, infile inquire (file=infile,exist=exists) if (.not. exists) then print 103 103 format (' **STOP. Pathname failure.') stop end if end if print 104 104 format (' TYPE a name for the Plouff output file:') read 565, outfile inquire (file=outfile,exist=exists) if (.not. exists) go to 1 print 105 105 format (' A file with that name exists. Do you want to ', 1 'overwrite it?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 1 print 104 inquire (file=outfile,exist=exists) if (.not. exists) go to 1 print 105 read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 1 print 103 stop 1 print 109 109 format (' Next, provide geographic boundaries of the area you ', 1 'want to select.',/,' TYPE the latitude of the south edge ', 2 '(degrees, minutes):') read (5,*,err=99) d,f c Convert to integer minutes to 0.01 lats=6000.0*abs(d)+100.0*abs(f) c Prepare fastest search for integer degrees part. ldegs=lats/6000 print 106 106 format (' TYPE the latitude at the north edge (deg, min):') read (5,*,err=99) d,f c 0.9 for user input to 0.001 minute. South was just a truncation. latn=6000.0*abs(d)+100.0*abs(f)+0.9 ldegn=latn/6000 print 108 108 format (' TYPE the longitude at the east edge (deg, min):') read (5,*,err=99) d,f lone=6000.0*abs(d)+100.0*abs(f) ldege=lone/6000 print 107 107 format (' TYPE the longitude at the west edge (deg, min):') read (5,*,err=99) d,f lonw=6000.0*abs(d)+100.0*abs(f)+0.9 ldegw=lonw/6000 c Plouff output (keyed to master terrain file) requires all positive c values for the conterminous U.S., and extraction here is faster by c reading longitudes without negative sign. c Program TESTDMA91 tested the longitude range -174D11.3M to -66D04.M c File latitudes range from 18D54.96M to 71D06.69M from editor if (lats .ge. latn .or. lone .ge. lonw .or. lats .lt. 108000 1 .or. latn .gt. 432000 .or. lone .lt. 396000 .or. lonw .gt. 2 1050000) then print 599 599 format (' **STOP. Your geographic range is not in the file ', 1 'or not in sequence. Maximum',/,' range is 18 to 72 ', 2 'degrees N latitude and 66 to 175 degrees W longitude.') stop end if open (10,file=infile,status='old',form='formatted',blank='zero') open (15,file=outfile,status='unknown',form='formatted') icnt=0 go to 10 99 print 188 188 format (' **STOP. Bad format for typing coordinates.') stop c Read DMA91.DAT format without coordinate signs in columns 4 and 12. 10 read (10,800,end=9,err=77) ltd,ltm,lnd,lnm,rest 800 format (bz,4x,i2,i4,2x,i3,i4,a53) c Get to south edge of data as fast as possible with next line if (ltd .lt. ldegs) go to 10 c Next fastest tests (longitude degrees are negative) if (lnd .gt. ldegw .or. lnd .lt. ldege) go to 10 c Stop reading if no stations were in northeastmost requested cell if (ltd .gt. ldegn) go to 9 lat=6000*ltd+ltm if (lat .lt. lats) go to 10 c Stop reading when past north minutes limit in NE-most degree cell if (lat .gt. latn) then if (lnd .eq. ldege) go to 9 go to 10 end if c Longitudes are treated as positive (both degrees and minutes are +). lon=6000*lnd+lnm if (lon .lt. lone .or. lon .gt. lonw) go to 10 icnt=icnt+1 read (rest,801,err=77) jel,iobs,source,sequence 801 format (bz,3x,i7,7x,i6,15x,a4,7x,a4) c Convert elevations to nearest foot. "iev" in column 22 (blank here) elft=0.1*float(jel)/0.3048 jel=elft+sign(0.5,elft) jel=10*jel iobs=iobs+7600000 c Convert to positive longitudes lnd=iabs(lnd) write (15,150) source,sequence,ltd,ltm,lnd,lnm,jel,iobs 150 format (2a4,i3,3i4,i6,i7) go to 10 9 print 900, icnt 900 format(i9,' data points were obtained. Sort output file, to ', 1 'gather by',/,' ascending station name, which optimizes later', 2 ' tests. Then run program',/,' cntsrc to list how many data ', 3 'data points were obtained from each source.') go to 999 77 n=n+1 print 577 577 format (' **STOP. Format error in DMA91 file.') print 578, icnt 578 format (' Your output file has',i5,' data points to there.') 999 close (10) close (15) stop end