PROGRAM GEOID *********************************************************************** * * * PROGRAM : GEOID * * * * PURPOSE: COMPUTATION PROGRAM TO INTERPOLOATE GEOID HEIGHT * * VALUES FROM A GRID OF REGULARLY-SPACED ESTIMATES. * * THE INPUT REQUIRES LATITUDE AND LONGITUDE VALUES * * REFERRENCED TO THE NORTH AMERICAN DATUM OF 1983 * * (NAD 83) AND THE WORLD GEODETIC SYSTEM OF 1980 (WGS 80).* * INPUT DATA REFERRENCED TO THE NORTH AMERICAN DATUM OF * * 1927 (NAD27) SHOULD BE CONVERTED TO NAD 83 VIA NADCON * * PRIOR TO USING THIS PROGRAM. * * * * THE ACTUAL COMPUTATION IS PERFORMED AS AN INTERPOLATION * * FROM A REGULARLY-SPACED GRID OF POINTS OBTAINED FROM * * TWO DIFFERENT SOURCES. ONE SOURCE IS THE GEOID93 MODEL * * COMPUTED BY DR. DENNIS G. MILBERT, NGS, NOAA. THE OTHER* * IS AN EVALUATION OF A GEOPOTENTIAL MODEL BASED ON A * * COMPLETE EXPANSION OF SPHERICAL HARMONIC COEFFICIENTS * * OBTAINED FROM OHIO STATE UNIVERSITY AND KNOWN AS OSU89B.* * THESE COEFFICIENTS WERE OBTAINED FROM DR. RICHARD RAPP * * OF THE GEODETIC SCIENCE AND SURVEYING DEPARTMENT. * * * * THE INTERPOLATION IS ACCOMPLISHED BY A LOCALLY FIT * * BIQUADRATIC FUNCTION. * * * * THE POLYNOMIAL SURFACE IS FIT TO NINE DATA POINTS * * DEFINING THE AREA SURROUNDING THE (X,Y) PAIR * * WHERE THE INTERPOLATION IS TO TAKE PLACE. * * * * THE PROGRAM REQUIRES THAT THE USER SPECIFY: * * * * 1) THE NAME OF AN OUTPUT FILE * * * * 2) THE NAME OF AN INPUT FILE (IF AVAILABLE). * * * * * * THIS PROGRAM ALLOWS FOR EITHER NGS STANDARD HORIZONTAL * * DATA FORMATS AS SPECIFIED IN THE FGCC PUBLICATION, * * COMMONLY KNOWN AS THE 'HORIZONTAL BLUE BOOK' (SEE * * SUBROUTINES TYPE3, TYPE4, TYPE5), OR IN A GENERIC FILE * * FORMAT (SEE SUBROUTINE TYPE1 OR SUBROUTINE TYPE2). * * * * THE CODE CAN BE EASILY MODIFIED TO ACCOMMODATE CUSTOM * * FILE SPECIFICATIONS BY MODIFYING SUBROUTINES: ENDREP, * * GETPT, IPARMS, WRTPT, AND (OPTIONALLY) FHELP. * * * * * * VERSION CODE: 1.00 (20-dec-90) * * VERSION CODE: 1.01 (27-feb-91) * * increaded MXAREA to 10 * * VERSION CODE: 2.00 (10-mar-93) * * * * VERSION DATE: march 10, 1993 * * * * DENNIS MILBERT, Ph.D. (GEOID93 model) * * DONALD SCHULTZ (GEOID93 model) * * ADVANCED GEODETIC SCIENCE BRANCH * * NATIONAL GEODETIC SURVEY, NOS, NOAA * * ROCKVILLE, MD 20852 * *********************************************************************** *********************************************************************** * * * DISCLAIMER * * * * THIS PROGRAM AND SUPPORTING INFORMATION IS FURNISHED BY THE * * GOVERNMENT OF THE UNITED STATES OF AMERICA, AND IS ACCEPTED AND * * USED BY THE RECIPIENT WITH THE UNDERSTANDING THAT THE UNITED STATES * * GOVERNMENT MAKES NO WARRANTIES, EXPRESS OR IMPLIED, CONCERNING THE * * ACCURACY, COMPLETENESS, RELIABILITY, OR SUITABILITY OF THIS * * PROGRAM, OF ITS CONSTITUENT PARTS, OR OF ANY SUPPORTING DATA. * * * * THE GOVERNMENT OF THE UNITED STATES OF AMERICA SHALL BE UNDER NO * * LIABILITY WHATSOEVER RESULTING FROM ANY USE OF THIS PROGRAM. THIS * * PROGRAM SHOULD NOT BE RELIED UPON AS THE SOLE BASIS FOR SOLVING A * * PROBLEM WHOSE INCORRECT SOLUTION COULD RESULT IN INJURY TO PERSON * * OR PROPERTY. * * * * THIS PROGRAM IS PROPERTY OF THE GOVERNMENT OF THE UNITED STATES * * OF AMERICA. THEREFORE, THE RECIPIENT FURTHER AGREES NOT TO ASSERT * * PROPRIETARY RIGHTS THEREIN AND NOT TO REPRESENT THIS PROGRAM TO * * ANYONE AS BEING OTHER THAN A GOVERNMENT PROGRAM. * * * *********************************************************************** * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) REAL VRSION INTEGER MXAREA PARAMETER (VRSION = 2.00E0, MXAREA =10) DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG DOUBLE PRECISION AGHT, VGHT DOUBLE PRECISION SDGHT, SDGHT2 REAL GHT, SMGHT, BGGHT INTEGER IPAGE, ITYPE, NCONV, IFILE CHARACTER*80 NAME LOGICAL PAGE, NODATA, SCREEN REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) ********************** * INITIALIZE VARIABLES ********************** CALL INITL (SCREEN, PAGE, NAME, IPAGE, ITYPE, + SMGHT, BGGHT, AGHT, VGHT, SDGHT, SDGHT2, + XSMALL, XBIG, YSMALL, YBIG) ************************** * PRINT HEADER INFORMATION ************************** CALL HEADR (VRSION) *********************** * OPEN GEOID DATA FILES *********************** CALL NGRIDS (NODATA) IF (NODATA) GOTO 9999 ****************************************************** * REQUEST FOR THE NEEDED VARIABLE VALUES FROM THE USER ****************************************************** CALL IPARMS (ITYPE, SCREEN) ******************************** * LOOP ONCE FOR EACH COMPUTATION ******************************** CALL MLOOP (NCONV, IPAGE, ITYPE, VRSION, GHT, + SDGHT, SDGHT2, SMGHT, BGGHT, + XSMALL, XBIG, YSMALL, YBIG, + PAGE, SCREEN, NAME) *********************************************** * FINISHED WITH ALL COMPUTATIONS - WRITE REPORT *********************************************** CALL ENDREP (IPAGE, NCONV, ITYPE, GHT, AGHT, VGHT, + SDGHT, SDGHT2, SMGHT, BGGHT, + XSMALL, XBIG, YSMALL, YBIG) ***************** * CLOSE ALL FILES ***************** DO 1010 IFILE = 1, NAREA CLOSE ( LUAREA(IFILE) ) 1010 CONTINUE CLOSE (NIN, STATUS='KEEP') CLOSE (NOUT, STATUS='KEEP') CLOSE (NAPAR, STATUS='KEEP') 9999 STOP END SUBROUTINE ASKPT (NCONV, NAME, IDLA, IMLA, SLA, + IDLO, IMLO, SLO, XPT, YPT, EOF, NOPT) * Interactively ask for the name and location of a point * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) CHARACTER*40 B40 PARAMETER (B40 = ' ') DOUBLE PRECISION XPT, YPT, RDLA, RDLO, DCARD REAL SLA, SLO, RMLA, RMLO, RCARD INTEGER NCONV INTEGER IDLA, IMLA, IDLO, IMLO INTEGER IFLAG1, IFLAG2, N1, N2 INTEGER LENG, IERR, IOS CHARACTER*80 NAME CHARACTER*40 ANS, DUM LOGICAL EOF, NOPT INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) DATA IFLAG1 /1/, IFLAG2 /2/ NAME = ' ' WRITE (LUOUT,*) ' What is the NAME for this station or', + ' point?' READ (LUIN,'(A80)') NAME ********** * LATITUDE ********** IF (NCONV .EQ. 0) THEN WRITE (LUOUT,110) 110 FORMAT (/, ' Latitudes and Longitudes may be entered', + ' in three formats:', /, + ' (1) degrees, minutes, and decimal seconds, OR', /, + ' (2) degrees, and decimal minutes OR', /, + ' (3) decimal degrees.', /, + ' Degrees, minutes and seconds may be separated', + ' by blanks or commas.', /, + ' A latitude or longitude of 0 will end data', + ' entry.', /) ENDIF WRITE (LUOUT,*) ' What is its NAD 83 latitude?' WRITE (LUOUT,*) ' ' READ (LUIN,170,ERR=9930,IOSTAT=IOS) ANS 170 FORMAT (A40) IF (ANS .EQ. B40) GOTO 9999 DUM = ANS CALL NBLANK (DUM, IFLAG2, N2) LENG = N2 RDLA = DCARD( DUM(1:N2), LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 IF (LENG .GT. 0) THEN RMLA = RCARD( DUM, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 IF (LENG .GT. 0) THEN SLA = RCARD( DUM, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 ELSE SLA = 0.E0 ENDIF ELSE RMLA = 0.E0 SLA = 0.E0 ENDIF IF ( (RDLA .EQ. 0.D0) .AND. (RMLA .EQ. 0.E0) .AND. + (SLA .EQ. 0.E0) ) GOTO 9999 * Check for illogical values IF (RDLA .LT. 0.D0) GOTO 9940 IF (RDLA .GT. 90.D0) GOTO 9950 IF (RMLA .LT. 0.E0 .OR. RMLA .GT. 60.E0) GOTO 9950 IF ( SLA .LT. 0.E0 .OR. SLA .GT. 60.E0) GOTO 9950 *********** * LONGITUDE *********** WRITE (LUOUT,*) ' What is its NAD 83 longitude?', + ' (Longitude is positive west.)' WRITE (LUOUT,*) ' ' READ (LUIN,170,ERR=9930,IOSTAT=IOS) ANS IF (ANS .EQ. B40) GOTO 9999 DUM = ANS CALL NBLANK (DUM, IFLAG2, N2) LENG = N2 RDLO = DCARD( DUM(1:N2), LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 IF (LENG .GT. 0) THEN RMLO = RCARD( DUM, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 IF (LENG .GT. 0) THEN SLO = RCARD( DUM, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 ELSE SLO = 0.E0 ENDIF ELSE RMLO = 0.E0 SLO = 0.E0 ENDIF IF ( (RDLO .EQ. 0.D0) .AND. (RMLO .EQ. 0.E0) .AND. + (SLO .EQ. 0.E0) ) GOTO 9999 * Check for illogical values IF (RDLO .LT. 0.D0) GOTO 9940 IF (RDLO .GT. 360.D0) GOTO 9960 IF (RMLO .LT. 0.E0 .OR. RMLO .GT. 60.E0) GOTO 9960 IF ( SLO .LT. 0.E0 .OR. SLO .GT. 60.E0) GOTO 9960 * Calculate decimal degrees YPT = RDLA + DBLE(RMLA)/60.D0 + DBLE(SLA)/3600.D0 XPT = RDLO + DBLE(RMLO)/60.D0 + DBLE(SLO)/3600.D0 * Get degrees, minutes, seconds CALL HMS (YPT, IDLA, IMLA, SLA) CALL HMS (XPT, IDLO, IMLO, SLO) 9000 RETURN * Error messages 9930 CONTINUE CALL NBLANK (ANS, IFLAG1, N1) CALL NBLANK (ANS, IFLAG2, N2) WRITE (LUOUT,9935) ANS(N1:N2) 9935 FORMAT (' ERROR - in the answer:', /, + 9X, '''', A, '''', /, + ' Must enter number in prescribed format!', /) NOPT = .TRUE. GOTO 9000 9940 WRITE (LUOUT,9945) CALL NBLANK (ANS, IFLAG1, N1) CALL NBLANK (ANS, IFLAG2, N2) WRITE (LUOUT,9945) ANS(N1:N2) 9945 FORMAT (' ERROR - in the answer:', /, + 9X, '''', A, '''', /, + ' Latitude and Longitudes must be positive!', /, + ' Longitude is positive west.', /) NOPT = .TRUE. GOTO 9000 9950 CONTINUE CALL NBLANK (ANS, IFLAG1, N1) CALL NBLANK (ANS, IFLAG2, N2) WRITE (LUOUT,9955) ANS(N1:N2) 9955 FORMAT (' ERROR - Illogical value for latitude in the answer:', /, + ' ''', A, '''', /, + ' Latitude must be between 0 and 90 degrees.', /, + ' Minutes and seconds must be between 0', + ' and 60.', /) NOPT = .TRUE. GOTO 9000 9960 CONTINUE CALL NBLANK (ANS, IFLAG1, N1) CALL NBLANK (ANS, IFLAG2, N2) WRITE (LUOUT,9965) ANS(N1:N2) 9965 FORMAT (' ERROR - Illogical value for longitude in the answer:',/, + ' ''', A, '''', /, + ' Longitude must be between 0 and 360 degrees.',/, + ' Minutes and seconds must be between 0', + ' and 60.', /) NOPT = .TRUE. GOTO 9000 9999 CONTINUE EOF = .TRUE. GOTO 9000 END SUBROUTINE DGRIDS * This subroutine opens the geoid grid(s) using the default grid * names and locations. The default names of the grid areas are * given in DAREAS and the default base file locations are in DFILES * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA, MXDEF PARAMETER (MXAREA =10, MXDEF = MXAREA) CHARACTER*80 B80 CHARACTER*20 B20 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) CHARACTER*65 B65 PARAMETER (B65 = B20//B20//B20//' ') REAL DX1, DY1, XMAX1, XMIN1, YMAX1, YMIN1 INTEGER IDEF, ITEMP, NR1, NC1 CHARACTER*80 DUM CHARACTER*65 AFILE, DFILES(MXAREA) CHARACTER*15 DAREAS(MXDEF) LOGICAL NOGO, GFLAG CHARACTER*15 AREAS COMMON /AREAS/ AREAS(MXAREA) INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA DATA DUM / B80 / * DFILES contains the default locations (pathname) of the grid files * without the .geo extension. (For example 'conus' would * indicate that the conus.geo grid file is in the * current working directory.) The length of each entry in DFILES may * be up to 65 characters. DAREAS contains the default names of these * areas. The names are used internally in the program and in the * program output. They may be no longer than 15 characters. They * must correspond on a one-to-one basis with the file locations in * the DFILES array. That is, the first area name in DAREAS must * be the name that you wish for the first data file set in the * DFILES array. You may, of course, have the arrays the same if * the location of the data file is no longer than 15 characters. * The locations of the grid files may be differ for each * installation. If the pathnames are not correct DFILES (and, possibly, * DAREAS) may be changed and the program recompiled. DATA DFILES /'conus', 'hawaii', 'prvi', 'alaska', + ' ', ' ', ' ', ' ', ' ', ' '/ DATA DAREAS /'Conus OSU89B', 'Hawaii OSU89B', + 'PR & VI OSU89B', 'Alaska OSU89B', + ' ', ' ', ' ', ' ', ' ', ' '/ GFLAG = .FALSE. WRITE (LUOUT, 80) 80 FORMAT (/, ' Default Data Grids', /, + ' # AREA NAME', /, 1X, 79('=') ) DO 140 IDEF = 1, MXDEF AFILE = DFILES(IDEF) IF (AFILE .EQ. B65) GOTO 999 * Try to open a set of default files. * Do not print error messages for non-existing files. ITEMP = NAREA + 1 CALL OPENFL (AFILE, ITEMP, NOGO, DX1, DY1, + XMAX1, XMIN1, YMAX1, YMIN1, NR1, NC1, DUM, GFLAG) IF (.NOT. NOGO) THEN * Set of files opened OK and variables read NAREA = ITEMP AREAS(NAREA) = DAREAS(IDEF) DX(NAREA) = DX1 DY(NAREA) = DY1 XMAX(NAREA) = XMAX1 XMIN(NAREA) = XMIN1 YMAX(NAREA) = YMAX1 YMIN(NAREA) = YMIN1 NR(NAREA) = NR1 NC(NAREA) = NC1 WRITE (LUOUT,120) NAREA, AREAS(NAREA) 120 FORMAT (2X, I2, 2X, A15) ENDIF 140 CONTINUE 999 RETURN END SUBROUTINE DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG) * This subroutine prints out a small diagram showing the area * that was computed. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION XTEMP, XSMALL, XBIG, YSMALL, YBIG REAL SLOMIN, SLOMAX, SLAMIN, SLAMAX INTEGER LODMIN, LOMMIN, LODMAX, LOMMAX INTEGER LADMIN, LAMMIN, LADMAX, LAMMAX INTEGER LU, NCONV XTEMP = XSMALL XSMALL = XBIG XBIG = XTEMP CALL HMS (XBIG, LODMIN, LOMMIN, SLOMIN) CALL HMS (XSMALL, LODMAX, LOMMAX, SLOMAX) CALL HMS (YSMALL, LADMIN, LAMMIN, SLAMIN) CALL HMS (YBIG, LADMAX, LAMMAX, SLAMAX) WRITE (LU,90) NCONV 90 FORMAT (//, ' The total number of geoid estimates: ', I8) WRITE (LU,100) 100 FORMAT (//, 30X, 'Region of Estimation') WRITE (LU,1000) LODMAX, LOMMAX, SLOMAX, LODMIN, LOMMIN, SLOMIN, + LADMAX, LAMMAX, SLAMAX, LADMAX, LAMMAX, SLAMAX, + LADMIN, LAMMIN, SLAMIN, LADMIN, LAMMIN, SLAMIN, + LODMAX, LOMMAX, SLOMAX, LODMIN, LOMMIN, SLOMIN 1000 FORMAT (5(/), T4, 'NAD 83', /, + T4, 'Longitude:', 8X, I4, 1X, I2.2, 1X, F6.3, 13X, I4, 1X, + I2.2, 1X, F6.3, /, + T4, 'Latitude:', 9X, I4, 1X, I2.2, 1X, F6.3, + ' ************', I4, 1X, I2.2, 1X, F6.3, /, + 5(T27, 3X, '*', T57, '*', /), + T27, 3X, '*', 9X, ' NAD 83', T57, '*', /, + T27, 3X, '*', 9X, ' data ', T57, '*', /, + T27, 3X, '*', 9X, ' points', T57, '*', /, + 5(T27, 3X, '*', T57, '*', /), + T4, 'Latitude:', 9X, I4, 1X, I2.2, 1X, F6.3, + ' ************', I4, 1X, I2.2, 1X, F6.3, /, + T4, 'Longitude:', 8X, I4, 1X, I2.2, 1X, F6.3, 13X, I4, 1X, + I2.2, 1X, F6.3, //) RETURN END SUBROUTINE ENDREP (IPAGE, NCONV, ITYPE, GHT, AGHT, VGHT, + SDGHT, SDGHT2, SMGHT, BGGHT, + XSMALL, XBIG, YSMALL, YBIG) *** Gather statistics and write the end-of-program report * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG DOUBLE PRECISION AGHT, VGHT DOUBLE PRECISION SDGHT, SDGHT2 REAL GHT, SMGHT, BGGHT INTEGER IPAGE, NCONV, ITYPE INTEGER LU, I LOGICAL PAGE REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) PAGE = .TRUE. IPAGE = IPAGE + 1 ******************* * DO THE STATISTICS ******************* IF (NCONV .GE. 2) THEN * calculate mean, variance, standard deviation for geoid height * estimates (meters). ******************* * GEOID HT (METERS) ******************* AGHT = SDGHT/DBLE(NCONV) VGHT = SDGHT2/DBLE(NCONV-1) - SDGHT**2/DBLE( NCONV*(NCONV-1 )) IF (VGHT .GT. 1.0D-6) THEN SDGHT = DSQRT(VGHT) ELSE VGHT = 0.0D0 SDGHT = 0.0D0 ENDIF ENDIF IF (NCONV .LT. 2) THEN AGHT = GHT VGHT = 0.0D0 SDGHT = 0.0D0 ENDIF *************************************** * PRINT OUT THE STATISTICS FOR THIS JOB *************************************** LU = LUOUT IF (NCONV .GT. 0) THEN ************************************************* * FIRST REPORT THE FINAL STATISTICS TO THE SCREEN ************************************************* CALL REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT, IPAGE, PAGE) CALL DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG) **************************************************** * NOW REPORT THE FINAL STATISTICS TO THE OUTPUT FILE **************************************************** IF (ITYPE .EQ. 0) THEN ************************************** * INTERACTIVE USE ONLY - NO INPUT FILE ************************************** LU = NOUT CALL REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT, + IPAGE, PAGE) CALL DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG) IF (NCONV .EQ. 0) THEN DO 1007 I = 1, 2 WRITE (LU,*) ' All of your NAD 83 LOCATIONS are out of', + ' bounds.' LU = NOUT 1007 CONTINUE ENDIF ELSEIF (ITYPE .EQ. 1) THEN * For file format types 1 and 2 LU = NOUT CALL REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT, + IPAGE, PAGE) CALL DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG) ELSEIF (ITYPE .EQ. 2) THEN * ITYPE = 2, (free format input, free format output) does not have * anything written to the output file in this subroutine ELSEIF (ITYPE .EQ. 3 .OR. ITYPE .EQ. 4 .OR. + ITYPE .EQ. 5 .or. ITYPE .EQ. 6) THEN * ITYPE=3, ITYPE=4, and ITYPE=5 & 6 the NGS Horizontal Blue Book * file formats do not have a * report written to the output file ENDIF ENDIF RETURN END SUBROUTINE FGRID (XPT, YPT, DX, DY, XMAX, XMIN, + YMAX, YMIN, XGRID, YGRID, IROW, JCOL, NOGO) ********************************************************************** ** SUBROUTINE FGRID: IDENTIFIES THE LOCAL GRID SQUARE FOR INTERP. * ********************************************************************** * This subroutine is designed to identify the grid square in which a * particular point is located and get the corner coordinates * converted into the index coordinate system. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION XPT, YPT REAL DX, DY, XMAX, XMIN, YMAX, YMIN, XGRID, YGRID INTEGER IROW, JCOL LOGICAL NOGO NOGO = .FALSE. * Check to see it the point is outside the area of the gridded data IF (XPT .GE. DBLE(XMAX) .OR. XPT .LE. DBLE(XMIN) .OR. + YPT .GE. DBLE(YMAX) .OR. YPT .LE. DBLE(YMIN) ) THEN NOGO = .TRUE. * WRITE (*,*) '***THE POINT IS OUT OF BOUNDS***' GOTO 200 ENDIF * Calculate the coordinate values for the point to be interpolated * in terms of grid indices XGRID = SNGL(XPT - DBLE(XMIN) )/DX + 1.E0 YGRID = SNGL(YPT - DBLE(YMIN) )/DY + 1.E0 * Find the I,J values for the SW corner of the local square IROW = IFIX(YGRID) JCOL = IFIX(XGRID) 200 RETURN END SUBROUTINE FHELP *** Print information about the formats of the input data *** file types used by GEOID. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA, ITYPE, IOS PARAMETER (MXAREA =10) CHARACTER*1 ANS INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) ***************************** * CHOSE THE INPUT FILE FORMAT ***************************** 9001 WRITE (LUOUT,*) ' ' WRITE (LUOUT,*) ' What format do you want information about?' WRITE (LUOUT,*) ' 1) Free Format Type 1' WRITE (LUOUT,*) ' 2) Free Format Type 2' WRITE (LUOUT,*) ' 3) NGS Blue Book (old) format (default)' WRITE (LUOUT,*) ' 4) NGS Blue Book (new *84*) format' WRITE (LUOUT,*) ' 5) NGS Blue Book (new *86*) format' WRITE (LUOUT,*) ' 6) NGS Blue Book (extended) format' READ (LUIN,'(A1)') ANS IF (ANS .EQ. ' ') THEN ITYPE = 3 ELSE READ (ANS,90,ERR=9940,IOSTAT=IOS) ITYPE 90 FORMAT (I1) IF (ITYPE .GT. 6 .OR. ITYPE .LT. 1) THEN * Not a good answer - Try again. WRITE (LUOUT,*) ' Gotta pick one of these -', + ' sorry try again.' GOTO 9001 ENDIF ENDIF * Print information IF (ITYPE .EQ. 1) THEN ******************************* * FOR FILE FORMAT ITYPE = 1 * FREE FORMAT TYPE1 INPUT FILE * PRETTY OUTPUT FORMAT ******************************* WRITE (LUOUT,2) * 2 FORMAT ('1') 2 FORMAT (' ') WRITE (LUOUT, 110) 110 FORMAT (' Free Format Type 1 - The first 40 characters of', + ' the input data record may', /, + ' contain the station name or be blank. The rest of', + ' the record (columns 41-80)', /, + ' must contain the latitude and longitude. They may', + ' be given in (1) decimal', /, + ' degrees; (2) integer degrees and decimal minutes,', + ' or (3) integer degrees,', /, + ' integer minutes, and decimal seconds. The decimal', + ' portion of the latitude', /, + ' MUST contain a decimal point as it is used to', + ' determine which is the last', /, + ' number forming part of the latitude. The output', + ' will be in "pretty" format.', /) WRITE (LUOUT, 120) 120 FORMAT (' The following three records are examples of valid', + ' input records:', /, + ' <------------ Columns 1-40 ------------>', + '<------------ Columns 41-80----------->', /, + ' AAA 34.', + '4444444 98.8888888', /, + ' BBB 25', + ' 55.55555 76 56.66666', /, + ' CCC 45 45', + ' 45.555 111 11 11.111', /) WRITE (LUOUT,931) 931 FORMAT (15X, ' (Hit RETURN to continue.)') READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) WRITE (LUOUT, 130) 130 FORMAT (' The following is an example of the output. Note', + ' that with Free Format Type 1', /, + ' data, the input latitude and', + ' longitude are expressed in', /, + ' degrees, minutes, and seconds regardless of the', + ' method of input.', /) WRITE (LUOUT, 140) 140 FORMAT (' Computation #: 1', + ' Region: Conus OSU89B', //, + ' Station Name: AAA', //, + ' Latitude ', + 'Longitude Geoid Height', /, + ' ', + ' (meters)', //, + ' 34 26 39.99984 98 53', + ' 19.99968 -26.868', /) ELSEIF (ITYPE .EQ. 2) THEN ******************************* * FOR FILE FORMAT ITYPE = 2 * FREE FORMAT TYPE2 INPUT FILE * FREE FORMAT OUTPUT ******************************* WRITE (LUOUT,2) WRITE (LUOUT, 210) 210 FORMAT (' Free Format Type 2 - The first 32 characters of', + ' the input data record must', /, + ' contain the latitude and longitude. They may be', + ' given in (1) decimal degrees;', /, + ' (2) integer degrees and decimal minutes, or (3)', + ' integer degrees, integer', /, + ' minutes, and decimal seconds. The decimal portion', + ' of the latitude MUST contain', /, + ' a decimal point as it is used to determine which', + ' is the last number forming', /, + ' part of the latitude. The next eight characters', + ' (columns 33-40) will have the', /, + ' calculated geoid height while the rest of the input', + ' record (columns 41-80) may', /, + ' contain the station name or be blank. The output', + ' will be in the same format as', /, + ' the input but will contain the geoid height.', /) WRITE (LUOUT,931) READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) WRITE (LUOUT, 220) 220 FORMAT (' The following three records are examples of', + ' valid input records:', //, + ' <------------ Columns 1-40 ------------>', + '<------------ Columns 41-80----------->', /, + ' 45 45 45.55555 111 11 11.11111 ', + 'one', /, + ' 25 55.5555555 76 56.6666666 ', + 'two', /, + ' 34.444444444 98.888888888 ', + 'three') WRITE (LUOUT, 230) 230 FORMAT (/, ' The following is an example of the output.', //, + ' GEOID Version 2.00 - NGS', + ' GEOID93 Model', /, + ' 45 45 45.55555 111 11 11.11111 -10.859one', /, + ' 25 55.5555555 76 56.6666666 -39.988two', /, + ' 34.444444444 98.888888888 -26.868three', /) ELSEIF (ITYPE .EQ. 3) THEN **************************************** * FOR INPUT FILE ITYPE = 3 * THE HORIZONTAL BLUE BOOK SPECIFICATION **************************************** * The following two format statements logically belong together * but are split in order to have fewer that 19 continuation lines WRITE (LUOUT,2) WRITE (LUOUT, 310) 310 FORMAT (' NGS Horizontal ''Old'' Blue Book format - *80*', + ' (Control Point) and *84* (Geoid', /, + ' Height) Records. GEOID uses information from the', + ' *80* and *84* records in Blue', /, + ' Book files. The Station Serial Number (SSN), the', + ' latitude and the longitude', /, + ' are read from the *80* records, the rest of the', + ' record is ignored. Only the', /, + ' *84* record is modified by GEOID - the *80*', + ' records and all other records are') WRITE (LUOUT, 311) 311 FORMAT (' passed through without change to the output file.', + ' On the *84* records, the SSN', /, + ' is used to determine which is the corresponding', + ' the *80* record. The', /, + ' information to the right of the SSN (columns', + ' 15-80) in the *84* record is', /, + ' created by GEOID. On the *80* records, the', + ' direction of the latitude must', /, + ' be north positive (''N'') and the direction of the', + ' longitude must be west', /, + ' positive (''W'').', /) WRITE (LUOUT, 320) 320 FORMAT (' This format is out of date and should no longer be', + ' used for submitting data', /, + ' to the National Geodetic Survey.', //, + ' This format was defined in:', /, + ' ''Input Formats and Specifications of the', + ' National Geodetic Survey Data Base''', /, + ' ''Volume 1. Horizontal Control Data''.', /, + ' Published by the Federal Geodetic Control Committee', + ' in 1980 and available', /, + ' from: the National Geodetic Survey, NOAA, Rockville,', + 'MD 20852.', /) WRITE (LUOUT,931) READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) WRITE (LUOUT, 330) 330 FORMAT (' The following input example is *80* and *84* records', + ' from an ''old'' Blue Book', /, + ' file:', //, + '004560*80*096 KNOXVILLE CT HSE', + ' 411906578 N0930548534 W 277 MIA33', /, + '004565*84*096 RAPP - OSU86G ', + ' -319 10') WRITE (LUOUT, 340) 340 FORMAT (/, ' The following example is of the output *84*', + ' record with the new, interpolated', /, + ' geoid height.', //, + '004565*84*096 NGS PROGRAM GEOID Conus OSU89B', + ' -316 5', /) ELSEIF (ITYPE .EQ. 4) THEN * The following two format statements logically belong together * but are split in order to have fewer that 19 continuation lines WRITE (LUOUT,2) WRITE (LUOUT,410) 410 FORMAT (' NGS Horizontal ''New'' Blue Book format - *80*', + ' (Control Point) and *84* (Geoid', /, + ' Height) Records. GEOID uses information from the', + ' *80* and *84* records in Blue', /, + ' Book files. The Station Serial Number (SSN), the', + ' latitude and the longitude', /, + ' are read from the *80* records, the rest of the', + ' record is ignored. Only the', /, + ' *84* record is modified by GEOID - the *80*', + ' records and all other records are') WRITE (LUOUT, 411) 411 FORMAT (' passed through without change to the output file.', + ' On the *84* records, the SSN', /, + ' is used to determine which is the corresponding', + ' the *80* record. The', /, + ' information to the right of the SSN (columns', + ' 15-80) in the *84* record is', /, + ' created by GEOID. On the *80* records, the', + ' direction of the latitude must', /, + ' be north positive (''N'') and the direction of the', + ' longitude must be west', /, + ' positive (''W'').', /) WRITE (LUOUT, 420) 420 FORMAT (' For more information on this format,', + ' please refer to:', /, + ' ''Input Formats and Specifications of the', + ' National Geodetic Survey Data Base''', /, + ' ''Volume 1. Horizontal Control Data''.', /, + ' Published by the Federal Geodetic Control Committee', + ' in January, 1989 and', /, + ' available from: the National Geodetic Survey, NOAA,', + ' Rockville, MD 20852.', /) WRITE (LUOUT,931) READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) WRITE (LUOUT, 430) 430 FORMAT (' The following input example is *80* and *84* records', + ' from a ''new'' Blue Book', /, + ' file:', //, + '004560*80*0096KNOXVILLE CT HSE ', + '411906578 N0930548534 W 277 MIA33', /, + '004565*84*0096RAPP ', + ' -3190 100') WRITE (LUOUT, 440) 440 FORMAT (/, ' The following example is of the output *84*', + ' record with the new, interpolated', /, + ' geoid height.', //, + '004565*84*0096NGS PROGRAM GEOID Conus OSU89B', + ' -3164 50', /) ********************************************************************* *** *86* records ********************************************************************* ELSEIF (ITYPE .EQ. 5) THEN * The following two format statements logically belong together * but are split in order to have fewer that 19 continuation lines WRITE (LUOUT,2) WRITE (LUOUT,510) 510 FORMAT (' NGS Horizontal ''New'' Blue Book format - *80*', + ' (Control Point) and *86* (Combo', /, + ' Height) Records. GEOID uses information from the', + ' *80* and *86* records in Blue', /, + ' Book files. The Station Serial Number (SSN), the', + ' latitude and the longitude', /, + ' are read from the *80* records, the rest of the', + ' record is ignored. Only the', /, + ' *86* record is modified by GEOID - the *80*', + ' records and all other records are') WRITE (LUOUT, 511) 511 FORMAT (' passed through without change to the output file.', + ' Any *84* records are ', /, + ' deleted. Information in columns', + ' 36-43) in the *86* record is', /, + ' created by GEOID. On the *80* records, the', + ' direction of the latitude must', /, + ' be north positive (''N'') and the direction of the', + ' longitude must be west', /, + ' positive (''W'').', /) WRITE (LUOUT, 520) 520 FORMAT (' For more information on this format,', + ' please refer to:', /, + ' ''Input Formats and Specifications of the', + ' National Geodetic Survey Data Base''', /, + ' ''Volume 1. Horizontal Control Data''.', /, + ' Published by the Federal Geodetic Control Committee', + ' in January, 1989 and', /, + ' available from: the National Geodetic Survey, NOAA,', + ' Rockville, MD 20852.', /) WRITE (LUOUT,931) READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) WRITE (LUOUT, 530) 530 FORMAT (' The following input example is *80* and *86* records', + ' from a ''new'' Blue Book', /, + ' file:', //, + '004560*80*0096KNOXVILLE CT HSE ', + '411906578 N0930548534 W 277 MIA33', /, + '004565*86*0096 any values in here ', + 'will be reproduced if ssn matches ') WRITE (LUOUT, 540) 540 FORMAT (/, ' The following example is of the output *86*', + ' record with the new, interpolated', /, + ' geoid height.', //, + '004565*86*0096 any values in here -31768D ', + 'will be reproduced if ssn matches ',/) ********************************************************************* ********************************************************************* ELSEIF (ITYPE .EQ. 6) THEN WRITE (LUOUT,2) WRITE (LUOUT,610) 610 FORMAT (' NGS Horizontal ''Extended'' Blue Book format', + ' (BIGADJUST format) - *80* (Control', /, + ' Point) and *84* (Geoid Height) Records. Only the', + ' *80* records are used from', /, + ' the input file. The output consists only of *84*' + ' records. The SSN''s for the', /, + ' *84* records are taken from the *80* records - that', + ' is, there is an *84*', /, + ' record for each *80* record, regardless of any', + ' *84* records in the input file.', //, + ' This format is for internal NGS use only!!!', /) ENDIF WRITE (LUOUT,*) ' Do you want more information (Y/N)?' WRITE (LUOUT,*) ' (Default is Y)' READ (LUIN,'(A1)') ANS IF (ANS .NE. 'N' .AND. ANS .NE. 'n') GOTO 9001 RETURN * Error message 9940 WRITE (LUOUT,*) ' Gotta pick ''1'' through ''6'' -', + ' sorry try again.' GOTO 9001 END SUBROUTINE GETPT (NCONV, ITYPE, NAME, IDLA, IMLA, SLA, IDLO, + IMLO, SLO, XPT, YPT, EOF, NOPT) * Get the name, latitude, and longitude of a point either interactively * or from an input data file * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) DOUBLE PRECISION XPT, YPT REAL SLA, SLO INTEGER NCONV, ITYPE INTEGER IDLA, IMLA, IDLO, IMLO CHARACTER*80 NAME CHARACTER*1 ANS LOGICAL EOF, NOPT INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) EOF = .FALSE. NOPT = .FALSE. IF (ITYPE .EQ. 0) THEN ************************************* * FOR INTERACTIVE USE - NO INPUT FILE ************************************* IF (NCONV .GE. 1) THEN WRITE (LUOUT,*) ' Do you want to do another geoid', + ' computation (Y/N)?' WRITE (LUOUT,*) ' (Default is Y)' READ (LUIN,'(A1)') ANS IF (ANS .EQ. 'n' .OR. ANS .EQ. 'N') GOTO 9999 ENDIF * Get a point (X,Y) to compute CALL ASKPT (NCONV, NAME, IDLA, IMLA, SLA, + IDLO, IMLO, SLO, XPT, YPT, EOF, NOPT) IF (NOPT) GOTO 9000 ELSEIF (ITYPE .EQ. 1) THEN * Free format type 1 CALL TYPE1 (NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO, + XPT, YPT, EOF, NOPT) ELSEIF (ITYPE .EQ. 2) THEN * Free format type 2 CALL TYPE2 (IDLA, IMLA, SLA, IDLO, IMLO, SLO, + XPT, YPT, EOF, NOPT) ELSEIF (ITYPE .EQ. 3 .OR. ITYPE .EQ. 4 .OR. + ITYPE .EQ. 5 .or. ITYPE .EQ. 6 ) THEN * ITYPE=3 is NGS old Horizontal Blue Book * ITYPE=4 is NGS new Horizontal Blue Book with *84* * ITYPE=5 is NGS new Horizontal Blue Book with *86* * ITYPE=6 is NGS extended Horizontal Blue Book CALL TYPE3 (ITYPE, NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO, + XPT, YPT, EOF, NOPT) ENDIF ********************************************************* * CHANGE THE LONGITUDE TO POSITIVE EAST FOR INTERPOLATION ********************************************************* XPT = -XPT 9000 RETURN * End of file 9999 CONTINUE EOF = .TRUE. GOTO 9000 END SUBROUTINE GRIDS * This subroutine opens the geoid grid(s) that are requested in * a file named 'AREA.PAR' (if it exists). * AREA.PAR (or area.par) will be read for the names * and locations of the gridded geoid height data set files. * The data in the AREA.PAR file has the following format: * Columns 1-15 contain the name of the area. This name may * contain blanks or any other characters. * Columns 16-80 (the rest of the record) contain the base name of the * grid files. That is the '.geo' extension is not * included - it is added by this subroutine before each file is * opened. * Comments are indicated by an '*' in column 1, blank records are * also treated as comments. Comment records are printed to the * output file but otherwise ignored. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) CHARACTER*80 B80 CHARACTER*20 B20 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) REAL DX1, DY1, XMAX1, XMIN1, YMAX1, YMIN1 INTEGER NR1, NC1 INTEGER IOS, I, IFLAG1, IFLAG2, N1, N2, LENG, IERR, ITEMP CHARACTER*80 CARD, CCARD, DUM CHARACTER*65 AFILE, GFILE CHARACTER*15 AAREA LOGICAL NOGO, GFLAG CHARACTER*15 AREAS COMMON /AREAS/ AREAS(MXAREA) INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA DATA DUM / B80 / DATA IFLAG1 /1/, IFLAG2 /2/ *********************************************************************** * TRY AND OPEN THE 'AREA.PAR' FILE, IF NOT THEN TRY AND OPEN 'area.par' *********************************************************************** GFILE = 'AREA.PAR' OPEN (NAPAR,FILE=GFILE,FORM='FORMATTED',STATUS='OLD', + ACCESS='SEQUENTIAL',ERR=9100,IOSTAT=IOS) * File containing grid names exits 10 GFLAG = .TRUE. WRITE (LUOUT, 80) 80 FORMAT (/, ' Data Grids named in the AREA.PAR file', /, + ' # AREA NAME LOCATION', /, 1X, 79('=') ) DO 120 I = 1, MXAREA 100 READ (NAPAR,110,ERR=9000,END=9000,IOSTAT=IOS) CARD 110 FORMAT (A80) * Check for comment records and blank records IF ( CARD(1:1) .EQ. '*' ) THEN CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,'(5X, A)') CARD(1:N2) GOTO 100 ELSEIF ( CARD .EQ. B80 ) THEN WRITE (LUOUT,*) ' ' GOTO 100 ENDIF * Get area name and basename of file (i.e. location without extensions) DUM = CARD AAREA = DUM(1:15) CALL NBLANK (CARD(16:80), IFLAG1, N1) DUM(1:65) = CARD(15+N1:80) LENG = 65 AFILE = CCARD(DUM, LENG, IERR) IF (IERR .NE. 0) STOP 'Coding error in GRIDS -- AFILE' IF (AFILE .EQ. ' ') GOTO 9000 * Try and open files ITEMP = NAREA + 1 CALL OPENFL (AFILE, ITEMP, NOGO, DX1, DY1, + XMAX1, XMIN1, YMAX1, YMIN1, NR1, NC1, CARD, GFLAG) IF (.NOT. NOGO) THEN * Files opened OK and variables read NAREA = ITEMP AREAS(NAREA) = AAREA DX(NAREA) = DX1 DY(NAREA) = DY1 XMAX(NAREA) = XMAX1 XMIN(NAREA) = XMIN1 YMAX(NAREA) = YMAX1 YMIN(NAREA) = YMIN1 NR(NAREA) = NR1 NC(NAREA) = NC1 CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,140) NAREA, CARD(1:N2) 140 FORMAT (2X, I2, 2X, A) ENDIF 120 CONTINUE 9000 RETURN * 'AREA.PAR' does not exist, try the name 'area.par' * If it exists, open it and continue, if it does not exist, return. 9100 CONTINUE GFILE = 'area.par' OPEN (NAPAR,FILE=GFILE,FORM='FORMATTED',STATUS='OLD', + ACCESS='SEQUENTIAL',ERR=9000,IOSTAT=IOS) GOTO 10 END SUBROUTINE HEADR (VRSION) *** This subroutine prints the header information and the disclaimer * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) REAL VRSION CHARACTER*1 ANS INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) WRITE (LUOUT,920) 920 FORMAT (10X, ' Welcome', /, + 10X, ' to the', /, + 10X, ' National Geodetic Survey', /, + 10X, ' Geoid Computation Program', //, + 10X, ' NGS GEOID93 Model', /) WRITE (LUOUT,930) VRSION 930 FORMAT (//,10X, ' (Version', F5.2, ')', /, + 10X, ' March 10, 1993', /, + 10X, ' Dennis Milbert, Ph.D.',/, + 10X, ' Donald Schultz',/, + 10X, ' Advanced Geodetic Science Branch', /) WRITE (LUOUT,931) 931 FORMAT (10X, ' (Hit RETURN to continue.)') READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) * 2 FORMAT ('1') 2 FORMAT (' ') WRITE (LUOUT,932) 932 FORMAT (/, ' DISCLAIMER', //, + ' This program and supporting information is furnished by', + ' the government of', /, + ' the United States of America, and is accepted/used by the', + ' recipient with', /, + ' the understanding that the U. S. government makes no', + ' warranties, express or', /, + ' implied, concerning the accuracy, completeness, reliability,', + ' or suitability', /, + ' of this program, of its constituent parts, or of any', + ' supporting data.', //, + ' The government of the United States of America shall be', + ' under no liability', /, + ' whatsoever resulting from any use of this program.', + ' This program should', /, + ' not be relied upon as the sole basis for solving a problem', + ' whose incorrect', /, + ' solution could result in injury to person or property.') WRITE (LUOUT,933) 933 FORMAT ( /, + ' This program is the property of the government of the', + ' United States of', /, + ' America. Therefore, the recipient further agrees not to', + ' assert proprietary', /, + ' rights therein and not to represent this program to anyone as', + ' being other', /, + ' than a government program.', /) WRITE (LUOUT,931) READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) RETURN END SUBROUTINE HMS (DD, ID, IM, S) * Use this to change from decimal degrees (double precision) * to integer degrees, integer minutes, and decimal seconds (real) * Seconds are assumed to have no more than 5 decimal places * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) REAL SMALL PARAMETER (SMALL = 1.E-5) DOUBLE PRECISION DD, TEMP REAL S INTEGER ID, IM ID = IDINT(DD) TEMP = ( DD - DBLE(ID) )*60.0D0 IM = IDINT(TEMP) S = SNGL( ( TEMP - DBLE(IM) )*60.0D0 ) IF (IM .EQ. 60) THEN IM = 0 ID = ID + 1 ENDIF IF (S .LT. SMALL) S = 0.E0 IF (S .GT. (60.E0-SMALL) ) THEN S = 0.E0 IM = IM + 1 ENDIF RETURN END SUBROUTINE INITL (SCREEN, PAGE, NAME, IPAGE, ITYPE, + SMGHT, BGGHT, AGHT, VGHT, SDGHT, SDGHT2, + XSMALL, XBIG, YSMALL, YBIG) *** This subroutine initializes all the variables needed in GEOID * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) CHARACTER*20 B20 CHARACTER*80 B80 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG DOUBLE PRECISION AGHT, VGHT DOUBLE PRECISION SDGHT, SDGHT2 REAL SMGHT, BGGHT INTEGER IPAGE, ITYPE CHARACTER*80 NAME LOGICAL PAGE, SCREEN CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) * Initialize card variable in common CURNT to blank CARD = B80 * Set the logical units for input/output common INOUT * Note that the unit numbers for the data grids (LUAREA) are defined * in subroutine OPENFL as the numbers 11 through MXAREA LUIN = 5 LUOUT = 6 NOUT = 101 NIN = 102 NAPAR = 103 ****************************************************************** * INITIALIZE ****************************************************************** * Defaults: SCREEN = .TRUE. => send results to screen * PAGE = .FALSE. => don't start a new page in the output file * NAME = ' ' => no station name * IPAGE = 0 => current output file page number is 0 * ITYPE = 0 => interactive input of points SCREEN = .TRUE. PAGE = .FALSE. NAME = ' ' IPAGE = 0 ITYPE = 0 SMGHT = 1.0E10 BGGHT = -1.0E10 AGHT = 0.0D0 VGHT = 0.0D0 SDGHT = 0.0D0 SDGHT2 = 0.0D0 XSMALL = 1.0D10 XBIG = -1.0D10 YSMALL = 1.0D10 YBIG = -1.0D10 RETURN END SUBROUTINE INTERP (NOGO, RESP, XPT, YPT, GHT, ITYPE) *** COMPUTE BIQUADRATIC INTERPOLATION, FLOATING POINT VERSION *** LOGIC NOT TESTED FOR '0/360 MERIDIAN WRAPAROUND' * THIS CODE IS ADAPTED FROM THAT OF DENNIS MILBERT - 6/5/90 * DENNIS(XPT,YPT,DX,DY,XMAX,XMIN,YMAX,YMIN, * & ZEE,NR,NC,LUAREA) * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MNCBIN, MXAREA PARAMETER (MNCBIN = 1000, MXAREA =10) CHARACTER*20 B20 CHARACTER*80 B80 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) DOUBLE PRECISION XPT, YPT REAL BUF, GHT REAL X, Y, FX0, FX1, FX2, QTERP2 REAL DX0, DY0, XMAX0, YMAX0, XMIN0, YMIN0 REAL XGRID, YGRID, VAL REAL TEE1, TEE2, TEE3, TEE4, TEE5, TEE6, TEE7, TEE8, TEE9 INTEGER NR0, NC0 INTEGER IAREA, IROW, JCOL, IFILE, IDUM, I, J, IDIM INTEGER IX, IX1, IX2, JY, JY1, JY2 INTEGER IFLAG1, IFLAG2, N1, N2, ITYPE CHARACTER*15 RESP LOGICAL NOGO, FLAG DIMENSION BUF(MNCBIN) CHARACTER*15 AREAS COMMON /AREAS/ AREAS(MXAREA) CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA DATA IFLAG1 /1/, IFLAG2 /2/, FLAG /.FALSE./ ****************************************************************** * INITIALIZE ****************************************************************** NOGO = .FALSE. **************************************************** * READ WHERE TO GET THE DATA AND HOW IT IS ORGANIZED **************************************************** * CHECK TO SEE WHICH SET OF GRIDDED FILES XPT,YPT is in. DO 100 IAREA = 1, NAREA DX0 = DX(IAREA) DY0 = DY(IAREA) XMAX0 = XMAX(IAREA) XMIN0 = XMIN(IAREA) YMAX0 = YMAX(IAREA) YMIN0 = YMIN(IAREA) NR0 = NR(IAREA) NC0 = NC(IAREA) CALL FGRID (XPT, YPT, DX0, DY0, XMAX0, XMIN0, + YMAX0, YMIN0, XGRID, YGRID, IROW, JCOL, NOGO) IF (.NOT. NOGO) GOTO 200 100 CONTINUE * Not in any of the grid areas NOGO = .TRUE. GOTO 950 200 CONTINUE * Point in area number IAREA and named AREAS(IAREA) IFILE = LUAREA(IAREA) RESP = AREAS(IAREA) IF (YPT .LT. YMIN0 .OR. YPT .GT. YMAX0 .OR. + XPT .LT. XMIN0 .OR. XPT .GT. XMAX0) THEN VAL=1.D30 ELSE *** Within grid boundaries, get indices in the grid IX = (XPT - XMIN0)/DX0 + 1.D0 IX1 = IX + 1 IX2 = IX + 2 *** Check if edge collision 1 IF (IX2 .GT. NC0) THEN IX = IX - 1 IX1 = IX1 - 1 IX2 = IX2 - 1 GO TO 1 ENDIF X = (XPT - DX0*(IX - 1) - XMIN0)/DX0 *** move grid to get nearer center of 3 X 3 IF (X .LT. 0.5 .AND. IX .GT. 1) THEN IX = IX - 1 IX1 = IX1 - 1 IX2 = IX2 - 1 X = X + 1. ENDIF IF ( X .LT. 0. .OR. X .GT. 2.) STOP 'INTERP - 1' *** Now do it for Y JY = (YPT - YMIN0)/DY0 + 1.D0 JY1 = JY + 1 JY2 = JY + 2 2 IF (JY2 .GT. NR0) THEN JY = JY - 1 JY1 = JY1 - 1 JY2 = JY2 - 1 GO TO 2 ENDIF Y = (YPT - DY0*(JY - 1) - YMIN0)/DY0 IF (Y .LT. 0.5 .AND. JY .GT. 1) THEN JY = JY - 1 JY1 = JY1 - 1 JY2 = JY2 - 1 Y = Y + 1. ENDIF IF (Y .LT. 0. .OR. Y .GT. 2.) STOP 'INTERP - 2' *** Get the nine data points needed from gridded geoid height file !!! *** lower boundary READ (IFILE,REC=JY+1) IDUM, (BUF(J), J=1,NC0) TEE1 = BUF(IX) TEE2 = BUF(IX+1) TEE3 = BUF(IX+2) *** middle row READ (IFILE,REC=JY1+1) IDUM, (BUF(J), J=1,NC0) TEE4 = BUF(IX) TEE5 = BUF(IX+1) TEE6 = BUF(IX+2) *** upper boundary READ (IFILE,REC=JY2+1) IDIM, (BUF(J), J=1,NC0) TEE7 = BUF(IX) TEE8 = BUF(IX+1) TEE9 = BUF(IX+2) *** Compute the value FX0 = QTERP2(X, TEE1, TEE2, TEE3) FX1 = QTERP2(X, TEE4, TEE5, TEE6) FX2 = QTERP2(X, TEE7, TEE8, TEE9) GHT = QTERP2(Y, FX0, FX1, FX2) ENDIF 9999 RETURN * Error Messages 950 CONTINUE IF (ITYPE .NE. 0) THEN CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,955) CARD(N1:N2) 955 FORMAT (' *** THIS POINT IS OUT OF BOUNDS ***', /, + 1X, '''', A, '''') ELSE WRITE (LUOUT,960) 960 FORMAT (' *** THE POINT IS OUT OF BOUNDS ***') ENDIF * Write out grid areas for the first out-of-bounds error message IF (.NOT.FLAG .OR. ITYPE .EQ. 0) THEN WRITE (LUOUT,*) ' It must be within one of the following grid', + ' areas;' WRITE (LUOUT,975) 975 FORMAT (18X, 7X, 'Latitude', 7X, 'Longitude', /, + 5X, 'Area Name', 5X, 2(5X, 'MIN', 4X, 'MAX', 1X), + '(degrees)' ) DO 970 I = 1, NAREA WRITE (LUOUT,965) AREAS(I), NINT( YMIN(I) ), NINT( YMAX(I) ), + NINT( -XMAX(I) ), NINT( -XMIN(I) ) 965 FORMAT (1X, '''', A15, '''', 2(2X, 2I7) ) 970 CONTINUE FLAG = .TRUE. ENDIF WRITE (LUOUT,*) ' ' GOTO 9999 END SUBROUTINE IPARMS (ITYPE, SCREEN) *** This subroutine interactively requests for information *** needed by GEOID. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) CHARACTER*80 B80 CHARACTER*20 B20 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) CHARACTER*80 INFILE, OFILE CHARACTER*1 ANS INTEGER ITYPE INTEGER IFLAG1, IFLAG2, N1, N2, IOS LOGICAL SCREEN, EFLAG INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) DATA IFLAG1 /1/, IFLAG2 /2/ ********************************** * GET THE NAME FOR THE OUTPUT FILE ********************************** 14 WRITE (LUOUT,*) ' What is the name of the file which will', + ' contain the program results?' WRITE (LUOUT,*) ' The default name is ''geoid.out''.' READ (LUIN,'(A80)') OFILE IF (OFILE .EQ. B80) OFILE = 'geoid.out' INQUIRE (FILE=OFILE, EXIST=EFLAG) IF (EFLAG) THEN CALL NBLANK (OFILE, IFLAG1, N1) CALL NBLANK (OFILE, IFLAG2, N2) WRITE (LUOUT,*) ' The file ''', OFILE(N1:N2), '''' WRITE (LUOUT,*) ' already exists. Do you want to overwrite', + ' it (Y/N)?' WRITE (LUOUT,*) ' (Default is Y)' READ (LUIN, '(A1)') ANS IF (ANS .EQ. 'n' .OR. ANS .EQ. 'N') GOTO 14 ENDIF OPEN (NOUT,FILE=OFILE,FORM='FORMATTED',STATUS='UNKNOWN', + ACCESS='SEQUENTIAL',ERR=9910,IOSTAT=IOS) ********************************** * GET THE NAME FOR THE INPUT FILE ********************************** WRITE (LUOUT,*) ' ' 13 WRITE (LUOUT,*) ' Do you have an input data file (Y/N)?' WRITE (LUOUT,*) ' (Default is N)' READ (LUIN,'(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') THEN WRITE (LUOUT,*) ' ' WRITE (LUOUT,*) ' What is the name of the input data file?' WRITE (LUOUT,*) ' The default name is ''BBOOK''.' READ (LUIN,'(A80)') INFILE IF (INFILE .EQ. B80) INFILE = 'BBOOK' OPEN (NIN,FILE=INFILE,FORM='FORMATTED',STATUS='OLD', + ACCESS='SEQUENTIAL',ERR=9920,IOSTAT=IOS) ***************************** * CHOSE THE INPUT FILE FORMAT ***************************** 9001 CONTINUE CALL NBLANK (INFILE, IFLAG2, N2) WRITE (LUOUT,*) ' ' WRITE (LUOUT,*) ' What is your file format?' WRITE (LUOUT,*) ' 0) Help - File format information' WRITE (LUOUT,*) ' 1) Free Format Type 1' WRITE (LUOUT,*) ' 2) Free Format Type 2' WRITE (LUOUT,*) ' 3) NGS Blue Book (old) format (default)' WRITE (LUOUT,*) ' 4) NGS Blue Book (new) format (make *84*)' WRITE (LUOUT,*) ' 5) NGS Blue Book (new) format (make *86*)' WRITE (LUOUT,*) ' 6) NGS Blue Book (extended) format' READ (LUIN,'(A1)') ANS IF (ANS .EQ. ' ') THEN * NGS Horizontal Blue Book format ITYPE = 3 ELSE READ (ANS,347,ERR=9940,IOSTAT=IOS) ITYPE 347 FORMAT (I1) IF (ITYPE .GT. 6 .OR. ITYPE .LT. 0) THEN * Not a good answer - Try again. WRITE (LUOUT,*) ' Gotta pick one of these -', + ' sorry try again.' GOTO 9001 ENDIF ENDIF * Get help information IF (ITYPE .EQ. 0) THEN CALL FHELP GOTO 9001 ENDIF ******************************** * CHECK FOR A SCREEN OUTPUT ALSO ******************************** WRITE (LUOUT,*) ' ' WRITE (LUOUT,*) ' Do you want the results written to the', + ' terminal screen as well as to' WRITE (LUOUT,*) ' the output file (Y/N)?' WRITE (LUOUT,*) ' (Default is N)' READ (LUIN,'(A1)') ANS IF (ANS .NE. 'y' .AND. ANS .NE. 'Y') SCREEN = .FALSE. GOTO 9002 * Error message 9940 WRITE (LUOUT,*) ' Gotta pick ''1'' or ''2'' or ''3'' -', + ' sorry try again.' GOTO 9001 9002 ENDIF RETURN * Error message 9910 CONTINUE CALL NBLANK (OFILE, IFLAG1, N1) CALL NBLANK (OFILE, IFLAG2, N2) WRITE (LUOUT,9915) IOS, OFILE(N1:N2) 9915 FORMAT (' ERROR (', I5, ') - The operating system could not', + ' open the file ', /, + ' ''', A, '''', /, + ' Try again.', /) GOTO 14 9920 CONTINUE CALL NBLANK (INFILE, IFLAG1, N1) CALL NBLANK (INFILE, IFLAG2, N2) WRITE (LUOUT,9915) IOS, INFILE(N1:N2) GOTO 13 END SUBROUTINE MLOOP (NCONV, IPAGE, ITYPE, VRSION, GHT, + SDGHT, SDGHT2, SMGHT, BGGHT, + XSMALL, XBIG, YSMALL, YBIG, + PAGE, SCREEN, NAME) ********************************************************************** * THIS SUBROUTINE LOOPS THROUGH THE INPUT DATA (EITHER AN INPUT DATA * * FILE OR INTERACTIVELY), CALCULATES THE GEOID HEIGHT VALUES, * * UPDATES THE MINIMUM, MAXIMUM, AND STATISTICAL SUMMATIONS, AND THEN * * PRINTS THE RESULTS TO THE OUTPUT FILE AND/OR THE SCREEN. * ********************************************************************** * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG, XPT, YPT DOUBLE PRECISION SDGHT, SDGHT2 REAL VRSION, GHT, SMGHT, BGGHT REAL SLA, SLO, GHT2 INTEGER NCONV, IPAGE, ITYPE INTEGER IDLA, IMLA, IDLO, IMLO CHARACTER*80 NAME CHARACTER*15 RESP LOGICAL PAGE, NOGO, SCREEN, NOPT, EOF REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) * set defaults for those variables not used by every format type ******************************************************************** * BEGIN THE COMPUTATION LOOP FOR EACH INTERPOLATION * DO UNTIL END OF FILE OR NO MORE COMPUTATIONS REQUESTED ******************************************************************** NCONV = 0 160 CONTINUE PAGE = .FALSE. ******************************************** * GET THE NAME AND LOCATION OF ANOTHER POINT ******************************************** CALL GETPT (NCONV, ITYPE, NAME, IDLA, IMLA, SLA, IDLO, + IMLO, SLO, XPT, YPT, EOF, NOPT) IF (NOPT) GOTO 155 IF (EOF) GOTO 9999 ********************** * DO THE INTERPOLATION ********************** NOGO = .FALSE. CALL INTERP (NOGO, RESP, XPT, YPT, GHT, ITYPE) **************************************************** * CHECK TO SEE IF THIS POINT IS VALID * IF NOGO IS TRUE THEN GET ANOTHER POINT AND DON'T * DO THE COMPUTATION - POINT IS OUT OF BOUNDS * IF NOGO IS NOT TRUE THEN PROCEED - ESTIMATE MADE **************************************************** IF (NOGO) GOTO 155 NCONV = NCONV + 1 ******************************************** * CHANGE THE LONGITUDE BACK TO POSITIVE WEST * AND CHANGE BACK TO D.M.S FORMAT ******************************************** XPT = -XPT ************************** * DO THE LITTLE STATISTICS ************************** * First, the basics * meters.... GHT2 = ( DBLE(GHT) )**2 SDGHT2 = SDGHT2 + GHT2 SDGHT = SDGHT + DBLE(GHT) * then the ranges XSMALL = DMIN1( XSMALL, XPT) XBIG = DMAX1( XBIG , XPT) YSMALL = DMIN1( YSMALL, YPT) YBIG = DMAX1( YBIG , YPT) SMGHT = MIN( SMGHT, GHT) BGGHT = MAX( BGGHT, GHT) ********************************** ** WRITE TO OUTPUT FILE AND SCREEN ********************************** CALL WRTPT (ITYPE, NCONV, VRSION, NAME, + IDLA, IMLA, SLA, IDLO, IMLO, SLO, + GHT, RESP, IPAGE, PAGE, SCREEN) ********************** * START THE LOOP AGAIN ********************** 155 GOTO 160 9999 RETURN END SUBROUTINE NBLANK (A, IFLAG, NBLK) *** Return position of last non-blank in string (IFLAG = 2) *** or position of first non-blank in string (IFLAG = 1) * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER IFLAG, NBLK, LENG, IBLK CHARACTER*(*) A LENG = LEN(A) IF (IFLAG .EQ. 2) THEN DO 1 IBLK = LENG, 1, -1 IF ( A(IBLK:IBLK) .NE. ' ' ) THEN NBLK = IBLK RETURN ENDIF 1 CONTINUE ELSEIF (IFLAG .EQ. 1) THEN DO 2 IBLK = 1, LENG, +1 IF ( A(IBLK:IBLK) .NE. ' ' ) THEN NBLK = IBLK RETURN ENDIF 2 CONTINUE ENDIF * String contains all blanks NBLK = 0 RETURN END SUBROUTINE NGRIDS (NODATA) * This subroutine opens the GEOID grids. * A total of one file is necessary for each area. * If a file named AREA.PAR exists it will be read for the names and * locations of the gridded data. The format of the data in * the AREA.PAR file is given in the GRIDS subroutine. * If the AREA.PAR file does not exist, or there is still room in the * arrays in the GDINFO common then the default area names used. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =10) CHARACTER*1 ANS LOGICAL NODATA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) REAL DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NR, NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NR(MXAREA), NC(MXAREA), NAREA * Initialize NODATA = .FALSE. NAREA = 0 WRITE (LUOUT,100) 100 FORMAT (' GEOID is now opening the files containing the', + ' gridded data.', /, + ' The areas listed below may be used for GEOID', + ' ESTIMATES.') * Try to open the 'AREA.PAR' file in the subroutine GRIDS CALL GRIDS * If NAREA>=MXAREA, then skip the section that opens the default files. * If NAREA