C SccsID = "@(#)usng.f 1.14 10/29/02" PROGRAM USNG * IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*10 CON_TYPE CHARACTER*13 IN_LAT CHARACTER*14 IN_LON CHARACTER*14 IN_NORTH CHARACTER*2 IN_ZONE CHARACTER*14 IN_EAST CHARACTER*2 IN_DATUM CHARACTER*4 IN_ACCURACY CHARACTER*25 USNG_VALUE CHARACTER*20 PARAM(10) INTEGER INUNIT, OUTUNIT INTEGER PCODE INTEGER NUM_ARGS LOGICAL FILFLAG,FILPRT,FIL81FL COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/SETS/SET1,SET2,SET3,SET4,SET5,SET6 COMMON /UNITS/INUNIT, OUTUNIT COMMON/PAR/PARAM COMMON/PRECISE/PCODE DATA PARAM/10*' '/ C ----------------------------------------- C Initialize UNITS to 0 C Will get reset in BY_FILE() if used C ----------------------------------------- INUNIT = 0 OUTUNIT = 0 FILFLAG=.FALSE. FILPRT=.FALSE. FIL81FL=.FALSE. PI=4.D0*DATAN(1.D0) RAD=180.D0/PI IN_LAT = ' ' IN_LON = ' ' IN_NORTH = ' ' IN_EAST = ' ' IN_ZONE = ' ' IN_DATUM = ' ' IN_ACCURACY = ' ' PCODE = 1 USNG_VALUE = ' ' CALL FILL_SET1() CALL FILL_SET2() CALL FILL_SET3() CALL FILL_SET4() CALL FILL_SET5() CALL FILL_SET6() CON_TYPE = ' ' IN_LAT = ' ' IN_LON = ' ' IN_DATUM = ' ' USNG_VALUE = ' ' NUM_ARGS = IARGC() C ---------------------------------------------------------- C Get the command line parameters and convert to UPPER case C Except FILE=myfile Parameter C ---------------------------------------------------------- DO I=1, NUM_ARGS CALL GETARG(I, PARAM(I)) IF (INDEX(PARAM(I), '=') .LE. 0) THEN CALL TO_UPPER(PARAM(I), LEN(PARAM(I))) ENDIF END DO C --------------------------------------------------------- C Give USAGE info if command line parameters are missing C --------------------------------------------------------- IF (PARAM(1) .EQ. 'CMD') THEN CALL PRINT_CMD_USAGE() STOP ELSEIF (PARAM(1) .EQ. 'FILE') THEN CALL PRINT_FILE_USAGE() STOP ENDIF IF (NUM_ARGS .LT. 2) THEN CALL PRINT_VERSION() CALL PRINT_DISCLAIMER() CALL BY_PROMPT() STOP ENDIF C ----------------------------- C Get the conversion type C GP2US, US2GP, UTM2US, US2UTM C ----------------------------- CON_TYPE = PARAM(1) IF ((CON_TYPE .NE. 'GP2US').AND. * (CON_TYPE .NE. 'US2GP').AND. * (CON_TYPE .NE. 'UTM2US').AND. * (CON_TYPE .NE. 'US2UTM')) THEN WRITE(6, 20) 20 FORMAT (/, * '----------------------------------------------------------',/, * 'ERROR - Parameter "GP2US","US2GP, "UTM2US", "US2UTM" ', * 'is missing', /, * ' TYPE "usng CMD" TO SEE EXAMPLE.', /, * '----------------------------------------------------------',/) STOP ENDIF C -------------------------------------------------------- C Get command line parameters and do the conversion C -------------------------------------------------------- IF ((INDEX(PARAM(2), "FILE") .GT. 0).OR. *(INDEX(PARAM(2), "file") .GT. 0)) THEN CALL BY_FILE() ELSEIF (CON_TYPE .EQ. 'GP2US') THEN IN_LAT = PARAM(2) IN_LON = PARAM(3) IN_DATUM = PARAM(4) CALL CHECK_IN_LAT(IN_LAT, 0, 0) CALL CHECK_IN_LON(IN_LON, 0, 0) CALL CHECK_IN_DATUM(IN_DATUM, 0, 0) CALL CHECK_IN_LAT(IN_LAT, 0, 0) CALL CHECK_IN_LON(IN_LON, 0, 0) CALL CHECK_IN_DATUM(IN_DATUM, 0, 0) CALL GP_TO_US(IN_LAT, IN_LON, IN_DATUM) ELSEIF (CON_TYPE .EQ. 'US2GP') THEN USNG_VALUE = PARAM(2) IN_DATUM = PARAM(3) CALL CHECK_IN_USNG(USNG_VALUE, 0, 0) CALL CHECK_IN_DATUM(IN_DATUM, 0, 0) CALL US_TO_GP (USNG_VALUE, IN_DATUM) ELSEIF (CON_TYPE .EQ. 'UTM2US') THEN IN_ZONE = PARAM(2) IN_EAST = PARAM(3) IN_NORTH = PARAM(4) IN_DATUM = PARAM(5) IN_ACCURACY = PARAM(6) CALL CHECK_IN_EAST(IN_EAST, 0, 0) CALL CHECK_IN_NORTH(IN_NORTH, 0, 0) CALL CHECK_IN_ZONE(IN_ZONE, 0, 0) CALL CHECK_IN_DATUM(IN_DATUM, 0, 0) CALL CHECK_IN_ACCURACY(IN_ACCURACY, 0, 0) IN_ACCURACY = PARAM(6) IF(IN_ACCURACY .EQ. 'A100') THEN PCODE = 3 ELSEIF (IN_ACCURACY .EQ. 'A10') THEN PCODE = 2 ELSE PCODE = 1 ENDIF CALL UTM_TO_US(IN_ZONE, IN_EAST, IN_NORTH, IN_DATUM) ELSEIF (CON_TYPE .EQ. 'US2UTM') THEN USNG_VALUE = PARAM(2) IN_DATUM = PARAM(3) CALL CHECK_IN_USNG(USNG_VALUE, 0, 0) CALL CHECK_IN_DATUM(IN_DATUM, 0, 0) CALL US_TO_UTM (USNG_VALUE, IN_DATUM) ELSE WRITE(6, 3000) 3000 FORMAT(//, * 'ERROR - Unable to determine direction of conversion',/) CALL PRINT_CMD_USAGE() ENDIF IF (OUTUNIT .GT. 6) THEN WRITE(6, 5000) 5000 FORMAT(//, * ' ----------------------------------------------------',/, * ' Output File = "usng.out"', /, * ' ----------------------------------------------------',/) ENDIF STOP END C---------------------------------------------- SUBROUTINE BY_FILE() INTEGER INUNIT, OUTUNIT CHARACTER*80 FILENAME INTEGER START_READ CHARACTER*10 CON_TYPE CHARACTER*70 LINE CHARACTER*13 IN_LAT CHARACTER*14 IN_LON CHARACTER*25 USNG_VALUE CHARACTER*2 IN_DATUM CHARACTER*14 IN_EAST CHARACTER*2 IN_ZONE CHARACTER*14 IN_NORTH CHARACTER*4 IN_ACCURACY INTEGER ICOUNT CHARACTER*20 PARAM(10) INTEGER PCODE COMMON /UNITS/INUNIT, OUTUNIT COMMON/PAR/PARAM COMMON/PRECISE/PCODE C ---------------------------------------------------------- C Extract the CON_TYPE and the 2nd parameter from the command line C US2GP FILE=myfile C ---------------------------------------------------------- CON_TYPE = PARAM(1) IF (INDEX(PARAM(2), "=") .LE. 0) then WRITE(6, 10) 10 FORMAT(/, * '----------------------------------------------------',/, * 'ERROR - Unable to determine file name.',/, * ' Filename should be passed in the format...',/, * ' FILE=myfile',/, * ' i.e. NO SPACES AROUND "=" ', /, * ' TYPE "usng FILE"TO SEE EXAMPLE', /, * '----------------------------------------------------',/) STOP ENDIF FILENAME = ' ' C ----------------------------------------------------- C Extract the input filename from the 2nd parameter C i.e. remove 'FILE=' text C ----------------------------------------------------- J=0 START_READ=0 DO I=1,80 IF ((START_READ .EQ. 1).AND.(PARAM(2)(I:I) .NE. ' ')) THEN J = J + 1 FILENAME(J:J) = PARAM(2)(I:I) ENDIF IF (PARAM(2)(I:I) .EQ. '=') THEN START_READ=1 ENDIF END DO WRITE(6, *) FILENAME C ----------------------- C Assign I/O Units C ----------------------- INUNIT = 11 CALL OPEN_OUTFILE() ICOUNT = 0 OPEN(INUNIT, FILE=FILENAME, STATUS='OLD', *ACCESS='SEQUENTIAL', ERR=70) GOTO 100 70 WRITE(6, 71) FILENAME 71 FORMAT ('Cannot open file ', A50) STOP C ---------------------------------------------------------- C READ THE INPUT FILE C ---------------------------------------------------------- 100 CONTINUE LINE = ' ' READ(INUNIT, 110, END=1000) LINE 110 FORMAT(A) ICOUNT = ICOUNT + 1 IF (CON_TYPE .EQ. 'GP2US') THEN IN_LAT = ' ' IN_LON = ' ' IN_DATUM = ' ' CALL X_GP_LINE(LINE, IN_LAT, IN_LON, IN_DATUM) CALL CHECK_IN_LAT(IN_LAT, ICOUNT, 0) CALL CHECK_IN_LON(IN_LON, ICOUNT, 0) CALL CHECK_IN_DATUM(IN_DATUM, ICOUNT, 0) CALL GP_TO_US(IN_LAT, IN_LON, IN_DATUM) ELSEIF (CON_TYPE .EQ. 'US2GP') THEN USNG_VALUE = ' ' IN_DATUM = ' ' CALL X_US_LINE(LINE, USNG_VALUE, IN_DATUM) CALL CHECK_IN_USNG(USNG_VALUE, ICOUNT, 0) CALL CHECK_IN_DATUM(IN_DATUM, ICOUNT, 0) CALL US_TO_GP (USNG_VALUE, IN_DATUM) ELSEIF (CON_TYPE .EQ. 'UTM2US') THEN IN_ZONE = ' ' IN_EAST = ' ' IN_NORTH = ' ' IN_DATUM = ' ' IN_ACCURACY = ' ' CALL X_UTM_LINE(LINE, IN_ZONE, IN_EAST, IN_NORTH, IN_DATUM, * IN_ACCURACY) CALL CHECK_IN_EAST(IN_EAST, ICOUNT, 0) CALL CHECK_IN_NORTH(IN_NORTH, ICOUNT, 0) CALL CHECK_IN_DATUM(IN_DATUM, ICOUNT, 0) CALL CHECK_IN_ACCURACY(IN_ACCURACY, ICOUNT, 0) IF(IN_ACCURACY .EQ. 'A100') THEN PCODE = 3 ELSEIF (IN_ACCURACY .EQ. 'A10') THEN PCODE = 2 ELSE PCODE = 1 ENDIF CALL UTM_TO_US(IN_ZONE,IN_EAST, IN_NORTH, IN_DATUM) ELSEIF (CON_TYPE .EQ. 'US2UTM') THEN USNG_VALUE = ' ' IN_DATUM = ' ' CALL X_US_LINE(LINE, USNG_VALUE, IN_DATUM) CALL CHECK_IN_USNG(USNG_VALUE, ICOUNT, 0) CALL CHECK_IN_DATUM(IN_DATUM, ICOUNT, 0) CALL US_TO_UTM (USNG_VALUE, IN_DATUM) ENDIF GOTO 100 1000 CONTINUE RETURN END C----------------------------------------------------------- SUBROUTINE X_GP_LINE(LINE, LAT, LON, DATUM) CHARACTER*70 LINE CHARACTER*13 LAT CHARACTER*14 LON CHARACTER*2 DATUM INTEGER L, O, D INTEGER COUNT COUNT = 0 L = 0 M = 0 P = 0 C --------------------------------------------------- C Extract LAT, LON, DATUM from the delimited line C --------------------------------------------------- DO I=1,70 IF (LINE(I:I) .EQ. ' ') THEN CONTINUE ELSEIF (LINE(I:I) .EQ. '|') THEN COUNT = COUNT + 1 ELSE IF (COUNT .EQ. 2) THEN L = L + 1 LAT(L:L) = LINE(I:I) ELSE IF (COUNT .EQ. 3) THEN M = M + 1 LON(M:M) = LINE(I:I) ELSE IF (COUNT .EQ. 4) THEN P = P + 1 DATUM(P:P) = LINE(I:I) ENDIF END DO RETURN END C----------------------------------------------------------- SUBROUTINE X_UTM_LINE(LINE, ZONE, EAST, NORTH, DATUM, * ACCURACY) CHARACTER*70 LINE CHARACTER*2 ZONE CHARACTER*14 EAST CHARACTER*14 NORTH CHARACTER*4 ACCURACY CHARACTER*2 DATUM INTEGER L, M, P, Q, R INTEGER COUNT COUNT = 0 L = 0 M = 0 P = 0 Q = 0 R = 0 C ------------------------------------------------------------- C Extract EAST, NORTH, ACCURACY, DATUM from the delimited line C ------------------------------------------------------------- DO I=1,70 IF (LINE(I:I) .EQ. ' ') THEN CONTINUE ELSEIF (LINE(I:I) .EQ. '|') THEN COUNT = COUNT + 1 ELSE IF (COUNT .EQ. 2) THEN L = L + 1 ZONE(L:L) = LINE(I:I) ELSE IF (COUNT .EQ. 3) THEN M = M + 1 EAST(M:M) = LINE(I:I) ELSE IF (COUNT .EQ. 4) THEN P = P + 1 NORTH(P:P) = LINE(I:I) ELSE IF (COUNT .EQ. 5) THEN Q = Q + 1 DATUM(Q:Q) = LINE(I:I) ELSE IF (COUNT .EQ. 6) THEN R = R + 1 ACCURACY(R:R) = LINE(I:I) ENDIF END DO RETURN END C----------------------------------------------------------- SUBROUTINE X_US_LINE(LINE, USNG_VALUE, DATUM) CHARACTER*70 LINE CHARACTER*25 USNG_VALUE CHARACTER*2 DATUM INTEGER L, O, D INTEGER COUNT COUNT = 0 L = 0 M = 0 P = 0 C --------------------------------------------------- C Extract USGN from the delimited line C --------------------------------------------------- DO I=1,70 IF (LINE(I:I) .EQ. ' ') THEN CONTINUE ELSEIF (LINE(I:I) .EQ. '|') THEN COUNT = COUNT + 1 ELSE IF (COUNT .EQ. 2) THEN L = L + 1 USNG_VALUE(L:L) = LINE(I:I) ELSE IF (COUNT .EQ. 3) THEN P = P + 1 DATUM(P:P) = LINE(I:I) ENDIF END DO RETURN END *********************************************************** SUBROUTINE GP_TO_US(IN_LAT, IN_LON, IN_DATUM) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 DATNUM CHARACTER*10 CON_TYPE CHARACTER*13 IN_LAT CHARACTER*14 IN_LON CHARACTER*2 IN_DATUM CHARACTER*25 USNG_VALUE CHARACTER*80 CARDR REAL*8 NORTH, EAST INTEGER PCOUNT INTEGER IZ INTEGER PCODE LOGICAL FILFLAG,FILPRT,FIL81FL COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/DATUM/DATNUM COMMON/XY/NORTH,EAST COMMON/SETS/SET1,SET2,SET3,SET4,SET5,SET6 COMMON/PRECISE/PCODE FILFLAG=.FALSE. FILPRT=.FALSE. FIL81FL=.FALSE. PI=4.D0*DATAN(1.D0) RAD=180.D0/PI CON_TYPE = 'GP2US' USNG_VALUE = ' ' IF (IN_DATUM .EQ. '27') THEN DATNUM = '1' ELSE DATNUM = '2' ENDIF CALL DATUMM(ER,RF,F,ESQ,DATNUM) CARDR = ' ' CALL BUILD_80(IN_LAT, IN_LON, CARDR) C ----------------------------------------------------- C PCODE is used to determine output precision C 1 = input 100th Sec, output = 1 Meter C 2 = input 10th Sec, output = 10 Meters C 3 = input 1 Sec, output = 100 Meters C ----------------------------------------------------- PCOUNT = 0 DO I=1, LEN(IN_LAT) IF (IN_LAT(I:I) .NE. ' ') THEN PCOUNT = PCOUNT + 1 ENDIF END DO IF (PCOUNT .GT. 9) THEN PCODE = 1 ELSEIF (PCOUNT .GT. 8) THEN PCODE = 2 ELSE PCODE = 3 ENDIF CALL DRGPUS(IN_LAT,IN_LON,IN_DATUM,CARDR,FILFLAG,FIL81FL) RETURN END C------------------------------------------------------- SUBROUTINE US_TO_GP (USNG_VALUE, IN_DATUM) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 DATNUM CHARACTER*10 CON_TYPE CHARACTER*13 IN_LAT CHARACTER*14 IN_LON CHARACTER*2 IN_DATUM CHARACTER*25 USNG_VALUE CHARACTER*80 CARDR REAL*8 NORTH, EAST INTEGER Z INTEGER PCODE LOGICAL FILFLAG,FILPRT,FIL81FL COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/DATUM/DATNUM COMMON/XY/NORTH,EAST COMMON/SETS/SET1,SET2,SET3,SET4,SET5,SET6 COMMON/PRECISE/PCODE FILFLAG=.FALSE. FILPRT=.FALSE. FIL81FL=.FALSE. PI=4.D0*DATAN(1.D0) RAD=180.D0/PI CON_TYPE = 'US2GP' IN_LAT = ' ' IN_LON = ' ' IF (IN_DATUM .EQ. '27') THEN DATNUM = '1' ELSE DATNUM = '2' ENDIF CALL DATUMM(ER,RF,F,ESQ,DATNUM) CARDR = ' ' Z=0 DO I=1, LEN(USNG_VALUE) IF (USNG_VALUE(I:I) .NE. ' ') THEN Z = Z + 1 ENDIF END DO USNG_VALUE = USNG_VALUE(1:Z) // '(NAD '// IN_DATUM // ')' NORTH = 0.0 EAST = 0.0 IZ = 0 C --------------------------------------------------- C Initialize PCODE, It should get reset in USNG2UTM() C --------------------------------------------------- PCODE = 1 CALL USNG2UTM(NORTH, EAST, IZ, IN_DATUM, * USNG_VALUE) IF ((INDEX("CDEFGHJKLM", USNG_VALUE(2:2)) .GT. 0) .OR. * (INDEX("CDEFGHJKLM", USNG_VALUE(3:3)) .GT. 0)) THEN NORTH = 20000000.D0 - NORTH ENDIF CARDR(69:69) = 'W' CALL DRUSGP(USNG_VALUE,CARDR,IN_DATUM, * IZ,FILFLAG,FILPRT) RETURN END *********************************************************** SUBROUTINE UTM_TO_US(IN_ZONE, IN_EAST, IN_NORTH, IN_DATUM) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 DATNUM CHARACTER*10 CON_TYPE CHARACTER*2 IN_ZONE CHARACTER*14 IN_EAST CHARACTER*14 IN_NORTH CHARACTER*2 IN_DATUM CHARACTER*25 USNG_VALUE REAL*8 NORTH, EAST REAL*8 DLAT REAL*8 ORIGIAN_DLAT INTEGER IZ CHARACTER*1 GZD CHARACTER*10 PTEXT INTEGER PCODE INTEGER INUNIT, OUTUNIT INTEGER OCOUNT CHARACTER*1 NORS SAVE OCOUNT LOGICAL FILFLAG,FILPRT,FIL81FL COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/DATUM/DATNUM COMMON/SETS/SET1,SET2,SET3,SET4,SET5,SET6 COMMON/PRECISE/PCODE COMMON /UNITS/INUNIT, OUTUNIT COMMON/XY/NORTH,EAST C -------------------------------------------- C Look for a minus sign in IN_NORTH C which indicates Southern Hemisphere C -------------------------------------------- IF (INDEX(IN_NORTH, '-') .GT. 0) THEN NORS= 'S' ELSE NORS='N' ENDIF FILFLAG=.FALSE. FILPRT=.FALSE. FIL81FL=.FALSE. PI=4.D0*DATAN(1.D0) RAD=180.D0/PI CON_TYPE = 'GP2US' USNG_VALUE = ' ' IF (IN_DATUM .EQ. '27') THEN DATNUM = '1' ELSE DATNUM = '2' ENDIF CALL DATUMM(ER,RF,F,ESQ,DATNUM) C ----------------------------------------------------- C PCODE is used to determine output precision C 1 = input 100th Sec, output = 1 Meter C 2 = input 10th Sec, output = 10 Meters C 3 = input 1 Sec, output = 100 Meters C and is obtained from the A1 parameter of command line C ----------------------------------------------------- IF (PCODE .EQ. 3) THEN PTEXT = '(100M)' ELSEIF (PCODE .EQ. 2) THEN PTEXT = '( 10M)' ELSE PTEXT = '( 1M)' ENDIF READ(IN_ZONE, 20) IZ 20 FORMAT(I2) READ(IN_EAST, 30) EAST 30 FORMAT(F11.3) READ(IN_NORTH, 40) NORTH 40 FORMAT(F11.3) C -------------------------------------------------- C Get the GZD C -------------------------------------------------- DLAT = 0.D0 CALL X_DLAT (IN_ZONE, DLAT) ORIGINAL_DLAT = DLAT CALL FIND_GZD(NORS, ORIGINAL_DLAT, GZD) C ---------------------------------------- C Compute the US National Grid C ---------------------------------------- USNG_VALUE = ' ' CALL UTM2USNG(GZD, NORTH, EAST, IZ, IN_DATUM, * USNG_VALUE) C USNG_VALUE = '18UJ2348306480' OCOUNT = OCOUNT + 1 IF (OCOUNT .EQ. 1) THEN CALL PRINT_UTM2US_HEADER() ENDIF WRITE(6,2000) IN_ZONE,IN_EAST,IN_NORTH,PTEXT,IN_DATUM, * USNG_VALUE 2000 FORMAT(A2,' ', A13,' ',A13,' ',A6, ' ','NAD ',A2,' -> ',A23) IF (OUTUNIT .GT. 6) THEN WRITE(OUTUNIT,2001) IN_ZONE,IN_EAST,IN_NORTH,PTEXT,IN_DATUM, * USNG_VALUE 2001 FORMAT(A2,' ', A13,' ',A13,' ',A6,' ','NAD ',A2,' -> ',A23) ENDIF RETURN END SUBROUTINE US_TO_UTM (USNG_VALUE, IN_DATUM) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*1 DATNUM CHARACTER*10 CON_TYPE CHARACTER*2 IN_DATUM CHARACTER*25 USNG_VALUE REAL*8 NORTH, EAST CHARACTER*1 SIGN INTEGER Z INTEGER PCODE CHARACTER*10 PTEXT INTEGER INUNIT, OUTUNIT INTEGER OCOUNT SAVE OCOUNT LOGICAL FILFLAG,FILPRT,FIL81FL COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/DATUM/DATNUM COMMON/XY/NORTH,EAST COMMON/SETS/SET1,SET2,SET3,SET4,SET5,SET6 COMMON/PRECISE/PCODE COMMON /UNITS/INUNIT, OUTUNIT FILFLAG=.FALSE. FILPRT=.FALSE. FIL81FL=.FALSE. PI=4.D0*DATAN(1.D0) RAD=180.D0/PI CON_TYPE = 'US2GP' IF (IN_DATUM .EQ. '27') THEN DATNUM = '1' ELSE DATNUM = '2' ENDIF CALL DATUMM(ER,RF,F,ESQ,DATNUM) Z=0 DO I=1, LEN(USNG_VALUE) IF (USNG_VALUE(I:I) .NE. ' ') THEN Z = Z + 1 ENDIF END DO USNG_VALUE = USNG_VALUE(1:Z) // '(NAD '// IN_DATUM // ')' NORTH = 0.0 EAST = 0.0 IZ = 0 C --------------------------------------------------- C Initialize PCODE, It should get reset in USNG2UTM() C --------------------------------------------------- PCODE = 1 CALL USNG2UTM(NORTH, EAST, IZ, IN_DATUM, USNG_VALUE) IF (PCODE .EQ. 3) THEN PTEXT = '(100M)' ELSEIF (PCODE .EQ. 2) THEN PTEXT = '( 10M)' ELSE PTEXT = '( 1M)' ENDIF OCOUNT = OCOUNT + 1 IF (OCOUNT .EQ. 1) THEN CALL PRINT_US2UTM_HEADER() ENDIF IF ((INDEX("CDEFGHJKLM", USNG_VALUE(2:2)) .GT. 0) .OR. * (INDEX("CDEFGHJKLM", USNG_VALUE(3:3)) .GT. 0)) THEN SIGN = '-' ELSE SIGN = ' ' ENDIF WRITE(6, 2000) USNG_VALUE, IZ, EAST, SIGN, NORTH, *PTEXT, IN_DATUM 2000 FORMAT(A23,' -> ',I2,' ',F10.1,' ',A1,F10.1,' ', *A6, ' ','NAD ',A2) IF (OUTUNIT .GT. 6) THEN WRITE(OUTUNIT, 2001) USNG_VALUE, IZ, EAST, SIGN, NORTH, * PTEXT, IN_DATUM 2001 FORMAT(A23,' -> ',I2,' ',F10.1,' ',A1,F10.1,' ', * A6, ' ','NAD ',A2) ENDIF RETURN END ************************************************************************ SUBROUTINE DRGPUS(IN_LAT,IN_LON,IN_DATUM,CARDR,FILFLAG,FIL81FL) ********************************************************************** * * * THIS IS THE DRIVER TO COMPUTE UTM NORTHINGS AND EASTINGS * FOR EACH PRIMARY ZONE AND THE AJACENT ZONE IF THE LONGITUDE * IS WITH 5 MINUTES OF THE ZONE BOUNDARIES * * THE OUTPUT IS FOR THE DATA SHEET PROGRAM * * VARIABLES * CARDR = A MODIFIED 80 RECORD CARD WITH A LENGTH OF 211 COLS * ER = EQUATORIAL RADIUS OF THE ELLIPSOID (SEMI-MAJOR AXIS) * RF = RECIPROCAL OF FLATTING OF THE ELLIPSOD * ESQ= E SQUARED * RAD = RADIAN CONVERSION FACTOR * CM = CENTRAL MERIDIAN ( COMPUTED USEING THE LONGITUDE) * SF = SCALE FACTOR OF CENTRAL MERIDIAN ( ALWAYS .9996 FOR UTM) * OR = SOUTHERNMOST PARALLEL OF LATITUDE ( ALWAYS ZERO FOR UTM) * R, A, B, C, U, V, W = ELLIPSOID CONSTANTS USED FOR COMPUTING * MERIDIONAL DISTANCE FROM LATITUDE * SO = MERIDIONAL DISTANCE (MULTIPLIED BY SCALE FACTOR ) * FROM THE EQUATOR TO THE SOUTHERNMOST PARALLEL OF LATITUDE * ( ALWAYS ZERO FOR UTM) * IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*13 IN_LAT CHARACTER*14 IN_LON CHARACTER*2 IN_DATUM CHARACTER*80 CARDR CHARACTER*1 EORW CHARACTER*1 NORS CHARACTER*4 ZONE INTEGER*4 LD,LM,LOD,LOM,FOUND INTEGER OCOUNT REAL*8 SLAT,SLON,FI,LAM,LCM,UCM REAL*8 KP REAL*8 TD,TM,ND,NM REAL*8 NORTH REAL*8 LOD1 REAL*8 DLAT REAL*8 ORIGINAL_DLAT CHARACTER*1 GZD CHARACTER*25 USNG_VALUE INTEGER INUNIT, OUTUNIT SAVE OCOUNT LOGICAL FILFLAG LOGICAL FIL81FL COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/SETS/SET1,SET2,SET3,SET4,SET5,SET6 COMMON /UNITS/INUNIT, OUTUNIT READ(CARDR,50)LD,LM,SLAT,NORS,LOD,LOM,SLON,EORW 50 FORMAT(T45,I2,I2,F7.5,A1,I3,I2,F7.5,A1) * * CONVERT THE LATITUDE AND LONGITUDE TO PI AND LAM * TD=DBLE(FLOAT(LD)) TM=DBLE(FLOAT(LM)) FI=(TD+(TM+SLAT/60.D0)/60.D0)/RAD ND=DBLE(FLOAT(LOD)) NM=DBLE(FLOAT(LOM)) IF((EORW.EQ.'E').OR.(EORW.EQ.'e')) THEN LAM=(360.D0-(ND+(NM+SLON/60.D0)/60.D0))/RAD LOD1=(360.D0-(ND+(NM+SLON/60.D0)/60.D0)) LOD=DINT(LOD1) ENDIF IF((EORW.EQ.'W').OR.(EORW.EQ.'w')) THEN LAM=(ND+(NM+SLON/60.D0)/60.D0)/RAD LOD=LOD ENDIF * ----------------------------------------- * Find the GZD * ----------------------------------------- IF ((NORS .EQ. 'S') .OR. (NORS .EQ. 's')) THEN DLAT = (-1.D0) * (TD+(TM+SLAT/60.D0)/60.D0) ELSE DLAT = (TD+(TM+SLAT/60.D0)/60.D0) ENDIF ORIGINAL_DLAT = DLAT CALL FIND_GZD(NORS, ORIGINAL_DLAT, GZD) * * FIND THE ZONE FOR LONGITUDE LESS THAN 180 DEGREES * IF(LOD.LT.180) THEN IZ=LOD/6 IZ= 30 -IZ ICM=(183-(6*IZ)) CM=DBLE(FLOAT(ICM))/RAD UCM=(ICM+3)/RAD LCM=(ICM-3)/RAD ENDIF * * FIND THE ZONE FOR LONGITUDE GREATER THAN 180 DEGREES * IF(LOD.GE.180) THEN IZ=(LOD)/6 IZ= 90 - IZ ICM=(543 - (6*IZ)) CM= DBLE(FLOAT(ICM))/RAD UCM=(ICM+3)/RAD LCM=(ICM-3)/RAD ENDIF TOL=(5.0D0/60.0D0)/RAD FN = 0.D0 IF((NORS.EQ.'S').OR.(NORS.EQ.'s')) THEN FN = 10000000.D0 ENDIF IF((NORS.EQ.'N').OR.(NORS.EQ.'n')) THEN FN = 0.D0 ENDIF FE=500000.0D0 SF=0.9996D0 OR=0.0D0 FOUND=0 CALL TCONST(ER,RF,SF,OR,ESQ,EPS,R,A,B,C,U,V,W,SO, & CM,FE,FN) * * COMPUTE THE NORTH AND EASTINGS * 200 CALL TMGRID(FI,LAM,NORTH,EAST,CONV,KP,ER,ESQ,EPS,CM, & FE,FN,SF,SO,R,A,B,C,U,V,W) * * WRITE THE ZONE NUMBER * IF (IZ.GT.9) THEN WRITE(ZONE,600) IZ 600 FORMAT(1X,I2) ELSE WRITE(ZONE,605) IZ 605 FORMAT(1X,I2.2) ENDIF * ----------------------------------------- * Find the GZD * ----------------------------------------- DLAT = (TD+(TM+SLAT/60.D0)/60.D0) CALL FIND_GZD(NORS, ORIGINAL_DLAT, GZD) C ---------------------------------------- C Compute the US National Grid C ---------------------------------------- USNG_VALUE = ' ' CALL UTM2USNG(GZD, NORTH, EAST, IZ, IN_DATUM, * USNG_VALUE) OCOUNT = OCOUNT + 1 IF (OCOUNT .EQ. 1) THEN CALL PRINT_GP2US_HEADER() ENDIF WRITE(6, 2000) IN_LAT, IN_LON, IN_DATUM, USNG_VALUE 2000 FORMAT(A13,' ',A14,' ','NAD ',A2,' -> ',A23) IF (OUTUNIT .GT. 6) THEN WRITE(OUTUNIT, 2001) IN_LAT, IN_LON, IN_DATUM, USNG_VALUE 2001 FORMAT(A13,' ',A14,' ','NAD ',A2,' -> ',A23) ENDIF RETURN END ******************************************************************* SUBROUTINE DRUSGP(USNG_VALUE,CARDR,IN_DATUM, * ICODE,FILFLAG,FILPRT) * * * * THIS IS THE DRIVER TO COMPUTE LATITUDES AND LONGITUDES FROM * THE UTMs FOR EACH ZONE * * * VARIABLES * CARDR = A MODIFIED 80 RECORD CARD WITH A LENGTH OF 211 COLS * ER = EQUATORIAL RADIUS OF THE ELLIPSOID (SEMI-MAJOR AXIS) * RF = RECIPROCAL OF FLATTING OF THE ELLIPSOD * ESQ= E SQUARED * RAD = RADIAN CONVERSION FACTOR * CM = CENTRAL MERIDIAN ( COMPUTED USEING THE LONGITUDE) * SF = SCALE FACTOR OF CENTRAL MERIDIAN ( ALWAYS .9996 FOR UTM) * OR = SOUTHERNMOST PARALLEL OF LATITUDE ( ALWAYS ZERO FOR UTM) * R, A, B, C, U, V, W = ELLIPSOID CONSTANTS USED FOR COMPUTING * MERIDIONAL DISTANCE FROM LATITUDE * SO = MERIDIONAL DISTANCE (MULTIPLIED BY SCALE FACTOR ) * FROM THE EQUATOR TO THE SOUTHERNMOST PARALLEL OF LATITUDE * ( ALWAYS ZERO FOR UTM) * IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*25 USNG_VALUE CHARACTER*80 CARDR CHARACTER*2 IN_DATUM CHARACTER*1 NORS CHARACTER*4 ZONE INTEGER*4 FOUND REAL*8 LAT,LON,LCM,UCM REAL*8 KP REAL*8 NORTH,EAST CHARACTER*13 OLAT CHARACTER*14 OLON INTEGER OCOUNT INTEGER INUNIT, OUTUNIT SAVE OCOUNT LOGICAL FILFLAG,FILPRT COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/XY/NORTH,EAST COMMON /UNITS/INUNIT, OUTUNIT * * FIND THE CENTRAL MERIDAIN IF THE ZONE NUMBER IS LESS THAN 30 * IF(ICODE.LT.30) THEN IZ=ICODE ICM=(183-(6*IZ)) CM=DBLE(FLOAT(ICM))/RAD UCM=(ICM+3)/RAD LCM=(ICM-3)/RAD ENDIF * * FIND THE CENTRAL MERIDAN IF THE ZONE NUMBER IS LARGER THAN 30 * IF(ICODE.GE.30) THEN IZ=ICODE ICM=(543 - (6*IZ)) CM= DBLE(FLOAT(ICM))/RAD UCM=(ICM+3)/RAD LCM=(ICM-3)/RAD ENDIF TOL=(5.0D0/60.0D0)/RAD IF(NORTH.GE.10000000.0) THEN FN= 10000000.0D0 NORS= 'S' ELSE FN=0.D0 NORS='N' ENDIF FE=500000.0D0 SF=0.9996D0 OR=0.0D0 FOUND=0 CALL TCONPC (SF,OR,EPS,R,SO,V0,V2,V4,V6,ER,ESQ,RF) * * COMPUTE THE LATITUDES AND LONGITUDES * 200 CALL TMGEOD (NORTH,EAST,LAT,LON,EPS,CM,FE,SF,SO,R,V0,V2, & V4,V6,FN,ER,ESQ,CONV,KP) * * WRITE THE ZONE NUMBER * IF (IZ.GT.9) THEN WRITE(ZONE,600) IZ 600 FORMAT(1X,I2) ELSE WRITE(ZONE,605) IZ 605 FORMAT(1X,I2.2) ENDIF CALL DATAGP(CARDR,LAT,LON,FILFLAG,FOUND,FILPRT,ZONE,CONV,KP,NORS, * OLAT, OLON) OCOUNT = OCOUNT + 1 IF (OCOUNT .EQ. 1) THEN CALL PRINT_US2GP_HEADER() ENDIF WRITE(6, 2000) USNG_VALUE, OLAT, OLON, IN_DATUM 2000 FORMAT(A23, ' -> ', A13,' ',A14,' ', 'NAD ',A2) IF (OUTUNIT .GT. 6) THEN WRITE(OUTUNIT, 2001) USNG_VALUE, OLAT, OLON, IN_DATUM 2001 FORMAT(A23, ' -> ', A13,' ',A14,' ', 'NAD ',A2) ENDIF RETURN END ********************************************************************* SUBROUTINE DATAGP(CARDR,LAT,LON,FILFLAG,J,FILPRT,ZONE,CONV,KP, & NORS, OLAT, OLON) ********************************************************************** *** IMPLICIT REAL*8 (A-H,O-Z) LOGICAL FILFLAG,FILPRT REAL*8 NORTH,EAST,KP,LAT,LON,CONV INTEGER*4 IDEG,IMIN CHARACTER*1 NORS CHARACTER*1 WORE CHARACTER*1 DATNUM CHARACTER*4 ZONE CHARACTER*80 CARDR CHARACTER*13 OLAT CHARACTER*14 OLON INTEGER PCODE INTEGER FOUND_DEC COMMON/XY/NORTH,EAST COMMON/DATUM/DATNUM COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/PRECISE/PCODE R360 = 360.D0/RAD C WORE = 'W' READ(CARDR,FMT='(T69,A1)') WORE IF((WORE.EQ.'E').OR.(WORE.EQ.'e')) THEN LON = R360 - LON ENDIF CALL TODMS(LAT,LD,LM,SLAT) CALL TODMS(LON,LOD,LOM,SLON) CALL TODMS(DABS(CONV),IDEG,IMIN,CSEC) IF ((PCODE .EQ. 3).OR.(PCODE .EQ. 4)) THEN WRITE(OLAT, 2000) NORS, LD, LM, SLAT 2000 FORMAT(A1,I2.2, I2.2, F3.0) WRITE(OLON, 2001) LOD, LOM, SLON 2001 FORMAT('W',I3.3, I2.2, F3.0) ELSEIF (PCODE .EQ. 2) THEN WRITE(OLAT, 2005) NORS, LD, LM, SLAT 2005 FORMAT(A1,I2.2, I2.2, F4.1) WRITE(OLON, 2006) LOD, LOM, SLON 2006 FORMAT('W',I3.3, I2.2, F4.1) ELSE WRITE(OLAT, 2010) NORS, LD, LM, SLAT 2010 FORMAT(A1,I2.2, I2.2, F5.2) WRITE(OLON, 2011) LOD, LOM, SLON 2011 FORMAT('W',I3.3, I2.2, F5.2) ENDIF FOUND_DEC = 0 DO I=1, LEN(OLAT) IF ((FOUND_DEC .EQ. 0).AND.(OLAT(I:I) .EQ. ' ')) THEN OLAT(I:I) = '0' ELSEIF (OLAT(I:I) .EQ. '.') THEN FOUND_DEC = 1 ENDIF END DO FOUND_DEC = 0 DO I=1, LEN(OLON) IF ((FOUND_DEC .EQ. 0).AND.(OLON(I:I) .EQ. ' ')) THEN OLON(I:I) = '0' ELSEIF (OLON(I:I) .EQ. '.') THEN FOUND_DEC = 1 ENDIF END DO CALL CHK_OLAT(OLAT) CALL CHK_OLON(OLON) RETURN END ********************************************************************* SUBROUTINE TODMS(RAD,IDG,MIN,SEC) ********************************************************************* C RADIANS TO DEGREES,MINUTES AND SECONDS C REAL*8 RAD,SEC,RHOSEC DATA RHOSEC/2.062648062471D05/ SEC=RAD*RHOSEC IDG=SEC/3600.D0 SEC=SEC-DBLE(IDG*3600) MIN=SEC/60.D0 SEC=SEC-DBLE(MIN*60) IF((60.D0-DABS(SEC)).GT.5.D-6) GO TO 100 SEC=SEC-DSIGN(60.D0,SEC) MIN=MIN+ISIGN(1,MIN) 100 IF(IABS(MIN).LT.60) GO TO 101 MIN=MIN-ISIGN(60,MIN) IDG=IDG+ISIGN(1,IDG) 101 MIN=IABS(MIN) SEC=DABS(SEC) IF(RAD.GE.0.D0) GO TO 102 IF(IDG.EQ.0) MIN=-MIN IF(IDG.EQ.0.AND.MIN.EQ.0)SEC=-SEC 102 RETURN END ******************************************************************** ******************************************************************* ********************************************************************** SUBROUTINE CHGDEC (NNN,MMM,SS,CHAR) C ----------------------------------------------------- CHARACTER*1 DASH,ZERO,DOL,BLK1,CHAR(*),IB(20),TT CHARACTER*20 JB INTEGER*4 IDG,MIN REAL*8 S,SS,DEC,W,SEC,TEN,TOL C EQUIVALENCE (IB,JB) C DATA BLK1,DOL,ZERO,DASH/' ','$','0','-'/,TEN/10.0D0/ C C CHAR 1-16 LENGTH OF CHARACTER ARRAY FIELD C NR 1-13 LENGTH OF FIELD TO BE USED FROM LEFT C NP 7-13 LOCATION OF DECIMAL POINT FROM RIGHT C NEG CHAR(1) PUT THE MINUS SIGN HERE C EXAMPLE: RANGE NR=14 AND WITH A POINT NP=6 C W= -3600.376541 C C CHAR FIELD 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 C BLANK FILLED B B B B B B B B B B B B B B B B BLANK C NR ADDED |B B B B B B B B B B B B B B| B B LIMIT C NP ADDED |B B B B B B B . B B B B B B| B B POINT C W ENTERED |B B - 3 6 0 0 . 3 7 6 5 4 1| 8 4 C DOL (SIGN) |B B - 3 6 0 0 . 3 7 6 5 4 1| $ 4 C NEG (SIGN) |B B - 3 6 0 0 . 3 7 6 5 4 1| $ 4 C C ------------------------------- C SETUP OUTSIDE CONSTANTS C NR=NNN NP=MMM DEC=SS C C CHECK TO SEE IF NR AND NF ARE WITHIN LIMITS C NNP=IABS(NP) KTST=0 M=IABS(NR) C C IF NR IS GREATER THAN ZERO -- DECIMAL NUMBER C IF(NR.GT.0)GOTO 1 C ------------------------------------- C THIS IS A DEG-MIN-SEC FORMAT C IF(NNP.GT.5)NNP=5 C C ENTRY IS DDD-MM-SS.SS C OR HH-MM-SS.SS C DEC=DABS(DEC) SEC=DEC*3600.0D0 IDG=SEC/3600.0D0 SEC=SEC-DBLE(FLOAT(IDG))*3600.0D0 MIN=SEC/60.0D0 SEC=SEC-DBLE(FLOAT(MIN))*60.0D0 DEC=DBLE(FLOAT(IDG))*1000000.0D0+DBLE(FLOAT(MIN))*1000.0D0+SEC C KTST=1 M=15 C C ROUND THE DECIMAL NUMBER C 1 IP=-1*(NNP+1) TOL=5.0D0*(TEN**IP) W=(DABS(DEC) + TOL) IF(DEC.LT.0.0D0) W=-W DEC=W C N=NNP IF(M.GT.15)M=15 IF(M.LT.4)M=4 IF(N.GE.M)N=M-1 W=DEC C C CONVERT THE DECIMAL NUMBER C WRITE(JB,100) W 100 FORMAT(F20.10) C C BLANK FILL THE ARRAY C DO 5 IQ=1,16 5 CHAR(IQ)=BLK1 C C LOOK FOR THE FIRST NON-BLANK CHARACTER C DO 6 I=1,10 TT=IB(I) 6 IF(TT.NE.BLK1)GOTO 7 C C COMPUTE THE PROPER NUMBER LENGTH WITH THE PROPER DECIMALS C 7 K=11-I K=K+N IF(K.GT.M)M=K L=10-(M-N) IF(L.LE.0)L=0 C MM=M+1 J=0 IDEC=0 DO 30 I=1,MM K=I+L IF(K.GT.20)K=20 J=J+1 CHAR(J)=IB(K) 30 IF(CHAR(J).EQ.'.') IDEC=J N=IDEC+NNP+1 C IF(KTST.EQ.0)GOTO 40 C C FILL-OUT THE DEG-MIN-SEC FIELD C CHAR(4)=DASH CHAR(7)=DASH C C ZERO-OUT ALL BLANK CHARACTERS C DO 10 I=1,9 10 IF(CHAR(I).EQ.BLK1)CHAR(I)=ZERO C 40 CHAR(N)=DOL C RETURN END ************************************************************ SUBROUTINE DATUMM(ER,RF,F,ESQ,DATNUM) CHARACTER*1 DATNUM REAL*8 ER REAL*8 RF REAL*8 F REAL*8 ESQ ** ** FIND THE RIGHT SEMI MAJOR AXIS AND FLATTING ** IF(DATNUM.EQ.'1') THEN ** FOR THE NAD 27 DATUM ** ER=6378206.4D0 RF=294.978698D0 F=1.D0/RF ESQ=(F+F-F*F) ELSEIF (DATNUM.EQ.'2') THEN ** FOR THE NAD83 DATUM ** ER=6378137.D0 RF=298.257222101D0 F=1.D0/RF ESQ=(F+F-F*F) ELSE PRINT *, ' INVALID DATUM ' PRINT *, ' ' STOP ENDIF RETURN END C------------------------------------------------------------ SUBROUTINE PRINT_CMD_USAGE() CALL PRINT_VERSION() WRITE(6, 10) 10 FORMAT (/, *' USAGE - usng ', /, *' EXAMPLE - usng GP2US N385322.08 W0770206.86 83', /, *' usng US2GP 18SUJ2348306480 83', /, *' usng UTM2US 18 323483.2 4306479.5 83 A10',/, *' ', /, *' USAGE2 - usng FILE=', /, *' EXAMPLE - usng GP2US FILE=pos.in', /, * ' ',/, * ' TYPE "usng FILE" to see file input formats.', /, * ' (SEE FILE "usng.doc" FOR MORE INFORMATION).',/, *' ') RETURN END C----------------------------------------------------------- SUBROUTINE PRINT_FILE_USAGE() CALL PRINT_VERSION() WRITE(6, 10) 10 FORMAT (/, *' USAGE2 - usng FILE=', /, *' EXAMPLE - usng GP2US FILE=pos.in', /, * ' ', /, * ' ** NOTE - THE INPUT FILE IS A "DELIMITED" FORMAT.', /, * ' ', /, * ' For "GP2US" the input file format is:',/, * ' -------------------------------------',/, * ' FIELD1 = ANYTHING (Recommend PID)',/, * ' FIELD2 = ANYTHING (Recommend STATION NAME)', /, * ' FIELD3 = LATITUDE', /, * ' FIELD4 = LONGITUDE', /, * ' FIELD5 = DATUM', /, * ' ',/, * ' EXAMPLE:', /, * ' ZZ0298|BENCH MARK 283|N384558.68|W0775532.78|83',/, * ' ZZ0299|BENCH MARK 284|N384558.7|W0775532.8|83',/, * ' ZZ0300|BENCH MARK 7|N384559|W0775533|83',/, * ' ') CALL WAIT_FOR_INPUT() WRITE(6, 11) 11 FORMAT(/, * ' For "US2GP" the input file format is:',/, * ' -------------------------------------',/, * ' FIELD1 = ANYTHING (Recommend PID)',/, * ' FIELD2 = ANYTHING (Recommend STATION NAME)', /, * ' FIELD3 = USNG', /, * ' FIELD4 = DATUM', /, * ' ',/, * ' EXAMPLE:', /, * ' ZZ0401|BENCH MARK 283|18SUJ2348306479|83',/, * ' ZZ0402|BENCH MARK 284|18SUJ23480647|83',/, * ' ZZ0403|BENCH MARK 55|18SUJ234064|83',/, * ' ', * ' ',/, * ' (SEE FILE "usng.doc" FOR MORE INFORMATION).',/, *' ') RETURN END C------------------------------------------------------------ SUBROUTINE BUILD_80(LAT, LON, CARDR) CHARACTER*13 LAT CHARACTER*14 LON CHARACTER*80 CARDR CARDR(7:10) = '*80*' CARDR(45:50) = LAT(2:7) CARDR(51:55) = LAT(9:13) CARDR(56:56) = LAT(1:1) CARDR(57:63) = LON(2:8) CARDR(64:68) = LON(10:14) CARDR(69:69) = LON(1:1) C ---------------------------------------------------- C Change any spaces in Lat, Lon to 0 C ---------------------------------------------------- DO I= 45, 55 IF (CARDR(I:I) .EQ. ' ') CARDR(I:I) = '0' END DO DO I= 57, 68 IF (CARDR(I:I) .EQ. ' ') CARDR(I:I) = '0' END DO RETURN END C----------------------------------------------------------- SUBROUTINE CHK_OLAT(OLAT) CHARACTER*13 OLAT CHARACTER*1 DIR INTEGER DEG, MIN, SEC CHARACTER*3 DEC READ(OLAT, 20) DIR, DEG, MIN, SEC, DEC 20 FORMAT(A1, I2.2, I2.2, I2.2, A3) C ------------------------------------------------------ C Change any type N384560.0 to N384600.0 C ------------------------------------------------------ IF (SEC .LT. 60) THEN RETURN ENDIF IF (SEC .GE. 60) THEN SEC = SEC - 60 MIN = MIN + 1 ENDIF IF (MIN .GE. 60) THEN MIN = MIN - 60 DEG = DEG + 1 ENDIF IF (DEG .GE. 90) THEN DEG = DEG - 90 ENDIF OLAT = ' ' WRITE(OLAT, 100) DIR, DEG, MIN, SEC, DEC 100 FORMAT(A1, I2.2, I2.2, I2.2, A3) RETURN END C-------------------------------------------------------- SUBROUTINE CHK_OLON(OLON) CHARACTER*13 OLON CHARACTER*1 DIR INTEGER DEG, MIN, SEC CHARACTER*3 DEC C ------------------------------------------------------ C Change any type W0774560.0 to W0774600.0 C ------------------------------------------------------ READ(OLON, 20) DIR, DEG, MIN, SEC, DEC 20 FORMAT(A1, I3.3, I2.2, I2.2, A3) IF (SEC .LT. 60) THEN RETURN ENDIF IF (SEC .GE. 60) THEN SEC = SEC - 60 MIN = MIN + 1 ENDIF IF (MIN .GE. 60) THEN MIN = MIN - 60 DEG = DEG + 1 ENDIF IF (DEG .GE. 360) THEN DEG = DEG - 360 ENDIF OLON = ' ' WRITE(OLON, 100) DIR, DEG, MIN, SEC, DEC 100 FORMAT(A1, I3.3, I2.2, I2.2, A3) RETURN END C-------------------------------------------------------------- SUBROUTINE CHECK_IN_LAT(LAT, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*13 LAT CHARACTER*13 LAT_NUM CHARACTER*1 DIR INTEGER DEG, MIN, SEC INTEGER ERROR ERROR=0 C ------------------------------------ C Convert to upper case for user C ------------------------------------ IF (LAT(1:1) .EQ. 'n') THEN LAT(1:1) = 'N' ELSEIF (LAT(1:1) .EQ. 's') THEN LAT(1:1) = 'S' ENDIF DEG=-1 MIN=-1 SEC=-1 READ (LAT, 100, ERR=200) DIR, DEG, MIN, SEC 100 FORMAT(A1, I2, I2, I2) GOTO 210 200 ERROR = 1 210 CONTINUE IF ((LAT(1:1) .NE. 'N').AND.(LAT(1:1) .NE. 'S')) THEN ERROR = 1 ENDIF J=0 DO I=2,13 IF ((LAT(I:I) .NE. ' ').AND.(LAT(I:I) .NE. '.')) THEN J = J + 1 LAT_NUM(J:J) = LAT(I:I) ENDIF END DO IF (J .LT. 6) THEN ERROR = 1 ENDIF DO I=1, J IF (INDEX('0123456789', LAT_NUM(I:I)) .LE. 0) THEN ERROR = 1 ENDIF END DO IF ((DEG .LT. 0).OR.(DEG .GT. 90)) THEN ERROR = 1 ELSEIF ((MIN .LT. 0).OR.(MIN .GT. 59)) THEN ERROR = 1 ELSEIF ((SEC .LT. 0).OR.(SEC .GT. 59)) THEN ERROR = 1 ENDIF IF (ERROR .EQ. 1) THEN WRITE(6, 1000) LAT 1000 FORMAT (/, * 'ERROR - Invalid Format for Latitude ', A13,/, * ' Format should be Hemisphere-Deg-Min-Sec',/, * ' for a length of char 7 or greater.', /, * ' EXAMPLE1 = N385945.01',/, * ' EXAMPLE2 = S095945.01',/, * ' ',/, * 'with value between:',/, * ' 0 < Deg < 90',/, * ' 0 < Min < 60',/, * ' 0 < Sec < 60',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C-------------------------------------------------------------- SUBROUTINE CHECK_IN_LON(LON, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*14 LON CHARACTER*14 LON_NUM CHARACTER*1 DIR INTEGER DEG, MIN, SEC INTEGER ERROR ERROR=0 C ------------------------------------ C Convert to upper case for user C ------------------------------------ IF (LON(1:1) .EQ. 'w') THEN LON(1:1) = 'W' ELSEIF (LON(1:1) .EQ. 'e') THEN LON(1:1) = 'E' ENDIF DEG=-1 MIN=-1 SEC=-1 READ (LON, 100, ERR=200) DIR, DEG, MIN, SEC 100 FORMAT(A1, I3, I2, I2) GOTO 210 200 ERROR = 1 210 CONTINUE IF ((LON(1:1) .NE. 'W').AND.(LON(1:1) .NE. 'E')) THEN ERROR = 1 ENDIF J=0 DO I=2,14 IF ((LON(I:I) .NE. ' ').AND.(LON(I:I) .NE. '.')) THEN J = J + 1 LON_NUM(J:J) = LON(I:I) ENDIF END DO IF (J .LT. 7) THEN ERROR = 1 ENDIF DO I=1, J IF (INDEX('0123456789', LON_NUM(I:I)) .LE. 0) THEN ERROR = 1 ENDIF END DO IF ((DEG .LT. 0).OR.(DEG .GT. 360)) THEN ERROR = 1 ELSEIF ((MIN .LT. 0).OR.(MIN .GT. 59)) THEN ERROR = 1 ELSEIF ((SEC .LT. 0).OR.(SEC .GT. 59)) THEN ERROR = 1 ENDIF IF (ERROR .EQ. 1) THEN WRITE(6, 1000) LON 1000 FORMAT (/, * 'ERROR - Invalid Format for Longitude ', A14,/, * ' Format should be Hemisphere-Deg-Min-Sec',/, * ' for a length of char 8 or greater.', /, * ' EXAMPLE1 = W0775945.01',/, * ' EXAMPLE2 = E2835945.01',/, * ' ',/, * 'with value between:',/, * ' 0 < Deg < 360',/, * ' 0 < Min < 60',/, * ' 0 < Sec < 60',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C---------------------------------------------------------- SUBROUTINE CHECK_IN_DATUM(DATUM, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*2 DATUM INTEGER ERROR ERROR=0 IF ((DATUM .NE. '83').AND.(DATUM .NE. '27')) THEN ERROR=1 ENDIF IF (ERROR .EQ. 1) THEN WRITE(6, 1000) DATUM 1000 FORMAT (/, * 'ERROR - Invalid Format for Datum ', A2,/, * ' Format should be 83 or 27',/, * ' for a length of char 2.', /, * ' EXAMPLE1 = 83',/, * ' EXAMPLE2 = 27',/, * ' ',/, * 'where:',/, * ' 83 means Datum=NAD 83',/, * ' 27 means Datum=NAD 27',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C-------------------------------------------------------------- SUBROUTINE CHECK_IN_EAST(EAST, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*14 EAST CHARACTER*14 EAST_NUM REAL*8 R_EAST INTEGER FOUND_DEC INTEGER FOUND_COMMA INTEGER ERROR ERROR=0 FOUND_DEC = 0 FOUND_COMMA = 0 READ (EAST, 100, ERR=200) R_EAST 100 FORMAT(F11.3) GOTO 210 200 ERROR = 1 210 CONTINUE J=0 DO I=1,13 IF (EAST(I:I) .EQ. ',') THEN FOUND_COMMA = 1 ELSE IF (EAST(I:I) .EQ. '.') THEN FOUND_DEC=1 ELSE IF ((EAST(I:I) .NE. ' ').AND.(EAST(I:I) .NE. '.')) THEN J = J + 1 EAST_NUM(J:J) = EAST(I:I) ENDIF END DO IF (J .LT. 2) THEN ERROR = 1 ENDIF DO I=1, J IF (INDEX('0123456789', EAST_NUM(I:I)) .LE. 0) THEN ERROR = 1 ENDIF END DO IF ((R_EAST .LT. 0).OR.(R_EAST .GT. 900000)) THEN ERROR = 1 ENDIF IF (FOUND_DEC .EQ. 0) THEN ERROR = 1 ENDIF IF (FOUND_COMMA .NE. 0) THEN ERROR = 1 ENDIF IF (ERROR .EQ. 1) THEN IF (FOUND_DEC .EQ. 0) THEN WRITE (6, 998) 998 FORMAT(/,'*** DECIMAL POINT IS MISSING') ENDIF IF (FOUND_COMMA .NE. 0) THEN WRITE (6, 999) 999 FORMAT(/,'*** COMMAS SHOULD NOT BE INCLUDED') ENDIF WRITE(6, 1000) EAST 1000 FORMAT (/, * 'ERROR - Invalid Format for East ', A14,/, * ' Format should be a decimal number',/, * ' for a length of char 3 or greater.', /, * ' with a decimal point',/, * ' EXAMPLE1 = 323483.193',/, * ' ',/, * 'with value between:',/, * ' 0.0 < East < 900000.000',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C-------------------------------------------------------------- SUBROUTINE CHECK_IN_NORTH(NORTH, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*14 NORTH CHARACTER*14 NORTH_NUM REAL*8 R_NORTH INTEGER FOUND_DEC INTEGER FOUND_COMMA INTEGER ERROR ERROR=0 FOUND_DEC = 0 FOUND_COMMA = 0 READ (NORTH, 100, ERR=200) R_NORTH 100 FORMAT(F11.3) GOTO 210 200 ERROR = 1 210 CONTINUE J=0 DO I=1,14 IF (NORTH(I:I) .EQ. ',') THEN FOUND_COMMA = 1 ELSE IF (NORTH(I:I) .EQ. '.') THEN FOUND_DEC=1 ELSE IF ((NORTH(I:I) .NE. ' ').AND.(NORTH(I:I) .NE. '.')) THEN J = J + 1 NORTH_NUM(J:J) = NORTH(I:I) ENDIF END DO IF (J .LT. 2) THEN ERROR = 1 ENDIF DO I=1, J IF (INDEX('-0123456789', NORTH_NUM(I:I)) .LE. 0) THEN ERROR = 1 ENDIF END DO IF ((R_NORTH .LT. -10000000).OR.(R_NORTH .GT. 9999999)) THEN ERROR = 1 ENDIF IF (FOUND_DEC .EQ. 0) THEN ERROR = 1 ENDIF IF (FOUND_COMMA .NE. 0) THEN ERROR = 1 ENDIF IF (ERROR .EQ. 1) THEN IF (FOUND_DEC .EQ. 0) THEN WRITE (6, 998) 998 FORMAT(/,'*** DECIMAL POINT IS MISSING') ENDIF IF (FOUND_COMMA .NE. 0) THEN WRITE (6, 999) 999 FORMAT(/,'*** COMMAS SHOULD NOT BE INCLUDED') ENDIF WRITE(6, 1000) NORTH 1000 FORMAT (/, * 'ERROR - Invalid Format for North ', A14,/, * ' Format should be a decimal number',/, * ' for a length of char 3 or greater.', /, * ' with a decimal point',/, * ' EXAMPLE1 = 4306479.542',/, * ' ',/, * 'with value between:',/, * ' -9999999.999 < North < 9999999.999',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C-------------------------------------------------------------- SUBROUTINE CHECK_IN_ZONE(ZONE, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*2 ZONE CHARACTER*2 ZONE_NUM INTEGER IZONE INTEGER ERROR ERROR=0 READ (ZONE, 100, ERR=200) IZONE 100 FORMAT(I2) GOTO 210 200 ERROR = 1 210 CONTINUE J=0 DO I=1,2 IF ((ZONE(I:I) .NE. ' ').AND.(ZONE(I:I) .NE. '.')) THEN J = J + 1 ZONE_NUM(J:J) = ZONE(I:I) ENDIF END DO DO I=1, J IF (INDEX('0123456789', ZONE_NUM(I:I)) .LE. 0) THEN ERROR = 1 ENDIF END DO IF ((IZONE .LT. 1).OR.(IZONE .GT. 60)) THEN ERROR = 1 ENDIF IF (ERROR .EQ. 1) THEN WRITE(6, 1000) ZONE 1000 FORMAT (/, * 'ERROR - Invalid Format for Zone ', A2,/, * ' Format should be a integer',/, * ' EXAMPLE1 = 18',/, * ' ',/, * 'with value between:',/, * ' 1 < Zone < 60',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C---------------------------------------------------------- SUBROUTINE CHECK_IN_ACCURACY(ACCURACY, LINE_COUNT, ERR_STATUS) INTEGER LINE_COUNT INTEGER ERR_STATUS CHARACTER*4 ACCURACY INTEGER ERROR ERROR=0 IF ((ACCURACY .NE. 'A1').AND. * (ACCURACY .NE. 'A10').AND. * (ACCURACY .NE. 'A100')) THEN ERROR=1 ENDIF IF (ERROR .EQ. 1) THEN WRITE(6, 1000) LON 1000 FORMAT (/, * 'ERROR - Invalid Format for UTM Accuracy ', A4,/, * ' Format should be',/, * ' A1, A10, or A100',/, * ' for a length not greater than char 4.', /, * ' ',/, * 'where:',/, * ' A1 means good to 1 meter',/, * ' A10 means good to 10 meters',/, * ' A100 means good to 10 meters',/) IF (LINE_COUNT .GT. 0) THEN WRITE(6, 1010) LINE_COUNT 1010 FORMAT(/, 'Point of error = line ', I6, //) ENDIF IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF RETURN END C--------------------------------------------------------- SUBROUTINE CHECK_IN_USNG(USNG_VALUE, LINE_COUNT, ERR_STATUS) CHARACTER*25 USNG_VALUE INTEGER LINE_COUNT INTEGER ERR_STATUS INTEGER IZ CHARACTER*2 GRID_ID INTEGER NORTH_1M, EAST_1M CHARACTER*10 CEAST_NORTH CHARACTER*5 CEAST_1M CHARACTER*5 CNORTH_1M INTEGER COUNT INTEGER IS_NUMBER INTEGER LOC CHARACTER*1 GZD INTEGER PCODE COMMON/PRECISE/PCODE IF (INDEX('0123456789', USNG_VALUE(2:2)) .GT. 0) THEN READ(USNG_VALUE, 100, ERR=999) IZ,GZD, GRID_ID 100 FORMAT(I2, A1, A2) ELSE READ(USNG_VALUE, 110, ERR=999) IZ,GZD, GRID_ID 110 FORMAT(I1, A1, A2) ENDIF C C------------------------------------------------------ C C Southern Hemisphere enabled 9/25/02 C C------------------------------------------------------ C IF (INDEX("NPQRSTUVWX", GZD) .LE. 0) THEN C WRITE(6, 120) GZD C 120 FORMAT(/, C * 'ERROR - Grid Zone Designation Letter = ', A2,/, C * ' will not compute',/, C * ' Current version of this program only works for', /, C * ' the Northern Hemisphere', /, C * ' GZD letters N thru X', /) C IF (ERR_STATUS .NE. 9) THEN C STOP C ELSE C ERR_STATUS = 1 C RETURN C ENDIF C ENDIF COUNT=1 LOC=5 IS_NUMBER=1 DO WHILE(IS_NUMBER .EQ. 1) IF ((USNG_VALUE(LOC:LOC) .EQ. '') .OR. * (USNG_VALUE(LOC:LOC) .EQ. ' ') .OR. * (USNG_VALUE(LOC:LOC) .EQ. '(')) THEN IS_NUMBER=0 ELSE IF (INDEX('0123456789', USNG_VALUE(LOC:LOC)) .GT. 0) THEN CEAST_NORTH(COUNT:COUNT) = USNG_VALUE(LOC:LOC) COUNT = COUNT + 1 LOC = LOC + 1 ELSE LOC = LOC + 1 ENDIF END DO COUNT = COUNT - 1 IF (COUNT .EQ. 10) THEN CEAST_1M(1:5) = CEAST_NORTH(1:5) CNORTH_1M(1:5) = CEAST_NORTH(6:10) READ(CEAST_1M, 200, ERR=999) EAST_1M READ(CNORTH_1M, 200, ERR=999) NORTH_1M 200 FORMAT(I5.5) PCODE = 1 ELSEIF (COUNT .EQ. 8) THEN CEAST_1M(1:4) = CEAST_NORTH(1:4) CEAST_1M(5:5) = '0' CNORTH_1M(1:4) = CEAST_NORTH(5:8) CNORTH_1M(5:5) = '0' READ(CEAST_1M, 201, ERR=999) EAST_1M READ(CNORTH_1M, 201, ERR=999) NORTH_1M 201 FORMAT(I5.5) PCODE = 2 ELSEIF (COUNT .EQ. 6) THEN CEAST_1M(1:3) = CEAST_NORTH(1:3) CEAST_1M(4:5) = '00' CNORTH_1M(1:3) = CEAST_NORTH(4:6) CNORTH_1M(4:5) = '00' READ(CEAST_1M, 202, ERR=999) EAST_1M READ(CNORTH_1M, 202, ERR=999) NORTH_1M 202 FORMAT(I5.5) PCODE = 3 ELSEIF (COUNT .EQ. 4) THEN WRITE(6, 300) CEAST_NORTH IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF CEAST_1M(1:2) = CEAST_NORTH(1:2) CEAST_1M(3:5) = '000' CNORTH_1M(1:2) = CEAST_NORTH(3:4) CNORTH_1M(3:5) = '000' READ(CEAST_1M, 204, ERR=999) EAST_1M READ(CNORTH_1M, 204, ERR=999) NORTH_1M 204 FORMAT(I5.5) PCODE = 4 ELSE WRITE(6, 300) CEAST_NORTH 300 FORMAT('ERROR - Cannot convert coordinate ', A10,/, * ' Number of digits must be 10, 8, or 6', /, * ' EXAMPLE:', /, *' 18SUJ2348306479 ( 1-meter accuracy)',/, *' 18SUJ23430647 ( 10-meter accuracy)',/, *' 18SUJ234064 (100-meter accuracy)',/) IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF ENDIF GOTO 1000 999 IF (ERR_STATUS .NE. 9) THEN STOP ELSE ERR_STATUS = 1 RETURN ENDIF 1000 CONTINUE RETURN END C------------------------------------------- SUBROUTINE TO_UPPER(MYTEXT, VSIZE) CHARACTER*80 MYTEXT INTEGER VSIZE INTEGER IVALUE IF ((VSIZE .LT. 0) .OR. (VSIZE .GT. 80)) THEN RETURN ENDIF DO I=1,VSIZE IVALUE = ICHAR(MYTEXT(I:I)) IF ((IVALUE .GE. 97).AND.(IVALUE .LE. 122)) THEN IVALUE = IVALUE - 32 MYTEXT(I:I) = CHAR(IVALUE) ENDIF END DO RETURN END C--------------------------------------------------------------- SUBROUTINE X_DLAT(IN_ZONE, DLAT) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*2 IN_ZONE REAL*8 DLAT CHARACTER*1 NORS REAL*8 LAT,LON,LCM,UCM REAL*8 KP REAL*8 NORTH,EAST REAL*8 ABS_NORTH COMMON/CONST/RAD,ER,RF,ESQ,PI COMMON/XY/NORTH,EAST ABS_NORTH = DABS(NORTH) IF (NORTH .LT. 0.0) THEN ABS_NORTH = 20000000.D0 - ABS_NORTH ENDIF READ(IN_ZONE, 20) ICODE 20 FORMAT(I2) * * FIND THE CENTRAL MERIDAIN IF THE ZONE NUMBER IS LESS THAN 30 * IF(ICODE.LT.30) THEN IZ=ICODE ICM=(183-(6*IZ)) CM=DBLE(FLOAT(ICM))/RAD UCM=(ICM+3)/RAD LCM=(ICM-3)/RAD ENDIF * * FIND THE CENTRAL MERIDAN IF THE ZONE NUMBER IS LARGER THAN 30 * IF(ICODE.GE.30) THEN IZ=ICODE ICM=(543 - (6*IZ)) CM= DBLE(FLOAT(ICM))/RAD UCM=(ICM+3)/RAD LCM=(ICM-3)/RAD ENDIF TOL=(5.0D0/60.0D0)/RAD IF(ABS_NORTH.GE.10000000.0) THEN FN= 10000000.0D0 NORS= 'S' ELSE FN=0.D0 NORS='N' ENDIF FE=500000.0D0 SF=0.9996D0 OR=0.0D0 CALL TCONPC (SF,OR,EPS,R,SO,V0,V2,V4,V6,ER,ESQ,RF) * * COMPUTE THE LATITUDES AND LONGITUDES * 200 CALL TMGEOD (ABS_NORTH,EAST,LAT,LON,EPS,CM,FE,SF,SO,R,V0,V2, & V4,V6,FN,ER,ESQ,CONV,KP) DLAT = LAT * RAD IF (NORS .EQ. 'S') THEN DLAT = -1.D0 * DLAT ENDIF RETURN END