U.S. DEPARTMENT OF THE INTERIOR

U.S. GEOLOGICAL SURVEY

Field Estimates of Gravity Terrain Corrections and Y2K-Compatible

Method to Convert from Gravity Readings with Multiple Base Stations

to Tide- and Long-Term Drift-Corrected Observations

By Donald Plouff1

Open-File Report OF 00-140

Version 1.0

2000

U.S. DEPARTMENT OF THE INTERIOR

BRUCE BABBITT, Secretary

U.S. GEOLOGICAL SURVEY

Charles G. Groat, Director

This report is preliminary and has not been reviewed for conformity with U.S. Geological Survey editorial standards or with

the North American Stratigraphic Code. Any use of trade, firm, or product names is for descriptive purposes only and does

not imply endorsement by the U.S. Government. Although the computer programs and supporting information have been

developed, tested, and used by the U.S. Geological Survey, no warranty (expressed or implied) is made by the author or the

U.S. Geological Survey as to the accuracy and functioning of the computer programs, constituent program parts, and related

program materials. The fact of distribution shall not constitute any warranty, and no responsibility is assumed by the U.S.

Government in connection therewith.

1 U.S. Geological Survey, 345 Middlefield Road, Menlo Park, CA 94025
i

CONTENTS

Page

Abstract 1

Introduction 1

Acknowledgments 3

Field data format 3

Steps of reduction process 6

Step 1—Test field data and output lists of data points and maps 6

Step 2—Obtain geographic coordinates 6

Step 3—Prepare gravity meter calibration table 7

Step 4—Convert readings to milligals and correct for effect of tide 8

Step 5—Utilize known or determined values of observed gravity 8

Step 6—Calculate gravity meter drift and values of observed gravity 9

Step 7—Plot gravity meter drift 11

Step 8—Delete redundant data 11

Additional comments about reduction process 13

Gravity terrain corrections 13

Programming considerations 17

References 17

Appendixes

1. Accuracy codes 19

2. Assembling maps, testing data, and listing data points (fieldg) 20

3. Obtaining geographic coordinates to merge with field data (extractg, geog5) 22

4. Creating a gravity meter calibration file (calibg) 24

5. Converting readings to milligals and correcting for tide effect (mergeg) 25

6. Converting tide-corrected values to observed gravity (setbaseg, driftg) 27

7. Plotting instrument drift (driftplt) 31

8. Deleting redundant data (extraobs) 32

9. Estimating gravity terrain correction 33

Terrain correction of a post (tcpost) 33

Terrain correction for a two-dimensional profile (fieldtc) 35

ILLUSTRATIONS

FIGURE 1. Organization of gravity reduction method 2

2. Field data sheet 4

3. Drift plots 12

4. Gravity terrain correction of a displaced slope (Hayford system) 15

5. Gravity terrain correction of a displaced slope (Hammer system) 16

6. Gravity terrain correction for reading on a concrete post 33
1

ABSTRACT

FORTRAN computer programs to convert gravity readings to absolute values of observed gravity with

minimum data entry and ample prompting are documented and released with this report. Inasmuch as

gravity meter drift is linked for up to one calendar year of field data, the first line of the data input file only

has the gravity meter name, the Y2K-compatible 4-digit year, the correction to Greenwich Mean Time

(Coordinated Universal Time), and a code for the distance to which field terrain corrections are estimated.

Only the 2-digit month and day are entered before each day of data. Each line of data has a 5-digit station

name, local time, gravity meter reading, elevation, elevation units, an optional accuracy code, the estimated

terrain correction, and an optional map code. Geographic locations and base station values are searched for

and extracted from other files. Calibration tables in digital format are interactively created only for the

needed span of gravity readings. Decisions for successively assigning absolute gravity values for gravity

stations with more than one observation are based on printed summaries of least-square values computed in

previous iterations and are confirmed by subsequent gravity meter drift plots. Field estimates of gravity

terrain corrections utilize lists and graphs or entry of distances, heights, and angular slopes into a program

based on models for thin vertical slabs. 

INTRODUCTION

Starting with a simple gravity reduction method in 1970 (Plouff, 1972) and a modular system in 1980

(Plouff, 1982), a comprehensive system of FORTRAN computer programs to convert field gravity readings

to values of tide- and drift-corrected observed gravity has been developed. Input field data files consist of

readings for one gravity meter during as much as one calendar year. Output data and print files include

values of observed gravity, lists of tide and drift corrections, standard deviations for repeated observations,

and plots of gravity meter drift. The system is still in a state of development to minimize the apparent

complexity of eight steps (fig. 1), to print gravity tide correction tables before fieldwork, to improve

accuracy to the microgal level, to utilize elevations from Global Positioning Systems, to correct horizontal

and vertical datum shifts, to accept data collected in east longitudes and south latitudes, to accept different

data formats including spreadsheets, and to utilize window environments. 

The main criteria for design of the reduction system are to link gravity meter drift for a calendar year, to

minimize the amount of data and program entry, to accept elevations in either meters or feet, to interactively

input instrument calibration factors only for the needed range of gravity readings, to identify data input

errors, to provide detailed information and a summary so that the user iteratively can assign best values of

observed gravity at "stations" (sites or locations) with multiple "readings" (observations), to explain the role

of the program at each step, and to suggest the next step in the process. Dates of readings and closest base

station in time are retained as elements of the final output for each reading. Retaining dates and abbreviated

map locations before including geographic coordinates is less cumbersome than unnecessarily embedding

these quantities into station names. The original raw field data file is not changed unless data fields are

incomplete or errors are detected. Vulnerability to typographical errors (written commun., D.L. Healey,

1980; written commun., Dolores Wilson, 1980; and Plouff, 1990, p. 7-8) in formats that require daily input

of base station values is minimized by utilizing a separate file with previously determined values of

observed gravity, which is expanded throughout reduction processes. If later readings yield more accurate

values of observed gravity, inspection of earlier print summaries suggests entry points to reprocess

concurrently and subsequently determined values of observed gravity. 

The user types the same short dataset name for each program step of the reduction system, and files

consequently are read or created with that root name followed by a decimal point and a 3-digit suffix for the

role of the file in that processing step. A file of unique station names and associated numbers of readings

is output from the first step of the process. Data entry errors and ranges of dates and gravity readings are

printed. After correcting errors found in step one, files are searched to obtain geographic coordinates for a

diminishing list of unique station names during step two. Unless previously available, the user can create

an easily proofread gravity meter calibration table needed to cover the range of field readings in step three.

Tide corrections and conversions from gravity meter readings to equivalent-milligal observations are

determined in step four. Values of observed gravity and other parameters such as defining the time of a

sudden datum shift ("tare"), which may be caused by change of instrument temperature or helicopter

vibration, for example, are edited or interactively input for the recursively repeated step five. Values of
Figure 1.--Organization of gravity reduction method. Program names are in shaded boxes. All

programs have an output print file with a .PNT suffix after the program name. Program input and

output file names have root name of data set followed by a decimal point and the indicated 3-digit

suffix. *, user- selected file name; **, name of a file with geographic coordinates.

2

FIELDTC 

(interactive input)

STEP 1 .FLD FIELDG .NAM

.MAP

.FLD .FLD* GMTCHNG 

.FLD .PNT FIELDMAP MAP*

STEP 2 EXTRACTG .NAM

. **

.LFT

.FCT

.REF

.NAM

.FCT

STEP 3 CALIBG 

(interactive input)

.CAL

STEP 4 MERGEG 

.FLD

.FCT

.CAL

.TID

STEP 5 SETBASEG .BAS

.TID

.BAS

.GIV .TID

STEP 6 .TID DRIFTG .DFT

.OBS

STEP 7 .DFT .PST DRIFTPLT

STEP 8 .OBS .SRT .PLF EXTRAOBS
3

observed gravity determined for reoccupied gravity stations are listed in detail, and their average values are

summarized in the recursively repeated step six. Output files from step six are a file that is input to a

program that creates a gravity meter drift plot in step seven and a file that is input to a program that retains

only the latest or best value of observed gravity in step eight. Examples of usage and data formats for each

step are included in following sections of this report. Validation of computer algorithms and highlights of

computer sessions for theoretical examples and reduction of data primarily collected during a short field

session in 1993 (Plouff, 1996) are included in appendixes. 

Copies of programs described in this report can be down-loaded from Web Site

http://geopubs.wr.usgs.gov/open-file/of00-140/

at the U.S. Geological Survey. Comments and suggestions about program performance are welcome. 

ACKNOWLEDGMENTS

David F. Barnes suggested a computer algorithm to compute tide corrections. Robert C. Jachens

provided an algorithm to correct for the effect of circular dial errors of gravity meters. Carter W. Roberts

calculated times and places of solar eclipses to determine times of local maxima for earth tides after 1999.

David A. Ponce supported implementation of short geographic-based names for labeling field maps and

files. Robert L. Morin suggested improvements to the text and supplied a test data set for Alaska. 

FIELD DATA FORMAT

Copies of field sheets (fig. 2) to be regularly mailed or faxed to the office includes information needed

to recover station locations (odometer readings, map names, and descriptions), gravity readings, elevations,

and closein gravity terrain corrections in case field maps are lost. The computer program fldform.f to

create the field sheet can be modified or the print file FLDFORM.PS can be edited to include the project

name and the gravity meter name, for example. The following sample of field data (file with suffix .FLD)

illustrates the simplicity of data format and optional comments (from Plouff, 1996):

G8B 1993 7 CHAMMER

 6-21 [including the year as 06-21-93 is optional]

WNBAS 0912 346591 13100MN22002WD7NW City Hall sidewalk; 13' SE of K18 Reset

WINNA 0943 344968 4298 WD8NE Winne airport phone booth. ACIC 0474-1

WNBAS 1012 346599 Later check, uncertain reading?

WINNA 1038 344965

WN 6R1551 3533.98 12579MN22001VB1NE 1ft S of BM STR. Near station V8604?

WNBAS 1809 346594

999999999999999999999999999999999999

 6-22

WNBAS 0746 346600

WN 6 0842 353412

WN 11R0902 352782 1257 MG55001VB1NE Rd CL, about max curve of abandoned rd to NE

WNBAS 1911 346592

9999999999999999999999999999999999999

The gravity meter name on line one is G8B. The previous version of this reduction system did not read

the year but arbitrarily assigned years after 1999 if the 2-digit years in lines 2 and 10 were less than 40. A

4-digit calendar year typed only in line one, however, saves typing a hyphen and a 2-digit year for every

day, treats the year-2000 problem in a straightforward way, and is carried forward with the gravity meter

name throughout the reduction process to be available at the stage of plotting gravity meter drift. The

number 7 in line one is the correction in hours to be added to local Pacific Daylight Time in order to obtain

Greenwich Mean Time (Coordinated Universal Time). Terrain corrections were estimated to a distance of

175 ft (53 m) in the field, which includes zone "C" of the Hammer (1939) system. 

The month and day are entered as five keystrokes in lines 2 and 10 before each day of readings. For

example, " 6-21" refers to June 21, 1993. In a previous program, eight keystrokes were needed for " 6-21-

93," which could have been replaced by a more universally accepted expression, "19930621" (followed by

hours and minutes), with same number of keystrokes (written commun., C.W. Roberts, 1997). Tests for

the two minus signs or diagonal slashes in the date were used in a preliminary version of the computer
START: 

LUNCH: 

 STOP: 

TOTAL: 

PROJECT DATE (1-8) DATA SHEET 

GRAVITY METER/CALIB LAST NEW STATION 

TIME ZONE OPERATOR RECORDER 

ODOMETER STATION R 

? TIME READING TIDE-TEMP 

CORR-READ 

ELEVATION /INCREMENT 

M 

? T.C. 29-31 

ACC. 26-28 STATION DESCRIPTION MAP 

1-5 6 7-10 11-18 19-24 25 T.C. RING 32-36 

LINE WITH ALL NINES AFTER LAST READING OF DAY D. PLOUFF 4-99 FORM

Figure 2.--Field data sheet. Reduced from full size of page.

4
5

program to flag the start of a new day before delimiter lines with nines were introduced to eliminate the

need to re-program a re-read function for different computers. 

The 5-digit "station name" (data point identifier) in columns 1-5 can consist of a one- or two-letter

prefix and numbers increasing to 999, for example. Inasmuch as the date and suggested map designation

are carried forward in the process and the five digits later are expanded to eight digits, the year or

abbreviation of a map name are not needed at this stage. If many new prefixes are assigned and numbers

are restarted often during a field season, the original field data are harder to locate after the data are

absorbed into larger sets. The letter, "R," in column 6 on lines 7 and 13 indicates that preliminary printout

of a computed absolute value of observed gravity is wanted for a gravity reading possibly near a station

from a previous gravity survey. Detailed printout for gravity stations with only one observation otherwise

is suppressed during preliminary evaluation until all values of observed gravity are established for

reoccupied gravity stations. The local time in hours and minutes is in columns 7-8 and 9-10, respectively. 

The time and effort of data entry are minimized by using a format without decimal points for gravity

readings, elevations, and gravity terrain corrections. Elevations, accuracy codes, terrain corrections, and map

designations are not re-typed for later readings at the same gravity station. Space is provided to express

gravity readings to the equivalent of microgals. Column 18 for the microgal digit was left blank in the

above example to denote readings recorded in units of 0.01 mGal. Unless field data collection methods

include additional readings with other gravity meters or frequent re-occupations, readings should only be

carried to the equivalent of 0.01 milligal (mGal) to minimize data recording errors. Elevations expressed to

0.1 unit are entered in units of meters (followed by "M") in column 25 or feet (followed by a blank space),

depending on the units expressed on the map or other source of elevations. An auxiliary computer

program might be needed to transfer elevations downloaded from Global Positioning Systems or surveying

instruments to the field data file. Space is provided on the data sheet so that "1310.2-0.2" originally was

written to retain the information that the reading was made at a location 0.2 m lower than the benchmark

elevation of 1310.2 m, but the elevation on line 3 was entered as "13100" in the digital file. Other data forms

may include a column for "height to reference," which is ambiguous. 

The three digits in columns 26-28 are optional codes for a description, elevation accuracy, and

horizontal location accuracy, respectively (appendix 1). The fourth digit of the accuracy code is the

accuracy of observed gravity, which is available after step 8 has been completed. To insure correct spacing,

leading zeroes may be typed for the value of the field terrain correction (TC) expressed in units of 0.01

mGal in columns 29-31. Complete TCs carried to a standard distance of 166.7 km probably will not have

an accuracy better than 0.1 mGal. Therefore, typing a negligible value of 0.01 for the field TC rather than a

true 0.00 can be used as a flag to indicate that the field TC has been estimated. Survey notes especially

should be recorded for terrain corrections greater than 0.05 mGal. 

The optional map name in columns 32-36 can be obtained from a system adopted by D.A. Ponce (oral

commun., 1992) and others for organizing 7.5-minute maps in study areas of 1 by 2 degrees. The first

