PROGRAM NADCON *********************************************************************** * * * PROGRAM : NADCON * * * * PURPOSE: COMPUTATION PROGRAM TO CONVERT (OR TRANSFORM) * * POSITIONAL DATA (E.G., LATITUDES AND LONGITUDES) FROM * * THE NORTH AMERICAN DATUM OF 1927 (NAD 27) TO THE * * NORTH AMERICAN DATUM OF 1983 (NAD 83). THIS PROGRAM * * CAN COMPUTE FROM FROM EITHER DATUM TO THE OTHER. * * * * THE ACTUAL COMPUTATION IS PERFORMED AS AN INTERPOLATION * * FROM A REGULARLY-SPACED GRID OF POINTS OBTAINED FROM THE* * FITTING OF A MINIMUM-CURVATURE SURFACE TO THE ACTUAL * * SHIFT DATA RESULTING FROM THE NAD 83 ADJUSTMENT. * * * * THE INTERPOLATION IS ACCOMPLISHED BY LOCALLY FITTING * * A CURVED POLYNOMIAL SURFACE TO THE FOUR DATA POINTS * * DEFINING THE SQUARE WHICH SURROUND THE (X,Y) PAIR * * WHERE THE INTERPOLATION IS TO TAKE PLACE. * * * * THE POLYNOMIAL SURFACE IS DEFINED BY: * * * * A+BX+CY+DXY=Z * * * * THE PROGRAM REQUIRES THAT THE USER SPECIFY: * * * * 1) THE NAME OF AN OUTPUT FILE * * * * 2) THE NAME OF AN INPUT FILE (IF AVAILABLE). * * * * * * * * ESTIMATES OF DATUM SHIFTS IN TERMS OF METERS ARE * * COMPUTED FROM THE SHIFT ESTIMATES USING ELLIPSOIDAL * * SCALING. * * * * THIS PROGRAM ALLOWS FOR EITHER NGS STANDARD HORIZONTAL * * DATA FORMATS AS SPECIFIED IN THE FGCC PUBLICATION, * * COMMONLY KNOWN AS THE 'HORIZONTAL BLUE BOOK' (SEE * * SUBROUTINE TYPE3), 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: 2.00 * * * * VERSION DATE: JULY 1, 1991 * * * * AUTHOR: WARREN T. DEWHURST, PH.D. * * LIEUTENANT COMMANDER, NOAA * * ALICE R. DREW * * SENIOR GEODESIST, HORIZONTAL NETWORK 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. * * * *********************************************************************** C C THE SUBROUTINES HAVE BEEN ALTERED AND REORDERED FOR C THE INTERGRATION WITH THE CORPSCON 3.0 PROGRAM..... C C LAST MODIFICATION: C 2 DECEMBER 1991 C Parameter Change 2 June 92 C Parameter Change 24 February 94 - changed C corpsd(22)->corpsd(24) C MODIFIED BY: C JEFF WALKER C U.S. ARMY TOPOGRAPHIC ENGINEERING CENTER C FORT BELVOIR, VA 22060-5546 C (703) 355-2766 C C IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION VRSION INTEGER MXAREA PARAMETER (VRSION = 2.00D0, MXAREA = 7) DOUBLE PRECISION DLAM, DLOM, DLAS, DLOS INTEGER IFILE LOGICAL NODATA, NOGO DOUBLE PRECISION DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) DOUBLE PRECISION CORPSD(29) INTEGER*2 FLAG, KEY CHARACTER*25 NAME CHARACTER*60 INFILE WRITE(*,*) 'STARTING CORPSNAD' ********************** * INITIALIZE VARIABLES ********************** CALL INITL ******************************************************* * OPEN NADCON DATA FILES (LATITUDE AND LONGITUDE GRIDS) ******************************************************* CALL NGRIDS (NODATA) IF (NODATA) GOTO 9999 C******************************** C READ DATA FROM CORPSCON PROGRAM C******************************** READ(5,10) INFILE 10 FORMAT(A60) READ(5,15) KEY 15 FORMAT(I1) WRITE(*,*) 'CORPSNAD INPUT PARAMETERS READ SUCCESSFULLY' C******************************** C OPEN BINARY INPUT/OUTPUT FILES C FROM THE CORPSCON PROGRAM C******************************** OPEN(91,FILE=INFILE,ACCESS='SEQUENTIAL',FORM='BINARY', &ERR=25,STATUS='OLD') WRITE(*,*) 'CORPSNAD INPUT FILE OPENED SUCCESSFULLY' OPEN(93,FILE='INT$OUT',ACCESS='SEQUENTIAL',FORM='BINARY', &ERR=25,STATUS='UNKNOWN') WRITE(*,*) 'CORPSNAD OUTPUT FILE OPENED SUCCESSFULLY' ********************************* * LOOP (ONCE FOR EACH CONVERSION) C UNTIL END OF BINARY FILE (91) ********************************* 20 READ (91,END=25) NAME,(CORPSD(II),II=1,29),FLAG IF(KEY.EQ.1) THEN YPT = CORPSD(9) XPT = -CORPSD(10) ELSE YPT = CORPSD(11) XPT = -CORPSD(12) ENDIF ************************ * DO THE TRANSFORMATION ************************ NOGO = .FALSE. CALL TRANSF (NOGO, XPT, YPT, XPT2, YPT2, + DLAM, DLOM, DLAS, DLOS, KEY) IF (NOGO) THEN FLAG = IOR(FLAG,4) IF(KEY.EQ.1) THEN CORPSD(11) = CORPSD(9) CORPSD(12) = CORPSD(10) ELSE CORPSD(9) = CORPSD(11) CORPSD(10) = CORPSD(12) ENDIF CORPSD(21) = 0 CORPSD(22) = 0 ELSE c DISTM = SQRT(DLAM**2 + DLOM**2) IF(KEY.EQ.1) THEN CORPSD(11) = YPT2 CORPSD(12) = -XPT2 ELSE CORPSD(9) = YPT CORPSD(10) = -XPT ENDIF CORPSD(21) = DLAM CORPSD(22) = DLOM ENDIF WRITE (93) NAME,(CORPSD(II),II=1,29),FLAG GOTO 20 ***************** * CLOSE ALL FILES ***************** 25 DO 1010 IFILE = 1, 2*NAREA CLOSE (LUAREA(IFILE), STATUS='KEEP') 1010 CONTINUE C WRITE(LUOUT,*) 'ERROR OPENING FILES' CLOSE(91) CLOSE(93) WRITE(*,*) 'CORPSNAD INPUT/OUTPUT FILES CLOSED' WRITE(*,*) 'CORPSNAD COMPLETE' 9999 STOP END SUBROUTINE COEFF (TEE1, TEE2, TEE3, TEE4, AY, BEE, CEE, DEE) ********************************************************************** ** SUBROUTINE COEFF: GENERATES COEFFICIENTS FOR SURFACE FUNCTION * ********************************************************************** IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION AY, BEE, CEE, DEE DOUBLE PRECISION TEE1, TEE2, TEE3, TEE4 AY = TEE1 BEE = TEE3 - TEE1 CEE = TEE2 - TEE1 DEE = TEE4 - TEE3 - TEE2 + TEE1 RETURN END ********************************************************************** SUBROUTINE DGRIDS ********************************************************************** * This subroutine opens the NADCON grids 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 DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA, MXDEF PARAMETER (MXAREA = 7, MXDEF = MXAREA) CHARACTER*80 B80 CHARACTER*65 B65 CHARACTER*20 B20 PARAMETER (B20 = ' ') PARAMETER (B80 = B20//B20//B20//B20) PARAMETER (B65 = B20//B20//B20//' ') DOUBLE PRECISION XMAX1, XMIN1, YMAX1, YMIN1 DOUBLE PRECISION DX1, DY1 INTEGER IDEF, ITEMP, NC1 CHARACTER*80 DUM CHARACTER*65 AFILE, DFILES(MXAREA) CHARACTER*15 DAREAS(MXDEF) LOGICAL NOGO, GFLAG CHARACTER*15 AREAS COMMON /AREAS/ AREAS(MXAREA) DOUBLE PRECISION DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) DATA DUM / B80 / * DFILES contains the default locations (pathname) of the grid files * without the .las and .los extensions. (For example 'conus' would * indicate that the conus.las and conus.los grid files are 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', + 'stlrnc', 'stgeorge', 'stpaul', 'alaska'/ DATA DAREAS /'Conus', 'Hawaii', 'P.R. and V.I.', + 'St. Lawrence I.', 'St. George I.', 'St. Paul I.' + ,'Alaska'/ 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, GFLAG, NOGO, DX1, DY1, + XMAX1, XMIN1, YMAX1, YMIN1, NC1, DUM) 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 NC(NAREA) = NC1 WRITE (LUOUT,120) NAREA,AREAS(NAREA),YMIN1,XMIN1,YMAX1,XMAX1 120 FORMAT (2X, I2, 2X, A15, 4(2X, F10.5)) ENDIF 140 CONTINUE 999 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 INTRP. * ********************************************************************** * 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 DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION XPT, YPT, XGRID, YGRID DOUBLE PRECISION XMAX, XMIN, YMAX, YMIN DOUBLE PRECISION DX, DY INTEGER IROW, JCOL LOGICAL NOGO NOGO = .FALSE. * Check to see it the point is outside the area of the gridded data IF (XPT .GE. XMAX .OR. XPT .LE. XMIN .OR. + YPT .GE. YMAX .OR. YPT .LE. YMIN ) THEN NOGO = .TRUE. c WRITE (LUOUT,*) '***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 = ( XPT - XMIN )/DX + 1.D0 YGRID = ( YPT - YMIN )/DY + 1.D0 * Find the I,J values for the SW corner of the local square IROW = IDINT(YGRID) JCOL = IDINT(XGRID) 200 RETURN END ********************************************************************** SUBROUTINE INITL ********************************************************************** *** This subroutine initializes all the variables needed in NADCON IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA = 7) CHARACTER*20 B20 CHARACTER*80 B80 PARAMETER (B20 = ' ') PARAMETER (B80 = B20//B20//B20//B20) CHARACTER*80 CARD COMMON /CURNT/ CARD INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*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 2*MXAREA LUIN = 5 LUOUT = 6 NOUT = 101 NIN = 102 NAPAR = 103 RETURN END ********************************************************************** SUBROUTINE INTRP (IAREA, IROW, NC, JCOL, XGRID, YGRID, + XPT, YPT, XPT2, YPT2, DLOS, DLAS, DLAM, DLOM) ********************************************************************** ** DETERMINE SURFACE FUNCTION FOR THIS GRID SQUARE * ** AND INTERPOLATE A VALUE, ZEE, FOR XPT, YPT * ********************************************************************** IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA, MAXCOL PARAMETER (MAXCOL = 600, MXAREA = 7) DOUBLE PRECISION XPT, YPT, XPT2, YPT2, XGRID, YGRID DOUBLE PRECISION DLOS, DLAS, DLAM, DLOM DOUBLE PRECISION TEE1, TEE2, TEE3, TEE4, ZEE INTEGER IROW, JCOL, NC, IAREA, IFILE, IDUM, J REAL BUF(MAXCOL) INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) DOUBLE PRECISION AY1, BEE1, CEE1, DEE1, AY2, BEE2, CEE2, DEE2 SAVE AY1, BEE1, CEE1, DEE1, AY2, BEE2, CEE2, DEE2 INTEGER IROWL, JCOLL, IAREAL SAVE IROWL, JCOLL, IAREAL DATA IROWL / 0 /, JCOLL / 0 /, IAREAL / 0 / ********** * LATITUDE ********** IF ( IROW .NE. IROWL .OR. JCOL .NE. JCOLL .OR. + IAREA .NE. IAREAL ) THEN * Lower boundary IFILE = LUAREA( 2*IAREA - 1 ) READ (IFILE,REC=IROW+1) IDUM, (BUF(J), J=1,NC) TEE1 = DBLE( BUF(JCOL) ) * TEE4 = DBLE( BUF(JCOL+1) ) TEE3 = DBLE( BUF(JCOL+1) ) * Upper boundary READ (IFILE,REC=IROW+2) IDUM, (BUF(J), J=1,NC) TEE2 = DBLE( BUF(JCOL) ) * TEE3 = DBLE( BUF(JCOL+1) ) TEE4 = DBLE( BUF(JCOL+1) ) CALL COEFF (TEE1, TEE2, TEE3, TEE4, AY1, BEE1, CEE1, DEE1) ENDIF CALL SURF (XGRID, YGRID, ZEE, AY1, BEE1, CEE1, DEE1, IROW, JCOL) DLAS = ZEE *********** * LONGITUDE *********** IF ( IROW .NE. IROWL .OR. JCOL .NE. JCOLL .OR. + IAREA .NE. IAREAL ) THEN * Lower boundary IFILE = LUAREA( 2*IAREA ) READ (IFILE,REC=IROW+1) IDUM, (BUF(J), J=1,NC) TEE1 = DBLE( BUF(JCOL) ) * TEE4 = DBLE( BUF(JCOL+1) ) TEE3 = DBLE( BUF(JCOL+1) ) * Upper boundary READ (IFILE,REC=IROW+2) IDUM, (BUF(J), J=1,NC) TEE2 = DBLE( BUF(JCOL) ) * TEE3 = DBLE( BUF(JCOL+1) ) TEE4 = DBLE( BUF(JCOL+1) ) CALL COEFF (TEE1, TEE2, TEE3, TEE4, AY2, BEE2, CEE2, DEE2) ENDIF CALL SURF (XGRID, YGRID, ZEE, AY2, BEE2, CEE2, DEE2, IROW, JCOL) DLOS = ZEE ************************** * COMPUTE THE NAD 83 VALUES ************************** YPT2 = YPT + DLAS/3600.D0 * Longitude is positive west in this subroutine XPT2 = XPT - DLOS/3600.D0 ********************************************************************* * USE THE NEW ELLIPSOIDAL VARIABLES TO COMPUTE THE SHIFTS IN METERS ********************************************************************* CALL METERS (YPT, XPT, YPT2, XPT2, DLAM, DLOM) * Update the last-value variables IROWL = IROW JCOLL = JCOL IAREAL = IAREA RETURN END ********************************************************************** SUBROUTINE METERS (LAT1, LONG1, LAT2, LONG2, LATMTR, LONMTR) ********************************************************************** * This subroutine computes the difference in two positions in meters. * * This method utilizes ellipsoidal rather than spherical * parameters. I believe that the original approach and code * for this came from Ed McKay. * The reference used by Ed McKay for this was: * 'A Course in Higher Geodesy' by P.W. Zakatov, Israel Program * for Scientific Translations, Jerusalem, 1962 * * Warren T. Dewhurst * 11/1/89 * Note that this subroutine is set up for +west longitude IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) * I think that these are GRS80 parameters DOUBLE PRECISION AXIS, E2, RHOSEC PARAMETER (AXIS = 6378137.0D0) PARAMETER (E2 = 0.0066943800229D0) PARAMETER (RHOSEC = 206264.806247D0) DOUBLE PRECISION W, LM, LP, AVLAT DOUBLE PRECISION LAT1S, LAT2S, LONG1S, LONG2S, LAT1, LAT2 DOUBLE PRECISION LONG1, LONG2, DLAT, DLONG DOUBLE PRECISION LATMTR, LONMTR * LAT1 = (LATSEC + 60.D0*( LATMIN + 60.D0*LATDEG) )/RHOSEC * LONG1 = (LONSEC + 60.D0*( LONMIN + 60.D0*LONDEG) )/RHOSEC * LAT2 = (LATSEC + 60.D0*( LATMIN + 60.D0*LATDEG) )/RHOSEC * LONG2 = (LONSEC + 60.D0*( LONMIN + 60.D0*LONDEG) )/RHOSEC * Change into sec.ddd and convert to +west longitude LAT1S = LAT1*60.D0*60.D0/RHOSEC LONG1S = -LONG1*60.D0*60.D0/RHOSEC LAT2S = LAT2*60.D0*60.D0/RHOSEC LONG2S = -LONG2*60.D0*60.D0/RHOSEC DLAT = ( LAT2S - LAT1S)*RHOSEC DLONG = (LONG2S - LONG1S)*RHOSEC AVLAT = (LAT1S + LAT2S)/2.0D0 W = DSQRT(1.0D0 - E2*DSIN(AVLAT)**2) LM = AXIS*(1.0D0 - E2)/(W**3*RHOSEC) LP = AXIS*DCOS(AVLAT)/(W*RHOSEC) LATMTR = LM*DLAT LONMTR = LP*DLONG RETURN END ********************************************************************** SUBROUTINE NGRIDS (NODATA) ********************************************************************** * This subroutine opens the NADCON grids which contain datum shifts. * A total of two files are necessary for each area; 1 for each latitude * and longitude shift table (gridded data set) expressed in arc seconds. * 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 DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA = 7) LOGICAL NODATA DOUBLE PRECISION DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) * Initialize NODATA = .FALSE. NAREA = 0 CALL DGRIDS 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, GFLAG, NOGO, DX, DY, + XMAX1, XMIN1, YMAX1, YMIN1, NC1, CARD) ********************************************************************** *** Given base name of gridded data files, open the two data files IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA = 7) CHARACTER*23 B23 PARAMETER (B23 = ' ') CHARACTER*69 B69 PARAMETER (B69 = B23//B23//B23) DOUBLE PRECISION XMAX1, XMIN1, YMAX1, YMIN1, DX, DY REAL DX1, DY1, DX2, DY2 REAL X01, Y01, ANGLE1, X02, Y02, ANGLE2 INTEGER IFLAG1, IFLAG2, N1, N2, N3, N4 INTEGER ITEMP, LRECL, ILA, ILO, IFILE, IOS INTEGER NC1, NR1, NZ1, NC2, NR2, NZ2 CHARACTER*80 CARD CHARACTER*69 ALAS, ALOS CHARACTER*65 AFILE CHARACTER*56 RIDENT CHARACTER*8 PGM LOGICAL GFLAG, NOGO, OFLAG, EFLAG1, EFLAG2 INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) DATA IFLAG1 /1/, IFLAG2 /2/ DATA OFLAG /.FALSE./, EFLAG1 /.FALSE./, EFLAG2 /.FALSE./ * Initialize NOGO = .FALSE. * Form complete names of grid files CALL NBLANK (AFILE, IFLAG2, N2) IF (N2 .EQ. 0) STOP 'Logical Coding Error in OPENF' ALAS = B69 ALAS(1:N2) = AFILE ALAS(N2+1:N2+4) = '.las' ALOS = B69 ALOS(1:N2) = AFILE ALOS(N2+1:N2+4) = '.los' ******************************************************* * 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 ******************************************************* * Seconds of latitude grid file LRECL = 256 ILA = 2*ITEMP - 1 IFILE = ILA + 10 LUAREA(ILA) = IFILE INQUIRE (FILE=ALAS, EXIST=EFLAG1, OPENED=OFLAG) IF (.NOT. EFLAG1) GOTO 100 IF (OFLAG) GOTO 980 OPEN (IFILE,FILE=ALAS,FORM='UNFORMATTED',STATUS='OLD', + ACCESS='DIRECT',RECL=LRECL,ERR=940,IOSTAT=IOS) READ (IFILE,REC=1) RIDENT, PGM, NC1, NR1, NZ1, X01, DX1, + Y01, DY1, ANGLE1 CLOSE (IFILE) LRECL = 4*(NC1+1) OPEN (IFILE,FILE=ALAS,FORM='UNFORMATTED',STATUS='OLD', + ACCESS='DIRECT',RECL=LRECL,ERR=940,IOSTAT=IOS) * Seconds of longitude grid file 100 LRECL = 256 ILO = 2*ITEMP IFILE = ILO + 10 LUAREA(ILO) = IFILE INQUIRE (FILE=ALOS, EXIST=EFLAG2, OPENED=OFLAG) IF (.NOT. EFLAG1) GOTO 910 IF (.NOT. EFLAG2) GOTO 920 IF (OFLAG) GOTO 980 OPEN (IFILE,FILE=ALOS,FORM='UNFORMATTED',STATUS='OLD', + ACCESS='DIRECT',RECL=LRECL,ERR=940,IOSTAT=IOS) READ (IFILE,REC=1) RIDENT, PGM, NC2, NR2, NZ2, X02, DX2, + Y02, DY2, ANGLE2 CLOSE (IFILE) LRECL = 4*(NC2+1) OPEN (IFILE,FILE=ALOS,FORM='UNFORMATTED',STATUS='OLD', + ACCESS='DIRECT',RECL=LRECL,ERR=940,IOSTAT=IOS) * Check to see if the two files have the same variables IF ( (NC2 .NE. NC1) .OR. (NR2 .NE. NR1) .OR. + (NZ2 .NE. NZ1) .OR. + (X02 .NE. X01) .OR. (DX2 .NE. DX1) .OR. + (Y02 .NE. Y01) .OR. (DY2 .NE. DY1) .OR. + (ANGLE2 .NE. ANGLE1) ) GOTO 960 * Calculate values used in this program XMIN1 = DBLE(X01) YMIN1 = DBLE(Y01) XMAX1 = DBLE(X01) + (NC1-1)*DBLE(DX1) YMAX1 = DBLE(Y01) + (NR1-1)*DBLE(DY1) DX = DBLE( ABS(DX1) ) DY = DBLE( 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 files do not exist 910 CONTINUE NOGO = .TRUE. IF (GFLAG) THEN CALL NBLANK (ALAS, IFLAG1, N1) CALL NBLANK (ALAS, IFLAG2, N2) CALL NBLANK (CARD, IFLAG1, N3) CALL NBLANK (CARD, IFLAG2, N4) IF (EFLAG2) THEN * .los exists, .las does not exist WRITE (LUOUT, 925) ALAS(N1:N2), CARD(N3:N4), ALOS(N1:N2) ELSE * neither .los nor .las exist WRITE (LUOUT, 915) ALAS(N1:N2), ALOS(N1:N2), CARD(N3:N4) 915 FORMAT (/, 5X, ' *********** WARNING ***********', /, + 5X, ' The grid files', /, + 6X, '''', A, '''', /, + 6X, '''', A, '''', /, + 5X, ' from record:', /, + 6X, '''', A, '''', /, + 5X, ' do not exist!', /, + 5X, ' *******************************', /) ENDIF ENDIF CLOSE ( LUAREA(ILA) ) CLOSE ( LUAREA(ILO) ) GOTO 9999 920 CONTINUE NOGO = .TRUE. IF (GFLAG) THEN * .las exists, .los does not exist CALL NBLANK (ALAS, IFLAG1, N1) CALL NBLANK (ALAS, IFLAG2, N2) CALL NBLANK (CARD, IFLAG1, N3) CALL NBLANK (CARD, IFLAG2, N4) WRITE (LUOUT, 925) ALOS(N1:N2), CARD(N3:N4), ALAS(N1:N2) 925 FORMAT (/, 5X, ' *********** WARNING ***********', /, + 5X, ' The grid file', /, + 6X, '''', A, '''', /, + 5X, ' from record:', /, + 6X, '''', A, '''', /, + 5X, ' does not exist! However, the grid file', /, + 6X, '''', A, '''', /, + 5X, ' does exist!', /, + 5X, ' *******************************', /) ENDIF CLOSE ( LUAREA(ILA) ) CLOSE ( LUAREA(ILO) ) GOTO 9999 * Grid file(s) already open 940 CONTINUE NOGO = .TRUE. IF (GFLAG) THEN CALL NBLANK (ALAS, IFLAG1, N1) CALL NBLANK (ALAS, IFLAG2, N2) CALL NBLANK (CARD, IFLAG1, N3) CALL NBLANK (CARD, IFLAG2, N4) WRITE (LUOUT, 945) ALAS(N1:N2), ALOS(N1:N2), CARD(N3:N4), IOS 945 FORMAT (/, 5X, ' *********** WARNING ***********', /, + 5X, ' The grid file', /, + 6X, '''', A, '''', /, + 5X, ' and/or grid file', /, + 6X, '''', A, '''', /, + 5X, ' from record:', /, + 6X, '''', A, '''', /, + 5X, ' cannot be opened!', + ' These files will be ignored.', 5X, I5, /, + 5X, ' *******************************', /) ENDIF CLOSE ( LUAREA(ILA) ) CLOSE ( LUAREA(ILO) ) GOTO 9999 * Grid files do not agree 960 CONTINUE NOGO = .TRUE. CALL NBLANK (ALAS, IFLAG1, N1) CALL NBLANK (ALAS, IFLAG2, N2) WRITE (LUOUT, 965) ALAS(N1:N2), ALOS(N1:N2) 965 FORMAT (/, 5X, ' *********** ERROR ***********', /, + 5X, ' The header information in grid files', /, + 6X, '''', A, '''', /, + 5X, ' and', /, + 6X, '''', A, '''', /, + 5X, ' do not agree! One or both of these files must', + ' be corrupted.', /, + 5X, ' These files will be ignored.', /, + 5X, ' *****************************', /) CLOSE ( LUAREA(ILA) ) CLOSE ( LUAREA(ILO) ) GOTO 9999 * Grid files already open 980 CONTINUE NOGO = .TRUE. IF (GFLAG) THEN CALL NBLANK (ALAS, IFLAG1, N1) CALL NBLANK (ALAS, IFLAG2, N2) WRITE (LUOUT, 985) ALAS(N1:N2), ALOS(N1:N2) 985 FORMAT (/, 5X, ' *********** ERROR ***********', /, + 5X, ' The grid file', /, + 6X, '''', A, '''', /, + 5X, ' and the grid file', /, + 6X, '''', A, '''', /, + 5X, ' have already been opened! These files', + ' will not be reopened.', /, + 5X, ' *****************************', /) ENDIF GOTO 9999 END ********************************************************************** SUBROUTINE TO83 (NOGO, XPT, YPT, XPT2, YPT2, + DLAM, DLOM, DLAS, DLOS) ********************************************************************** * This subroutine predicts the NAD 83 latitude and longitude values * given the NAD 27 latitude and longitude values in degree decimal * format. In addition, the program returns the shift values between * The datums in both arc secs and meters. * All of the predictions are based upon a straight-forward interpolation * of a gridded data set of datum shifts. The datum shifts are assumed * to be provided in the files opened in the NGRIDS subroutine. The * common AREAS contains the names of the valid areas while the common * GDINFO contains the grid variables. NAREA is the number of areas * which had data files opened. A total of two files are necessary for * each area: one latitude and one longitude shift table (gridded data * set) expressed in arc seconds. * For this subroutine, it is important to remember that the * input longitude is assumed to be positive east and the * output longitude will be positive east. * Author: Warren T. Dewhurst, PH. D. * National Geodetic Survey * November 1, 1989 IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA PARAMETER (MXAREA = 7) DOUBLE PRECISION XPT, YPT, XPT2, YPT2 DOUBLE PRECISION XGRID, YGRID DOUBLE PRECISION DLAM, DLOM, DLAS, DLOS DOUBLE PRECISION DX0, DY0, XMAX0, XMIN0, YMAX0, YMIN0 INTEGER IROW, JCOL, IAREA, I, NC0, ITYPE INTEGER IFLAG1, IFLAG2, N1, N2 LOGICAL NOGO, FLAG CHARACTER*15 AREAS COMMON /AREAS/ AREAS(MXAREA) DOUBLE PRECISION DX, DY, XMAX, XMIN, YMAX, YMIN INTEGER NC, NAREA COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA), + XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA), + NC(MXAREA), NAREA INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) CHARACTER*80 CARD COMMON /CURNT/ CARD SAVE FLAG DATA IFLAG1 /1/, IFLAG2 /2/, FLAG /.FALSE./ ****************************************************************** * INITIALIZE ****************************************************************** NOGO = .FALSE. * ITYPE = 0 **************************************************** * 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) NC0 = NC(IAREA) CALL FGRID (XPT, YPT, DX0, DY0, XMAX0, XMIN0, + YMAX0, YMIN0, XGRID, YGRID, IROW, JCOL, NOGO) IF (.NOT. NOGO) GOTO 200 WRITE (LUOUT,28) YPT,XPT,IAREA 28 FORMAT ('THE POINT ',F10.5,2X,F10.5,' IS NOT IN AREA ',I2) 100 CONTINUE * Not in any of the grid areas NOGO = .TRUE. GOTO 950 200 CONTINUE * Point in area number IAREA and named AREAS(IAREA) CALL INTRP (IAREA, IROW, NC0, JCOL, XGRID, YGRID, + XPT, YPT, XPT2, YPT2, DLOS, DLAS, DLAM, DLOM) 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 GOTO 9999 END ********************************************************************** SUBROUTINE TRANSF (NOGO, XPT, YPT, XPT2, YPT2, + DLAM, DLOM, DLAS, DLOS, KEY) ********************************************************************** * This subroutine computes either the forward or inverse coordinate * transformation depending upon the value of the integer variable 'key' IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) INTEGER MXAREA, ITMAX DOUBLE PRECISION SMALL PARAMETER (MXAREA = 7, ITMAX = 10, SMALL = 1.0D-9 ) DOUBLE PRECISION XPT, YPT, XPT2, YPT2 DOUBLE PRECISION XTEMP, YTEMP, XDIF, YDIF DOUBLE PRECISION DLAM, DLOM, DLAS, DLOS DOUBLE PRECISION DXLAST, DYLAST INTEGER*2 KEY, NUM LOGICAL NOGO INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA) IF (KEY .EQ. 1) THEN ********************** * FOR NAD 27 TO NAD 83 ********************** CALL TO83 (NOGO, XPT, YPT, XPT2, YPT2, + DLAM, DLOM, DLAS, DLOS) ELSEIF (KEY .EQ. 2) THEN *************************** * FOR NAD 83 TO NAD 27) * THIS IS DONE BY ITERATION *************************** NUM = 0 ************************************************** * SET THE XPT,YPT TO TEMPORARY VALUES * (REMEMBER, XPT AND YPT ARE REALLY NAD 83 VALUES) ************************************************** XTEMP = XPT YTEMP = YPT ************************************************************** * PRETEND THAT THESE TEMPORARY VALUES ARE REALLY NAD 27 VALUES * FOR A FIRST GUESS AND COMPUTE PSEUDO-NAD 83 COORDINATES ************************************************************** 200 CONTINUE NUM = NUM + 1 ******************************** * CHECK THE NUMBER OF ITERATIONS ******************************** IF (NUM .GE. ITMAX) THEN WRITE (LUOUT,*) ' *** MAXIMUM ITERATIONS EXCEEDED!! ***' WRITE (LUOUT,*) ' *** CALL PROGRAMMER FOR HELP ***' WRITE (LUOUT,*) ' LATITUDE =', YTEMP, ' LONGITUDE =', XTEMP NOGO = .TRUE. GOTO 1000 ENDIF CALL TO83 (NOGO, XTEMP, YTEMP, XPT2, YPT2, + DLAM, DLOM, DLAS, DLOS) DXLAST = DLOS DYLAST = DLAS ************************************** * COMPARE TO ACTUAL NAD 83 COORDINATES ************************************** XDIF = XPT - XPT2 YDIF = YPT - YPT2 **************************************************************** * COMPUTE A NEW GUESS UNLESS THE DIFFERENCES ARE LESS THAN SMALL * WHERE SMALL IS DEFINED (ABOVE) TO BE; SMALL = 1.0D-9 **************************************************************** IF (NUM .EQ. 1) THEN IF (DABS(XDIF) .GT. SMALL) THEN XTEMP = XPT - DLOS/3600.D0 ENDIF IF (DABS(YDIF) .GT. SMALL) THEN YTEMP = YPT - DLAS/3600.D0 ENDIF ELSE IF (DABS(XDIF) .GT. SMALL) THEN XTEMP = XTEMP - (XPT2 - XPT) ENDIF IF (DABS(YDIF) .GT. SMALL) THEN YTEMP = YTEMP - (YPT2 - YPT) ENDIF ENDIF IF (DABS(YDIF) .LE. SMALL .AND. DABS(XDIF) .LE. SMALL) THEN ****************************** * IF CONVERGED THEN LEAVE LOOP ****************************** XPT = XTEMP YPT = YTEMP GOTO 1000 ENDIF ******************************* * IF NOT CONVERGED THEN ITERATE ******************************* GOTO 200 ENDIF 1000 RETURN END ********************************************************************** SUBROUTINE SURF (XGRID, YGRID, ZEE, AY, BEE, CEE, DEE, IROW, JCOL) ********************************************************************** ** SUBROUTINE SURF: INTERPOLATES THE Z VALUE * ********************************************************************** * Calculated the value of the grid at the point XPT, YPT. The * interpolation is done in the index coordinate system for convenience. IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) * IMPLICIT UNDEFINED (A-Z) DOUBLE PRECISION XGRID, YGRID DOUBLE PRECISION AY, BEE, CEE, DEE DOUBLE PRECISION ZEE, ZEE1, ZEE2, ZEE3, ZEE4 INTEGER IROW, JCOL ZEE1 = AY ZEE2 = BEE*(XGRID - DBLE(JCOL) ) ZEE3 = CEE*(YGRID - DBLE(IROW) ) ZEE4 = DEE*(XGRID - DBLE(JCOL) )*(YGRID - DBLE(IROW) ) ZEE = ZEE1 + ZEE2 + ZEE3 + ZEE4 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 DOUBLE PRECISION (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