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 (NAD83) * * FOR THE CONTERMINOUS UNITED STATES. * * INPUT DATA REFERRENCED TO THE NORTH AMERICAN DATUM OF * * 1927 (NAD27) SHOULD BE CONVERTED TO NAD 83 VIA NADCON * * PRIOR TO USING THIS PROGRAM, FOR USE IN THE * * CONTERMINOUS UNITED STATES. FOR MEXICO AND THE * * CARIBBEAN, NAD27 COORDINATES ARE ACCEPTABLE. * * * * VERSION 3.10: * * NOTE THAT THE MODELS G96SSS, MEXICO97, AND CARIB97 ARE * * ALL GEOID HEIGHTS ABOVE THE GRS-80 ELLIPSOID, CENTERED * * AT THE ORIGIN OF THE ITRF94(1996.0) REFERENCE FRAME. * * YOU CAN *STILL* INPUT LAT/LON IN NAD83 AND EXPECT TO * * GET THE SAME VALUE AS IF YOU HAD USED GRS-80/ITRF94 * * LAT/LON VALUES, HOWEVER IT *MUST* BE POINTED OUT THAT * * THE GEOID HEIGHTS OF THOSE MODELS ARE *NOT* GEOID * * HEIGHTS ABOVE NAD83! (THE REASON THIS WORKS IS BECAUSE * * THE HORIZONTAL SHIFT FROM NAD83 TO ITRF94 IS ABOUT 1 OR * * 2 METERS, WHICH WILL HAVE *NO* IMPACT ON THE GEOID * * HEIGHT VALUE; HOWEVER THE VERTICAL SHIFT FROM ITRF94 * * TO NAD83 IS ALSO 1-2 METERS, SO IT IS IMPORTANT TO * * REMEMBER WHICH ELLIPSOID THE GEOIDS REFER TO. * * * * * * THE ACTUAL COMPUTATION IS PERFORMED AS AN INTERPOLATION * * FROM A REGULARLY-SPACED GRID OF POINTS. * * 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 CODE: 3.00 (20-oct-96) * * added support for inexact header in GRIDS() * * increaded MXAREA to 15 * * VERSION CODE: 3.10 ( 6-may-97) * * * * VERSION DATE: may 6, 1997 * * july 16, 1997 -- still retained version as 3.10 * * (just allowed NAD 27 for MEX/CAR stuff) * * * * DRU SMITH, Ph.D. * * DENNIS MILBERT, Ph.D. * * SYSTEMS DEVELOPMENT DIVISION * * NATIONAL GEODETIC SURVEY, NOS, NOAA * * SILVER SPRING, MD 20910 * *********************************************************************** *********************************************************************** * * * 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 = 3.10E0, MXAREA =15) 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) logical lsss,lmex,lcar common/scienc/lsss,lmex,lcar ********************** * 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 *********************** lsss=.false. lmex=.false. lcar=.false. 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 =15) 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.', /) WRITE (LUOUT,1111) 1111 FORMAT (' *************** Which latitude/longitude to use? ', + ' ***************',/, + ' You may enter latitudes and longitudes in either the', + ' NAD 83 or the',/, + ' ITRF94/GRS-80 system, with no change in the interpolated', + ' geoid value.',/, + ' Do not use NAD 27 in the United States.',/, + ' Transform to NAD 83 first,', + ' using NADCON software.',/) ENDIF WRITE (LUOUT,*) ' What is its geodetic latitude?' c WRITE (LUOUT,*) ' What is its latitude?' c WRITE (LUOUT,*) ' (Use NAD 83 for GEOID96 and' c WRITE (LUOUT,*) ' geocentric GRS-80 for other models)' c 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 geodetic 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 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, ' ', /, + 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, ' ', 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 =15) 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 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 =15) 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 3.00 - NGS', + ' GEOID96 Model', /, + ' 45 45 45.55555 111 11 11.11111 -11.237one', /, + ' 25 55.5555555 76 56.6666666 -38.758two', /, + ' 34.444444444 98.888888888 -26.850three', /) 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, ', + 'Silver Spring, MD 20910.', /) 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,', + ' Silver Spring, MD 20910.', /) 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,', + ' Silver Spring, MD 20910.', /) 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 -31814E ', + '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 =15) 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 =15) 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 *** patch for inexact headers (due to 2' spacing) 8-oct-96 dgm idx1=idnint(DX1*3600.d0) DX(NAREA) = dble(idx1)/3600.d0 idy1=idnint(DY1*3600.d0) DY(NAREA) = dble(idy1)/3600.d0 ***** DX(NAREA) = DX1 **** old code ***** DY(NAREA) = DY1 **** old code 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 =15) 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 GEOID96 Model',//, + 10X, ' (can also interpolate G96SSS,', / + 10X, ' MEXICO97, and CARIB97 models)',/) WRITE (LUOUT,930) VRSION 930 FORMAT (/,10X, ' (Version', F5.2, ')', /, + 10X, ' July 16, 1997', /, + 10X, ' Dru Smith, Ph.D.',/, + 10X, ' Dennis Milbert, Ph.D.',/, + 10X, ' Systems Development Division', /) 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 =15) 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 10+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 =15) 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. write(luout,*)' ' c pause ' Press return to proceed...' 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 =15) 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 =15) 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 =15) 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 *** g96sss flag *** mexico97,carib97 flags logical lsss,lmex,lcar common/scienc/lsss,lmex,lcar * 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 *** extracted default file name support 23-sep-96 dgm *** added MEXICO97, CARIB97 8-apr-97 das if(lsss) write(luout,*) ' ***************************************' if(lsss) write(luout,*) ' G96SSS models encountered *************' if(lsss) write(luout,*) ' ***************************************' if(lsss) write(luout,*) ' ' if(lmex) write(luout,*) ' ***************************************' if(lmex) write(luout,*) ' MEXICO97 models encountered ***********' if(lmex) write(luout,*) ' ***************************************' if(lmex) write(luout,*) ' ' if(lcar) write(luout,*) ' ***************************************' if(lcar) write(luout,*) ' CARIB97 models encountered ************' if(lcar) write(luout,*) ' ***************************************' if(lcar) write(luout,*) ' ' WRITE (LUOUT,975) 975 FORMAT (/, ' (Hit RETURN to continue.)') READ (LUIN,'(A1)') ANS WRITE (LUOUT,2) * 2 FORMAT ('1') 2 FORMAT (' ') IF (NAREA .EQ. 0) THEN NODATA = .TRUE. WRITE (LUOUT, 970) 970 FORMAT (/, ' ******* ERROR *********', /, + ' No grid files were opened -- program ending!') ENDIF RETURN END SUBROUTINE OPENFL (AFILE, ITEMP, NOGO, DX1, DY1, XMAX1, + XMIN1, YMAX1, YMIN1, NR1, NC1, CARD, GFLAG) *** Given base name of gridded data files, open the data file * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =15) CHARACTER*23 B23 PARAMETER (B23 = ' ') CHARACTER*69 B69 PARAMETER (B69 = B23//B23//B23) REAL DX1, DY1, XMAX1, XMIN1, YMAX1, YMIN1 REAL X01, Y01, ANGLE1 INTEGER LRECL, IFILE, IOS, NZ1 INTEGER IFLAG1, IFLAG2 INTEGER ITEMP, NR1, NC1 INTEGER N1, N2, N3, N4 CHARACTER*80 CARD CHARACTER*69 AGEO CHARACTER*65 AFILE CHARACTER*56 RIDENT CHARACTER*8 PGM LOGICAL NOGO, OFLAG, EFLAG1, GFLAG INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) logical lsss,lmex,lcar common/scienc/lsss,lmex,lcar DATA OFLAG /.FALSE./, EFLAG1 /.FALSE./ DATA IFLAG1 /1/, IFLAG2 /2/ * Initialize NOGO = .FALSE. * Form complete names of grid files CALL NBLANK (AFILE, IFLAG2, N2) IF (N2 .EQ. 0) STOP 'Logical Coding Error in OPENF' AGEO = B69 AGEO(1:N2) = AFILE AGEO(N2+1:N2+4) = '.geo' ******************************************************* * DIRECT ACCESS GRID FILES * Each file is opened once to get the grid variables. * The file is then closed and reopened to ensure that * the record length is correct ******************************************************* LRECL = 256 IFILE = 10 + ITEMP LUAREA(ITEMP) = IFILE INQUIRE (FILE=AGEO, EXIST=EFLAG1, OPENED=OFLAG) IF (.NOT. EFLAG1) GOTO 910 IF (OFLAG) GOTO 980 OPEN (IFILE,FILE=AGEO,FORM='UNFORMATTED',STATUS='OLD', + ACCESS='DIRECT',RECL=LRECL,ERR=910,IOSTAT=IOS) READ (IFILE,REC=1) RIDENT, PGM, NC1, NR1, NZ1, X01, DX1, + Y01, DY1, ANGLE1 *** set flag if a g96sss model was encountered, dgm *** similarly for mexico97 or carib97, das if(ageo(1:3).eq.'sss'.or.ageo(1:3).eq.'SSS') lsss=.true. if(ageo(1:3).eq.'mex'.or.ageo(1:3).eq.'MEX') lmex=.true. if(ageo(1:3).eq.'car'.or.ageo(1:3).eq.'CAR') lcar=.true. CLOSE (IFILE) LRECL = 4*(NC1+1) OPEN (IFILE,FILE=AGEO,FORM='UNFORMATTED',STATUS='OLD', + ACCESS='DIRECT',RECL=LRECL,ERR=910,IOSTAT=IOS) * Calculate values used in this program XMIN1 = X01 YMIN1 = Y01 XMAX1 = X01 + (NC1-1)*DX1 YMAX1 = Y01 + (NR1-1)*DY1 DX1 = ABS(DX1) DY1 = ABS(DY1) ***************************************** * REPORT SOMETHING ABOUT THE GRIDDED DATA ***************************************** * WRITE (LUOUT,4050) RIDENT, PGM, NC1, NR *4050 FORMAT (1X, A56, /, 1X, A8, /, I5, I5) * WRITE (LUOUT,*) 'DX,DY,NR,NC', DX1, DY1, NR1, NC1 * WRITE (LUOUT,4055) -XMAX1, -XMIN1, YMIN1, YMAX1 *4055 FORMAT (' MIN Longitude = ', F10.4, ' MAX Longitude = ', F10.4, /, * + ' MIN Latitude = ', F10.4, ' MAX Latitude = ', F10.4, /) ***************************************** 9999 RETURN **************************** * WARNING AND ERROR MESSAGES **************************** * Grid file does not exist 910 CONTINUE NOGO = .TRUE. IF (GFLAG) THEN CALL NBLANK (AGEO, IFLAG1, N1) CALL NBLANK (CARD, IFLAG1, N3) CALL NBLANK (AGEO, IFLAG2, N2) CALL NBLANK (CARD, IFLAG2, N4) * .geo does not exist WRITE (LUOUT, 915) AGEO(N1:N2), CARD(N3:N4) 915 FORMAT (' File "',a,'" from record "',a,'" does not exist!') ENDIF CLOSE (IFILE) GOTO 9999 * Grid file already open 980 CONTINUE NOGO = .TRUE. CALL NBLANK (AGEO, IFLAG1, N1) CALL NBLANK (CARD, IFLAG1, N3) CALL NBLANK (AGEO, IFLAG2, N2) CALL NBLANK (CARD, IFLAG2, N4) WRITE (LUOUT, 985) AGEO(N1:N2) 985 FORMAT (/, 5X, ' *********** ERROR ***********', /, + 5X, ' The grid file', /, + 6X, '''', A, '''', /, + 5X, ' has already been opened! This file', + ' will not be reopened.', /, + 5X, ' *****************************', /) GOTO 9999 END SUBROUTINE PRINT (LU, NCONV, NAME, VRSION, IDLA, IMLA, SLA, + IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE) * This subroutine prints out the actual results using * a pretty format - not the same as the input file format (if there * is one). This subroutine is used by type-1 format input and * interactive input. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) REAL VRSION REAL SLA, SLO, GHT INTEGER LU, NCONV, IPAGE INTEGER IDLA, IMLA, IDLO, IMLO CHARACTER*80 NAME CHARACTER*15 RESP LOGICAL PAGE IF (NCONV .EQ. 1) THEN ******************** * FIRST PAGE HEADING ******************** WRITE (LU,10) IPAGE WRITE (LU,5) WRITE (LU,26) WRITE (LU,7) VRSION WRITE (LU,8) ENDIF IF (PAGE) THEN IF (IPAGE .GT. 1) THEN WRITE (LU,2) WRITE (LU,10) IPAGE ENDIF ENDIF 10 FORMAT (70(' '), 'Page ', I4, /) 5 FORMAT (20X, ' Geoid Height Computation') 26 FORMAT (15X, ' GEOID96/G96SSS/MEXICO97/CARIB97 Interpolation' ) 7 FORMAT (20X, ' GEOID Program ', /, 20X, + ' Version ', F4.2 ) 8 FORMAT (20X, ' ', /, 1X, 79('=') ) * 2 FORMAT ('1') 2 FORMAT (' ') WRITE (LU,921) NCONV, RESP IF (NAME .NE. ' ') WRITE (LU,922) NAME WRITE (LU,900) WRITE (LU,927) WRITE (LU,923) IDLA, IMLA, SLA, IDLO, IMLO, SLO, GHT 921 FORMAT (//, 27X, 'Computation #: ', I4, 8X, 'Region: ', A15, /) 922 FORMAT (2X, 'Station Name: ', A80, /) 900 FORMAT (22X, 'Latitude', 17X, 'Longitude', 7X, 'Geoid Height') 927 FORMAT (65X, '(meters)', /) 923 FORMAT (19X, I2, 1X, I2.2, F9.5, 10X, I3, 1X, I2.2, F9.5, + 5X, F9.3, /) RETURN END SUBROUTINE PRINT2 (LU, NCONV, VRSION, GHT) * This subroutine prints out the actual results using * a free format - the same as the input file format. This is used * for type 2 format. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) REAL VRSION, GHT INTEGER LU, NCONV CHARACTER*80 CARD COMMON /CURNT/ CARD * Write header record to identify source of geoid height and value IF (NCONV .EQ. 1) THEN WRITE (LU, 10) VRSION 10 FORMAT ('GEOID Version', F5.2, + ' - NGS GEOID96/G96SSS/MEXICO97/CARIB97 Model') ENDIF * In this format, the variable CARD contains the image of the input * card. This is written to the output file instead of using the * latitude, longitude, and name variables. The geoid height * overwrites whatever is in columns 33-40 WRITE (LU,100) CARD(1:32), GHT, CARD(41:80) 100 FORMAT (A32, F8.3, A40) RETURN END SUBROUTINE PRINT3 (ITYPE, LU, GHT, RESP) * This subroutine prints out the actual results into * a *84* or *86* record. This is used for type 3 format. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) CHARACTER*15 B15 PARAMETER (B15 = ' ') REAL GHT INTEGER LU, IGHT, ITYPE CHARACTER*125 CRD125 CHARACTER*80 CRD80,card2 CHARACTER*15 RESP CHARACTER*10 ASTD2 CHARACTER*4 ASTD CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER MXAREA PARAMETER (MXAREA =15) INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) *** g96sss flag logical lsss,lmex,lcar common/scienc/lsss,lmex,lcar *** Standard deviation is assumed to be 0.5 meters DATA ASTD /' 5'/, ASTD2 /' 0.5000'/ * In this format, the variable CARD contains the image of the input * card. This is written to the output file instead of using the * latitude, longitude, and name variables. The geoid height * overwrites whatever is in columns 33-40 IF (ITYPE .EQ. 3) THEN * Old Blue Book CRD80(1:7) = CARD(1:7) CRD80(8:9) = '84' CRD80(10:14) = CARD(10:14) CRD80(15:20) = 'NGS ' CRD80(21:71) = 'PROGRAM GEOID '//RESP//B15//' ' IGHT = NINT(GHT*10.) WRITE (CRD80(72:76),'(I5)') IGHT CRD80(77:80) = ASTD WRITE (LU,'(A80)') CRD80 ELSEIF (ITYPE .EQ. 4) THEN * New Blue Book *84* CRD80(1:7) = CARD(1:7) CRD80(8:9) = '84' CRD80(10:14) = CARD(10:14) CRD80(15:20) = 'NGS ' CRD80(21:69) = 'PROGRAM GEOID '//RESP//B15 IGHT = NINT(GHT*100.) WRITE (CRD80(70:75),'(I6)') IGHT CRD80(76:80) = ASTD//'0' WRITE (LU,'(A80)') CRD80 ******************************************************************* ELSEIF (ITYPE .EQ. 5) THEN * New Blue Book *86* *** look ahead for an exisiting 86 record *** kludgy because no clean way in standard fortran to check if at eof read(nin,'(a80)',end=66565) card2 if(card2(8:9).eq.'86'.and.card2(11:14).eq.card(11:14)) then crd80=card2 go to 66566 else crd80=' ' backspace(nin) go to 66566 endif stop 34287 *************** eof error trap *************************** 66565 crd80=' ' 66566 CRD80(1:7) = CARD(1:7) CRD80(8:9) = '86' CRD80(10:14) = CARD(10:14) IGHT = NINT(GHT*1000.0) WRITE (CRD80(36:42),'(I7)') IGHT *** code for geoid96 CRD80(43:43) = 'E' *** code for g96sss if(lsss) CRD80(43:43) = 'F' *** code for mexico97 if(lmex) CRD80(43:43) = 'J' *** code for carib97 if(lcar) CRD80(43:43) = 'H' WRITE (LU,'(A80)') CRD80 ******************************************************************* ELSEIF (ITYPE .EQ. 6) THEN * Extended Blue Book CRD125(1:7) = CARD(1:7) CRD125(8:9) = '84' CRD125(10:14) = CARD(10:14) CRD125(15:20) = 'NGS ' CRD125(21:71) = 'PROGRAM GEOID '//RESP//B15//' ' IGHT = NINT(GHT*10.) WRITE (CRD125(72:76),'(I5)') IGHT CRD125(77:80) = ASTD CRD125(81:101) = B15//' ' WRITE (CRD125(102:115),'(F14.4)') GHT CRD125(116:125) = ASTD2 WRITE (LU,'(A125)') CRD125 ENDIF RETURN END SUBROUTINE REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT, + IPAGE, PAGE) * This subroutine prints out the statistics * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION AGHT, VGHT, SDGHT REAL SMGHT, BGGHT INTEGER LU, IPAGE LOGICAL PAGE ************ * NEW PAGE ************ WRITE (LU,2) * 2 FORMAT ('1') 2 FORMAT (' ') IF (PAGE) THEN ****************** * NUMBER THIS PAGE ****************** WRITE (LU,15) IPAGE 15 FORMAT (70(' '), 'Page ', I4, /) ENDIF WRITE (LU,90) 90 FORMAT (28X, 'Geoid Height Computation') WRITE (LU,100) 100 FORMAT (//, 30X, 'Statistics for Region', /, 1X, 79('='), /) WRITE (LU,910) 910 FORMAT (/,' Geoid Height ', + 8X, 'MIN', 6X, 'MAX' ) WRITE (LU,110) SMGHT, BGGHT 110 FORMAT (' Range of Height (meters) ', 2F9.3) WRITE (LU,130) AGHT 130 FORMAT (/, 10X, 'Mean Geoid Height ', F9.3) WRITE (LU,131) VGHT 131 FORMAT ( 10X, 'Variance of Geoid Height ', F9.3) WRITE (LU,132) SDGHT 132 FORMAT ( 10X, 'Std. Dev. of Geoid Height', F9.3) RETURN END SUBROUTINE TYPE1 (NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO, + XPT, YPT, EOF, NOPT) * Read a record from a file of type 1. In this type there is a station * name (or blanks) in columns 1-40, and free-format latitude and * longitude values in columns 41-80. By free format we mean that the * numbers making up the degrees, minutes and seconds of latitude, * degrees, minutes, seconds of longitude must appear in that order in * columns 41 through 80 but are not restricted to any specific columns. * The latitude and longitude may be either (1) integer degrees, integer * minutes, decimal seconds, or (2) integer degrees, decimal minutes, or * (3) decimal degrees. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =15) CHARACTER*1 DOT, BLK PARAMETER ( DOT = '.', BLK = ' ' ) CHARACTER*80 B80 CHARACTER*20 B20 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) DOUBLE PRECISION XPT, YPT, RDLA, RDLO, DCARD REAL SLA, SLO, RCARD, RMLA, RMLO INTEGER IFLAG1, IFLAG2, N1, N2 INTEGER IDOT, IBLK, LENG, IERR INTEGER IDLA, IMLA, IDLO, IMLO CHARACTER*80 NAME CHARACTER*40 DUMLA, DUMLO LOGICAL EOF, NOPT CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) DATA IFLAG1 /1/, IFLAG2 /2/ *********************************** * FOR INPUT FILE OF ITYPE = 1 *********************************** 1 READ (NIN,'(A80)',END=9999) CARD READ (CARD(1:40), '(A40)') NAME * Check for blank line IF (CARD .EQ. B80) GOTO 1 * Find position of the first decimal point (to indicate the last * number in the latitude) IDOT = INDEX(CARD(41:80), DOT) * Error - no decimal point IF (IDOT .EQ. 0) GOTO 9980 * find position of the first blank after the first decimal point (to * indicate the blank after the last number in the latitude) IDOT = IDOT + 40 IBLK = INDEX(CARD(IDOT+1:80), BLK) IBLK = IBLK + IDOT DUMLA = CARD(41:IBLK) LENG = IBLK - 41 RDLA = DCARD( DUMLA, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 IF (LENG .GT. 0) THEN RMLA = RCARD( DUMLA, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 IF (LENG .GT. 0) THEN SLA = RCARD( DUMLA, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 ELSE SLA = 0.E0 ENDIF ELSE RMLA = 0.E0 SLA = 0.E0 ENDIF * 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 *********** DUMLO = CARD(IBLK+1:80) CALL NBLANK (DUMLO, IFLAG2, N2) LENG = N2 RDLO = DCARD( DUMLO, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 IF (LENG .GT. 0) THEN RMLO = RCARD( DUMLO, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 IF (LENG .GT. 0) THEN SLO = RCARD( DUMLO, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 ELSE SLO = 0.E0 ENDIF ELSE RMLO = 0.E0 SLO = 0.E0 ENDIF * 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 9940 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9945) CARD(N1:N2) 9945 FORMAT (' ERROR - in the following record:', /, + 9X, '''', A, '''', /, + ' Latitude and Longitudes must be positive!', /, + ' Longitude is positive west.', /) NOPT = .TRUE. GOTO 9000 9950 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9955) CARD(N1:N2) 9955 FORMAT (' ERROR - Illogical values for latitude', + ' in the following record:', /, + 9X, '''', 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 (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9965) CARD(N1:N2) 9965 FORMAT (' ERROR - Illogical values for longitude', + ' in the following record:', /, + 9X, '''', A, '''', /, + ' Longitude must be between 0 and 360 degrees.',/, + ' Minutes and seconds must be between 0', + ' and 60.', /) NOPT = .TRUE. GOTO 9000 9980 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9985) CARD(N1:N2) 9985 FORMAT (' ERROR - The following record does not have a decimal', + ' point in the latitude.', /, + 9X, '''', A, '''', /, + ' In the free format a decimal point is used', + ' to determine what is', /, + ' the last number in the latitude. Please', + ' correct this record', /, + ' and check all of the data in this file to', + ' ensure that it follows', /, + ' the correct format.', /) NOPT = .TRUE. GOTO 9000 9999 CONTINUE EOF = .TRUE. GOTO 9000 END SUBROUTINE TYPE2 (IDLA, IMLA, SLA, IDLO, IMLO, SLO, + XPT, YPT, EOF, NOPT) * Read a record from a file of type 2. In this type there is free-format * latitude and longitude values in columns 1-32, and a station name * (or blanks) in columns 41-80. Columns 33-40 will contain the geoid * height in the output record and are ignored in the input record. * By free format we mean that the numbers making up the degrees, * minutes and seconds of latitude, degrees, minutes, seconds of * longitude must appear in that order in columns 1 through 32 but are * not restricted to any specific columns. The latitude and longitude * may be either (1) integer degrees, integer minutes, decimal seconds, * or (2) integer degrees, decimal minutes, or (3) decimal degrees. * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =15) CHARACTER*1 DOT, BLK PARAMETER ( DOT = '.', BLK = ' ' ) CHARACTER*20 B20 CHARACTER*80 B80 PARAMETER (B20 = ' ', B80 = B20//B20//B20//B20) DOUBLE PRECISION XPT, YPT, RDLA, RDLO, DCARD REAL SLA, SLO, RCARD, RMLA, RMLO INTEGER IFLAG1, IFLAG2, N1, N2 INTEGER IDOT, IBLK, LENG, IERR INTEGER IDLA, IMLA, IDLO, IMLO CHARACTER*32 DUMLA, DUMLO LOGICAL EOF, NOPT CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) DATA IFLAG1 /1/, IFLAG2 /2/ *********************************** * FOR INPUT FILE OF ITYPE = 2 *********************************** 1 READ (NIN,'(A80)',END=9999) CARD * Check for blank line IF (CARD .EQ. B80) GOTO 1 * Find position of the first decimal point (to indicate the last * number in the latitude) IDOT = INDEX(CARD(1:32), DOT) * Error - no decimal point IF (IDOT .EQ. 0) GOTO 9980 * find position of the first blank after the first decimal point (to * indicate the blank after the last number in the latitude) IBLK = INDEX(CARD(IDOT+1:32), BLK) IBLK = IBLK + IDOT DUMLA = CARD(1:IBLK) LENG = IBLK - 1 RDLA = DCARD( DUMLA, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 IF (LENG .GT. 0) THEN RMLA = RCARD( DUMLA, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 IF (LENG .GT. 0) THEN SLA = RCARD( DUMLA, LENG, IERR ) IF (IERR .NE. 0) GOTO 9950 ELSE SLA = 0.E0 ENDIF ELSE RMLA = 0.E0 SLA = 0.E0 ENDIF * 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 *********** DUMLO = CARD(IBLK+1:32) CALL NBLANK (DUMLO, IFLAG2, N2) LENG = N2 RDLO = DCARD( DUMLO, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 IF (LENG .GT. 0) THEN RMLO = RCARD( DUMLO, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 IF (LENG .GT. 0) THEN SLO = RCARD( DUMLO, LENG, IERR ) IF (IERR .NE. 0) GOTO 9960 ELSE SLO = 0.E0 ENDIF ELSE RMLO = 0.E0 SLO = 0.E0 ENDIF * 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 9940 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9945) CARD(N1:N2) 9945 FORMAT (' ERROR - in the following record:', /, + 9X, '''', A, '''', /, + ' Latitude and Longitudes must be positive!', /, + ' Longitude is positive west.', /) NOPT = .TRUE. GOTO 9000 9950 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9955) CARD(N1:N2) 9955 FORMAT (' ERROR - Illogical values for latitude', + ' in the following record:', /, + 9X, '''', 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 (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9965) CARD(N1:N2) 9965 FORMAT (' ERROR - Illogical values for longitude', + ' in the following record:', /, + 9X, '''', A, '''', /, + ' Longitude must be between 0 and 360 degrees.',/, + ' Minutes and seconds must be between 0', + ' and 60.', /) NOPT = .TRUE. GOTO 9000 9980 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9985) CARD(N1:N2) 9985 FORMAT (' ERROR - The following record does not have a decimal', + ' point in the latitude.', /, + 9X, '''', A, '''', /, + ' In the free format a decimal point is used', + ' to determine what is', /, + ' the last number in the latitude. Please', + ' correct this record', /, + ' and check all of the data in this file to', + ' ensure that it follows', /, + ' the correct format.', /) NOPT = .TRUE. GOTO 9000 9999 CONTINUE EOF = .TRUE. GOTO 9000 END SUBROUTINE TYPE3 (ITYPE, NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO, + XPT, YPT, EOF, NOPT) * Read a record from a file of type 3 ('old' Blue Book) or * of type 4 or 5 ('new' Blue Book) or type 6 (extended Blue Book). * These formats are defined in 'Input Formats and Specifications of * the' National Geodetic Survey Data Base', Volume 1. Horizontal Control * Data. For type 3 the publishing date was 1980. For type 4 the * publishing date was January, 1989. *** Type 5 uses the new *86* record * Type 6 is for internal NGS use. * This publication is available from the National Geodetic Survey for a * fee by calling (301) 443-8631 * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =15) DOUBLE PRECISION XPT, YPT REAL SLA, SLO INTEGER ITYPE INTEGER IDLA, IMLA, IDLO, IMLO INTEGER IFLAG1, IFLAG2, N1, N2 CHARACTER*80 NAME CHARACTER*1 DIRLAT, DIRLON LOGICAL EOF, NOPT CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA) DATA IFLAG1 /1/, IFLAG2 /2/ *********************************** * FOR INPUT FILE OF ITYPE = 3 *********************************** READ (NIN,6000,END=9999) CARD 6000 FORMAT (A80) *** strip any preexisting 84 and 86 (retained in extended format) IF((CARD(8:9).NE.'84'.AND.CARD(8:9).NE.'86').AND.ITYPE.NE.6)THEN * For ITYPE = 3 and ITYPE 4 copy input non-*80*, non-*84* record * to output WRITE (NOUT,6000) CARD ENDIF IF (CARD(8:9) .EQ. '80') THEN * The station names and locations are in the *80* records IF (ITYPE .EQ. 3) THEN READ (CARD,6003) NAME, IDLA, IMLA, SLA, DIRLAT, + IDLO, IMLO, SLO, DIRLON 6003 FORMAT (BZ, 10X, 3X, 1X, A30, I2, I2, F7.5, A1, + I3, I2, F7.5, A1) ELSEIF (ITYPE .EQ. 4 ) THEN READ (CARD,6004) NAME, IDLA, IMLA, SLA, DIRLAT, + IDLO, IMLO, SLO, DIRLON 6004 FORMAT (BZ, 10X, 4X, A30, I2, I2, F7.5, A1, + I3, I2, F7.5, A1) ELSEIF (itype.eq.5) THEN READ (CARD,6004) NAME, IDLA, IMLA, SLA, DIRLAT, + IDLO, IMLO, SLO, DIRLON ELSEIF (ITYPE .EQ. 6) THEN READ (CARD,6005) NAME, IDLA, IMLA, SLA, DIRLAT, + IDLO, IMLO, SLO, DIRLON 6005 FORMAT (BZ, 9X, 5X, A30, I2, I2, F7.5, A1, + I3, I2, F7.5, A1) ENDIF * Check for illogical values IF (DIRLAT .NE. 'N' .AND. DIRLAT .NE. 'n') GOTO 9940 IF (IDLA .LT. 0) GOTO 9940 IF (IDLA .GT. 90) GOTO 9950 IF (IMLA .LT. 0 .OR. IMLA .GT. 60 ) GOTO 9950 IF ( SLA .LT. 0.E0 .OR. SLA .GT. 60.E0) GOTO 9950 IF (DIRLON .NE. 'W' .AND. DIRLON .NE. 'w') GOTO 9940 IF (IDLO .LT. 0) GOTO 9940 IF (IDLO .GT. 360) GOTO 9950 IF (IMLO .LT. 0 .OR. IMLO .GT. 60 ) GOTO 9950 IF ( SLO .LT. 0.E0 .OR. SLO .GT. 60.E0) GOTO 9950 ELSE NOPT = .TRUE. GOTO 9000 ENDIF YPT = DBLE(IDLA) + DBLE(IMLA)/60.D0 + DBLE(SLA)/3600.D0 XPT = DBLE(IDLO) + DBLE(IMLO)/60.D0 + DBLE(SLO)/3600.D0 9000 RETURN * Error messages 9940 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9945) CARD(N1:N2) 9945 FORMAT (' ERROR - in the following record:', /, + 9X, '''', A, '''', /, + ' Latitude and Longitudes must be positive!', /, + ' Longitude is positive west.', /) NOPT = .TRUE. GOTO 9000 9950 CONTINUE CALL NBLANK (CARD, IFLAG1, N1) CALL NBLANK (CARD, IFLAG2, N2) WRITE (LUOUT,9955) CARD(N1:N2) 9955 FORMAT (' ERROR - Illogical values for latitude or longitude', + ' in the following record:', /, + 9X, '''', A, '''', /, + ' This record will be skipped.', /) NOPT = .TRUE. GOTO 9000 9999 CONTINUE EOF = .TRUE. GOTO 9000 END SUBROUTINE WRTPT (ITYPE, NCONV, VRSION, NAME, + IDLA, IMLA, SLA, IDLO, IMLO, SLO, + GHT, RESP, IPAGE, PAGE, SCREEN) *** Write the computations to output file ***(and screen). * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA =15) REAL VRSION, GHT REAL SLA, SLO INTEGER ITYPE, NCONV, IPAGE INTEGER IDLA, IMLA, IDLO, IMLO CHARACTER*80 NAME CHARACTER*15 RESP 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) ********************* * PAGE NUMBER COUNTER ********************* * this is where you change how many on a page IF ( MOD(NCONV,5) .EQ. 0 .OR. NCONV .EQ. 1) THEN PAGE = .TRUE. IPAGE = IPAGE + 1 ENDIF *********************** ** WRITE TO OUTPUT FILE *********************** IF (ITYPE .EQ. 0) THEN ************************************** * ONLY INTERACTIVE USE - NO INPUT FILE ************************************** CALL PRINT (NOUT, NCONV, NAME, VRSION, IDLA, IMLA, SLA, + IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE) ELSEIF (ITYPE .EQ. 1) THEN ************************** * FOR FREE FORMAT TYPE 1 ************************** CALL PRINT (NOUT, NCONV, NAME, VRSION, IDLA, IMLA, SLA, + IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE) ELSEIF (ITYPE .EQ. 2) THEN ************************** * FOR FREE FORMAT TYPE 2 ************************** CALL PRINT2 (NOUT, NCONV, VRSION, GHT) ELSEIF (ITYPE .EQ. 3 .OR. ITYPE .EQ. 4 .OR. + ITYPE .EQ. 5 .or. ITYPE .EQ. 6) THEN ********************************************************** * FOR ITYPE = 3 THE HORIZONTAL BLUE BOOK OLD SPECIFICATION * FOR ITYPE = 4 THE HORIZONTAL BLUE BOOK NEW SPECIFICATION *84* * FOR ITYPE = 5 THE HORIZONTAL BLUE BOOK NEW SPECIFICATION *86* * FOR ITYPE = 6 - THE EXTENDED BLUE BOOK SPECIFICATION ********************************************************** CALL PRINT3 (ITYPE, NOUT, GHT, RESP) ENDIF ******************* * FOR SCREEN OUTPUT ******************* IF (SCREEN) THEN CALL PRINT (LUOUT, NCONV, NAME, VRSION, IDLA, IMLA, SLA, + IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE) ENDIF RETURN END CHARACTER*(*) FUNCTION CCARD (CHLINE, LENG, IERR) *** Read a character variable from a line of card image. *** LENG is the length of the card *** blanks are the delimiters of the character variable * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER LENG, IERR, I, J, ILENG CHARACTER*80 CHLINE IERR = 0 * Find first non-blank character * DO WHILE line character is blank, I is first non-blank character I = 1 10 IF ( CHLINE(I:I) .EQ. ' ' .OR. CHLINE(I:I) .EQ. ',' ) THEN I = I + 1 * Check for totally blank card (assume length of 2) IF ( I .GE. LENG) THEN CCARD = ' ' RETURN ENDIF GOTO 10 ENDIF * Find first blank character (or end of line) * DO WHILE line character is not a blank J = I + 1 20 IF ( CHLINE(J:J) .NE. ' ' .AND. CHLINE(J:J) .NE. ',' ) THEN J = J + 1 * Check for totally filed card IF ( J .GT. LENG) THEN GOTO 40 ENDIF GOTO 20 ENDIF * J is now 1 more than the position of the last non-blank character 40 J = J - 1 * ILENG is the length of the character string, it can be any length * up to the length of the line ILENG = J - I + 1 IF (ILENG .GT. LENG) THEN STOP 'CCARD' ENDIF * Read the char variable from the line, and set the return VAR to it READ (CHLINE(I:J), 55, ERR=9999) CCARD 55 FORMAT (A) * Now reset the values of LENG and CHLINE to the rest of the card CHLINE( 1 : LENG ) = CHLINE( (J+1) : LENG ) LENG = LENG - J RETURN * Read error 9999 IERR = 1 RETURN END DOUBLE PRECISION FUNCTION DCARD (CHLINE, LENG, IERR) *** Read a double precision number from a line of card image. *** LENG is the length of the card *** blanks are the delimiters of the REAL*8 variable * IMPLICIT DOUBLE PRECISION (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION VAR INTEGER LENG, IERR, I, J, ILENG CHARACTER*80 CHLINE IERR = 0 * Find first non-blank character * DO WHILE line character is blank, I is first non-blank character I = 1 10 IF ( CHLINE(I:I) .EQ. ' ' .OR. CHLINE(I:I) .EQ. ',' ) THEN I = I + 1 * Check for totally blank card IF ( I .GE. LENG) THEN DCARD = 0.0D0 LENG = 0 RETURN ENDIF GOTO 10 ENDIF * Find first blank character (or end of line) * DO WHILE line character is not a blank J = I + 1 20 IF ( CHLINE(J:J) .NE. ' ' .AND. CHLINE(J:J) .NE. ',' ) THEN J = J + 1 * Check for totally filed card IF ( J .GT. LENG) THEN GOTO 40 ENDIF GOTO 20 ENDIF * J is now 1 more than the position of the last non-blank character 40 J = J - 1 * ILENG is the length of the real string, it cannot be greater * than 15 characters ILENG = J - I + 1 IF (ILENG .GT. 20) THEN STOP 'DCARD' ENDIF * Read the real variable from the line, and set the return VAR to it READ (CHLINE(I:J), 55, ERR=9999) VAR 55 FORMAT (F20.0) DCARD = VAR * Now reset the values of LENG and CHLINE to the rest of the card CHLINE( 1 : LENG ) = CHLINE( (J+1) : LENG ) LENG = LENG - J RETURN * Read error 9999 IERR = 1 RETURN END INTEGER FUNCTION ICARD (CHLINE, LENG, IERR) *** Read an integer from a line of card image. *** LENG is the length of the card *** blanks are the delimiters of the integer * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER LENG, IERR, I, J, ILENG, IVAR CHARACTER*80 CHLINE IERR = 0 * Find first non-blank character * DO WHILE line character is bland, I is first non-blank character I = 1 10 IF ( CHLINE(I:I) .EQ. ' ' .OR. CHLINE(I:I) .EQ. ',' ) THEN I = I + 1 * Check for totally blank card IF ( I .GE. LENG) THEN ICARD = 0 LENG = 0 RETURN ENDIF GOTO 10 ENDIF * Find first blank character (or end of line) * DO WHILE line character is not a blank J = I + 1 20 IF ( CHLINE(J:J) .NE. ' ' .AND. CHLINE(J:J) .NE. ',' ) THEN J = J + 1 * Check for totally filed card IF ( J .GT. LENG) THEN GOTO 40 ENDIF GOTO 20 ENDIF * J is now 1 more than the position of the last non-blank character 40 J = J - 1 * ILENG is the length of the integer string, it cannot be greater * than 13 characters ILENG = J - I + 1 IF (ILENG .GT. 13) THEN STOP 'ICARD' ENDIF * Read the integer variable from the line, and set the return VAR to it READ (CHLINE(I:J), 55, ERR=9999) IVAR 55 FORMAT (I13) ICARD = IVAR * Now reset the values for LENG and CHLINE to the rest of the card CHLINE( 1 : LENG ) = CHLINE( (J+1) : LENG ) LENG = LENG - J RETURN * Read error 9999 IERR = 1 RETURN END REAL FUNCTION QTERP2 (X, F0, F1, F2) *** Linear Quadratic interpolation from equally spaced (real) values *** Uses Newton-Gregory forward polynomial *** X ranged from 0 through 2. (thus S = X) * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) REAL X, F0, F1, F2 REAL DF0, DF1, D2F0 *** forward differences DF0 = F1 - F0 DF1 = F2 - F1 D2F0 = DF1 - DF0 QTERP2 = F0 + X*DF0 + 0.5*X*(X-1.)*D2F0 RETURN END REAL FUNCTION RCARD (CHLINE, LENG, IERR) *** Read a real number from a line of card image. *** LENG is the length of the card *** blanks are the delimiters of the REAL*4 variable * IMPLICIT REAL (A-H, O-Z) * IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER LENG, IERR, I, J, ILENG REAL VAR CHARACTER*80 CHLINE IERR = 0 * Find first non-blank character * DO WHILE line character is blank, I is first non-blank character I = 1 10 IF ( CHLINE(I:I) .EQ. ' ' .OR. CHLINE(I:I) .EQ. ',' ) THEN I = I + 1 * Check for totally blank card IF ( I .GE. LENG) THEN RCARD = 0.0E0 LENG = 0 RETURN ENDIF GOTO 10 ENDIF * Find first blank character (or end of line) * DO WHILE line character is not a blank J = I + 1 20 IF ( CHLINE(J:J) .NE. ' ' .AND. CHLINE(J:J) .NE. ',' ) THEN J = J + 1 * Check for totally filed card IF ( J .GT. LENG) THEN GOTO 40 ENDIF GOTO 20 ENDIF * J is now 1 more than the position of the last non-blank character 40 J = J - 1 * ILENG is the length of the real string, it cannot be greater * than 15 characters ILENG = J - I + 1 IF (ILENG .GT. 15) THEN STOP 'RCARD' ENDIF * Read the real variable from the line, and set the return VAR to it READ (CHLINE(I:J), 55, ERR=9999) VAR 55 FORMAT (F15.0) RCARD = VAR * Now reset the values of LENG and CHLINE to the rest of the card CHLINE( 1 : LENG ) = CHLINE( (J+1) : LENG ) LENG = LENG - J RETURN * Read error 9999 IERR = 1 RETURN END