digit in the 5-digit map name is a symbol for the 1° by 2° quadrangle, in which the gravity station is located.

The second digit reflects the northward sequence, A-B-C-D , of 15-minute quadrangles in a 1° by 2°

quadrangle, and the number in column 34 reflects the westward sequence 1-2-3-4-5-6-7-8. Letters NE, NW, SE,

or SW in columns 35-36 refer to the location of a 7.5-minute quadrangle within a 15-minute quadrangle.

The map name WD7NW in line 3 is the 7.5-minute quadrangle with a southeast corner of lat 40° 52' 30" N.,

lon 117° 37' 30" W. in the Winnemucca ("W"), Nevada 1° by 2° quadrangle. The map designation of the

U.S. Geological Survey is "40117-H8" for that 7.5-minute quadrangle, based on a northward sequence of

A-B-C-D-E-F-G-H and a westward sequence of 1-2-3-4-5-6-7-8 relative to the southeast corner of a 1°

quadrangle. Therefore, a more universal 4-digit designation "07H8" based on a 10-degree system can

replace "WD7NW" based on a 1- by 2-degree system for arranging and identifying field maps. 

Field data recorders should not hesitate to write numerous comments about the location, gravity meter

behavior, and topography, for example, on field sheets even if the comments are not entered into a digital

file. Geographic coordinates acquired by a Global Positioning System should be recorded. Odometer

readings should be recorded for each gravity station location to speed future re-occupations, to insure

against loss of field maps, and to recalculate locations if the odometer calibration is determined to be

erroneous. Odometer readings and notes about road conditions such as locked gates, difficult crossings,

and unmapped roads also should be written on field maps. Abbreviations for compass directions

(N, S, E, W), road (rd), highway (hwy), centerline (CL), Vee, Tee, River (R), crossing (Xing), and, benchmark
6

(BM), for example, save space on data sheets. Ambiguous comments such as "at spot elevation" are not

useful without specific descriptions such as "rd CLs" or "hillcrest." The first field sheet of a field session

should include comments such as confirmation that level bubbles and thermostat temperature, have been

checked. For recording numbers, the number four should include a horizontal line that distinctly crosses

the vertical line (unlike a 9), the top parts of fives should have square corners (unlike a 6), and diagonal

slashes are not needed. Observers should call out gravity readings, and recorders should repeat the

numbers with different enunciated groupings. 

STEPS OF REDUCTION PROCESS

Step 1—Test field data and output lists of data points and maps

Program fieldg.f (appendix 2) prints a detailed summary of the contents of the input .FLD file,

including the number of days, the range of gravity readings, number of gravity stations, and number of

maps, and lists the names of the gravity stations and maps, including the number of occurrences of each.

Typographical and other errors, including an unexpected range of gravity readings and wrong time

chronology, are printed so that errors can be corrected and program fieldg.f can be rerun before further

processing continues. In addition to the printout, the program outputs a .NAM file of station names and

number of occurrences, for example,

WNBAS 29

WINNA 3

WN 1 1

WN 6 9

and a .MAP file of map names and number of occurrences, for example,

 Maps listed in file blm938.FLD

MB8NW 5

VB1NE 33

MB8SW 24

Alternatively, program fieldmap.f reads file .FLD to prepare a list of maps to be digitized or to show data

coverage, for example. 

Program fieldg.f allows only one Greenwich Time correction. In case field data are collected during

both Daylight Savings and Standard Times, and the gravity meter drift correction is to be preserved without

interruption, the readings of one of the two data sets can be converted to the other Greenwich time

correction by executing program gmtchng.f, and the two data sets can be merged before executing

program fieldg.f. 

Step 2—Obtain geographic coordinates

Program extractg.f searches files to find geographic coordinates for stations with the same 5-digit

name as listed in the .NAM file (appendix 3). The first file to be searched is a file with a suffix .LOC,

which includes most of the new data. The format of the .LOC file is either the 5-digit station name

followed by manually digitized geographic coordinates in compressed form, for example,

WINNA 405423 1174822

WN 1 412388 1174945

WN 6 412458 118 117

or digitizer output with three columns before the 5-digit station name, for example,

 WN 1 412388 1174945

 WN 6 412458 118 117

Latitudes and longitudes without algebraic signs are expressed in integer degrees (40 and 41; 117 and 118)

and minutes to 0.01 without decimal points. Program extractg.f prints and outputs (in a .LFT file in the

same format as the .NAM file) a list of stations not found. For example, two of 263 stations were not

found in a simulation (appendix 3) of the first run for data collected in 1993 by Plouff (1996): 
7

WNA11 2

 P163 4

A .FCT file also was created by program extractg.f, which includes stations

WINNA4054231174822 0 3 2 0

WN 14123881174945 0 1 1 0

WN 64124581180117 0 9 2 0

The column with numbers 3, 1, and 9 refers to the number of observations carried forward from the .NAM

file. 

Station P163 was found in a file of previously established stations with the final processing format of

 P163 405500 1182281 416307984851G533 -1088 71284 2 7S-15409 13WINNA

The first .NAM file was discarded or re-named, the .LFT file was re-named to become the second .NAM

file, and program extractg.f was executed to acquire station P163. Station WNA11 remained in a .LFT file.

Executions of program extractg.f to acquire other locations after .LOC input create .REF files rather than

.FCT files. The .REF files are appended to the .FCT file in preparation for step 4. Inasmuch as station

WNA11 is another reading near station WN 11, a text editor was used to duplicate the coordinates of WN 11 to

create WNA11 in the final .FCT file needed for step 4. If only a few stations remain in the .LFT file, a

program to prompt for manually-derived (Plouff, 1998) geographic coordinates might be useful. 

Step 3—Prepare gravity meter calibration table

Program calibg.f creates a calibration table by prompting the user for values expressed in units of

milligals per dial reading unit to include the range of readings printed by program fieldg.f (appendix 4).

Inasmuch as the output .CAL file is easy to prepare and proofread without programmer assistance, the

alternatives of attempting to embed many calibration tables into one computer program has proved to be

vulnerable to programmer error (Plouff, 1990, p. 8). Programming to accept calibration tables as

"PARAMETER" arrays also is more complicated. Another .CAL file easily can be created or an initial file

can be expanded in a text editor if a larger range of gravity meter readings is needed for another gravity

survey. 

Calibration values are obtained from tables provided by the manufacturer of the gravity meter. Program

calibg.f first requests a constant factor—based on empirical observations (Barnes and others, 1969)—to

be multiplied by the manufacturer's values. Next, the "counter reading" interval, which refers to the interval

between values in the manufacturer's table, for example, 100 dial units, is requested. The lowest counter

reading needed then is requested. The following interactive prompting requests the milligal-equivalent of

the first and successive counter readings, for example, to create a G8B.CAL calibration file. 

 TYPE value in milligals after counter-reading prompt.

 Dial 3200 milligal-equivalent (zero to STOP): 3343.072

The milligal-equivalent of each counter reading can be expressed in units of microgals. Inasmuch as

counter readings may only be expressed by the manufacturer in units of 0.01 milligal, the interval factors

may be used to convert manufacturer's counter readings to units of microgals. 

The first three lines of file G8B.CAL created to reduce data collected in 1993 (Plouff, 1996) are:

G8B INTERVAL= 100 CIRC ERR= 0 FACTOR= 1.000612

 3200 3343.072 1.045900

 3300 3447.662 1.045960

The first line indicates the gravity meter name (also in line one of the .FLD file), the counter interval for

linear interpolation, the lack of a "circular error" correction, and the constant multiplier to be applied.

Editing can change the constant multiplier if the file is to be used after another constant is determined. The

circular error correction, which requires typing five constants and a subsequent additional line in the

calibration table, was developed by R.C. Jachens (written commun., 1981) to correct for assumed harmonic

deviations within dial revolutions of some gravity meters (appendixes 4 and 5). Lines of the calibration

table essentially reproduce the manufacturer's table by listing the counter reading, milligal-equivalent to that

reading, and the linear interpolation factor to multiply readings in the following interval. 
8

Step 4—Convert readings to milligals and correct for effect of tide

Program mergeg.f merges dates, times, readings, elevations, accuracy codes, and field terrain

corrections in the .FLD file with the geographic coordinates in the .FCT file and then calculates tidecorrected milligal-equivalent readings by applying file .CAL to create a .TID file (fig. 1 and appendix 5).

The program confirms which calibration factors were applied by printing, for example,

 Table of readings ranging from 3200 to 3600 units.

 MANUFACTURER'S REVISION

 COUNTER MGAL FACTOR MGAL FACTOR

 3200 3343.072 1.045901 3345.118 1.046541

 3300 3447.662 1.045959 3449.772 1.046600

Lines from a .TID file, with time and date as supplementary notes in italics, are, for example,

G8B 1993

WNBAS4058351174402 42979 62193 912 7 023623385N22 2C [AT 9:12, 6-21-93]

WINNA4054231174822 42980 62193 943 7 023606421 0 [AT 9:43, 6-21-93]

WN 6412458118 117 41270 621931551 7 023694775N22 1C [AT 15:51, 6-21-93]

WN 6412458118 117 41270 62293 842 7 023694741N22 1C [AT 8:42, 6-22-93]

WNBAS4058351174402 42979 622931912 7 023623401N22 2C [AT 19:11, 6-22-93]

After the first line with the gravity meter name and the year, lines include the station name, geographic

coordinates (degrees, minutes to 0.01 minute), elevation (to 0.1 foot) , date (month-day-year), local time

(hours, minutes), Greenwich Time correction (7), a flag ("2") for re-occupations, tide-corrected milligalequivalent readings (3623.385 to microgals for the first reading), optional accuracy code (“N22”), and field

terrain correction (to 0.01 milligal) with code ("C"). The gap after the Greenwich Time correction followed

by a zero previously included a value of observed gravity copied by program extractg.f from existing

processed data. Components of the tide-corrected milligal-equivalent readings are printed in file

MERGEG.PNT. 

Step 5—Utilize known or determined values of observed gravity

A key feature of this gravity reduction system is to create and recursively modify a file with assumed

absolute values of observed gravity and other options either by creating a .BAS file in an editor or by

executing program setbaseg.f (appendix 6). Inasmuch as a reading was observed at a new base station

WNBAS before a reading was made at a previously established base station WINNA, and program driftg.f

calculates values of observed gravity by linear interpolation between previously determined values of

observed gravity, a temporary value of 979,800.000 mGal for observed gravity at station WNBAS was

inserted into the .BAS file and consequently into a .TID file during the first execution by program

setbaseg.f. The first .BAS file consisted of one line:

WNBAS 3 79800.000 BASE

The status flag "3" indicates that station WNBAS is to be treated as a base station, to which other stations can

be tied. The value of 79800.000 (979,800.000 mGal) was substituted by program setbaseg.f for the value of

observed gravity at all occurrences of station WNBAS in a short .TID file that consisted of the first day of

readings at the two base stations. As is the case for programs and formats in this system, a leading "9" is

assumed for observed gravity to maximize computational accuracy . Program driftg.f next calculated the

best value in the least-squares sense for the observed gravity of station WINNA relative to WNBAS so that the

next execution of program setbaseg.f could substitute the consequently determined value for WNBAS as

well as for WINNA into the complete .TID data set. A value of 1.0 for observed gravity indicates that the

existing value of observed gravity is to be retained, but the status flag is to be changed. 

Other flags in addition to "3" for a base station are: "0" or "1" for a gravity station with only one

observation; "2," a station with re-occupations or a station identified in the .FLD file as station with a

calculated value of observed gravity to be printed; "4," a reoccupied station without output values; "5," a time

with an instantaneous change of gravity readings (a "tare"); "6," a duplicate base station observation to

obtain output; and "7," a reading to be deleted. Inasmuch as flags 1 and 2 occur as default flags in the first

.TID file, their stations are not in the .BAS file unless changes are needed. Flags 5, 6, and 7 in the .BAS

file require a following line with a date and time to be typed in an editor. For example, 
9

TONOJ 5

06-20-80 1535 TARE

 BELC 6

05-12-83 1644 BASE

After substituting values of observed gravity for the two base stations, a .GIV file is created by program

setbaseg.f, for example,

G8B 1993

WNBAS4058351174402 42979 62193 912 7798275133623385N22 2C

WINNA4054231174822 42980 62193 943 7798104833606421 0

WN 64124581180117 41270 621931551 7 023694775N22 1C

WN 64124581180117 41270 62293 842 7 023694741N22 1C

WNBAS4058351174402 42979 622931911 7798275133623401N22 2C

Program setbaseg.f creates a .GIV file, which is a duplicate of the input .TID file with new flags and

values of observed gravity inserted at all occurrences of stations identified in the .BAS file. The .GIV file is

renamed (at command level or in a window) to be the next .TID file needed as input to program driftg.f. 

Step 6—Calculate gravity meter drift and values of observed gravity

Program driftg.f is executed in two modes. The first mode is to print values of gravity meter drift and

observed gravity only for reoccupied stations or stations flagged for printed output in the .FLD file. When

the last values of observed gravity are established for all reoccupied stations with recursive execution of

programs setbaseg.f and program driftg.f, the second mode outputs a file with values of observed gravity

for all readings except at base stations in addition to a printout for all readings. A .DFT file of values

needed for plotting gravity meter drift in program driftplt.f is created in both modes. 

A detailed list in file DRIFTG.PNT includes the input line number, station name, geographic

coordinates, elevation, time, date, reading, base station reference, gravity meter drift, observed gravity, "*" to

flag an assumed gravity value, "P" to flag output for file .OBS, and the line number of the previous reading

at reoccupied stations. For example,

LINE NAME LATIT LONGIT ELEV TIME DATE READING BASE GRAVITY LAST

 DEG-MIN DEG-MIN (FEET) H-M M-D (MGAL) STA DRIFT (MGAL) LINE

 2 WNBAS 4058.35 11744.02 4297.9 912 621 3623.385 WNBAS 0.000 79827.51* 2

 3 WINNA 4054.23 11748.22 4298.0 943 621 3606.421 WINNA 0.066 79810.48* 3

 4 WNBAS 4058.35 11744.02 4297.9 1012 621 3623.517 WNBAS 0.132 79827.51* 2

 13 WN 4 4122.66 118 1.84 4171.9 15 4 621 3688.248 WNBAS 0.051 79892.32 P

 14 WN 5 4121.67 118 2.28 4166.7 1535 621 3689.416 WNBAS 0.050 79893.49 P

 15 WN 6 4124.58 118 1.17 4127.0 1551 621 3694.775 WNBAS 0.050 79898.85*P 15

The above detailed list is viewed to focus on doubtful values and may be retained as documentation but

seldom would need to be printed. The summary after the detailed list, which shows calculated average

values of observed gravity is useful to print as documentation. 

Program driftg.f was executed once to establish a value for the local base station occupied in 1993

(Plouff, 1996) and five additional executions—generally in the combined sequence of most-to-least

observations at reoccupied stations and chronologically outer-to-inner ties of reoccupied stations—to

obtain values of observed gravity for all gravity readings (appendix 6). With closer scrutiny of the detailed

printout from program driftg.f, less executions of the combination of programs setbaseg.f and driftg.f

may have been needed, but more user time would have been required. No sudden changes of gravity meter

readings ("tares"), which exceed the expected variation of re-occupations, were noticed in the gravity meter

drift curves for the 1993 data. The method of correction for tares would be to insert a flag at the dubious

reading and assume no gravity meter drift occurred since the last reoccupation or base station with a

previously determined value of observed gravity. 

The printed summary at the end of the first execution of driftg.f for an abbreviated file listed an

average of linearly interpolated value of observed gravity for base station WINNA:
10

 PAST SEQ NAME LAST TIMES GRAVITY STD

 NAMES LINE READ (MGAL) DEV

 1 WNBAS 8 5 979800.00 base

 2 WINNA 7 3 979782.97 0.02

Therefore, the value of observed gravity at local base station WNBAS is 17.03±0.02 mGal higher than at

regional base station WINNA . Inasmuch as the accepted value of observed gravity at station WINNA is

979,810.48 mGal, the value of observed gravity thus determined for newly established base station WNBAS

is 979,827.51 mGal. After inserting the updated values of observed gravity at the two base stations into the

.TID file for all the 1993 data by using program setbaseg.f, the first six lines of the printed summary from

program driftg.f for the 1993 data set were:

 PAST SEQ NAME LAST TIMES GRAVITY STD

 NAMES LINE READ (MGAL) DEV

 1 WNBAS 317 29 979827.51 base

 2 WINNA 8 3 979810.48 base

 3 WN 2 62 3 979888.25 0.03

 4 WN 6 294 8 979898.85 0.02

Accepted values of observed gravity for stations WN 2, WN 6, and other stations were entered into the

.BAS and .TID files by executing program setbaseg.f. Program driftg.f subsequently averaged linearly

interpolated values of observed gravity for reoccupied stations. If the standard deviation is not acceptable,

the detailed printout must be examined to eliminate typographical errors or misreadings. The standard

deviation and number of re-occupations are useful to identify a station with re-occupations as a new base

station or to assign an accuracy code for the value of observed gravity. 

Each execution of program driftg.f creates a .DFT file with calculated values of gravity meter drift for

base and other stations reoccupied in 1993 as input to plot program driftplt.f, for example,

 1-Apr-98 16:00:43 {The date later is expressed as "1998:0401"}

blm938

G8B 1993

 7 6 21 93 912

 7 6 21 912WNBASWNBAS+ 0.000

 7 6 21 943WINNAWINNA+ 0.066

 7 6 21 1012WNBASWNBAS+ 0.132

 7 6 21 1244WNBASWNBAS+ 0.048

 7 6 21 14 2WN 2WNBAS. 0.053

 7 6 21 1551WN 6WNBAS. 0.050

 7 6 21 1642WN 2WNBAS. 0.030

 7 6 21 18 9WNBASWNBAS+ 0.061

The first line of output file blm938.DFT has the date and time of execution of program driftg.f. The

second line has the root name of the data set, to which suffixes are appended. The third line has the gravity

meter name and the year of the gravity survey. The fourth line has the Greenwich Time correction and the

date and time of the first reading. The remaining lines include the Greenwich Time correction, the date

(without the year) and time of the reading for which the value of observed gravity has been determined, the

station name, the name of the base station to which the station is referred, a plot symbol, and a value of

gravity meter drift in milligals relative to the first reading of the data set. The final line of the .DFT file has

the Greenwich Time correction and the date (without the year) and time of the last reading of the data set.

Plot symbols are "+" for a base station, a square for a reoccupied station, and "X" for a tare or for the first

or last station of a survey without a value of observed gravity. 

After the last values of observed gravity have been entered into the .BAS file and subsequently into the

.TID file by program setbaseg.f, the second mode of program driftg.f is selected to output a .OBS file

with values of observed gravity for all readings except for unselected base stations, for example,

 WN 1 412388 1174945 421927987771N22 61993 1 C WNBAS

 WN 2 412445 1175781 412277988825N22 61993 1 C WNBAS
11

The output .OBS file includes station names, geographic coordinates (degrees and minutes to 0.1),

elevations (to 0.1 foot), values of observed gravity (to 0.01 milligal without the leading 9), the first three

digits of accuracy codes, months and years of the readings, close-in terrain corrections (to 0.01 milligal),

ring codes of the terrain corrections, and base stations to which stations are referenced. The base station to

which a reading at a station is referenced is assigned by program driftg.f to be the closest base station in

time read either before or after the reading of the station. 

Step 7—Plot gravity meter drift

Plots of gravity meter drift created by program driftplt.f (fig. 3; appendix 7) are useful to highlight

reading errors or errors of values of the assumed observed gravity, show the general error level for

re-occupations, and document the quality of gravity observations for a report. For example, Plouff (1964,

fig. 2) plotted the difference of gravity meter drift between two gravity meters read at the same isolated

location for five months to display error caused by a malfunctioning relay in one gravity meter for two

months, and to show that although the drift for one gravity meter was about 0.4 mGal/day, the long term

difference of drift was linear. Plouff (1982, figs. 4 and 5) plotted gravity meter drift to identify tares and to

document the accuracy of field data acquired by the U.S. Geological Survey and other groups. 

Drift plot files with suffixes .PST created by program driftplt.f are in the format of the commonly

used PostScript language (Adobe Systems Incorporated, 1990). Therefore, the plot files can be printed on

a laser printer or viewed on a monitor. Erroneous points on a drift plot can be identified in a .DFT file by

interpolating the time and amplitude of their occurrence. Large gaps in time and shifts in amplitudes are

interactively accommodated in program driftplt.f by splitting drift plots into two user-selected adjacent

panes. The ability to "click" on a point displayed on the monitor so that the station name, time, date, and

amplitude of the point appears as an annotation would be useful to implement in the future. 

Step 8—Delete redundant data

The .OBS output file from program driftg.f includes a line for each gravity reading except base

stations. After the values of observed gravity for each reoccupied station have been entered into the .BAS

file and the final execution of programs setbaseg.f and driftg.f have been made, only the last reading or

the reading with the closest tie to a base station should be retained for computation of gravity anomalies.

Obtaining a file with only one line per station was accomplished by first sorting file .OBS so that lines with

consecutive readings for each station are adjacent rather than mixed throughout the file. The UNIX (trade

name of American Telephone and Telegraph Co.) system "sort" command

sort +0.0 -0.8 +0.48 -0.52 +0.46 -0.48 blm938.OBS -o blm938.SRT

rearranged lines in file blm938.OBS into a sequence of lines arranged in the descending order of station

name, year, and month, respectively, to create file blm938.SRT that starts with lines

 P163 405500 1182281 416307984850G53 62893 1 C WNBAS

 P163 405500 1182281 416307984850G53 62893 1 C WNBAS

 P163 405500 1182281 416307984850G53 63093 1 C WNBAS

 AUX53 411469 1175988 550527980310G75 62393 2 C WNBAS

 AUX53 411469 1175988 550527980310G75 62593 2 C WNBAS

 WN 1 412388 1174945 421927987771N22 62193 1 C WNBAS

 WN 2 412445 1175781 412277988825N22 62193 1 C WNBAS

 WN 2 412445 1175781 412277988825N22 62193 1 C WNBAS

(appendix 8). For this simple data set, execution of program extraobs.f, which did not need interactive

decisions for assignment of alternative base stations, printed an EXTRAOBS.PNT file with lines

 3 stations named P163 Dates range from 61993 to 61993. Bases are:

 WNBAS ( 3)

 2 stations named AUX53 Dates range from 61993 to 61993. Bases are:

 WNBAS ( 2)

 3 stations named WN 2 Dates range from 61993 to 61993. Bases are:

 WNBAS ( 3)

and a user-named output file blm938.plf that starts with lines
Figure 3.--Drift plots. Top plot, drift after values for base stations were determined for 1993 data.

Middle plot, drift after values for all repeat stations were determined for 1993 data. Bottom plot,

inverted drift for constant gravity readings at hypothetical gravity station located at lat 12 degrees

N. and lon 80 degrees W. during a solar eclipse at local time 12:22 on August 21 in the year 2017.

12

BASE 

REPEAT 

END OR TARE COMPILATION: 1-Apr-98 16:00:43

GRAVITY METER: G8B 

DRIFT START: JUNE 1993

-0.2

0.0

0.2

MILLIGALS 

21 22 23 24 25 26 27 28 29 30

J U N E 

1 2

J U L Y 

BASE 

REPEAT 

END OR TARE COMPILATION: 2-Apr-98 08:14:39

GRAVITY METER: G8B 

DRIFT START: JUNE 1993

-0.2

0.0

0.2

MILLIGALS 

21 22 23 24 25 26 27 28 29 30

J U N E 

1 2

J U L Y 

BASE 

REPEAT 

END OR TARE COMPILATION: 27-Mar-98 13:16:04

GRAVITY METER: FAKE 

DRIFT START: AUGUST 2017

-0.3

0.0

0.1

MILLIGALS 

20 21 22

A U G U S T 
13

 P163 405500 1182281 416307984850G53 0 61993 1 0C 0 0WNBAS

 AUX53 411469 1175988 550527980310G75 0 61993 2 0C 0 0WNBAS

 WN 1 412388 1174945 421927987771N22 0 61993 1 0C 0 0WNBAS

 WN 2 412445 1175781 412277988825N22 0 61993 1 0C 0 0WNBAS

If rearrangement is acceptable, a combination of sorting and execution of program extraobs.f can identify

and delete superfluous data points from a file in the same format with more than one year of data. 

ADDITIONAL COMMENTS ABOUT THE REDUCTION PROCESS

Computer programs described in this report should not be merged into one program because each

program may identify errors to be corrected or may request additional information before the next

processing step is begun. Inasmuch as the gravity reduction system easily creates reproducible results, no

files except those with suffixes .FLD, .BAS, .CAL , and .OBS and a record of processing steps need to be

retained. For example, Plouff (1996, table 1 and Part B) published those files except for the .OBS file,

which was replaced by a new file with terrain corrections and gravity anomalies. If the value of observed

gravity assumed for a base station is changed, a program addog.f provides an option to increment the value

of observed gravity for gravity stations referred to that base station. If values of observed gravity are

determined with an inaccurate gravity meter and some of the gravity stations are re-occupied with a more

accurate gravity meter, best fits between repeated readings can yield a revised gravity meter constant for the

less accurate gravity meter. Therefore, improved values of observed gravity can be obtained for data

collected with the less accurate gravity meter. For example, Plouff (1972, p. 7) utilized 72 re-occupations to

re-compute values of observed gravity for 231 gravity stations. 

Program driftg.f would need to be modified in order to improve accuracy by utilizing parabolic drift,

for example (written commun., R.C. Jachens and C.W. Roberts, 1992), rather than the presently assumed

linear instrument drift. Gravity meter drift is not directly linked across the time interval between New Years

Eve and New Years Day, but tests utilizing output from programs mergeg.f and driftg.f confirmed that

there are no discontinuities of tide-corrected instrument drift in that interval. 

GRAVITY TERRAIN CORRECTIONS

Unless reflectorless laser range finding equipment and software (Lyman and others, 1997) are utilized,

notes about topography should be recorded in the field so that gravity terrain corrections (TCs), the gravity

effect of departures of topography from an assumed flat surface, can be determined. Inasmuch as TCs are

determined to farther distances by overlaying topographic maps with "templates" or by utilizing digital map

elevations (Plouff, 1977; Cogbill, 1990), field TCs commonly are determined to distances of 53.3 m (175

ft) from the gravity station location to include the C ring of the Hammer (1939) TC system or 68 m (223 ft)

to include the B ring of the Hayford TC system (Hayford and Bowie, 1912, p. 18; Swick, 1942). If a point

with a known elevation is located at a site where the local terrain effect is high and a description of

topography is complicated, the gravity station preferably should be located nearby where the effect of

topography is minimized, and the elevation is adjusted to the station location by leveling. TCs for the

special case of gravity readings taken on a post are computed by program tcpost.f (appendix 9). 

TCs first were calculated at the U.S. Geological Survey by recording elevation estimates on standard

forms, subtracting the station elevation, looking up the TC for each compartment in tables, for example,

from Hammer (1939) or Swick (1942, p. 68), and summing the values. TCs for compartments designed to

match local topography can be calculated in the field with a programmable calculator (oral commun., C.W.

Roberts, 1985). Spielman and Ponce (1984 compute TCs directly from table entries in digital format. The

computations are based on the assumption that the topography within a compartment is horizontal. These

TCs commonly are too large because topography near a station tends to slope toward the station—with

elevations closer to that of the station neart the closer edges of compartments—rather than having a

stonehengelike form. In view of knowledge of the effect of the inverse square of distance, experienced

estimators subjectively recorded elevations closer to that of the station than either average elevations or

elevations near the centroid of compartments (oral commun., Kenneth Holden, 1975). This error can be

minimized by using smaller compartment sizes or by utilizing formulas for sloping compartments

(Campbell, 1980). 
14

Close to a gravity station, the ground surface more likely slopes through the gravity station as a first

approximation rather than forming a flat-topped stonehengelike horizontal slope throughout a TC

compartment. Therefore, Hammer (1939) derived a formula to calculate the two-dimensional gravity effect

of a sloping surface (two-dimensional wedge bounded by a circular ring), which passes through the gravity

station. The formula was implemented by Sandberg's (1958) tables for Hammer and Hayford rings with

TCs determined at slope intervals of one degree to a maximum slope of 30 degrees . Entries for

Sandberg's (1958) tables either are angles of slope (measured in the field) or differences of elevation

between the station and the interpolated contour elevation tangent to the outer edge of a semicircular area

centered at a gravity station. Campbell (1980) extended Hammer's (1939) slope formula to determine TCs

bounded by individual cylindrical compartments rather than semicircular rings. 

S.L. Robbins and H.W. Oliver (written commun., 1970) computed TCs for displaced slopes, in which

the ground surface is horizontal to a measured distance from a gravity station and, thence, an assumed

planar ground surface slopes upward or downward at a measured angle of slope. Their page-size plot for

TCs of displaced slopes still is useful to estimate field TCs for most gravity stations. Inasmuch as the

unpublished computer algorithm used by Robbins and Oliver was not found, program dspslope.f was

written to compute and plot curves for the TCs of displaced slopes (figs. 4 and 5). Families of curves show,

for example, that TCs to a distance of 68 m are 0.005 mGal if slopes on both sides of a station are 3

degrees or if the ground is flat to 9 ft (3 m) and thence the slope on one side of a gravity station is 5

degrees (fig. 4). The algorithm to reproduce the logarithmic plot of Robbins and Oliver is based on a

summation of 44 or less vertical rectangular prisms (Plouff, 1976), which are bounded by a vertical halfcylinder with a user-selected radius of 68 m (223 ft) or less. The exact TC for an inclined plane that passes

through the station location also is computed by program dspslope.f (figs. 4 and 5). Computation of the

exact TC is based on evaluation of a complete elliptical integral of the first kind (Bulirsch, 1965, p. 86). 

Program fieldtc.f computes TCs for more complex forms of two-dimensional topography (appendix

9). The TCs are computed as the gravity effect of 87 or less vertical rectangular prisms that start at a

distance of 0.3 ft (0.09 m) from the gravity station. TCs commonly are calculated within semicircles with

outer radii of 175 ft (53.3 m) for Hammer's (1939) C ring or 223 ft (68 m) for Hayford's (Swick, 1942) B

ring. Interactive program input consists of three components measured from a succession of points

starting at the gravity station: 1) angle of slope in degrees viewed at or parallel to ground level; 2) tapemeasured or usually paced distance along the slope in feet; and 3) height difference in feet. One of the

three entries always is zero to signify that the component was not measured or is a redundant component .

Field notes to describe nearly two-dimensional topography first should be written in the direction of

maximum slope, then in the opposite direction, and finally in other directions if needed to estimate the

maximum effect normal to the steepest slope or for defining three-dimensional topography. Program

entries can be made to obtain TCs for simple cases of slopes or displaced slopes. For a 3-degree slope in

one direction from a station, the entry to the program is the triplet of numbers "3 0 0" (regardless of

upward or downward direction of slope), which yields a TC of 0.003 mGal to a distance of 68-m. For a 5-

degree slope displaced at a distance of 9 ft, the entries are "0 9 0" followed by "5 0 0," which yields a TC

of 0.005 mGal to a distance of 68 m. Field notes for a station may be "flat to 10' NW; down 90° for 6'; flat

to 15'; up 90° for 6'; flat 8'; up 15°." The corresponding triplets of interactive entries to program fieldtc.f

are "0 10 0," "-90 0 -6," "0 15 0," "90 0 6," "0 8 0," and "15 0 0," which yields a TC of 0.034 mGal to a

distance of 68-m for the semicircular area to the northwest of the station. The TC for the contribution

associated with each triplet is printed in file FIELDTC.PNT. Using either 44 or 87 prisms yielded a TC of

0.118 mGal for an interrupted 20° slope ("20 5 0," "20 0 0") compared to an exact solution of 0.117 mGal

("20 0 0"). 

Barrows and Fett (1991) described a TC method for field or map application with a three-dimensional

model that consists of a series of sloping "wedges" bounded in plan view by radial lines. The base of the

wedges is a horizontal plane at the elevation of the gravity station, and the upper surface of each wedge is

conical (Olivier and Simard, 1981). Wedges join along the upper surfaces at changes of slope with heights

and distances measured along wedge bisectors in the field or obtained from maps. Lyman and others

(1997) described a method to rapidly "shoot" a reflectorless laser range finding gun at key points such as

slope changes, maxima, and minima. The input data are automatically recorded, and TCs, for example, to

50 m can be computed in the field by using apparently proprietary software partly based on algorithms

referenced by Cogbill (1990). 
Figure 4.--Gravity terrain correction of a displaced slope (Hayford system). Slope is

horizontal from gravity station to beginning of an inclined plane. Plus sign is located at

terrain correction for zero distance to base of slope.

15

DISTANCE TO BASE OF SLOPE, IN FEET

ANGLES OF SLOPE ALONG CURVES ARE IN DEGREES. OUTER RADIUS= 223.1 FEET, HAYFORD B

1

0.002

2 3 4 5 6 7 8 9 10 15 20 30 40 50 60 70 8090100 150 200 250

0.005

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.10

0.15

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.0

1.5

2.0

3.0

4.0

1 2 3 4 5 6 7 8 9 10 15 20 30 40 50 60 70 8090100 150 200 250

 3

 4

 5

 6

 7

 8

 9

10

12

15

20

25

30

35

40

45

50

55

60

65

70

75

80

85

TERRAIN CORRECTION, IN MILLIGALS
Figure 5.--Gravity terrain correction of a displaced slope (Hammer system). Slope is

horizontal from gravity station to beginning of an inclined plane. Plus sign is located at

terrain correction for zero distance to base of slope.

16

DISTANCE TO BASE OF SLOPE, IN FEET

ANGLES OF SLOPE ALONG CURVES ARE IN DEGREES. OUTER RADIUS= 175.0 FEET, HAMMER C

1

0.002

2 3 4 5 6 7 8 9 10 15 20 30 40 50 60 70 8090100 150 200 250

0.005

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.10

0.15

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.0

1.5

2.0

3.0

4.0

1 2 3 4 5 6 7 8 9 10 15 20 30 40 50 60 70 8090100 150 200 250

 3

 4

 5

 6

 7

 8

 9

10

12

15

20

25

30

35

40

45

50

55

60

65

70

75

80

85

TERRAIN CORRECTION, IN MILLIGALS
17

PROGRAMMING CONSIDERATIONS

The computer programs discussed in this report are written in the FORTRAN-77 language and are

compiled and executed in a UNIX (trade name of American Telephone and Telegraph Co.) system.

FORTRAN-77 is a scientific programming language, which is commonly referenced, for example, in

algorithms published by the Association for Computing Machinery (ACM) Transactions on Mathematical

Software. The FORTRAN language also is used for many other USGS geophysical applications (for

example, Cordell and others, 1992; Blakely, 1995). The generic level of FORTRAN coding and selfdocumentation by extensive prompting simplifies adaptation of programs discussed in this report to other

platforms and conversion of the codes to other languages. Computer files for plots and forms were

prepared in the commonly used and available PostScript language (Adobe Systems Incorporated, 1990) so

that the software would not be dependent on the availability of specialized plot equipment and software. 

Current dates and times recorded during the reduction process were derived from VAX/VMS (Digital

Equipment Corporation) date and time functions. The VMS format for a date, for example, expressed as

"14-Apr-98" later was replaced by a standard sequence of year-month-day expressed as "1998:0414,"

which was derived from a Sun Microsystems, Inc., Solaris date_and_time function. The colon was added

for readability. A subroutine, datetim.f, which calls and converts arguments from the Solaris Y2Kcompatible date_and_time function, replaced calls to the VMS date and time functions. 

The following features reflect the programming approach used in developing the computer programs

discussed in this report. Standard FORTRAN-77 functions called by the programs include the INQUIRE

command, which searches a directory to find if a requested file is present or absent, and an option of the

READ and WRITE commands to re-format variables. A "bz" command in READ statements protects against

carriage returns embedded in lines before all possible variables are read. A carriage return (ENTER key) in

response to prompting is a “yes.” Files generally are not overwritten during program execution except

with user permission in response to prompting. Accuracy of computations are maximized throughout

processes by utilizing a format with integer degrees separated from minutes. Although tide corrections

expressed in microgals require only three significant figures, tide corrections are calculated in double

precision because of the need to express time in minutes since the beginning of the year 1900. "Print" files

generally are only viewed on a computer screen. Print files with line lengths greater than the common 80-

column page width require rotation on a printed page to be retained as a "hard copy" for quality control.

Programs in this report perform essential tasks for which they are designed, but comprehensive testing

has not been done for a full range of input data. 

REFERENCES

Adobe Systems Incorporated, 1990, PostScript Language reference manual, 2nd ed.: New York, AddisonWesley Publishing Company, Inc., 764 p.

Barnes, D.F., Oliver, H.W., and Robbins, S.L., 1969, Standardization of gravimeter calibrations in the

Geological Survey: Eos (American Geophysical Union), v. 50, no. 10, p. 526-527.

Barrows, L.J., and Fett, J.D., 1991, A sloping wedge technique for calculating gravity terrain corrections:

Geophysics, v. 56, no. 7, p. 1061-1063.

Blakely, R.J., 1995, Potential theory in gravity and magnetic applications: New York, Cambridge University

Press, 441 p.

Broucke, R.A., Zürn, W.E., and Slichter, L.B., 1972, Lunar tidal acceleration on a rigid earth, in Flow and

Fracture of Rocks, H.C. Heard, I.Y. Borg, N.L. Carter, and C.B. Raleigh (eds.): American Geophysical

Union Geophysical Monograph 16, p. 319-324.

Bulirsch, Roland, 1965, Numerical calculation of elliptic integrals and elliptic functions: Numerische

Mathematik, Handbook series special functions, v.7, no. 1, p. 78-90.

Campbell, D.L., 1980, Gravity terrain corrections for stations on a uniform slope—A power law

approximation: Geophysics, v. 45, no. 1, p. 109-112.

Cogbill, Allen, 1990, Gravity terrain corrections using Digital Elevation Models: Geophysics, v. 55, no. 1, p.

102-106.
18

Cordell, Lindrith, Phillips, J.D., and Godson, R.H., 1992, USGS potential-fields geophysics software,

version 2.0: U.S. Geological Survey Open-File Report 92-018, Parts A-G, 18 p., 6 diskettes.

Goguel, Jean, 1954, A universal table for the prediction of the lunar-solar correction in gravimetry (tidal

gravity correction): Geophysical Prospecting (Journal of the European Association of Exploration

Geophysicists), v. 2, Supplement, p. 2-5.

Hammer, Sigmund, 1939, Terrain corrections for gravimeter stations: Geophysics, v. 4, no. 3,

p. 184-194.

Hayford, J.F., and Bowie, William, 1912, The effect of topography and isostatic compensation upon the

intensity of gravity: U.S. Coast and Geodetic Survey Special Publication no. 10, 132 p.

Longman, I.M., 1959, Formulas for computing the tidal accelerations due to the moon and the sun: Journal

of Geophysical Research, v. 64, no. 12, p. 2351-2355.

Lyman, G.D., Aiken, C.L.V., Cogbill, Allen, Balde, Mamadou, and Lide, Chet, 1997, Terrain mapping by

reflectorless laser rangefinding systems for inner zone gravity terrain corrections: 67th Annual

International Meeting, Society of Exploration Geophysicists, Expanded Abstracts, p. 482-485.

Olivier, R.J., and Simard, R.G., 1981, Improvement of the conic prism model for terrain correction in rugged

topography: Geophysics, v. 46, no. 7, p. 1054-1056.

Plouff, Donald, 1964, Gravity measurements in the Beaufort Sea area: Arctic, v. 17, no. 3, p. 150-161.

_____1972, Gravity observations near the San Juan Mountains, Colorado: National Technical Information

Service, U.S. Department of Commerce, Springfield, VA 22161, PB-72-211430, 77 p.

_____1976, Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections:

Geophysics, v. 41, no. 4, p. 727-741.

_____1977, Preliminary documentation for a FORTRAN program to compute gravity terrain corrections

based on topography digitized on a geographic grid: U.S. Geological Survey Open-File Report

77-535, 45 p.

_____1982, Gravity observations in the Walker Lake 1° x 2° quadrangle, California-Nevada: U.S.

Geological Survey Open-File Report 82-405, 105 p.

_____1990, Principal facts for gravity data compiled for the Conterminous United States Mineral Appraisal

Program, Tonopah 1- by 2-degree quadrangle, Nevada: U.S. Geological Survey Open-File Report

90-457; Part A, documentation, 22 p.; Part B, tables, 61 p.; Part C, diskette.

_____1996, Principal facts and field observations for gravity data in and adjacent to the Bureau of Land

Management's Winnemucca District and Surprise Resource Area, northwest Nevada and northeast

California: U.S. Geological Survey Open-File Report 96-290; Part A, 26 p.; Part B, diskette.

_____1998, Computer programs to display and modify data in geographic coordinates and methods to

transfer positions to and from maps, with applications to gravity data processing, global positioning

systems, and 30-meter digital elevation models: U.S. Geological Survey Open-File Report 98-233,

43 p.

Sandberg, C.H., 1958, Terrain corrections for an inclined plane in gravity computation: Geophysics,

v. 23, no. 4, p. 701-711.

Service Hydrographique de la Marine and Compagnie Générale de Géophysique, 1969, Tidal gravity

corrections for 1970: Geophysical Prospecting (Journal of the European Association of Exploration

Geophysicists), v. 17, Supplement No. 1, p. 1-53.

Spielman, J.B., and Ponce, D.A., 1984, HANDTC, A FORTRAN program to calculate inner-zone terrain

corrections: U.S. Geological Survey Open-File Report 84-777, 20 p.

Swick, C.H., 1942, Pendulum gravity measurements and isostatic reductions: U.S. Coast and Geodetic

Survey Special Publication 232, 82 p.
19

______________________________________________________________________________________________________________________________

APPENDIX 1

ACCURACY CODES ______________________________________________________________________________________________________________________________

A 4-digit code was defined by group consensus in Menlo Park, Calif., in the 1960's to describe the

accuracy of individual gravity stations. Slight modifications were made by data collectors and compilers to

tailor the usage to the requirements of specific data sets. The first digit serves to describe the location,

usually referred to a point of known elevation. 

Code

B On level-line bench mark or other permanent mark incorporated into USGS vertical control system,

including National Geodetic Survey bench marks.

N Near or possibly on level-line bench mark of USGS system.

V On vertical-angle bench mark.

H Near or possibly on vertical-angle bench mark.

D Near assumed location of any of the above marks that was destroyed or not found.

P Near surveyed elevation not printed on topographic map.

X Near well-defined marks such as wells, windmills, microwave towers, or section corners.

F Near location with or without a mark, at which a surveyed elevation is shown on a published map.

G Near location (on a manuscript or published map) at which a spot elevation is determined by

photogrammetry; near doubtful F-location. Printed in brown or with letter "T" appended.

T Elevation based on photogrammetry done by USGS; or spot elevation not printed on published map.

K Spot elevation based on photogrammetry determined by other organizations.

W Edge of lake, canal, or reservoir; interpolated elevation or elevation given for water or dam at unknown

height relative to present level.

A Elevation determined by using altimetry with unknown accuracy.

C Elevation determined by topographic contour interpolation.

Q Elevation determined by topographic contour interpolation along drainages.

2 Along regularly spaced profile of stations; elevation source is unknown.

3 Near road or section line; elevation source is unknown.

U Not near road or other landmark; elevation source is unknown.

(blank space) Code not assigned .

The second digit serves to estimate the effect of elevation error and the consequent effect on the

accuracy of the value of gravity anomalies. The second column is the elevation accuracy in meters. The

third column is the consequent error in calculating gravity anomalies, in milligals, assuming 0.02 mGal/m. 

Code

1 0.05 0.01 On or tied to level-line bench mark by surveying.

2 0.15 0.03 Elevation difference hand-leveled to nearby bench mark; elevation recorded to nearest foot.

3 0.3 0.06 Near bench mark.

4 0.6 0.12 On or near vertical-angle bench mark; flat area near level-line bench mark that was not found.

5 1.5 0.3 Near surveyed elevation on topographic map; elevation from map with 10-ft contour interval

6 3.0 0.6 Photogrammetric elevation or contour interpolation on map with 20-ft contour interval.

7 6.0 1.2 Uncertain location of photogrammetric spot elevation.

8 15.0 3.0 Contour interpolation along road or stream on map with 80-ft contour interval; doubtful

combination of elevation and horizontal location.

9 30.0 6.0 Location and elevation doubtful but no other data points nearby.

(blank space) Code not assigned.

The third digit serves to estimate the effect of location error and the consequent effect on the accuracy

of the value of gravity anomalies. The second column is the location accuracy in meters. The third column

is the consequent error in calculating gravity anomalies, in milligals, assuming that all the location error is

along the north component of direction. 

Code

1 13 0.01 Near vertical-angle bench mark; established by Global Positioning System..

2 25 0.02 Near permanent mark on map such as bench mark, section corner, or well.

3 65 0.05 Road intersection or stream fork.

4 130 0.1 Broad road curve or gentle hill crest.

5 250 0.2 Location depends on odometer for interval greater than 1 mi, or other estimate.
20

6 650 0.5 Location consistent with elevation but located 0.5 mi from a definitive landmark.

7 1,300 1.0 No likely landmark within 1 mile.

8 2,500 2.0 Original location changed by 1 to 2 mi, to agree with location of the given elevation.

9 5,000 5.0 Doubtful location but no other data points nearby.

(blank space) Code not assigned.

The fourth digit, the accuracy of the value of observed gravity, is not available until the field data are

processed. The accuracy estimate, in milligals, indicated in the second column, usually is relative to the

closest available regional gravity base station. 

Code

1 0.01 Base station established with LaCoste and Romberg or equivalent gravity meter.

2 0.02 Station established with multiple ties of high precision.

3 0.05 Repeated readings with LaCoste and Romberg gravity meter.

4 0.1 One reading with LaCoste and Romberg gravity meter.

5 0.2 Single tide-corrected reading with Worden gravity meter with good drift characteristics.

6 0.5 Value is average of two observations that are 1.0 mGal apart; uncertain datum shift.

7 1.0 Single reading of Worden gravity meter with doubtful drift characteristics.

8 2.0 Value is average of two observations that are 4 mGal apart.

9 5.0 Doubtful value but no other data available.

(blank space) Code not assigned.

______________________________________________________________________________________________________________________________

APPENDIX 2

ASSEMBLING MAPS , TESTING DATA, AND LISTING DATA POINTS ______________________________________________________________________________________________________________________________

The following computer session for program fieldg.f simulated a test of file blm938.FLD for gravity

data collected in 1993 (Plouff, 1996) and created files blm938.NAM and blm938.MAP. In practice,

previous executions of fieldg.f may have indicated typographical errors to be corrected. Typed responses

are in bold print. Note that the root name for files must have eight or less characters. Additional comments

not in the original session are bracketed in italics. After datetim.f was appended to program fieldg.f and

other programs, the date in the format "18-Mar-98" was replaced by "1998:0318 ."

fieldg [COMMAND TO EXECUTE PROGRAM]

 FIELDG, Plouff, 3-98. 18-Mar-98 16:42:56 Program to test field readings and

 to create files of unique station names and maps (west longitudes).

 A data file consists of a lead line with gravity meter name (1-4),

 year (5-9), positive Greenwich correction (10-11), terrain correction

 ring (13), and method (14-20). Days of field data follow. A day consists of

 a line with the date (month-day-year, 1-8), reading lines, and a line with

 all nines as a delimiter. The reading format is: station name (1-5), a code

 letter, time (hours-minutes, 7-10), microgal reading (11-18), elevation

 (19-24), elevation units, accuracy code (26-28), field terrain correction

 (29-31), and optional map name (32-36).

 TYPE the root name associated with your files (8-digit max): blm938

 Field data are in a file called blm938.FLD

 Station names are output in a file called blm938.NAM (EXTRACTG input)

 Map names are output in a file called blm938.MAP

 Print output is in a file called: FIELDG.PNT

 Gravity meter: G8B ; Year: 1993; Greenwich correction (+): 7 hours.

 Outer ring for field TC estimate is HAMMER -C

 Next, sort and print file blm938.MAP that lists maps to be digitized.

 Then run program EXTRACTG to find locations and existing observed gravity

 values for file: blm938.NAM You need to provide a calibration file named

 G8B.CAL when you run program MERGEG. Run program CALIBG if you need to

 create that file.

 7- 2-93 is last day read. Last line is number 341

 Readings range from 3397.570 at station WN115 to 3547.440 at station WN249

 Total days= 12 Total unique stations= 263

 Total input lines= 341 Total stations repeated from another survey= 18

 Total maps named= 38
21

[PRINT OF FILE FIELDG.PNT]

 FIELDG, Plouff, 3-98. 18-Mar-98 16:42:56 Program to test field readings and

 to create files of unique station names and maps (west longitudes).

 Field data are in a file called blm938.FLD

 Station names are output in a file called blm938.NAM (EXTRACTG input)

 Map names are output in a file called blm938.MAP

 Gravity meter: G8B; Year: 1993; Greenwich correction (+): 7 hours.

 Outer ring for field TC estimate is HAMMER -C

 Next, sort and print file blm938.MAP that lists maps to be digitized.

 Then run program EXTRACTG to find locations and existing observed gravity

 values for file: blm938.NAM You need to provide a calibration file named

 G8B.CAL when you run program MERGEG. Run program CALIBG if you need to

 create that file.

 7- 2-93 is last day read. Last line is number 341

 Readings range from 3397.570 at station WN115 to 3547.440 at station WN249

 Station list:

 NUMBER NAME NUMBER NAME NUMBER NAME NUMBER NAME NUMBER NAME NUMBER NAME

 29 WNBAS 3 WINNA 1 WN 1 3 WN 2 1 WN 3 1 WN 4

 1 WN 5 9 WN 6 1 WN 7 2 WN 8 1 WN 9 1 WN 10

 2 WN 11 1 WN 12 2 WN 13 1 WN 14 1 WN 15 1 WN 16

 1 WN 17 1 WN 18 1 WN 19 2 WN 20 1 WN 21 1 WN 22

[OTHER LINES]

 4 P163 1 WN154 1 WN155 1 WN156 1 WN157 1 WN158

 2 WN243 1 WN244 1 WN245 1 WN246 1 WN247 1 WN248

 2 WN249 1 WN250 1 WN251 1 WN252 1 WN253 1 WN254

 1 WN255 1 WN256 1 WN257 1 WN258 2 OROV

 Map list: ["NUMBER" BELOW REFERS TO NUMBER OF STATIONS ON THE MAP]

 NUMBER NAME NUMBER NAME NUMBER NAME NUMBER NAME NUMBER NAME NUMBER NAME

 1 WD7NW 2 WD8NE 1 MB8NE 5 MB8NW 33 VB1NE 13 VB1SE

 5 VB1NW 4 VB2NE 3 VB2SE 3 VB1SW 24 VA1NW 15 VA1NE

 24 MB8SW 12 MA8NW 3 MB8SE 7 WD8NW 10 LD1NW 13 VA1SW

 5 VA1SE 15 VA2SE 3 LD2NE 3 MA8SW 3 LD1NE 3 LD2NW

 14 VA2NE 14 VC1SE 1 MC8SW 1 VA2NW 1 MA7NW 2 MA7SW

 2 MA8SE 6 MA8NE 1 VC1SW 1 VC1NE 5 VC2NE 1 VC2SE

 3 VC2SW 1 MC8SE

 Total days= 12 Total unique stations= 263

 Total input lines= 341 Total stations repeated from another survey= 18

 Total maps named= 38

sort blm938.MAP -o blm938.SRT [COMMAND TO SORT FILE]

[PRINT OF FILE blm938.SRT WITH LINES COMBINED]

 Maps listed in file blm938.FLD

LD1NE 3 LD1NW 10 LD2NE 3 LD2NW 3 MA7NW 1 MA7SW 2 MA8NE 6 MA8NW 12

MA8SE 2 MA8SW 3 MB8NE 1 MB8NW 5 MB8SE 3 MB8SW 24 MC8SE 1 MC8SW 1

VA1NE 15 VA1NW 24 VA1SE 5 VA1SW 13 VA2NE 14 VA2NW 1 VA2SE 15 VB1NE 33

VB1NW 5 VB1SE 13 VB1SW 3 VB2NE 4 VB2SE 3 VC1NE 1 VC1SE 14 VC1SW 1

VC2NE 5 VC2SE 1 VC2SW 3 WD7NW 1 WD8NE 2 WD8NW 7

The above list can be used to extract maps from the field files so that locations can be digitized and

manual gravity terrain corrections can be made. If there are less than four stations on a map, hand

digitization of geographic coordinates with templates (Plouff, 1998) is faster than machine digitization. 
22

______________________________________________________________________________________________________________________________

APPENDIX 3

OBTAINING GEOGRAPHIC COORDINATES TO MERGE WITH FIELD DATA ______________________________________________________________________________________________________________________________

The following computer session for program extractg.f simulated step 2 (fig. 1) for finding geographic

coordinates for the 263 stations listed in file blm938.NAM from program fieldg.f. In the original

execution of this step, locations were found in files obtained from hand and machine digitization. The

digitization files, however, were no longer needed after merging was accomplished. Therefore, a previously

processed output file for the gravity survey (Plouff, 1996) was temporarily renamed blm938.LOC required

to simulate a first execution of program extractg.f. The last five digits of an 8-digit station name in the file

to be searched form the 5- digit station name carried throughout programs described in this report until the

output of step 8. Typed responses are in bold print. Additional comments not in the original session are

bracketed in italics.

extractg [COMMAND TO EXECUTE PROGRAM]

 EXTRACTG (Plouff, 8-94). Program to extract locations to be merged with

 field data. Files with newly digitized data (.LOC) or old files are searched

 for data points with matching unique 5-digit station names (.NAM)

 generated by the FIELDG program.

 TYPE the root name associated with your files (8-digit max): blm938

 Field readings are in blm938.FLD ; Station names are in blm938.NAM

 Is this a first run to extract data for blm938.FLD ? ["y" NOT NEEDED]

 Do you want to obtain digitized locations from file blm938.LOC now?

 One digitizer format is: 3 blank spaces, 5-digit name, 7 spaces for latitude,

 and 8 for longitude.

 Is this your format?

 Total of 263 unique station names in blm938.NAM

 261 station names were matched.

 2 stations now in blm938.LFT still need to be found.

 The next successful run of EXTRACTG will create a file blm938.REF , to be

 appended to file blm938.FCT

 Replace file blm938.NAM with blm938.LFT , and re-run EXTRACTG.

 The print file from this program is: EXTRACTG.PNT

[PRINT OF FILE blm938.LFT--LIST OF STATIONS NOT FOUND]

WNA11 2

 P163 4

[PRINT OF FILE EXTRACTG.PNT]

 EXTRACTG (Plouff, 8-94) 16-Jul-98 16:23:21 {Date later "1998:0716"}

 Field readings are in blm938.FLD ; Station names are in blm938.NAM

 You are creating file blm938.FCT for the first output file.

 You specified a digitizer format with the station name in the last 5 of

 8-digit names in file blm938.LOC

 Total of 263 unique station names in blm938.NAM

 261 station names were matched.

 2 stations now in blm938.LFT still need to be found.

 The next successful run of EXTRACTG will create a file blm938.REF , to be

 appended to file blm938.FCT

 Replace file blm938.NAM with blm938.LFT , and re-run EXTRACTG.

 List of missing stations in file blm938.LFT :

 WNA11 P163

[PRINT OF FIRST LINES OF OUTPUT FILE blm938.FCT]

 P1634055001182281 0 4 2 0

WN 664057241180755 0 1 1 0

WN 674058691180836 0 1 1 0

WN 684059131180752 0 1 1 0

WN 694059561180833 0 1 1 0
23

Each line of output file blm938.FCT has a 5-digit station name, latitude and longitude in degrees and

minutes, a position for elevation (zeroes), the number of readings (from file blm938.NAM), a flag "2" for

reoccupied locations or station that will need later printed output or a "1," a previously used position for

observed gravity, and a position for a terrain correction (zeroes). Only the two data points listed in file

blm938.LFT, which has the same format as file blm938.NAM, still need to be found. Station P163, as

indicated on the field sheet, is a reoccupation of a previously established gravity station. Before re-running

program extractg.f to find station P163, the name of file blm938.LFT needs to be changed to

blm938.NAM so that file blm938.LFT will have the same role as file blm938.NAM had as input to

program extractg.f. If desired, first file blm938.NAM can be saved. 

mv blm938.NAM blm938.NAM1 [SAVES THE FIRST .NAM FILE BY RE-NAMING]

mv blm938.LFT blm938.NAM [RE-NAMES THE .LFT FILE TO .NAM]

The "plouff" format for station P163 found in an existing file is:

 P163 405500 1182281 416307984851G533 -1088 71284 2 7S-15409 13WINNA

Only key lines from the following computer session are printed.

extractg [COMMAND TO EXECUTE PROGRAM ]

 A print file EXTRACTG.PNT already exists.

 Do you want to stop to save it by renaming it? n

 TYPE the root name associated with your files (8-digit max): blm938

 Field readings are in blm938.FLD ; Station names are in blm938.NAM

 Is this a first run to extract data for blm938.FLD ? n

 TYPE the name of your file that provides station locations: (<66 cols):

/gp/plouff/blmwinne/lovelock/plouflov.isn

 Two data formats provide locations only, and one also supplies elevations

 and hand terrain corrections. Their formats are:

 1. Location appended to 5-digit station name.

 2. Plouff format with station name in columns 4-8.

 3. Same plouff format with elevation, handTC, and ring code carried forward.

 TYPE number 1, 2, or 3: 2

 Total of 2 unique station names in blm938.NAM

 1 station names were matched.

 1 stations now in blm938.LFT still need to be found.

 Append file blm938.REF to file blm938.FCT , to update blm938.FCT

 Replace file blm938.NAM with blm938.LFT , and re-run EXTRACTG.

[PRINT OF OUTPUT FILE blm938.REF]

 P1634055001182281 0 4 2 0

Step 4 of the reduction process requires a file with a .FCT suffix (fig. 1). Therefore, later output files

with .REF suffixes must be appended to the .FCT file. This can be accomplished with text editing or with

system level commands to append or concatenate files, for example,

cat blm938.FCT blm938.REF > tempfile [CONCATENATE .FCT AND .REF FILES]

rm blm938.FCT [DELETE FIRST .FCT NAME]

mv tempfile blm938.FCT [RENAME MERGED FILES TO .FCT FOR STEP 4]

Station WNA11 still was not found. Inasmuch as the reading at station WNA11 was an attempted

reoccupation of station WN 11 in a flat area with an unexpectedly high gravity gradient, station WNA11 was

inserted into file blm938.FCT by copying station WN 11 in a text editor and assigning flags "2" for

reoocupied locations. The geographic range of the data can be obtained at the end of this reduction step by

executing program geog5.f, with file blm938.FCT as input. 
24

______________________________________________________________________________________________________________________________

APPENDIX 4

CREATING A GRAVITY METER CALIBRATION FILE ______________________________________________________________________________________________________________________________

The following computer session for program calibg.f created a calibration table needed for the gravity

data collected in 1993 (Plouff, 1996). The required range of gravity meter readings was printed as output

from the fieldg.f program (" Readings range from 3397.570 at station WN115 to 3547.440 at

station WN249"). In the following example, a reading range beyond what was needed was input, and

values were converted to an accuracy of 0.001 mGal by using the manufacturer's interval factors. Typed

responses are in bold print. Additional comments not in the original session are bracketed in italics.

calibg [COMMAND TO EXECUTE PROGRAM]

 CALIBG, Plouff, 3-98. Create a gravity meter calibration file from

 manufacturer's table.

 TYPE a 1- to 4-digit name for the gravity meter (left-adjusted): G8B

 Your calibration file will be named: G8B.CAL

 IF TYPING ERRORS FOR NUMBERS OCCUR, QUESTIONS WILL BE REPEATED.

 An experimentally determined constant may be multiplied times the instrument

 calibration factors (in milligals per dial unit) provided by the manufacturer.

 TYPE constant correction factor (about 1.): 1.000612

 The counter reading interval is 100 dial units for LaCoste-Romberg G-meters

 and 10 or 50 for D-meters. An interval of zero signifies gravity meters

 where a single constant applies.

 TYPE counter reading interval (integer) (0 for single-constant meters): 100

 Do you have constants to apply a circular error correction? N

 You will be prompted with counter readings. You only need to type enough table

 entries to cover the range indicated by the printout from program FIELDG.

 The table is constructed in 2 ways. One way is to type milligal-equivalents

 (preferably to microgals) of enough counter readings to cover the range.

 The second way is to type the interval factors, starting at counter 000.

 Do you want the preferred first way? Y

 TYPE lowest counter reading needed (integer): 3300

 TYPE value in milligals after counter-reading prompt.

 Dial 3300 milligal-equivalent (zero to STOP): 3447.662

 Dial 3400 milligal-equivalent (zero to STOP): 3552.258

 Dial 3500 milligal-equivalent (zero to STOP): 3656.859

 Dial 3600 milligal-equivalent (zero to STOP): 3761.464

 Dial 3700 milligal-equivalent (zero to STOP): 3866.071

 Dial 3800 milligal-equivalent (zero to STOP): 0

 Proofread file: G8B.CAL before running MERGEG.

[PRINT OF OUTPUT FILE G8B.CAL]

G8B INTERVAL= 100 CIRC ERR= 0 FACTOR= 1.000612

 3300 3447.662 1.045960

 3400 3552.258 1.046010

 3500 3656.859 1.046050

 3600 3761.464 1.046070 {FORMAT REPRODUCES MANUFACTURER'S TABLES}

R.C. Jachens (written commun., 1981) developed a method to correct approximately harmonic

deviations within the period of one dial revolution for gravity meters. Initially, the circular correction, c0, for

a dial reading, r, is c0=r-R, where R is a reading datum. The corrected reading, ri is calculated recursively

by the sequence

ai = 2πci−1

Fi

, ci = ci−1 − Sisin ai − Ci cosai , and ri = ci + R

with F1=70.9412, F2=35.4706, F3=23.6471, and F4=17.7353 for the kind of gravity meters studied. Subscripts

refer to the number of the harmonic. R, the sine (S), and cosine (C) terms for the harmonics are determined

by comparing gravity meter readings at intervals of about 5 units (100 units is one dial revolution) with

readings of other gravity meters. In practice, as indicated in the following example, contributions to the

circular error from terms after the first two harmonic terms is negligible.
25

calibg [COMMAND TO EXECUTE PROGRAM]

 CALIBG, Plouff, 3-98. Create a gravity meter calibration file from

 manufacturer's table.

 TYPE a 1- to 4-digit name for the gravity meter (left-adjusted): 17G

 Your calibration file will be named: 17G.CAL

 IF TYPING ERRORS FOR NUMBERS OCCUR, QUESTIONS WILL BE REPEATED.

 An experimentally determined constant may be multiplied times the instrument

 calibration factors (in milligals per dial unit) provided by the manufacturer.

 TYPE constant correction factor (about 1.): 1.00252

 The counter reading interval is 100 dial units for LaCoste-Romberg G-meters

 and 10 or 50 for D-meters. An interval of zero signifies gravity meters

 where a single constant applies.

 TYPE counter reading interval (integer) (0 for single-constant meters): 100

 Do you have constants to apply a circular error correction? y

 TYPE zero reading datum: 3093.918

 Sine term of first harmonic: 0.01287

 Cosine term of first harmonic: 0.03091

 Sine term of second harmonic: -0.00298

 Cosine term of second harmonic: 0.01279

 You will be prompted with counter readings. You only need to type enough table

 entries to cover the range indicated by the printout from program FIELDG.

 The table is constructed in 2 ways. One way is to type milligal-equivalents

 (preferably to microgals) of enough counter readings to cover the range.

 The second way is to type the interval factors, starting at counter 000.

 Do you want the preferred first way? y

 TYPE lowest counter reading needed (integer): 2700

 TYPE value in milligals after counter-reading prompt.

 Dial 2700 milligal-equivalent (zero to STOP): 2819.379

 Dial 2800 milligal-equivalent (zero to STOP): 2923.887

 Dial 2900 milligal-equivalent (zero to STOP): 3028.405

 Dial 3000 milligal-equivalent (zero to STOP): 3132.935

 Dial 3100 milligal-equivalent (zero to STOP): 0

 Proofread file: 17G.CAL before running MERGEG.

[PRINT OF OUTPUT FILE 17G.CAL]

17G INTERVAL= 100 CIRC ERR= 1 FACTOR= 1.002520

R0= 3093.918 S1= 0.012870 C1= 0.030910 S2=-0.002980 C2= 0.012790

 2700 2819.379 1.045080

 2800 2923.887 1.045180

 2900 3028.405 1.045300

______________________________________________________________________________________________________________________________

APPENDIX 5

 CONVERTING READINGS TO MILLIGALS AND CORRECTING FOR TIDE EFFECT ______________________________________________________________________________________________________________________________

Program mergeg.f merges dates, times, readings, elevations, accuracy codes, and field terrain

corrections in the .FLD file with the geographic coordinates and the flag for multiple readings in the .FCT

file and then calculates tide-corrected milligal-equivalent readings by applying file .CAL to create a .TID

file (step 4 in fig. 1). Fictitious data sets, including constant readings to span midnight, December 31,

1999, first were created so that conversion of readings to milligals and calculation of tide corrections

printed in output print file MERGEG.PNT could be verified. Conversions of readings to milligals easily

were verified for fictitious data sets by manual interpolation of gravity meter calibration tables reproduced

by file .CAL and multiplication by a constant. 

Broucke and others (1972, tables 1 and 4) listed separate components for the vertical acceleration due to

the attractions of the moon and the sun. Two calculations were made for a time of 22h 25m 37.128s

Coordinated Universal Time on November 11, 1970 at 89° 49.693' lon W. One point was located at 60°

56.4' and another at 15° 56.4' lat N. After the effects of the moon and the sun were added and the effect of

a rigid Earth was multiplied by an assumed elastic constant of 1.16, Earth tide corrections obtained by

Broucke and others (1972) were -89.7 and -96.9 microgals, respectively, which was stated as an

improvement to equivalent values of -89.2 and -96.6 microgals obtained from Longman's (1959) formulas. 
26

Fictitious files first were created, in which gravity readings were constant for field data and the instrument

calibration was a constant 1.0 mGal/dial unit. Using a local time of 16h 25m 37.128s and a correction to

Greenwich Mean Time (Coordinated Universal Time) of 6 hours, values of -89 and -96, respectively, were

obtained by program mergeg.f, which uses an algorithm originally derived from Longman's (1959)

formulas (oral commun., David F. Barnes, 1998). 

Earth tide corrections computed by program mergeg.f also were compared with values interpolated

from tables and graphs published by the European Association of Exploration Geophysicists (EAEG)

(Goguel, 1954; Service Hydrographique de la Marine and Compagnie Générale de Géophysique 1969)

until 1972. Tide corrections (in units of 0.01 mGal) in the following table were compared between values

from program mergeg.f and interpolated values from EAEG tables (from Service Hydrographique de la

Marine and Compagnie Générale de Géophysique, 1969, p. 1, 2, and 50) for the location at 60° 56.4' in the

above example from Broucke and others (1972). 

24-hour time 0200 0600 1000 1400 1800 2000 2200

mergeg.f value +1.9 -6.9 +10.6 -1.2 -6.2 +6.7 +18.4

EAEG value +1 1/2 -8 +11 -1 1/2 -6 +7 +18

Tide corrections require the time interval between the date of a gravity observation and noon December

31, 1899 in order to utilize Longman's (1959) algorithm. The mergeg.f program takes into account

that—unlike the year 1900—the year 2000 is a leap year. No attempt has been made to determine if

computer algorithms are available, which may supersede Longman's (1959) algorithm. Carter W. Roberts

(written commun., 1998) suggested that tide corrections plotted at close time intervals on the date and at a

location of an eclipse in the year 2017, for example, should indicate an Earth tide maximum at a predicted

time and place along the path of the eclipse (fig. 3, bottom plot) and, therefore, tends to verify that tide

corrections after the year 2000 are correctly calculated by program mergeg.f. 

The following computer session for program mergeg.f reproduces step 4 (fig. 1) of the gravity

reduction process needed to prepare file blm938.TID as input to the following steps of determining

absolute values of observed gravity for gravity data collected in 1993 (Plouff, 1996). Typed responses are

in bold print. Additional comments not in the original session are bracketed in italics.

mergeg [COMMAND TO EXECUTE PROGRAM]

 MERGEG, Plouff, 3-98. Merges file of locations (.FCT) with file of gravity

 readings (.FLD) and converts tide-corrected readings to milligals.

 A print file MERGE.PNT already exists in your directory.

 Do you want to STOP to save (rename) it? n

 TYPE the root name associated with your files (8-digit max): blm938

 Gravity meter: G8B ; Year: 1993; Greenwich correction: 7 hours.

 Outer ring for field TC estimate is HAMMER C

 Readings for a total of 12 days. If no errors, next run SETBASEG program,

 which creates or needs file blm938.BAS .

 Print file MERGEG.PNT requires a landscape page orientation for 93-columns.

[MODIFIED PRINT OF BEGINNING OF FILE MERGEG.PNT]

 24-Mar-98 09:22:31 MERGEG, Plouff, 3-98. Merges locations in file blm938.FCT

 with readings in file blm938.FLD for output in file blm938.TID

 Gravity meter: G8B ; Year: 1993; Greenwich correction: 7 hours.

 Outer ring for field TC estimate is HAMMER C

 Calibration table (interval= 100 units) for gravity meter G8B

 and values using a multiplication factor of 1.000612

 Table of readings ranging from 3200 to 3600 units.

 MANUFACTURER'S REVISION

 COUNTER MGAL FACTOR MGAL FACTOR

 3200 3343.072 1.045901 3345.118 1.046541

----------------------------------------------

 3600 3761.464 0.000000 3763.766 0.000000

[PRINT OF END OF FILE MERGEG.PNT, WITH SHORTENED DATA LINES]

 338 314 WN258 4135.56 11823.89 4132.0 1638 3532.980 0.000 -0.006 3693.611 G64 0.01C

 339 315 OROV 4134.17 11747.06 4337.0 1851 3509.310 0.000 -0.083 3668.759 0.00

 340 316 WNBAS 4058.35 11744.02 4297.9 1942 3465.900 0.000 -0.091 3623.315 N22 0.02C

 Readings for a total of 12 days. If no errors, next run SETBASEG program,

 which creates or needs file blm938.BAS .
27

[BEGINNING OF OUTPUT FILE blm938.TID]

G8B 1993

WNBAS4058351174402 42979 62193 912 7 023623385N22 2C

WINNA4054231174822 42980 62193 943 7 023606421 0

WNBAS4058351174402 42979 621931012 7 023623517N22 2C

WINNA4054231174822 42980 621931038 7 023606437 0

[END OF OUTPUT FILE blm938.TID]

WN2564137331182281 43040 7 29316 4 7 013697899G64 1C

WN2574136331182444 41540 7 2931621 7 013697163G63 1C

WN2584135561182389 41320 7 2931638 7 013693611G64 1C

 OROV4134171174706 43370 7 2931851 7 023668759 0

WNBAS4058351174402 42979 7 2931942 7 023623316N22 2C

The first line of output file blm938.TID retains the gravity meter name and the year of data

observations for later processing steps. Lines in the body of the file include (without decimal points for

units): 5-digit station name; latitude in degrees and minutes to hundredths; longitude in degrees and

minutes to hundredths (without negative sign); elevation to 0.1 ft; date abbreviated as month-day-year; 24-

hour time in hours and minutes; Greenwich correction in hours; an existing value of observed gravity

expressed to 0.01 mGal without the leading 9—usually zero at this stage; code "2" for reoccupied locations

or code "1" for stations with only one observation; tide-corrected gravity reading, converted to units of

microgals; 3-digit accuracy code; position for previously established observed-gravity accuracy code; field

terrain correction to 0.01 mGal; and ring code for field terrain correction. Although the reduction process

yields observed gravity expressed to 0.01 mGal, readings are carried forward in units of microgals (0.001

mGal) to minimize rounding errors. 

______________________________________________________________________________________________________________________________

APPENDIX 6

CONVERTING TIDE-CORRECTED READINGS TO OBSERVED GRAVITY ______________________________________________________________________________________________________________________________

Converting tide-corrected, milligal-equivalent readings in .TID files obtained in step 4 (fig. 1) to

absolute values of observed gravity in .OBS files requires a repetition of inserting previously established

values of observed gravity in step 5, determining the best linear drift-corrected values of observed gravity in

step 6, and inspecting drift plots in step 7 (appendix 7). After the value of observed gravity was established

at the principal base station, the following program executions ("runs") were used to correct for instrument

drift and to establish values of observed gravity at the remaining stations collected in 1993 (Plouff, 1996).

Runs of programs setbaseg.f and driftg.f were applied in a sequence of first determining values for the

most occupied stations, stations that had the least standard deviations, or stations closest in time to

previously established stations. After final values were determined, program extraobs.f was used to

remove redundant observations (appendix 8). Typed responses in the following sessions are in bold print.

Additional comments not in the original session are bracketed in italics. To minimize redundancy, print

format is shortened and sessions are abbreviated after the first run. 

RUN 1

setbaseg [COMMAND TO EXECUTE PROGRAM]

 SETBASEG, D. Plouff, 3-98. Direct editing of .TID files is impractical.

 This program provides station option code and absolute OG's to DRIFTG program.

 TYPE the root name of your files (8-digit max): blm938

 Input file blm938.BAS that provides initial OG's is not in your directory.

 Do you want to create it interactively? y

 Writing interactively into file blm938.BAS

 First, base stations will be entered. Next, repeat stations. Last, tares

 (instantaneous reading shift at place with known OG). Some cases are not

 covered by this list (safer to insert by editing blm938.BAS ): station

 converted back to code 1 (need to assure in FIELDG printout that it is the

 only reading at that station); code 6 to create a fictitious base reading

 for file output in DRIFTG (need to pick a date and time); and code 7 to

 delete a reading.

 TYPING ERRORS WILL CAUSE A QUESTION TO BE REPEATED.
28

 Do you have a base station to be identified now? y

 TYPE the base station name (5 digits): WINNA

 Typing 1.0 next indicates that this station is to be upgraded from a repeat

 station, which already has an OG.

 TYPE a value for OG (mGal with decimal point, without leading 9): 79810.48

 Do you have more base stations to update? Y

 TYPE the base station name (5 digits): WNBAS

 Typing 1.0 next indicates that this station is to be upgraded from a repeat

 station, which already has an OG.

 TYPE a value for OG (mGal with decimal point, without leading 9): 79827.51

 Do you have more base stations to update? N

 Do you have a repeat station to be identified now? N

 Do you have a tare to identify now at a station with a known OG? N

 Do you want to print option code definitions in your print file? N

 Tide-corrected (but not drift-corrected) input file is blm938.TID

 Print file is: SETBASEG.PNT

 Temporary output file, blm938.GIV , is to be renamed blm938.TID

 2 OG's and option codes are in file blm938.BAS

 Gravity meter: G8B ; year of field data: 1993

 List of stations updated and number of occurrences. WINNA 3 WNBAS 29

 Next, inspect print file SETBASEG.PNT

 Rename blm938.GIV to blm938.TID and run DRIFTG program to find values of

 observed gravity.

[OUTPUT FILE blm938.BAS]

WINNA 3 79810.480 BASE

WNBAS 3 79827.510 BASE

mv blm938.TID blm938.TID1 [RENAME ORIGINAL FILE TO SAVE RECORD]

mv blm938.GIV blm938.TID [UPDATE .TID FILE FOR INPUT TO NEXT STEP]

driftg [COMMAND TO EXECUTE PROGRAM]

 DRIFTG (Plouff, 7-98)

 To convert tide-corrected gravity readings in microgals to absolute OG.

 TYPE the root name associated with your files (8-digit max): blm938

 Do you want to evaluate preliminary drift without station output? y

 DRIFTG (Plouff, 7-98). Your input file is blm938.TID

 Run program DRIFTPLT, to get a drift plot from file blm938.DFT

 Gravity meter: G8B ; year of field data: 1993

 Run program DRIFTPLT, to get a drift plot from file blm938.DFT

 If further refinement is needed, next, prepare a new file blm938.BAS

 and again run SETBASEG program.

 Print file (112 columns, landscape) is DRIFTG.PNT

[PRINT FILE DRIFTG.PNT]

 DRIFTG (Plouff, 7-98). Your input file is blm938.TID

[SELECTED LINES OF SUMMARY AT END OF FILE DRIFTG.PNT]

 List of base and repeat stations for 317 readings.

 PAST SEQ NAME LAST TIMES GRAVITY STD

 NAMES LINE READ (MGAL) DEV

 1 WNBAS 317 29 979827.51 base

 2 WINNA 8 3 979810.48 base

 3 WN 2 62 3 979888.25 0.03 [ACCEPT THIS VALUE]

 4 WN 6 294 8 979898.85 0.02 [ACCEPT THIS VALUE]

 5 WN 8 126 2 979898.88 0.03

 9 WN 32 248 3 979850.34 0.01

 12 AUX53 146 2 979803.09 0.01

 16 WN 60 291 2 979886.23 0.02

 19 WN 85 176 3 979847.09 0.02 [ACCEPT THIS VALUE]

 23 WN131 169 2 979872.12 0.01

 25 P163 232 3 979848.49 0.00 [AVERAGE WITH PREVIOUS .51]

Before continuing, a drift plot was inspected to check for possible errors (fig. 3, top plot; appendix 7). 
29

RUN 2

setbaseg [COMMAND TO EXECUTE PROGRAM]

 TYPE the root name of your files (8-digit max): blm938

 Input file blm938.BAS that provides initial OG's is in your directory.

 Do you want to append new options interactively? y

 2 lines in file to be appended blm938.BAS

 Writing interactively into file blm38.BAS

[SUBSTITUTED THE FOLLOWING VALUES FOUND IN RUN 1 FOR REPEAT STATIONS]:

[ P163, 79848.50; WN 2, 79888.25; WN 6, 79898.85; WN 85, 79847.09]

 Do you have more repeat stations to update? N

 Do you have a tare to identify now at a station with a known OG? N

 Do you want to print option code definitions in your print file? N

 Temporary output file, blm938.GIV , is to be renamed blm938.TID

 6 OG's and option codes are in file blm938.BAS

 List of stations updated and number of occurrences.

 WINNA 3 WNBAS 29 P163 3 WN 2 3 WN 6 8 WN 85 3

 Rename blm938.GIV to blm938.TID and run DRIFTG program.

[FILE blm938.BAS, WITH LINES PAIRED]

WINNA 3 79810.480 BASE WNBAS 3 79827.510 BASE

 P163 2 79848.500 REPEAT WN 2 2 79888.250

[LAST LINES OF FILE DRIFTG.PNT]

 PAST SEQ NAME LAST TIMES GRAVITY STD

 NAMES LINE READ (MGAL) DEV

 5 WN 8 126 2 979898.90 0.02 [ACCEPT THIS VALUE]

 7 WN 13 29 2 979883.28 0.01

 9 WN 32 248 3 979850.35 0.01 [ACCEPT THIS VALUE]

 11 WN 42 67 2 979883.68 0.02

 15 WN 57 279 2 979880.12 0.03

 16 WN 60 291 2 979886.23 0.02 [ACCEPT THIS VALUE]

 18 WN 71 98 2 979849.18 0.02

 23 WN131 169 2 979872.13 0.01 [ACCEPT THIS VALUE]

RUN 3

setbaseg [COMMAND TO EXECUTE PROGRAM]

[SUBSTITUTED THE FOLLOWING VALUES FOUND IN RUN 2 FOR REPEAT STATIONS]:

[ WN 8, 79898.90; WN 32, 79850.35; WN 60, 79886.23; WN131, 79872.13]

 List of stations updated and number of occurrences.

 WINNA 3 WNBAS 29 P163 3 WN 2 3 WN 6 8 WN 85 3 WN 8 2

 WN 32 3 WN 60 2 WN131 2

mv blm938.TID blm938.TID3 [RENAME FILE TO SAVE RECORD

mv blm938.GIV blm938.TID [UPDATE .TIDFILE FOR INPUT TO NEXT STEP]

mv SETBASEG.PNT SETBASEG.PNT3 [RENAME FILE TO SAVE RECORD]

driftg

[EXCERPTS FROM FILE DRIFTG.PNT]

 PAST SEQ NAME LAST TIMES GRAVITY STD

 NAMES LINE READ (MGAL) DEV

 6 WN 11 24 1 979892.26 [AVERAGE WITH WNA11]

 7 WN 13 29 2 979883.28 0.01 [ACCEPT THIS VALUE]

 11 WN 42 67 2 979883.68 0.01 [ACCEPT THIS VALUE]

 12 AUX53 146 2 979803.10 0.02 [ACCEPT THIS VALUE]

 15 WN 57 279 2 979880.12 0.01 [.13 FITS DETAIL BETTER]

 18 WN 71 98 2 979849.18 0.02 [ACCEPT THIS VALUE]

 21 WNA11 127 1 979892.28 [APPROXIMATE REOCCUPATION WN 11]

 30 OROV 316 1 979872.96 [WAS 873.05 FROM ANOTHER BASE]

mv blm938.DFT blm938.DFT3 [RENAME PLOT INPUT TO SAVE RECORD]

mv DRIFTG.PNT DRIFTG.PNT3 [RENAME FILE TO SAVE RECORD]
30

RUN 4—FINAL VALUES WITH OUTPUT OF VALUES FOR OBSERVED GRAVITY

setbaseg [COMMAND TO EXECUTE PROGRAM]

[SUBSTITUTED THE FOLLOWING REPEAT STATIONS AND VALUES: WN 11, 79892.27;

 WN11A, 79892.27; WN 13, 79883.28; WN 42, 79883.68; WN 57, 79880.13;

 WN 71, 79849.18; WN249, 79913.13; AUX53, 79803.10; OROV, 79872.98]

[FILE blm938.BAS FINAL VALUES, PRINTED HERE AS PAIRED LINES]

WINNA 3 79810.480 BASE WNBAS 3 79827.510 BASE

 P163 2 79848.500 REPEAT WN 2 2 79888.250 REPEAT

WN 6 2 79898.850 REPEAT WN 85 2 79847.090 REPEAT

WN 8 2 79898.900 REPEAT WN 32 2 79850.350 REPEAT

WN 60 2 79886.230 REPEAT WN131 2 79872.130 REPEAT

WN 11 2 79892.270 REPEAT WN11A 2 79892.270 REPEAT

WN 13 2 79883.280 REPEAT WN 42 2 79883.680 REPEAT

WN 57 2 79880.130 REPEAT WN 71 2 79849.180 REPEAT

WN249 2 79913.130 REPEAT AUX53 2 79803.100 REPEAT

 OROV 2 79872.980 REPEAT

mv blm938.TID blm938.TID4 [RENAME ORIGINAL FILE TO SAVE RECORD]

mv blm938.GIV blm938.TID [UPDATE .TID FILE FOR NEXT STEP]

driftg [COMMAND TO EXECUTE PROGRAM]

 TYPE the root name associated with your files (8-digit max): blm938

 Do you want to evaluate preliminary drift without station output? n

 Your station output file is blm938.OBS

[EXCERPTS FROM DRIFTG.PNT PRINT FILE; DATE LATER IS "1998:0722"]

 DRIFTG (Plouff, 7-98). Your input file is blm938.TID

 Run program DRIFTPLT, to get a drift plot from file blm938.DFT

 Your station output file is blm938.OBS

 Asterisk on the following list refers to an assumed gravity value.

 Base stations are indicated by blanks in the second NAME-column.

 Trailing letter P indicates that .OBS file has that output data point.

 22-Jul-98 14:53:25 Gravity meter: G8B ; year of field data: 1993.

 LINE NAME LATITUDE LONGITUDE ELEVATION TIME DATE READING BASE ALL GRAVITY LAST TC

 DEG MIN DEG MIN (FEET) H-M GMT (MGAL) STA DRIFT (MGAL) LINE

 8 WINNA 40 54.23 117 48.22 4298.0 1227 7 6-21-93 3606.416 WINNA 0.061 979810.48* 5

 9 WNBAS 40 58.35 117 44.02 4297.9 1244 7 6-21-93 3623.433 WNBAS 0.048 979827.51* 7

 10 WN 1 41 23.88 117 49.45 4219.2 1342 7 6-21-93 3673.637 WNBAS 0.052 979877.71 P

 11 WN 2 41 24.45 117 57.81 4122.7 14 2 7 6-21-93 3684.178 WNBAS 0.053 979888.25* P 11 C

[OMITTED THE COLUMN WITH REPEATED STATION NAME]

 314 WN257 41 36.33 118 24.44 4154.0 1621 7 7- 2-93 3697.163 WNBAS -0.105 979901.39 P

 315 WN258 41 35.56 118 23.89 4132.0 1638 7 7- 2-93 3693.611 WNBAS -0.104 979897.84 P

 316 OROV 41 34.17 117 47.06 4337.0 1851 7 7- 2-93 3668.759 WNBAS -0.096 979872.98* P 316

 317 WNBAS 40 58.35 117 44.02 4297.9 1942 7 7- 2-93 3623.316 WNBAS -0.069 979827.51* 293

 22-Jul-98 14:53:25 List of base and repeat stations for 317 readings.

[BEGINNING AND END OF OUTPUT FILE blm938.OBS]

[ NAME LAT-D/M LON-D/M EL-FT OG-MGAL ACC DATE HTC R BASE]

 WN 1 412388 1174945 421927987771N22 61993 1 C WNBAS

 WN 2 412445 1175781 412277988825N22 61993 1 C WNBAS

 WN 3 412375 1180150 418047989360N22 61993 1 C WNBAS

--------------------------------------------------------------------------------

 WN256 413733 1182281 430407990213G64 71993 1 C WNBAS

 WN257 413633 1182444 415407990139G63 71993 1 C WNBAS

 WN258 413556 1182389 413207989784G64 71993 1 C WNBAS

 OROV 413417 1174706 433707987298 71993 0 WNBAS

Lines in file blm938.OBS include (without decimal points for units): 3 blank digits as prefixes to the

following 5-digit station name; latitude in degrees and minutes to hundredths; longitude in degrees and

minutes to hundredths (without negative sign); elevation to 0.1 ft; observed gravity to 0.01 mGal without

the leading 9; a 3-digit accuracy code; date expressed as month and 4-digit year; field terrain correction to

0.01 mGal; outer ring code for field terrain correction; and base station to which station is referenced. 
31

______________________________________________________________________________________________________________________________

APPENDIX 7

PLOTTING INSTRUMENT DRIFT ______________________________________________________________________________________________________________________________

Program driftg.f outputs a file to be used for creating an instrument drift plot. For example, the file

used to create the final drift plot (fig. 3, middle plot) for data collected in 1993 (Plouff, 1996) follows. If

datetime.f is appended to program driftg.f, the compilation date "1-Apr-98" is replaced by "19980401" in

figure 3. The following explanatory comments not in the original file are bracketed in italics. 

[BEGINNING AND END OF driftg.f OUTPUT FILE blm938.DFT]

 2-Apr-98 08:14:39 [DATE LATER IS "1998:0402"]

blm938 [ROOT NAME OF DATA SET] ]

G8B 1993 [GRAVITY METER; YEAR OF DATA COLLECTION]

 7 6 21 93 912 [BEGINNING FLAG A REPEAT OF FIRST TIME]

 7 6 21 912WNBASWNBAS+ 0.000 [GMT CORRECTION; MONTH-DAY; STATION; BASE;

 7 6 21 943WINNAWINNA+ 0.066 PLOT SYMBOL; DRIFT IN MILLIGALS FROM FIRST]

 7 6 21 1012WNBASWNBAS+ 0.132

 7 6 21 1038WINNAWINNA+ 0.082

--------------------------------------------------------------------------------

 7 7 2 1320WN249WNBAS. -0.116

 7 7 2 1851 OROVWNBAS. -0.096

 7 7 2 1942WNBASWNBAS+ -0.069

 7 7 2 1942 [ENDING FLAG A REPEAT OF FINAL TIME]

The following computer session for program driftplt.f created a plot file blm938.PST in the

PostScript language (Adobe Systems Incorporated, 1990) from input file blm938.DFT. Typed responses

are in bold print. Explanatory comments not in the original session are bracketed in italics. A carriage

return (ENTER key) signifies a response of “yes” (or “y” or “Y”). 

driftplt [COMMAND TO EXECUTE PROGRAM]

 Plouff, 3-98.

 Program to plot drift of observed gravity output from program DRIFTG.

 TYPE the root name associated with your files (8-digit max): blm938

 72 drift values range from -0.134 to 0.132 mGal.

 Your data extend for 12 days between JUNE 21 and JULY 2

 The largest inter-day drift change is 0.05 between JUNE 22 and JUNE 23 (gap 23)

 You may delete either the largest time or drift gap on your plot.

 TYPE the gap number (integer) with zero for none: 0

 TYPE the time scale in days/inch: 2

 Plot width is 7.0 inches. OKAY?

 You may select drift plot boundaries to exclude extreme values.

 Your all inclusive boundaries are -0.2 to 0.2 mGal. Acceptable?

 TYPE the drift scale in mGal/inch: 0.35

 Plot height is 2.1 inches. OKAY?

 A total of 72 points were plotted.

 Rename the plot file if it is to be saved or other specs are wanted.

 Your print file is DRIFTPLT.PNT

 Your plot file is blm938.PST . Edit BoundingBox?

[BEGINNING OF driftplt.f OUTPUT FILE blm938.PST]

%!PS-Adobe [COMMAND LINE REQUIRED FOR POSTSCRIPT PLOTS]

%%BoundingBox 30 30 800 800 [FOR 8.5 BY 11 INCH PAPER; OPTIONAL?]

/M {moveto} def /S {show} def

/LCSM {lineto currentpoint stroke M} def

72 72 scale 90 rotate 0.25 -8.2 translate [CHANGE TO ROTATE AND TRANSLATE PLOT]

0.004 setlinewidth

 0.8500 0.1600 M

 0.8900 0.2000 LCSM

 0.8500 0.2000 M

Program driftplt.f splits a plot into two adjacent parts with a small gap in response to prompting that

prints gap numbers for greatest number of missing days or the largest inter-day drift value change. 
32

______________________________________________________________________________________________________________________________

APPENDIX 8

DELETING REDUNDANT DATA ______________________________________________________________________________________________________________________________

Unless flags are set or the data points are defined as base stations, program driftg.f outputs one line in

the .OBS file for every gravity reading (See example at end of appendix 6). Inasmuch as the value of

observed gravity at this stage is the same for all readings at a given station, output lines that correspond to

readings before the last station occupation, for example, are superfluous. File .OBS first needs to be

rearranged by sorting in the sequence of station name, year, and month, respectively, so that program

extraobs.f can delete superfluous lines. The following command in the VAX/VMS operating system

(Digital Equipment Corporation), for example, can rearrange lines in the .OBS file in this way.

SORT OBS:==SORT/KEY=(POS:1,SIZ:8)/KEY=(POS:49,SIZ:4)/KEY=(POS:47,SIZ:2)

In the UNIX (trade name of American Telephone and Telegraph Co.) system, the following command

converted file blm938.OBS for data collected in 1993 (Plouff, 1996) into rearranged file blm938.SRT. 

sort +0.0 -0.8 +0.48 -0.52 +0.46 -0.48 blm938.OBS -o blm938.SRT

(BEGINNING AND END OF FILE blm938.SRT]

 P163 405500 1182281 416307984850G53 61993 1 C WNBAS

 P163 405500 1182281 416307984850G53 61993 1 C WNBAS

 P163 405500 1182281 416307984850G53 61993 1 C WNBAS

 AUX53 411469 1175988 550527980310G75 61993 2 C WNBAS

 AUX53 411469 1175988 550527980310G75 61993 2 C WNBAS

 WN 1 412388 1174945 421927987771N22 61993 1 C WNBAS

 WN 2 412445 1175781 412277988825N22 61993 1 C WNBAS

-----------------------------------------------------------------------------

 WN257 413633 1182444 415407990139G63 71993 1 C WNBAS

 WN258 413556 1182389 413207989784G64 71993 1 C WNBAS

 WNA11 412468 1180426 412407989228G55 61993 1 C WNBAS

The following computer session for program extraobs.f converted the sorted output file blm938.OBS

from program driftg.f to an arbitrarily named output file, blm938.plf, without superfluous stations. Typed

responses are in bold print. Explanatory comments not in the original session are bracketed in italics. The

date later is expressed as "1998:0510." 

extraobs [COMMAND TO EXECUTE PROGRAM]

 EXTRAOBS, Plouff, 5-98. 10-May-98 15:28:26 Interactive program to delete

 repeated station names in an extended plouff file. The file first must be

 sorted by station name, year, and month, respectively:

 sort +0.0 -0.8 +0.48 -0.52 +0.46 -0.48 file.OBS -o file.srt

 Dates (col. 47-52) and base names (76-80) are included. Therefore, the

 station with the latest date will be selected. The name of the most

 commonly used base is user selected. Other variables must be identical.

 TYPE the name of your plouff file (<66 cols): blm938.SRT

 TYPE the name of a plouff output file (<66 cols): blm938.plf

 A total of 284 input lines were read and 261 output lines were written.

[FILE EXTRAOBS.PNT]

 EXTRAOBS, Plouff, 5-98. 10-May-98 15:28:26 Interactive program to delete

 repeated station names in an extended plouff file. The file first must be

 sorted by station name, year, and month, respectively:

 sort +0.0 -0.8 +0.48 -0.52 +0.46 -0.48 file.OBS -o file.srt

 Dates (col. 47-52) and base names (76-80) are included. Therefore, the

 station with the latest date will be selected. The name of the most

 commonly used base is user selected. Other variables must be identical.
33

[FILE EXTRAOBS.PNT CONTINUED, WITH ABBREVIATED LINES]

 Input file: blm938.SRT Output file: blm938.plf

 3 stations named P163 Dates range from 61993 to 61993. Bases are: WNBAS ( 3)

 2 stations named AUX53 Dates range from 61993 to 61993. Bases are: WNBAS ( 2)

 3 stations named WN 2 Dates range from 61993 to 61993. Bases are: WNBAS ( 3)

 8 stations named WN 6 Dates range from 61993 to 71993. Bases are: WNBAS ( 8)

[NOT PRINTED HERE: WN 8; WN 13; WN 32; WN 42; WN 57; WN 60; WN 71]

 3 stations named WN 85 Dates range from 62493 to 62693. Bases are WNBAS ( 3)

 2 stations named WN131 Dates range from 62693 to 62693. Bases are: WNBAS ( 2)

 A total of 284 input lines were read and 261 output lines were written.

Recalling that base station names are automatically assigned by program driftg.f based on closeness in

time between readings of the station and the base station, the user would have been offered a choice of lines

to delete if a station was tied to more than one base station. 

[BEGINNING AND END OF FILE blm938.plf]

 OROV 413417 1174706 433707987298 0 71993 0 0 0 0WNBAS

 P163 405500 1182281 416307984850G53 0 61993 1 0C 0 0WNBAS

 AUX53 411469 1175988 550527980310G75 0 61993 2 0C 0 0WNBAS

 WN 1 412388 1174945 421927987771N22 0 61993 1 0C 0 0WNBAS

 WN 2 412445 1175781 412277988825N22 0 61993 1 0C 0 0WNBAS

--------------------------------------------------------------------------------

 WN257 413633 1182444 415407990139G63 0 71993 1 0C 0 0WNBAS

 WN258 413556 1182389 413207989784G64 0 71993 1 0C 0 0WNBAS

 WNA11 412468 1180426 412407989228G55 0 61993 1 0C 0 0WNBAS

Program extraobs.f also can be used to delete superfluous data if data from different years are merged

into one file. The step of sorting, however, will rearrange the combined file into a new sequence in

descending order of station names. If a new sequence is not wanted, program substute.f (Plouff, 1998,

p. 43) can extract stations with the same 8-digit name from data files for further examination before the

files are merged, for example, to compare values of observed gravity,. 

______________________________________________________________________________________________________________________________

APPENDIX 9

ESTIMATING GRAVITY TERRAIN CORRECTIONS ______________________________________________________________________________________________________________________________

Terrain correction for a post

The closest part of the TC is the effect of the post, upon which a base plate and a gravity meter sits.

Program tcpost.f creates a PostScript plot (Adobe Systems Incorporated, 1990) of the TC of a concrete

post (fig. 6). The following computer session created a plot for the TC of a post if the center of mass of

the gravity meter is 0.4 ft (0.12 m) above the base of the gravity meter (fig. 6). Successive responses of n

and 0.4 rather than y or a carriage return was to test input of a variable height. 

tcpost [COMMAND TO EXECUTE PROGRAM]

 TCPOST, Plouff. 5-98. Terrain corrections for gravity readings on posts with

 square and cylindrical cross sections. Assume standard rock density is

 2.67 g/cm^3, and density of the post is 2.0 g/cm^3

 Is the height of gravity meter mass above bottom of baseplate 0.4 ft? n

 TYPE the height of mass in feet (less than 10): 0.4

 Print PostScript file TCPOST.PST for plots of terrain corrections.
Figure 6.--Gravity terrain correction for reading on a concrete post. Numbers inside drawings are the

height of the post in feet; interval, 0.2 ft. Terrain corrections are in addition to the constant effect 

between the gravity meter mass and the bottom of the base plate. Right edge of drawings indicates

terrain correction for an infinitely extended concrete slab with a density of 2.0 g/cm^3.

34

WIDTH OF SQUARE POST, IN FEET

HEIGHT OF GRAVITY METER MASS ABOVE BOTTOM OF BASE PLATE=0.40 FOOT

AIR CONCRETE

0.0 0.5 1.0 1.5 2.0 2.5

0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.10

0.11

TERRAIN CORRECTION, IN MILLIGALS

AIR CONCRETE

0.0 0.5 1.0 1.5 2.0 2.5

0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.10

0.11

TERRAIN CORRECTION, IN MILLIGALS

0.4

0.8

1.2

1.6

2.0

2.4

2.8

DIAMETER OF CYLINDRICAL POST, IN FEET

0.4

0.8

1.2

1.6

2.0

2.4

2.8
35

Terrain corrections for two-dimensional profiles

Most terrain corrections (TCs) can be estimated in the field by utilizing sets of curves based on a model

of a displaced slope (figs. 4 and 5). Program fieldtc.f computes a TC for two-dimensional topography

more complicated than a displaced slope. The following computer session computes a two-dimensional TC

for field notes: "flat for 10 ft NW; vertically downward (into canal, for example) for 6 ft; flat for 15 ft;

vertically upward for 6 ft; flat for 8 ft; and 15 degrees upward (beyond ring limit)." Negative signs only are

entered in cases such as this, in which the direction of slope reverses along the profile. 

fieldtc [COMMAND TO EXECUTE PROGRAM]

 FIELDTC, 9-99 Plouff, estimate inner terrain corrections from field notes,

 assuming two-dimensional topography. Notes usually are carried to distances

 of 175 ft (Hammer C) or 223 ft (Hayford B) from the gravity station.

 TYPE the radial distance (not greater than 223) in feet: 223

 A record of your results will be stored in a file named FIELDTC.PNT.

 You will be asked to type a station name and then triplets of values (A, DS, and

 DH). At least one of the three is zero, (redundant or unmeasured).

 A is the angle of slope in degrees. DS is the slope distance in feet. DH is

 the height difference in feet. Positive values are upward, as viewed from the

 station, and negative downward. DH=DS=0 for the last triplet in a direction.

 TYPE station name/direction (carriage return to stop): test NW

 TYPE the angle of slope (degrees) A: 0 [FLAT]

 TYPE the slope distance (feet) DS: 10 [PACE 10 FEET]

 TYPE the height difference (feet) DH: 0 [REDUNDANT]

 TYPE the angle of slope (degrees) A: -90 [VERTICALLY DOWN]

 TYPE the slope distance (feet) DS: 0 [REDUNDANT]

 TYPE the height difference (feet) DH: -6 [TAPE 6 FEET DOWN]

 TYPE the angle of slope (degrees) A: 0 [FLAT]

 TYPE the slope distance (feet) DS: 15 [PACE 15 FEET]

 TYPE the height difference (feet) DH: 0 [REDUNDANT]

 TYPE the angle of slope (degrees) A: 90 [VERTICALLY UP]

 TYPE the slope distance (feet) DS: 0 [REDUNDANT]

 TYPE the height difference (feet) DH: 6 [ALREADY TAPED]

 TYPE the angle of slope (degrees) A: 0 [FLAT]

 TYPE the slope distance (feet) DS: 8 [PACE 8 FEET]

 TYPE the height difference (feet) DH: 0 [REDUNDANT]

 TYPE the angle of slope (degrees) A: 15 [MEASURED ANGLE UPWARD]

 TYPE the slope distance (feet) DS: 0 [CLOSURE PAIR]

 TYPE the height difference (feet) DH: 0 [CLOSURE PAIR]

 Total 2D TC= 0.034 mGal

 TYPE station name/direction (carriage return to stop): [NO MORE MEASUREMENTS]

[File FIELDTC.PNT, abbreviated]

 Terrain corrections from fieldnotes: angles; distances; height differences.

 Maximum distance=223.0 feet

 Station/direction: test NW

 Angle= 0.0; SLOPE DIST= 10.0; HT DIFF= 0.0

 TC= 0.000 mGal; x= 0.0 to 10.0; ht= 0.0 to 0.0 ft; indices: 1 1

 Angle=-90.0; SLOPE DIST= 0.0; HT DIFF= -6.0

 TC= 0.000 mGal; x= 10.0 to 10.0; ht= 0.0 to -6.0 ft; indices: 1 1

 Angle= 0.0; SLOPE DIST= 15.0; HT DIFF= 0.0

 TC= 0.011 mGal; x= 10.0 to 25.0; ht= -6.0 to -6.0 ft; indices: 1 25

 Angle= 90.0; SLOPE DIST= 0.0; HT DIFF= 6.0

 TC= 0.000 mGal; x= 25.0 to 25.0; ht= -6.0 to 0.0 ft; indices: 25 25

 Angle= 0.0; SLOPE DIST= 8.0; HT DIFF= 0.0

 TC= 0.000 mGal; x= 25.0 to 33.0; ht= 0.0 to 0.0 ft; indices: 25 25

 Angle= 15.0; SLOPE DIST= 0.0; HT DIFF= 0.0

 TC= 0.024 mGal; x= 33.0 to 223.0; ht= 0.0 to 50.1 ft; indices: 25 87

 Total 2D TC= 0.034 mGal

File FIELDTC.PNT lists the input data, resulting TC components, indices for program de-bugging, and

total TC for each set of measurements.