From 00251cf0a3e9f80622013cbcbdbb96fc758bff82 Mon Sep 17 00:00:00 2001 From: SriReddy-NOAA Date: Fri, 9 Apr 2021 00:20:37 -0400 Subject: version 3.3.0 --- htdp.f | 921 ++++++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 599 insertions(+), 322 deletions(-) (limited to 'htdp.f') diff --git a/htdp.f b/htdp.f index e793cee..910f8e2 100644 --- a/htdp.f +++ b/htdp.f @@ -24,8 +24,8 @@ C You must change HTDP version and date here, if necessary - HTDP_version = '3.2.9' - Version_date = 'May 17, 2020' + HTDP_version = '3.3.0' + Version_date = 'January 15, 2021' *** Introduce variables for file id's @@ -73,8 +73,9 @@ C 1 FORM='UNFORMATTED') WRITE(LUOUT,501) 501 FORMAT( 1 ' AUTHORS: Richard Snay, Chris Pearson, Jarir Saleh, '/ - 1 ' and Michael Dennis '/ - 1 ' Email: ngs.cors.htdp@noaa.gov '/ + 1 ' and Michael Dennis '// + 1 ' Web: https://geodesy.noaa.gov/TOOLS/Htdp/Htdp.shtml '/ + 1 ' Email: ngs.cors.htdp@noaa.gov '/ 1 '********************************************************'/) WRITE(LUOUT,10) 10 FORMAT( @@ -84,13 +85,13 @@ C 1 FORM='UNFORMATTED') WRITE(LUOUT,11) 11 FORMAT( 5 ' The User Guide contains additional information and a set'/ - 5 ' of exercises to familiarize users with the software.'// + 5 ' of exercises to familiarize users with the software.'/ 5 ' Hit ENTER or RETURN to continue. ') read(luin, '(a5)',err=51,iostat=ios) cont if (ios /= 0) goto 51 25 WRITE(LUOUT,26) - 26 FORMAT(' ***************************************'/ + 26 FORMAT('********************************************************'/ 1 ' MAIN MENU:',/ 6 ' 0... Exit software.',/ 7 ' 1... Estimate horizontal displacements between two dates.'/ @@ -129,7 +130,7 @@ C CLOSE(I5, STATUS = 'DELETE') write (*,*) 'You did not hit enter : ios = ',ios write (*,*) "ABNORMAL TERMINATION" STOP - + 52 write (*,'(/)') write (*,*) 'Problem with reading OPTION: ios = ',ios write (*,*) "ABNORMAL TERMINATION" @@ -175,7 +176,7 @@ C*** Set default reference epoch to Jan. 1, 2010 *** that correspond to the boundaries for the regions. *** Region 1 is the San Andreas fault in central California *** Region 2 is southern California -*** Region 3 is Northern California +*** Region 3 is northern California *** Region 4 is the Pacific Northwest *** Region 5 is western CONUS *** Region 6 is CONUS @@ -183,19 +184,38 @@ C*** Set default reference epoch to Jan. 1, 2010 *** Region 8 is south-central Alaska *** Region 9 is southeast Alaska *** Region 10 is All Mainland Alaska -*** Region 11 is the North American plate -*** Region 12 is the Caribbean plate -*** Region 13 is the Pacific plate -*** Region 14 is the Juan de Fuca plate -*** Region 15 is the Cocos plate -*** Region 16 is the Mariana plate -*** Region 17 is the Philippine Sea plate + +* Placeholder regions inserted to make compatible with new initbd.f +* that contains tectonic plates added in December 2020 (v3.3.0). +*** Region 11 is a placeholder region +*** Region 12 is a placeholder region +*** Region 13 is a placeholder region + +*** Region 14 is the North America plate +*** Region 15 is the Caribbean plate +*** Region 16 is the Pacific plate +*** Region 17 is the Juan de Fuca plate +*** Region 18 is the Cocos plate +*** Region 19 is the Mariana plate +*** Region 20 is the Philippine Sea plate + +* Plates added in December 2020 (v3.3.0). +*** Region 21 is the South America plate +*** Region 22 is the Nazca plate +*** Region 23 is the Panama plate +*** Region 24 is the North Andes plate +*** Region 25 is the Africa plate +*** Region 26 is the Eurasia plate +*** Region 27 is the Rivera plate +*** Region 28 is the Galapagos plate +*** Region 29 is the Tonga plate +*** Region 30 is the Niuafo'ou plate IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) - parameter (NMREGN = 17) + parameter (NMREGN = 30) COMMON /CONST/ A,F,E2,EPS,AF,PI,TWOPI,RHOSEC - COMMON /BNDRY/ X(4000), Y(4000), NPOINT(30) + COMMON /BNDRY/ X(7000), Y(7000), NPOINT(40) IEND = NPOINT(NMREGN + 1) - 1 DO 10 J = 1, IEND @@ -211,8 +231,8 @@ C*** Set default reference epoch to Jan. 1, 2010 IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) - parameter (NMREGN = 17) - COMMON /BNDRY/ X(4000), Y(4000), NPOINT(30) + parameter (NMREGN = 30) + COMMON /BNDRY/ X(7000), Y(7000), NPOINT(40) COMMON /FILES/ LUIN, LUOUT, I1, I2, I3, I4, I5, I6 COMMON /CONST/ A, F, E2, EPS, AF, PI, TWOPI, RHOSEC @@ -232,8 +252,11 @@ C*** Set default reference epoch to Jan. 1, 2010 RETURN END + ******************************************************** - SUBROUTINE POLYIN (X0,Y0,X,Y,N,NPC) + SUBROUTINE POLYIN (X0,Y0,X,Y,N,NPC) +* This version of subroutine is from TRANS4D v0.3.1 and replaces +* version in HTDP v3.2.9. C SUBROUTINE TO DETERMINE IF A POINT AT (X0,Y0) IS INSIDE OR C OUTSIDE OF A CLOSED FIGURE DESCRIBED BY A SEQUENCE OF CONNECTED C STRAIGHT LINE SEGMENTS WITH VERTICES AT X, Y. @@ -266,20 +289,18 @@ C MAINTAINS FULL ACCURACY OF INPUT COORDINATES. C IMPLICIT INTEGER*4 (I-N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) + COMMON /FILES/ LUIN, LUOUT, I1, I2, I3, I4, I5, I6 DIMENSION X(N),Y(N) - DATA I6/6/ + IS=0 NPC=0 C C FIND STARTING POINT WHERE X(I).NE.X0 IP=0 10 IP=IP+1 - IF ((X(IP)-X0) == 0) goto 12 - IF ((X(IP)-X0) < 0) goto 15 - IF ((X(IP)-X0) > 0) goto 16 -c IF (X(IP)-X0) 15,12,16 + IF(X(IP)-X0) 15,12,16 12 IF(IP.LE.N) GO TO 10 - WRITE(I6,6001) + WRITE(LUOUT,6001) 6001 FORMAT('0 POLYGON INPUT ERROR - ALL POINTS ON LINE X = X0') STOP 15 IL=-1 @@ -295,67 +316,29 @@ C DO 100 II=IP1,IPN I=II IF(I.GT.N) I=I-N - IF(IL == 0) goto 50 - IF(IL < 0) goto 30 - IF(IL > 0) goto 40 -c IF(IL) 30,50,40 - - 30 IF(X(I)-X0 == 0) goto 32 - IF(X(I)-X0 < 0) goto 90 - IF(X(I)-X0 > 0) goto 34 -c 30 IF(X(I)-X0) 90,32,34 + IF(IL) 30,50,40 + 30 IF(X(I)-X0) 90,32,34 32 IS=-1 GO TO 60 34 IL=1 GO TO 80 - 40 IF(X(I)-X0 == 0) goto 44 - IF(X(I)-X0 < 0) goto 42 - IF(X(I)-X0 > 0) goto 90 -c 40 IF(X(I)-X0) 42,44,90 + 40 IF(X(I)-X0) 42,44,90 42 IL=-1 GO TO 80 44 IS=1 GO TO 60 - 50 IF(X(I)-X0 == 0) goto 55 - IF(X(I)-X0 < 0) goto 52 - IF(X(I)-X0 > 0) goto 54 -c 50 IF(X(I)-X0) 52,55,54 + 50 IF(X(I)-X0) 52,55,54 52 IL=-1 - IF(IS == 0) goto 140 - IF(IS < 0) goto 90 - IF(IS > 0) goto 80 -c IF(IS) 90,140,80 + IF(IS) 90,140,80 54 IL=1 - IF(IS == 0) goto 140 - IF(IS < 0) goto 80 - IF(IS > 0) goto 90 -c IF(IS) 80,140,90 - - 55 IF(Y(I)-Y0 == 0) goto 120 - IF(Y(I)-Y0 < 0) goto 57 - IF(Y(I)-Y0 > 0) goto 58 -c 55 IF(Y(I)-Y0) 57,120,58 - - 57 IF(YL-Y0 == 0) goto 120 - IF(YL-Y0 < 0) goto 90 - IF(YL-Y0 > 0) goto 120 -c 57 IF(YL-Y0) 90,120,120 - - 58 IF(YL-Y0 == 0) goto 120 - IF(YL-Y0 < 0) goto 120 - IF(YL-Y0 > 0) goto 90 -c 58 IF(YL-Y0) 120,120,90 + IF(IS) 80,140,90 + 55 IF(Y(I)-Y0) 57,120,58 + 57 IF(YL-Y0) 90,120,120 + 58 IF(YL-Y0) 120,120,90 C 60 IL=0 - IF(Y(I)-Y0 == 0) goto 120 - IF(Y(I)-Y0 < 0) goto 90 - IF(Y(I)-Y0 > 0) goto 90 -c IF(Y(I)-Y0) 90,120,90 - - 80 IF(YL-Y0+(Y(I)-YL)*(X0-XL)/(X(I)-XL) == 0) goto 120 - IF(YL-Y0+(Y(I)-YL)*(X0-XL)/(X(I)-XL) < 0) goto 90 - IF(YL-Y0+(Y(I)-YL)*(X0-XL)/(X(I)-XL) > 0) goto 85 -c 80 IF(YL-Y0+(Y(I)-YL)*(X0-XL)/(X(I)-XL)) 90,120,85 + IF(Y(I)-Y0) 90,120,90 + 80 IF(YL-Y0+(Y(I)-YL)*(X0-XL)/(X(I)-XL)) 90,120,85 85 NPC=NPC+1 90 XL=X(I) YL=Y(I) @@ -364,11 +347,12 @@ c 80 IF(YL-Y0+(Y(I)-YL)*(X0-XL)/(X(I)-XL)) 90,120,85 RETURN 120 NPC=2 RETURN - 140 WRITE(I6,6002) + 140 WRITE(LUOUT,6002) 6002 FORMAT('0 POLYGON LOGIC ERROR - PROGRAM SHOULD NOT REACH THIS', . ' POINT') RETURN END + ***************************************************************** SUBROUTINE RADR8T (YLAT,VN,VE,VNR,VER) @@ -393,7 +377,7 @@ C Compute the NAD_83(CORS96) velocity at a point in mm/yr !Not anymore since IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) parameter (NUMGRD = 10) - parameter (NMREGN = 17) + parameter (NMREGN = 30) COMMON /CONST/ A,F,E2,EPS,AF,PI,TWOPI,RHOSEC COMMON /FILES/ LUIN, LUOUT, I1,I2,I3,I4,I5,I6 COMMON /VGRID/ B(210000) @@ -406,6 +390,10 @@ c1001 FORMAT( 'JREGN = ', I6) *** Use tectonic plate model to compute velocity relative *** to ITRF2008 IPLATE = JREGN - NUMGRD + +*** Subtract the number of placeholder regions + IPLATE = IPLATE - 3 + ELON = - YLON HT = 0.D0 CALL TOXYZ(YLAT, ELON, HT, X, Y, Z) @@ -470,33 +458,62 @@ c ENDIF **************************************** SUBROUTINE PLATVL(IPLATE, X, Y, Z, VX, VY, VZ) *** Compute the ITRF2008 velocity at point on plate = IPLATE -*** with coordinates X, Y, Z (in meters) -*** The resulting velocities--VX, VY, and VZ--will be in meters/yr -*** References -*** Altamimi et al. 2012 = JGR (Paper on ITRF2008 plate motion) -*** DeMets et al. 2010 = Geophysical Journal Int'l, vol 181, -*** Snay 2003 = SALIS, Vol 63, No 1 (Paper on Frames for Pacific) +*** with coordinates X, Y, Z (in meters) +*** The resulting velocities--VX, VY, and VZ--will be in meters/yr +*** References: +*** Altamimi et al. 2012 = JGR (Paper on ITRF2008 plate motion). +*** Altamimi et al. 2017. "ITRF2014 plate motion model." Geophys. J. Int, 209. +*** Bird 2003. "An updated digital model of plate boundaries." Geochem. Geophys. Geosyst., 4(3), 1027 +*** DeMets et al. 2010. "Geologically current plate motions." Geophys. J. Int., 181. +*** Kreemer et al. 2014. "A geodetic plate motion and global strain rate model." Geochem. Geophys. Geosyst., 15. +*** Mora-Páez et al. 2018. "Crustal deformation in the northern Andes – A new GPS velocity field." J. S. Amer. Earth Sci, 2018. +*** Snay 2003 = SALIS, Vol 63, No 1 (Paper on Frames for Pacific). IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER*4 (I-N) - DIMENSION WX(7), WY(7), WZ(7) + PARAMETER (NUMPLATE = 17) + + DIMENSION WX(NUMPLATE), WY(NUMPLATE), WZ(NUMPLATE) COMMON /FILES/ LUIN, LUOUT, I1, I2, I3, I4, I5, I6 + +*** All plate boundaries are from Bird 2003 with minor topological corrections. *** IPLATE = 1 --> North America (from Altamimi et al. 2012) *** 2 --> Caribbean (from Altamimi et al. 2012) *** 3 --> Pacific (from Altamimi et al. 2012) *** 4 --> Juan de Fuca (from DeMets et al. 2010) *** 5 --> Cocos (from DeMets et al. 2010) -*** 6 --> Mariana (from Snay, 2003) +*** 6 --> Mariana (from Snay 2003) *** 7 --> Philippine Sea (from DeMets et al. 2010) - DATA WX /0.170D-9, 0.238D-9, -1.993D-9, - 1 6.626D-9, -10.390D-9, -.097D-9, -0.841D-9/ - DATA WY /-3.209D-9, -5.275D-9, 5.023D-9, - 1 11.708D-9, -14.954D-9, .509D-9, 3.989D-9/ - DATA WZ /-0.485D-9, 3.219D-9,-10.501D-9, - 1 -10.615D-9, 9.148D-9,-1.682D-9, -10.626D-9/ - - IF (IPLATE .LE. 0 .OR. IPLATE .GT. 7) THEN +*** New plates added in December 2020 (v3.3.0): +*** 8 --> South America (from Altamimi et al. 2017) +*** 9 --> Nazca (from Altamimi et al. 2017) +*** 10 --> Panama (from Kreemer et al. 2014) +*** 11 --> North Andes (from Mora-Paez, 2018) +*** 12 --> Africa (from Altamimi et al. 2017) +*** 13 --> Eurasia (from Altamimi et al. 2017) +*** 14 --> Rivera (from DeMets et al. 2010) +*** 15 --> Galapagos (from Bird 2003) +*** 16 --> Tonga (from Kreemer et al. 2014) +*** 17 --> Niuafo'ou (from Bird 2003) + +* Plate rotation rates in radians per year. Existing plates (1-7) have ITRF2008 rates +* (except for plate 6, the Mariana plate referenced to ITRF2000). The new plates added +* in December 2020 have ITRF2014 rates (plates 8 thorugh 17). + DATA WX / 0.170D-9, 0.238D-9, -1.993D-9, 6.626D-9, -10.390D-9, + & -0.097D-9, -0.841D-9, -1.309D-9, -1.614D-9, 2.088D-9, + & -1.964D-9, 0.480D-9, -0.412D-9, -21.933D-9, 14.273D-9, + & 137.841D-9, -57.324D-9 / + DATA WY /-3.209D-9, -5.275D-9, 5.023D-9, 11.708D-9, -14.954D-9, + & 0.509D-9, 3.989D-9, -1.459D-9,-7.486D-9, -23.037D-9, + & -1.518D-9, -2.977D-9, -2.574D-9, -70.432D-9, 94.440D-9, + & 10.447D-9, -5.814D-9 / + DATA WZ /-0.485D-9, 3.219D-9,-10.501D-9, -10.615D-9, 9.148D-9, + & -1.682D-9, -10.626D-9, -0.679D-9, 7.869D-9, 6.729D-9, + & 0.400D-9, 3.554D-9, 3.733D-9, 27.071D-9, 4.520D-9, + & 68.458D-9, -3.722D-9 / + + IF (IPLATE .LE. 0 .OR. IPLATE .GT. NUMPLATE) THEN WRITE (LUOUT, 1) IPLATE 1 FORMAT(' Improper plate ID in PLATVL = ', I6) STOP @@ -519,8 +536,21 @@ c ENDIF VX = VX/1000.d0 VY = VY/1000.d0 VZ = VZ/1000.d0 + +*** The rotation rates for plates added in Dec 2020 refer to ITRF2014 +*** (plates 8 through 17). The following code converts the rates to +*** ITRF2008 velocities for these plates. + ELSEIF (IPLATE .GE. 8) THEN + VX = VX*1000.d0 + VY = VY*1000.d0 + VZ = VZ*1000.d0 + CALL VTRANF(X, Y, Z, VX, VY, VZ, 16, 15) + VX = VX/1000.d0 + VY = VY/1000.d0 + VZ = VZ/1000.d0 + *** The following translations rates are added per Altamimi et al. (2012) -*** for the other six plates +*** for the other six plates referenced to ITRF2008. ELSE VX = 0.00041d0 + VX VY = 0.00022d0 + VY @@ -1424,9 +1454,11 @@ C U3DS = U3DS + C3DS RETURN END -************************************************************* - SUBROUTINE OKADAW(PSI,ETA,Q,SDIP,CDIP,RATIO,TWOPI, - 1 VERT,U1SS,U2SS,U3SS,U1DS,U2DS,U3DS) + +C********1*********2*********3*********4*********5*********6*********7** + SUBROUTINE OKADAW(PSI,ETA,Q,SDIP,CDIP,RATIO,TWOPI,VERT,U1SS,U2SS, + & U3SS,U1DS,U2DS,U3DS) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) LOGICAL VERT @@ -1444,68 +1476,59 @@ C IF(VERT) THEN F5 = -RATIO*PSI*SDIP/(R + DBAR) F4 = -RATIO*Q/(R + DBAR) - F3 = 0.5D0*RATIO*(ETA/(R + DBAR) - 1 + YBAR*Q/((R + DBAR)*(R + DBAR)) - 2 - DLOG(R + ETA)) - F1 = -0.5D0*RATIO*PSI*Q/ - 1 ((R + DBAR)*(R + DBAR)) + F3 = 0.5D0*RATIO*(ETA/(R + DBAR) + YBAR*Q/ + & ((R + DBAR)*(R + DBAR)) - DLOG(R + ETA)) + F1 = -0.5D0*RATIO*PSI*Q/((R + DBAR)*(R + DBAR)) ELSE IF(DABS(PSI) .LE. 0.1D0) then F5 = 0.d0 ELSE - F5 = 2.D0*RATIO* - 1 DATAN((ETA*(X+Q*CDIP)+X*(R+X)*SDIP)/(PSI*(R+X)*CDIP)) - 2 /CDIP + F5 = 2.D0*RATIO*DATAN((ETA*(X+Q*CDIP)+X*(R+X)*SDIP)/ + & (PSI*(R+X)*CDIP))/CDIP ENDIF F4 = RATIO*(DLOG(R+DBAR)-SDIP*DLOG(R+ETA))/CDIP - F3 = RATIO*(YBAR/(CDIP*(R+DBAR)) - DLOG(R+ETA)) - 1 + SDIP*F4/CDIP + F3 = RATIO*(YBAR/(CDIP*(R+DBAR)) - DLOG(R+ETA)) + SDIP*F4/CDIP F1 = -RATIO*(PSI/(CDIP*(R+DBAR))) - SDIP*F5/CDIP ENDIF - F2 = -RATIO*DLOG(R+ETA) - F3 - - U1SS = -(PSI*Q/(R*(R+ETA)) - 1 + TERM + F1*SDIP)/TWOPI - U2SS = -(YBAR*Q/(R*(R+ETA)) - 1 + Q*CDIP/(R+ETA) - 2 + F2*SDIP)/TWOPI - U3SS = -(DBAR*Q/(R*(R+ETA)) - 1 + Q*SDIP/(R+ETA) - 2 + F4*SDIP)/TWOPI - U1DS = -(Q/R - F3*SDIP*CDIP)/TWOPI - U2DS = -(YBAR*Q/(R*(R+PSI)) - 1 + CDIP*TERM - F1*SDIP*CDIP)/TWOPI - U3DS = -(DBAR*Q/(R*(R+PSI)) - 1 + SDIP*TERM - F5*SDIP*CDIP)/TWOPI - RETURN - END -******************************************************************* - -C SUBROUTINE GRDCHK (POSX, POSY, INSIDE) - -C -C ROUTINE CHECKS IF THE POINT HAVING COORDINATES (POSX, POSY) -C IS WITHIN THE REGION SPANNED BY THE GRID -C -C IMPLICIT DOUBLE PRECISION (A-H,O-Z) -C IMPLICIT INTEGER*4 (I-N) -C LOGICAL INSIDE -C COMMON /CDGRID/ GRDLX, GRDUX, GRDLY, GRDUY, ICNTX, ICNTY + F2 = -RATIO*DLOG(R+ETA) - F3 -C INSIDE = .TRUE. - -C IF (POSX .LT. GRDLX .OR. POSX .GT. GRDUX) THEN -C INSIDE = .FALSE. -C ENDIF -C IF (POSY .LT. GRDLY .OR. POSY .GT. GRDUY) THEN -C INSIDE = .FALSE. -C ENDIF + U1SS = -(PSI*Q/(R*(R+ETA)) + TERM + F1*SDIP)/TWOPI + U2SS = -(YBAR*Q/(R*(R+ETA)) + Q*CDIP/(R+ETA) + F2*SDIP)/TWOPI + U3SS = -(DBAR*Q/(R*(R+ETA)) + Q*SDIP/(R+ETA) + F4*SDIP)/TWOPI + U1DS = -(Q/R - F3*SDIP*CDIP)/TWOPI -C RETURN -C END +* Following part of subroutine modified based on changes provided by +* Jarir Saleh on 11/19/2020 to avoid floating point exceptions when a +* point is located on the extension of a fault plane. +* Two conditional statements added to avoid floating point exception +* using variables TERM1 and TERM2: + TERM1 = YBAR*Q + IF (ABS(TERM1) < 1D-3) THEN + U2DS1 = 0.D0 + ELSE + U2DS1 = -YBAR*Q/R/(R+PSI)/TWOPI + ENDIF -******************************************************************* + U2DS2 = -CDIP*TERM/TWOPI + U2DS3 = F1*SDIP*CDIP/TWOPI + U2DS = U2DS1 + U2DS2 + U2DS3 + + TERM2 = DBAR*Q + IF (ABS(TERM2) < 1D-3) THEN + U3DS1 = 0.D0 + ELSE + U3DS1 = -DBAR*Q/R/(R+PSI)/TWOPI + ENDIF + + U3DS2 = -SDIP*TERM/TWOPI + U3DS3 = F5*SDIP*CDIP/TWOPI + U3DS = U3DS1 + U3DS2 + U3DS3 +* End of changes provided by Jarir + RETURN + END + +C********1*********2*********3*********4*********5*********6*********7** SUBROUTINE GRDWEI (YLON, YLAT, JREGN, I, J, WEI) C @@ -1713,6 +1736,7 @@ C*************************************************** C********************************************************* C SUBROUTINE TOMNT( IYR, IMON, IDAY, IHR, IMN, MINS ) +*** This subroutine is not called; apparently replaced by IYMDMJ. C C********1*********2*********3*********4*********5*********6*********7** C @@ -1785,13 +1809,22 @@ C C C........ 1.0 CALCULATION C - A= IYRP*0.01D0 - B= 2 - A + DINT( A*0.25D0 ) - C= 365.25D0*IYRP - D= 30.6001D0*(IMOP + 1) +*** Following 4 lines edited in v3.3.0 to address compiler issue with +*** real to integer conversion. +C A= IYRP*0.01D0 +C B= 2 - A + DINT( A*0.25D0 ) +C C= 365.25D0*IYRP +C D= 30.6001D0*(IMOP + 1) + + A = IDINT(DBLE(IYRP)*0.01D0) + B = 2 - A + IDINT(DBLE(A)*0.25D0) + C = IDINT(365.25D0*DBLE(IYRP)) + D = IDINT(30.6001D0*(DBLE(IMOP + 1))) + MINS = (B + C + D + IDAY - 679006) * (24 * 60) & + (60 * IHR) + IMN -C + + RETURN END ***************************************************** @@ -1842,7 +1875,7 @@ C COMMON /FILES/ LUIN,LUOUT, I1, I2, I3, I4, I5, I6 CHARACTER CARD*80 CHARACTER record*120 - CHARACTER NAMEF*30,NAME*30,NAMEBB*30, NAMEIF*30 + CHARACTER NAMEF*80, NAME*80, NAMEBB*80, NAMEIF*80 CHARACTER NAME24*24 CHARACTER BLAB*17 CHARACTER NAMEG*10 @@ -1884,7 +1917,7 @@ C 2 ' the predicted displacements. ') READ(LUIN,50,err=601,iostat=ios) NAMEF if (ios /= 0) goto 601 - 50 FORMAT(A30) + 50 FORMAT(A80) OPEN(I2,FILE=NAMEF,STATUS='UNKNOWN') CALL HEADER @@ -1912,7 +1945,7 @@ C 2 ' FROM ',F8.3, ' TO ',F8.3, ' (decimal years)'// 3 'NAME OF SITE LATITUDE LONGITUDE ', 4 ' NORTH EAST UP ') - 75 WRITE(LUOUT,80) + WRITE(LUOUT,80) !Remove label 75 (v3.3.0) 80 FORMAT(' ********************************'/ 1 ' Displacements will be predicted at each point whose',/ 2 ' horizontal position is specified.',/ @@ -1939,7 +1972,7 @@ C ELSEIF(OPTION .EQ. '1') THEN CALL GETPNT(LATD,LATM,SLAT,LATDIR,LOND,LONM,SLON, 1 LONDIR,NAME24, X,Y,Z,YLAT,YLON,EHT) - ELON = - YLON + ELON = -YLON call GETVLY(YLAT,ELON,VX,VY,VZ,VN,VE,VU,VOPT,210) if (vopt .eq. '0') then call PREDV( ylat, ylon, eht, date1, iopt, @@ -1957,11 +1990,11 @@ C 2 ' For additional displacements, please indicate how'/ 3 ' you wish to supply the horizontal coordinates.'/) WRITE(LUOUT,150) DN, DE,DU - 150 FORMAT(' *************************************'/ - 1 ' Northward displacement = ',F7.3,' meters.'/ - 1 ' Eastward displacement = ',F7.3,' meters.'/ - 1 ' Upward displacement = ',F7.3,' meters.'/ - 1 ' ****************************************'// + 150 FORMAT(' *****************************************'/ + 1 ' Northward displacement = ',F7.3,' meters'/ + 1 ' Eastward displacement = ',F7.3,' meters'/ + 1 ' Upward displacement = ',F7.3,' meters'/ + 1 ' *****************************************'// 2 ' For additional displacements, please indicate how'/ 3 ' you wish to supply the horizontal coordinates.'/) WRITE(I2,160) NAME24,LATD,LATM,SLAT,LATDIR, @@ -2015,7 +2048,7 @@ C 300 FORMAT(' Enter name of blue-book file ') READ(LUIN,310,err=603,iostat=ios) NAMEBB if (ios /= 0) goto 603 - 310 FORMAT(A30) + 310 FORMAT(A80) OPEN(I1,FILE=NAMEBB,STATUS='OLD') 320 READ(I1,330,END=350,err=604,iostat=ios) CARD if (ios /= 0) goto 604 @@ -2107,16 +2140,29 @@ C IF(JN.EQ.'S' .OR. JW.EQ.'E')GO TO 320 450 format(' Enter name of batch file ') read(luin, 451,err=606,iostat=ios) NAMEIF if (ios /= 0) goto 606 - 451 format(a30) + 451 format(A80) + open(I1,FILE=NAMEIF,STATUS='OLD') + LINE = 0 !Inititalize line number of input file. 455 read(I1,'(a)',END=460,err=607,iostat=ios) record if (ios /= 0) goto 607 -c write(i2, 457) record call interprate_latlon_record (record,XLAT,XLON,name24) -c write (i2,456) xlat, xlon -c write(i2,457) record -c 457 format (a50) -c 456 format(1x,' xlat = ', f15.3, 'xlon = ', f15.3) + LINE = LINE + 1 !Increment line number of input file. + +*** If latitude magnitudes > 90 degrees or longitude magnitude > 360 degrees, +*** write error message and terminate program (added for v3.3.0). + IF ((DABS(XLAT) .GT. 90).OR.(DABS(XLON) .GT. 360)) THEN + WRITE(*,*)'***********************************************' + WRITE(*,*)'Invalid latitude or longitude in input file' + WRITE(*,4551)'on line', LINE, ':' + WRITE(*,4552) XLAT, XLON, name24 + WRITE(*,*)'Please check your input file and try again.' + WRITE(*,*)'***********************************************' + 4551 FORMAT (1X, A, I8, A/) + 4552 FORMAT (1X, 2F17.10, 2X, A/) + STOP + ENDIF + YLAT = (XLAT*3600.D0)/RHOSEC YLON = (XLON*3600.D0)/RHOSEC CALL TODMSS(YLAT,LATD,LATM,SLAT,ISIGN) @@ -2154,6 +2200,7 @@ c 456 format(1x,' xlat = ', f15.3, 'xlon = ', f15.3) WRITE(LUOUT,500) 500 FORMAT(' Improper entry--select again. ') ENDIF + GO TO 85 510 CONTINUE CLOSE(I2, STATUS = 'KEEP') @@ -2220,7 +2267,7 @@ c 456 format(1x,' xlat = ', f15.3, 'xlon = ', f15.3) COMMON /CONST/ A,F,E2,EPS,AF,PI,TWOPI,RHOSEC COMMON /FILES/ LUIN,LUOUT, I1, I2, I3, I4, I5, I6 CHARACTER CARD*80 - CHARACTER NAMEF*30,NAME*30,NAMEBB*30,PVFILE*30 + CHARACTER NAMEF*80,NAME*80,NAMEBB*80,PVFILE*80 CHARACTER NAME24*24 CHARACTER BLAB*17 CHARACTER NAMEG*10 @@ -2230,11 +2277,12 @@ c 456 format(1x,' xlat = ', f15.3, 'xlon = ', f15.3) character frame1*24 character record*120 +*** Commented out following temporary code (v3.3.0). *** temporary code to plot results **************** - character TEMPNA - - TEMPNA = ' ' - DUMMY = 0.0D0 +C character TEMPNA*8 !Added string length to match string below. +C +C TEMPNA = ' ' +C DUMMY = 0.0D0 *** end of temporary code ************************* BLAB = 'OUTSIDE OF REGION' @@ -2244,7 +2292,7 @@ c 456 format(1x,' xlat = ', f15.3, 'xlon = ', f15.3) 1 ' predicted velocities. ') READ(LUIN,20,err=700,iostat=ios) NAMEF if (ios /= 0) goto 700 - 20 FORMAT(A30) + 20 FORMAT(A80) OPEN(I2,FILE=NAMEF,STATUS='UNKNOWN') CALL HEADER @@ -2260,7 +2308,7 @@ c 456 format(1x,' xlat = ', f15.3, 'xlon = ', f15.3) WRITE(LUOUT,22) 22 FORMAT(' Enter name for the file that is to contain'/ 1 ' the PV records. ') - READ(LUIN,'(A30)',err=702,iostat=ios) PVFILE + READ(LUIN,'(A80)',err=702,iostat=ios) PVFILE if (ios /= 0) goto 702 OPEN(I3,FILE = PVFILE, STATUS = 'UNKNOWN') ENDIF @@ -2434,27 +2482,28 @@ c 6 ' a specified point having a specified velocity.') 1 1X,F8.5,1X,A1, 1X,A17) ELSE +*** Commented out following temporary code (v3.3.0). *** temporary code to plot results - ZLAT = DBLE(LATD) + DBLE(LATM)/60.D0 + SLAT/3600.D0 - ZLON = DBLE(LOND) + DBLE(LONM)/60.D0 + SLON/3600.D0 - ZLON = 360.D0 - ZLON - TOTVEL = DSQRT(VN*VN + VE*VE) - +C ZLAT = DBLE(LATD) + DBLE(LATM)/60.D0 + SLAT/3600.D0 +C ZLON = DBLE(LOND) + DBLE(LONM)/60.D0 + SLON/3600.D0 +C ZLON = 360.D0 - ZLON +C TOTVEL = DSQRT(VN*VN + VE*VE) +C *** code to create vectors to be plotted -c IF (TOTVEL .GE. 5.0D0 .AND. TOTVEL .LE. 50.D0) THEN -c TVE = VE/10.D0 -c TVN = VN/10.D0 -c WRITE(I2,129) ZLON,ZLAT,TVE,TVN,DUMMY,DUMMY,DUMMY, -c 1 TEMPNA -c 129 FORMAT(F10.6, 2x, F9.6, 2X, 5(F7.2, 2x), A8) -c ENDIF - +C IF (TOTVEL .GE. 5.0D0 .AND. TOTVEL .LE. 50.D0) THEN +C TVE = VE/10.D0 +C TVN = VN/10.D0 +C WRITE(I2,129) ZLON,ZLAT,TVE,TVN,DUMMY,DUMMY,DUMMY, +C 1 TEMPNA +C 129 FORMAT(F10.6, 2x, F9.6, 2X, 5(F7.2, 2x), A8) +C ENDIF +C *** code to create contour plots C WRITE(I2,129) ZLON, ZLAT, TOTVEL C 129 FORMAT(F6.2, 1X, F6.2, 1X, F5.1) - *** end of temporary code + *** beginning of original code WRITE(I2,130)NAMEG,I,J,LATD,LATM,SLAT,LATDIR, 1 LOND,LONM,SLON,LONDIR,VN,VE,VU @@ -2479,7 +2528,7 @@ C 129 FORMAT(F6.2, 1X, F6.2, 1X, F5.1) 200 FORMAT(' Enter name of blue-book file. ') READ(LUIN,210,err=705,iostat=ios) NAMEBB if (ios /= 0) goto 705 - 210 FORMAT(A30) + 210 FORMAT(A80) OPEN(I1,FILE=NAMEBB,STATUS='OLD') 220 READ(I1,230,END=250,err=706,iostat=ios) CARD if (ios /= 0) goto 706 @@ -2571,11 +2620,29 @@ C 129 FORMAT(F6.2, 1X, F6.2, 1X, F5.1) 400 FORMAT(' Enter name of input file. ') READ(LUIN,410,err=708,iostat=ios) NAMEBB if (ios /= 0) goto 708 - 410 FORMAT(A30) + 410 FORMAT(A80) + OPEN(I1,FILE=NAMEBB,STATUS='OLD') + LINE = 0 !Inititalize line number of input file. 420 READ(I1,'(a)',END=450,err=709,iostat=ios) record call interprate_latlon_record (record,XLAT,XLON,NAME24) if (ios /= 0) goto 709 + LINE = LINE + 1 !Increment line number of input file. + +*** If latitude magnitudes > 90 degrees or longitude magnitude > 360 degrees, +*** write error message and terminate program (added for v3.3.0). + IF ((DABS(XLAT) .GT. 90).OR.(DABS(XLON) .GT. 360)) THEN + WRITE(*,*)'***********************************************' + WRITE(*,*)'Invalid latitude or longitude in input file' + WRITE(*,4201)'on line', LINE, ':' + WRITE(*,4202) XLAT, XLON, name24 + WRITE(*,*)'Please check your input file and try again.' + WRITE(*,*)'***********************************************' + 4201 FORMAT (1X, A, I8, A/) + 4202 FORMAT (1X, 2F17.10, 2X, A/) + STOP + ENDIF + YLAT = (XLAT*3600.D0)/RHOSEC YLON = XLON*3600.d0/RHOSEC CALL TODMSS(YLAT, LATD, LATM, SLAT, ISIGN) @@ -2692,14 +2759,21 @@ C 129 FORMAT(F6.2, 1X, F6.2, 1X, F5.1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) logical Is_iopt_NAD83 + COMMON /CONST/ A,F,E2,EPS,AF,PI,TWOPI,RHOSEC Is_iopt_NAD83 = (IOPT == 1) *** Get reference latitude RLAT and reference longitude RLON in NAD 83 -c IF(IOPT .EQ. 0 .OR. IOPT .EQ. 1) THEN !Not NAD83 antmore +c IF(IOPT .EQ. 0 .OR. IOPT .EQ. 1) THEN !Not NAD83 anymore IF(IOPT .EQ. 15) THEN RLAT = YLAT - RLON = YLON +*** Added conditional statement to get rid of out-of-region error when +*** YLON is negative for reference frame (v3.3.0): + IF(YLON .GE. 0.D0) THEN + RLON = YLON + ELSE + RLON = YLON + TWOPI + ENDIF ELSE ELON = -YLON CALL TOXYZ(YLAT,ELON,EHT,X,Y,Z) @@ -2819,7 +2893,7 @@ c write (*,*) "FROM XT08 ",RLAT*180.d0/pi,WLON*180.d0/pi,EHT08 character card*80,namebb*80,nameif*80,name24*80 character record*120 - character namef*30 + character NAMEF*80 character frame1*24, frame2*24 character jn*1,jw*1,LATDIR*1,LONDIR*1 character option*1, answer*1, vopt*1 @@ -2840,7 +2914,7 @@ c write (*,*) "FROM XT08 ",RLAT*180.d0/pi,WLON*180.d0/pi,EHT08 80 format( 1 ' Please enter the name for the file to contain'/ 2 ' the transformed positions. ') - read(luin,'(a30)',err=600,iostat=ios) namef + read(luin,'(A80)',err=600,iostat=ios) namef if (ios /= 0) goto 600 open(i2, file = namef, status = 'unknown') CALL HEADER @@ -2864,7 +2938,7 @@ c write (*,*) "FROM XT08 ",RLAT*180.d0/pi,WLON*180.d0/pi,EHT08 write(luout,*) ' Improper selection -- try again. ' go to 105 endif - Is_out_NAD83 = (iopt2 == 1 .and. frame1(1:3) /= "WGS") + Is_out_NAD83 = (iopt2 == 1 .and. frame2(1:3) /= "WGS") !Changed from frame1 in v3.3.0. Is_out_NAD83PAC = (iopt2 == 12) Is_out_NAD83MAR = (iopt2 == 13) @@ -3070,10 +3144,10 @@ c open (i6,file='transformed_'//namebb) C write some comments in the transformed files call extract_name (namebb,iii) -c write (I2,1309) namebb(1:iii) -c write (I2,1309) trim(namebb) - 1309 format (/'The file, transformed_',a, - & ', contains the transformed input file.'/) + +*** Commented out due to unused label (v3.3.0). +*1309 format (/'The file, transformed_',a, +* & ', contains the transformed input file.'/) write (i2,1310) HTDP_version 1310 format (' ***CAUTION: This file was processed using HTDP', @@ -3093,15 +3167,16 @@ c write (I2,1309) trim(namebb) 310 read(i1, '(a80)', end = 390,err=604,iostat=ios) card if (ios /= 0) goto 604 if (card(7:10) .eq. '*80*') then -c write (*,'(a80)' ) card read(card,320,err=605,iostat=ios) isn,name24(1:30),latd, & latm,slat,jn,lond, lonm, slon, jw -c write (*,321) latd, latm, slat, jn, lond, lonm, slon, jw if (ios /= 0) goto 605 320 format(10x, i4, a30, i2, i2, f7.5, a1, & i3, i2, f7.5, a1) - 321 format(i2, i2, f7.5, a1,3x, i3, i2, f7.5, a1) + +*** Commented out due to unused label (v3.3.0). +* 321 format(i2, i2, f7.5, a1,3x, i3, i2, f7.5, a1) + ylat = (dble((latd*60 + latm)*60) + slat) / rhosec ylon = (dble((lond*60 + lonm)*60) + slon) / rhosec if (jn .eq. 'S') ylat = -ylat @@ -3223,16 +3298,30 @@ C open (i6,file='transformed_'//nameif) C write some comments in the transformed files call extract_name (nameif,iii) -c write (I2,1309) nameif(1:iii) -c write (I2,1309) trim(nameif) write (i2,1310) HTDP_version write (i2,1311) frame2 WRITE (I2, 1312) MONTH2, IDAY2, IYEAR2, DATE2 -c const = 180.d0/PI + LINE = 0 !Inititalize line number of input file. 410 read (i1,'(a)',end=450,err=607,iostat=ios) record if (ios /= 0) goto 607 call interprate_XYZ_record (record,xlat,xlon,eht,name24) + LINE = LINE + 1 !Increment line number of input file. + +*** If latitude magnitudes > 90 degrees or longitude magnitude > 360 degrees, +*** write error message and terminate program (added for v3.3.0). + IF ((DABS(xlat) .GT. 90).OR.(DABS(xlon) .GT. 360)) THEN + WRITE(*,*)'***********************************************' + WRITE(*,*)'Invalid latitude or longitude in input file' + WRITE(*,4101)'on line', LINE, ':' + WRITE(*,4102) xlat, xlon, eht, name24 + WRITE(*,*)'Please check your input file and try again.' + WRITE(*,*)'***********************************************' + 4101 FORMAT (1X, A, I8, A/) + 4102 FORMAT (1X, 2F17.10, F11.3, 2X, A) + STOP + ENDIF + ylat = (xlat*3600.d0) / rhosec ylon = (xlon*3600.d0) / rhosec elon = -ylon @@ -3299,7 +3388,7 @@ c 1 name24,0,vxsave,vysave,vzsave,vnsave,vesave,vusave) write (i2,449) outlat,outlon,ehtnew,name24(1:iii) c write (i6,449) outlat,outlon,ehtnew,trim(name24) - 449 format (2f16.10,f10.3,4x,a) + 449 format (2f16.10,f10.3,4x,a) go to 410 450 close(i1, status = 'keep') @@ -3314,7 +3403,7 @@ c close(i6, status = 'keep') vesave = 0.d0 vusave = 0.d0 write (luout, 4001) -4001 format (' Enter name of input file: ') + 4001 format (' Enter name of input file: ') read (luin, '(a)',err=600,iostat=ios) nameif if (ios /= 0) goto 600 open (i1, file = nameif, status = 'old') @@ -3323,12 +3412,6 @@ c open (i6,file='transformed_'//nameif) C write some comments in the transformed files call extract_name (nameif,iii) -c write (i2,1309) nameif(1:iii) -c write (i2,1309) trim(nameif) -c write (i6,1310) HTDP_version -c write (i6,1311) frame2 -c WRITE (i6, 1312) MONTH2, IDAY2, IYEAR2, DATE2 -c const = 180.d0/PI 411 read (i1,'(a)',end=451,err=607,iostat=ios) record if (ios /= 0) goto 607 call interprate_XYZ_record (record,x,y,z,name24) @@ -3420,12 +3503,6 @@ c open (i6,file='transformed_'//nameif) C write some comments in the transformed files call extract_name (nameif,iii) -c write (i2,1309) nameif(1:iii) -c write (i2,1309) trim(nameif) -c write (i6,1310) HTDP_version -c write (i6,1311) frame2 -c WRITE (i6, 1312) MONTH2, IDAY2, IYEAR2, DATE2 -c const = 180.d0/PI 511 read (i1,'(a)',end=551,err=607,iostat=ios) record if (ios /= 0) goto 607 call interprate_XYZVxVyVz_record (record,x,y,z,Vx,Vy,Vz,name24) @@ -3630,7 +3707,7 @@ C Important note: C The parameters in common block tranpa are computed using the IGS convention ITRF96 <> ITRF97. C The parameters in common block tranpa1 are computed using the IERS convention ITRF96 = ITRF97. C The latter parameters were added to HTDP in 09/2014. They will be used to transform between -C ITRF systems. They will not be used if the transformation involves NAD83 or WGS84 (transit). +C ITRF systems. They will not be used if the transformation involves NAD83 or WGS84 (original/transit). C However, the NAD 83 Pacific frames use the IERS convention ITRF96 = ITRF97. @@ -3670,6 +3747,7 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 refepc(1) = 1997.0d0 *** From ITRF94 to ITRF88 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx(2) = 0.018d0 ty(2) = 0.000d0 tz(2) = -0.092d0 @@ -3682,11 +3760,12 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 drx(2) = 0.0d0 dry(2) = 0.0d0 drz(2) = 0.0d0 - scale(2) = 0.74d-8 + scale(2) = 0.749d-8 !Previous scale = 0.74d-8 (prior to v3.3.0) dscale(2) = 0.0d0 refepc(2) = 1988.0d0 *** From ITRF94 to ITRF89 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx(3) = 0.023d0 ty(3) = 0.036d0 tz(3) = -0.068d0 @@ -3699,11 +3778,12 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 drx(3) = 0.0d0 dry(3) = 0.0d0 drz(3) = 0.0d0 - scale(3) = 0.43d-8 + scale(3) = 0.439d-8 !Previous scale = 0.43d-8 (prior to v3.3.0) dscale(3) = 0.0d0 refepc(3) = 1988.0d0 *** From ITRF94 to ITRF90 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx(4) = 0.018d0 ty(4) = 0.012d0 tz(4) = -0.030d0 @@ -3716,11 +3796,12 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 drx(4) = 0.0d0 dry(4) = 0.0d0 drz(4) = 0.0d0 - scale(4) = 0.09d-8 + scale(4) = 0.099d-8 !Previous scale = 0.09d-8 (prior to v3.3.0) dscale(4) = 0.0d0 refepc(4) = 1988.0d0 *** From ITRF94 to ITRF91 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx(5) = 0.020d0 ty(5) = 0.016d0 tz(5) = -0.014d0 @@ -3733,11 +3814,12 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 drx(5) = 0.0d0 dry(5) = 0.0d0 drz(5) = 0.0d0 - scale(5) = 0.06d-8 + scale(5) = 0.069d-8 !Previous scale = 0.06d-8 (prior to v3.3.0) dscale(5) = 0.0d0 refepc(5) = 1988.0d0 *** From ITRF94 to ITRF92 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx(6) = 0.008d0 ty(6) = 0.002d0 tz(6) = -0.008d0 @@ -3750,11 +3832,20 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 drx(6) = 0.0d0 dry(6) = 0.0d0 drz(6) = 0.0d0 - scale(6) = -0.08d-8 + scale(6) = -0.071d-8 !Previous scale = -0.08d-8 (prior to v3.3.0) dscale(6) = 0.0d0 refepc(6) = 1988.0d0 *** From ITRF94 to ITRF93 +*** Update scale in v3.3.0 to match current value published by IERS +*** (https://itrf.ign.fr/trans_para.php). Supersedes scale in file +*** "ITRF94.TX" at ftp://itrf-ftp.ign.fr/pub/itrf/itrf94/. +*** Previous link under “Other transformation parameters” at +*** https://itrf.ign.fr/ITRF_solutions/94/ITRF94.php?page=2 currently +*** does not work (as of 1/15/2021). +*** This update also affects the scale for transformations from ITRF94 +*** to all earlier ITRFs. +*** Note that this transformation is not published in IERS TN20 (1996). tx(7) = 0.006d0 ty(7) = -0.005d0 tz(7) = -0.015d0 @@ -3767,7 +3858,7 @@ C Parameters computed with the IGS convention ITRF96 <> ITRF97 drx(7) = 0.00011d0 / rhosec dry(7) = 0.00019d0 / rhosec drz(7) = -0.00005d0 / rhosec - scale(7) = 0.04d-8 + scale(7) = 0.049d-8 !Previous scale = 0.040d-8 (prior to v3.3.0) dscale(7) = 0.0d0 refepc(7) = 1988.0d0 @@ -3968,6 +4059,7 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 refepc1(1) = 1997.0d0 *** From ITRF94 to ITRF88 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx1(2) = 0.018d0 ty1(2) = 0.000d0 tz1(2) = -0.092d0 @@ -3980,11 +4072,12 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 drx1(2) = 0.0d0 dry1(2) = 0.0d0 drz1(2) = 0.0d0 - scale1(2) = 0.74d-8 + scale1(2) = 0.749d-8 !Previous scale = 0.74d-8 (prior to v3.3.0) dscale1(2) = 0.0d0 refepc1(2) = 1988.0d0 *** From ITRF94 to ITRF89 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx1(3) = 0.023d0 ty1(3) = 0.036d0 tz1(3) = -0.068d0 @@ -3997,11 +4090,12 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 drx1(3) = 0.0d0 dry1(3) = 0.0d0 drz1(3) = 0.0d0 - scale1(3) = 0.43d-8 + scale1(3) = 0.439d-8 !Previous scale = 0.43d-8 (prior to v3.3.0) dscale1(3) = 0.0d0 refepc1(3) = 1988.0d0 *** From ITRF94 to ITRF90 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx1(4) = 0.018d0 ty1(4) = 0.012d0 tz1(4) = -0.030d0 @@ -4014,11 +4108,12 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 drx1(4) = 0.0d0 dry1(4) = 0.0d0 drz1(4) = 0.0d0 - scale1(4) = 0.09d-8 + scale1(4) = 0.099d-8 !Previous scale = 0.09d-8 (prior to v3.3.0) dscale1(4) = 0.0d0 refepc1(4) = 1988.0d0 *** From ITRF94 to ITRF91 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx1(5) = 0.020d0 ty1(5) = 0.016d0 tz1(5) = -0.014d0 @@ -4031,11 +4126,12 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 drx1(5) = 0.0d0 dry1(5) = 0.0d0 drz1(5) = 0.0d0 - scale1(5) = 0.06d-8 + scale1(5) = 0.069d-8 !Previous scale = 0.06d-8 (prior to v3.3.0) dscale1(5) = 0.0d0 refepc1(5) = 1988.0d0 *** From ITRF94 to ITRF92 +*** Update scale in v3.3.0 due to ITRF94 to ITRF93 transformation update. tx1(6) = 0.008d0 ty1(6) = 0.002d0 tz1(6) = -0.008d0 @@ -4048,11 +4144,20 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 drx1(6) = 0.0d0 dry1(6) = 0.0d0 drz1(6) = 0.0d0 - scale1(6) = -0.08d-8 + scale1(6) = -0.071d-8 !Previous scale = -0.08d-8 (prior to v3.3.0) dscale1(6) = 0.0d0 refepc1(6) = 1988.0d0 *** From ITRF94 to ITRF93 +*** Update scale in v3.3.0 to match current value published by IERS +*** (https://itrf.ign.fr/trans_para.php). Supersedes scale in file +*** "ITRF94.TX" at ftp://itrf-ftp.ign.fr/pub/itrf/itrf94/. +*** Previous link under “Other transformation parameters” at +*** https://itrf.ign.fr/ITRF_solutions/94/ITRF94.php?page=2 currently +*** does not work (as of 1/15/2021). +*** This update also affects the scale for transformations from ITRF94 +*** to all earlier ITRFs. +*** Note that this transformation is not published in IERS TN20 (1996). tx1(7) = 0.006d0 ty1(7) = -0.005d0 tz1(7) = -0.015d0 @@ -4065,7 +4170,7 @@ C Parameters computed with the IERS convention ITRF96 = ITRF97 drx1(7) = 0.00011d0 / rhosec dry1(7) = 0.00019d0 / rhosec drz1(7) = -0.00005d0 / rhosec - scale1(7) = 0.04d-8 + scale1(7) = 0.049d-8 !Previous scale = 0.04d-8 (prior to v3.3.0) dscale1(7) = 0.0d0 refepc1(7) = 1988.0d0 @@ -4489,7 +4594,7 @@ c write (*,*) "TO TIIT94_IERS ",x2,y2,z2 iframe(4) = 10 nframe(4) = 'WGS_72 ' iframe(5) = 1 - nframe(5) = 'WGS_84(transit) ' + nframe(5) = 'WGS_84(original) ' c iframe(6) = 6 (This was incorrect in all versions of HTDP) iframe(6) = 5 nframe(6) = 'WGS_84(G730) ' @@ -4532,14 +4637,14 @@ c nframe(10)= 'PNEOS_90 or NEOS_90 ' iframe(24)= 16 nframe(24)= 'ITRF2014 or IGS14/IGb14 ' - write(luout, 100) + write(luout, 100) 100 format( 1' 1...NAD_83(2011/CORS96/2007) (North American plate fixed) '/ 1' 2...NAD_83(PA11/PACP00) (Pacific plate fixed) '/ 1' 3...NAD_83(MA11/MARP00) (Mariana plate fixed) '/ 1' '/ c 1' 4...WGS_72 '/ - 1' 5...WGS_84(transit) (NAD_83(2011) used) 15...ITRF91 '/ + 1' 5...WGS_84(original) (NAD_83(2011) used) 15...ITRF91 '/ 1' 6...WGS_84(G730) (ITRF91 used) 16...ITRF92 '/ 1' 7...WGS_84(G873) (ITRF94 used) 17...ITRF93 '/ 1' 8...WGS_84(G1150) (ITRF2000 used) 18...ITRF94 '/ @@ -4561,6 +4666,7 @@ c 1' 14...ITRF90 '/ ) mframe = ' ' kopt = iopt endif + return 50 write (*,*) 'Failed to read option in MENU1:ios=',ios @@ -4579,9 +4685,9 @@ c 1' 14...ITRF90 '/ ) IMPLICIT INTEGER*4 (I-N) parameter (numref = 16) parameter (nbbdim = 10000) - CHARACTER OLDBB*30,NEWBB*30, NAMEIF*30 + CHARACTER OLDBB*80, NEWBB*80, NAMEIF*80 CHARACTER NAME24*24 - CHARACTER OPT*1,ANSWER*1,BBTYPE*1,VOPT*1 + CHARACTER OPT*1, ANSWER*1, BBTYPE*1, VOPT*1 CHARACTER LATDIR*1, LONDIR*1, LATDR*1, LONDR*1 character frame1*24, frame2*24 character HTDP_version*8 @@ -4589,7 +4695,6 @@ c 1' 14...ITRF90 '/ ) LOGICAL Is_iopt_NAD83 COMMON /CONST/ A,F,E2,EPS,AF,PI,TWOPI,RHOSEC COMMON /FILES/ LUIN, LUOUT, I1, I2, I3, I4, I5, I6 - COMMON /CAUTION/ MONTH2,IDAY2,IYEAR2,DATE2,frame1 WRITE(LUOUT,20) 20 FORMAT(' ********************************************'/ @@ -4656,7 +4761,7 @@ c 1' 14...ITRF90 '/ ) ENDIF WRITE(LUOUT,1030) 1030 FORMAT(' ENTER file to save updated positions. ') - READ(LUIN,'(A30)',err=502,iostat=ios) NEWBB + READ(LUIN,'(A80)',err=502,iostat=ios) NEWBB if (ios /= 0) goto 502 OPEN(I2,FILE=NEWBB,STATUS='UNKNOWN') CALL HEADER @@ -4729,6 +4834,7 @@ c 1' 14...ITRF90 '/ ) 1 f8.2, ' mm/yr',f8.3, ' m'/ 1 1X, 'Z ',2(6X,F13.3), 1 f8.2, ' mm/yr',F8.3,' m'/) + 1075 WRITE(LUOUT,1080) 1080 FORMAT(' Update more positions? (y/n) ') READ(LUIN,'(A1)',err=500,iostat=ios) ANSWER @@ -4758,7 +4864,7 @@ c 1' 14...ITRF90 '/ ) 100 FORMAT(' Enter name of the blue book file to be updated.'/) READ(LUIN,110,err=502,iostat=ios) OLDBB if (ios /= 0) goto 502 - 110 FORMAT(A30) + 110 FORMAT(A80) WRITE(LUOUT,120) 120 FORMAT(' Enter name for the new blue book file that is to '/ 1 'contain the updated information.'/) @@ -4815,9 +4921,9 @@ c * 'updated to ',I2,'-',I2.2,'-',I4, ' = (',F8.3,') ***') c ENDIF IF(BBTYPE .EQ. '1') THEN - CALL UPBB4(MIN1,MIN2,OPT,IOPT) + CALL UPBB4(MIN1,MIN2,OPT) ELSEIF(BBTYPE .EQ. '2') THEN - CALL UPBB5(MIN1,MIN2,OPT,IOPT) + CALL UPBB5(MIN1,MIN2,OPT) ENDIF CLOSE(I1, STATUS = 'KEEP') CLOSE(I2, STATUS = 'KEEP') @@ -4834,11 +4940,11 @@ c ENDIF ELSEIF(ANSWER .eq. 'Y' .or. ANSWER .eq. 'y')THEN WRITE(LUOUT,620) 620 FORMAT(' Enter name of old G-FILE to be updated. ') - READ(LUIN,'(A30)',err=502,iostat=ios) OLDBB + READ(LUIN,'(A80)',err=502,iostat=ios) OLDBB if (ios /= 0) goto 502 WRITE(LUOUT,630) 630 FORMAT(' Enter name for the new updated G-FILE. ') - READ(LUIN,'(A30)',err=502,iostat=ios) NEWBB + READ(LUIN,'(A80)',err=502,iostat=ios) NEWBB if (ios /= 0) goto 502 OPEN(I1,FILE=OLDBB,STATUS='OLD') OPEN(I2,FILE=NEWBB,STATUS='UNKNOWN') @@ -4897,7 +5003,7 @@ c ENDIF RETURN endif write(luout, 1030) - read(luin, '(a30)',err=502,iostat=ios)NEWBB + read(luin, '(A80)',err=502,iostat=ios)NEWBB if (ios /= 0) goto 502 open(I2, FILE=NEWBB, STATUS='UNKNOWN') CALL HEADER @@ -5433,12 +5539,8 @@ C 1 YLATT,YLONT,EHTNEW,DN,DE,DU,IOPT) ENDIF *** Unrecognized blue book record +*** Associated ELSE statement with label 500 deleted in v3.3.0. - ELSE -c WRITE(LUOUT,500) TYPE - 500 FORMAT(' This software does not recognize an ',A4,/ - 1 ' record. The record will be copied to the new'/ - 2 ' file without change.') ENDIF WRITE(I2,175) CARD GO TO 170 @@ -5617,7 +5719,7 @@ C 200 FORMAT(BZ, 9X,I5,T40,A6,T51,I5,T64,I3,I2,F3.1) ENDIF ANGLE = DABS(DA) CALL TODMSS(ANGLE,IDEG,MIN,SEC,ISIGN) - ISEC1 = SEC + ISEC1 = IDINT(SEC) !Convert SEC to integer (v3.3.0) ISEC2 = IDNINT((SEC - DBLE(ISEC1)) * 100.D0) WRITE(CARD(152:162),211)SIGN,IDEG,MIN,ISEC1,ISEC2 211 FORMAT(A1,I3.3,I2.2,I2.2,'.',I2.2) @@ -5679,7 +5781,7 @@ C 220 FORMAT(BZ, 9X,I5,I2,T40,A6,T51,I5) ENDIF ANGLE = DABS(DA - DA0) CALL TODMSS(ANGLE,IDEG,MIN,SEC,ISIGN) - ISEC1 = SEC + ISEC1 = IDINT(SEC) !Convert SEC to integer (v3.3.0) ISEC2 = IDNINT((SEC-DBLE(ISEC1))*100.D0) WRITE(CARD(132:142),251) SIGN,IDEG,MIN,ISEC1,ISEC2 251 FORMAT(A1,I3.3,I2.2,I2.2,'.',I2.2) @@ -5762,12 +5864,6 @@ C 1 YLATT,YLONT,EHTNEW,DN,DE,DU,IOPT) write (*,*) "PLEASE CHECK YOUR INPUT FILE AND TRY AGAIN" stop - 701 write (*,'(/)') - write (*,*) "Failed to read card *12* in UPBB4:ios=",ios - write (*,*) "ABNORMAL TERMINATION" - write (*,*) "PLEASE CHECK YOUR INPUT FILE AND TRY AGAIN" - stop - 702 write (*,'(/)') write (*,*) "Failed to read *50,51,52* in UPBB4:ios=",ios write (*,*) "ABNORMAL TERMINATION" @@ -5798,12 +5894,6 @@ C 1 YLATT,YLONT,EHTNEW,DN,DE,DU,IOPT) write (*,*) "PLEASE CHECK YOUR INPUT FILE AND TRY AGAIN" stop - 707 write (*,'(/)') - write (*,*) "Failed to read *30,32* in UPBB4:ios=",ios - write (*,*) "ABNORMAL TERMINATION" - write (*,*) "PLEASE CHECK YOUR INPUT FILE AND TRY AGAIN" - stop - 708 write (*,'(/)') write (*,*) "Failed to read *80* in UPBB4:ios=",ios write (*,*) "ABNORMAL TERMINATION" @@ -6188,12 +6278,6 @@ C DECYR2 = DBLE(IYEAR2) + DBLE(MINO2 - MIN00)/525600.D0 write (*,*) "PLEASE CHECK YOUR INPUT FILE AND TRY AGAIN" stop - 302 write (*,'(/)') - write (*,*) "Failed to read C card in UPGFIG:ios=",ios - write (*,*) "ABNORMAL TERMINATION" - write (*,*) "PLEASE CHECK YOUR INPUT FILE AND TRY AGAIN" - stop - 303 write (*,'(/)') write (*,*) "Failed to read F card in UPGFIG:ios=",ios write (*,*) "ABNORMAL TERMINATION" @@ -6379,8 +6463,8 @@ C coordinates as if they were NAD 83(2011)coordinates. IRFCON(37) = 16 -*** From HTDP identifier to blue book identifier -*** NAD 83 (set equal to WGS 84 (transit)) +*** From HTDP identifier to blue book identifier. +*** NAD 83 (2011/2007/CORS96/...) referenced to North America plate. c JRFCON(1) = 2 JRFCON(1) = 34 @@ -6418,17 +6502,11 @@ c JRFCON(1) = 2 c JRFCON(12) = 0 c JRFCON(12) = 2 JRFCON(12) = 35 -C Michael Dennis requested that HTDP identify -C NAD 83 (PA11) coordinates as WGS 84(transit) -C coordinates in the Bluebook context *** NAD 83 (MARP00) or NAD 83 (MA11) c JRFCON(13) = 0 c JRFCON(13) = 2 JRFCON(13) = 36 -C Michael Dennis requested that HTDP identify -C NAD 83 (MA11) coordinates as WGS 84(transit) -C coordinates in the Bluebook context *** ITRF2005 or IGS05 JRFCON(14) = 26 @@ -6675,8 +6753,88 @@ C The parameters in common block tranpa1 are computed using the IERS values of 2 ' with north being positive. For example, 35,17,28.3 '/ 2 ' For a point in the southern hemisphere, enter a minus sign'/ 3 ' before each value. For example, -35,-17,-28.3') + READ(LUIN,*,err=202,iostat=ios) LATD, LATM, SLAT if (ios /= 0) goto 202 + +*** START insertion of new code and comments for v3.3.0 to handle +*** incorrect input latitude and longitude. +*** +*** Test whether latitude magnitude is > 90 degrees. If it is, prompt +*** for re-entry. + TEST = DBLE(ABS(LATD)) + DBLE(ABS(LATM))/60.D0 + & + DABS(SLAT)/3600.D0 + DO WHILE (DABS(TEST) .GT. 90.D0) + WRITE(*,*)' Latitude magnitude (degrees) =', TEST + WRITE(*,*) + & ' Magnitude cannot exceed 90. Please re-enter all values.' + READ(LUIN,*,err=202,iostat=ios) LATD, LATM, SLAT + if (ios /= 0) goto 202 + TEST = DBLE(ABS(LATD)) + DBLE(ABS(LATM))/60.D0 + & + DABS(SLAT)/3600.D0 + END DO + +*** Non-zero latitude minutes and seconds must be the same sign, and +*** if latitude degrees is non-zero, minutes and seconds must be the +*** same sign as latitude degrees. Otherwise prompt for re-entry. + IF (LATD .EQ. 0) THEN + IF (DBLE(LATM)*SLAT .LT. 0.D0) THEN + WRITE(*,*)' Minutes and seconds must be the same sign.', + & ' Please re-enter minutes and seconds.' + READ(LUIN,*,err=202,iostat=ios) LATM, SLAT + if (ios /= 0) goto 202 + ENDIF + ENDIF + IF (LATD .GT. 0) THEN + IF ((LATM .LT. 0).OR.(SLAT .LT. 0.D0)) + & WRITE(*,*)' Latitude degrees are positive.' + DO WHILE (LATM .LT. 0) + WRITE(*,*) + & ' Minutes must also be positive. Please re-enter minutes.' + READ(LUIN,*,err=202,iostat=ios) LATM + if (ios /= 0) goto 202 + END DO + DO WHILE (SLAT .LT. 0.D0) + WRITE(*,*) + & ' Seconds must also be positive. Please re-enter seconds.' + READ(LUIN,*,err=202,iostat=ios) SLAT + if (ios /= 0) goto 202 + END DO + ENDIF + + IF (LATD .LT. 0) THEN + IF ((LATM .GT. 0).OR.(SLAT .GT. 0.D0)) + & WRITE(*,*)' Latitude degrees are negative.' + DO WHILE (LATM .GT. 0) + WRITE(*,*) + & ' Minutes must also be negative. Please re-enter minutes.' + READ(LUIN,*,err=202,iostat=ios) LATM + if (ios /= 0) goto 202 + END DO + DO WHILE (SLAT .GT. 0.D0) + WRITE(*,*) + & ' Seconds must also be negative. Please re-enter seconds.' + READ(LUIN,*,err=202,iostat=ios) SLAT + if (ios /= 0) goto 202 + END DO + ENDIF + +*** Test whether latitude minute or second magnitudes are >= 60. +*** If either is, prompt for re-entry. + DO WHILE (ABS(LATM) .GE. 60) + WRITE(*,*)' Minutes magnitude must be < 60.', + & ' Please re-enter latitude minutes.' + READ(LUIN,*,err=202,iostat=ios) LATM + if (ios /= 0) goto 202 + END DO + + DO WHILE (DABS(SLAT) .GE. 60.D0) + WRITE(*,*)' Seconds magnitude must be < 60.', + & ' Please re-enter latitude seconds.' + READ(LUIN,*,err=202,iostat=ios) SLAT + if (ios /= 0) goto 202 + END DO + WRITE(LUOUT,120) 120 FORMAT( 1 ' Enter longitude degrees-minutes-seconds in free format'/, @@ -6684,30 +6842,112 @@ C The parameters in common block tranpa1 are computed using the IERS values of 3 ' eastward, enter a minus sign before each value.') READ(LUIN,*,err=203,iostat=ios) LOND, LONM, SLON if (ios /= 0) goto 203 - WRITE(LUOUT, 125) - 125 FORMAT( - 1 ' Enter ellipsoid height in meters. (Note that'/, - 1 ' predicted motions are independent of this height.) ') - READ(LUIN,*,err=204,iostat=ios) EHT - if (ios /= 0) goto 204 - XLAT = (DBLE((LATD*60 + LATM)*60) + SLAT)/RHOSEC - LATDIR = 'N' - IF (XLAT .lt. 0.0D0) then - LATD = - LATD - LATM = - LATM - SLAT = - SLAT - LATDIR = 'S' - ENDIF - XLON = (DBLE((LOND*60 + LONM)*60) + SLON)/RHOSEC - ELON = -XLON - CALL TOXYZ(XLAT,ELON,EHT,X,Y,Z) - LONDIR = 'W' - IF (XLON .lt. 0.0D0) then - LOND = -LOND - LONM = -LONM - SLON = -SLON - LONDIR = 'E' + +*** Test whether longitude magnitude is > 360 degrees. If it is, prompt +*** for re-entry. + TEST = DBLE(ABS(LOND)) + DBLE(ABS(LONM))/60.D0 + & + DABS(SLON)/3600.D0 + DO WHILE (DABS(TEST) .GT. 360.D0) + WRITE(*,*)' Longitude magnitude (degrees) =', TEST + WRITE(*,*) + & ' Magnitude cannot exceed 360. Please re-enter all values.' + READ(LUIN,*,err=203,iostat=ios) LOND, LONM, SLON + if (ios /= 0) goto 203 + TEST = DBLE(ABS(LOND)) + DBLE(ABS(LONM))/60.D0 + & + DABS(SLON)/3600.D0 + END DO + +*** Non-zero longitude minutes and seconds must be the same sign, and +*** if longitude degrees is non-zero, minutes and seconds must be the +*** same sign as longitude degrees. Otherwise prompt for re-entry. + IF (LOND .EQ. 0) THEN + IF (DBLE(LONM)*SLON .LT. 0.D0) THEN + WRITE(*,*)' Minutes and seconds must be the same sign.', + & ' Please re-enter minutes and seconds.' + READ(LUIN,*,err=203,iostat=ios) LONM, SLON + if (ios /= 0) goto 203 ENDIF + ENDIF + + IF (LOND .GT. 0) THEN + IF ((LONM .LT. 0).OR.(SLON .LT. 0.D0)) + & WRITE(*,*)' Longitude degrees are negative.' + DO WHILE (LONM .LT. 0) + WRITE(*,*) + & ' Minutes must also be positive. Please re-enter minutes.' + READ(LUIN,*,err=203,iostat=ios) LONM + if (ios /= 0) goto 203 + END DO + DO WHILE (SLON .LT. 0.D0) + WRITE(*,*) + & ' Seconds must also be positive. Please re-enter seconds.' + READ(LUIN,*,err=203,iostat=ios) SLON + if (ios /= 0) goto 203 + END DO + ENDIF + + IF (LOND .LT. 0) THEN + IF ((LONM .GT. 0).OR.(SLON .GT. 0.D0)) + & WRITE(*,*)' Longitude degrees are negative.' + DO WHILE (LONM .GT. 0) + WRITE(*,*) + & ' Minutes must also be negative. Please re-enter minutes.' + READ(LUIN,*,err=203,iostat=ios) LONM + if (ios /= 0) goto 203 + END DO + DO WHILE (SLON .GT. 0.D0) + WRITE(*,*) + & ' Seconds must also be negative. Please re-enter seconds.' + READ(LUIN,*,err=203,iostat=ios) SLON + if (ios /= 0) goto 203 + END DO + ENDIF + +*** Test whether longitude minute or second magnitudes are >= 60. +*** If either is, prompt for re-entry. + DO WHILE (ABS(LONM) .GE. 60) + WRITE(*,*)' Minutes magnitude must be < 60.', + & ' Please re-enter longitude minutes.' + READ(LUIN,*,err=203,iostat=ios) LONM + if (ios /= 0) goto 203 + END DO + + DO WHILE (DABS(SLON) .GE. 60.D0) + WRITE(*,*)' Seconds magnitude must be < 60.', + & ' Please re-enter longitude seconds.' + READ(LUIN,*,err=203,iostat=ios) SLON + if (ios /= 0) goto 203 + END DO + +*** END insertion of new code and comments for v3.3.0 to handle +*** incorrect input latitude and longitude. + + WRITE(LUOUT, 125) + 125 FORMAT( + 1 ' Enter ellipsoid height in meters. (Note that'/, + 1 ' predicted motions are independent of this height.) ') + READ(LUIN,*,err=204,iostat=ios) EHT + if (ios /= 0) goto 204 + + XLAT = (DBLE((LATD*60 + LATM)*60) + SLAT)/RHOSEC + LATDIR = 'N' + IF (XLAT .lt. 0.0D0) then + LATD = - LATD + LATM = - LATM + SLAT = - SLAT + LATDIR = 'S' + ENDIF + + XLON = (DBLE((LOND*60 + LONM)*60) + SLON)/RHOSEC + ELON = -XLON + CALL TOXYZ(XLAT,ELON,EHT,X,Y,Z) + LONDIR = 'W' + IF (XLON .lt. 0.0D0) then + LOND = -LOND + LONM = -LONM + SLON = -SLON + LONDIR = 'E' + ENDIF ELSEIF(COPT .EQ. '2') THEN WRITE(LUOUT,130) @@ -7027,9 +7267,9 @@ c write (*,*) "FROM NEWCOR ",YLAT3*180.d0/pi,YLON3*180.d0/pi,HTNEW IMPLICIT DOUBLE PRECISION (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) logical Is_iopt_NAD83 + COMMON /CONST/ A,F,E2,EPS,AF,PI,TWOPI,RHOSEC ** Get reference latitude (RLAT) and reference longitude (RLON) - C The following 2 lines were added on 07/22/2015 after Rich found this bug elon = -ylon call TOXYZ(ylat, elon, eht, x, y, z) @@ -7037,7 +7277,13 @@ C The following 2 lines were added on 07/22/2015 after Rich found this bug c IF(IOPT .EQ. 0 .OR. IOPT .EQ. 1) THEN !Velocity grids are No longer in NAD83 IF(IOPT .EQ. 15) THEN !They are in ITRF2008 RLAT = YLAT - RLON = YLON +*** Added conditional statement to get rid of out-of-region error when +*** YLON is negative for reference frame (v3.3.0): + IF(YLON .GE. 0.D0) THEN + RLON = YLON + ELSE + RLON = YLON + TWOPI + ENDIF ELSE c elon = -ylon !Was a bug, commented out on 07222015 c call TOXYZ(ylat, elon, eht, x, y, z) !Was a bug, commented out on 07222015 @@ -7085,7 +7331,7 @@ c write (*,*) "From PREDV ",VN, VE, VU implicit integer*4 (i-n) parameter (numref = 16) character nameif*80,name24*80 - character namef*30 + character NAMEF*80 character frame1*24, frame2*24 character option*1 character vopt*1, LATDIR*1,LONDIR*1 @@ -7098,7 +7344,7 @@ c write (*,*) "From PREDV ",VN, VE, VU 100 format( 1 ' Please enter the name of the file to contain '/ 1 ' the transformed velocities. ') - read( luin, '(a30)',err=600,iostat=ios) namef + read( luin, '(A80)',err=600,iostat=ios) namef if (ios /= 0) goto 600 open( i2, file = namef, status = 'unknown') @@ -7176,11 +7422,27 @@ c write (*,*) "From PREDV ",VN, VE, VU if (ios /= 0) goto 602 open (i1, file = nameif, status = 'old') + LINE = 0 !Inititalize line number of input file. 210 read(i1,'(a)',end = 220,err=603,iostat=ios) record -c 210 read(i1, *, end = 220,err=603,iostat=ios) xlat, xlon,vn,ve,vu, name24 - call interprate_velocity_record (record,xlat,xlon,vn,ve, - & vu,name24) + call interprate_velocity_record (record,xlat,xlon,vn,ve,vu, + & name24) if (ios /= 0) goto 603 + LINE = LINE + 1 !Increment line number of input file. + +*** If latitude magnitudes > 90 degrees or longitude magnitude > 360 degrees, +*** write error message and terminate program (added for v3.3.0). + IF ((DABS(xlat) .GT. 90).OR.(DABS(xlon) .GT. 360)) THEN + WRITE(*,*)'***********************************************' + WRITE(*,*)'Invalid latitude or longitude in input file' + WRITE(*,2101)'on line', LINE, ':' + WRITE(*,2102) xlat, xlon, vn, ve, vu, name24 + WRITE(*,*)'Please check your input file and try again.' + WRITE(*,*)'***********************************************' + 2101 FORMAT (1X, A, I8, A/) + 2102 FORMAT (1X, 2F17.10, 3F8.2, 2X, A) + STOP + ENDIF + ylat = (xlat*3600.d0) / rhosec ylon = (xlon*3600.d0) / rhosec elon = -ylon @@ -7499,10 +7761,8 @@ C The parameters in common block tranpa1 are computed using the IERS values of TEST = .TRUE. RETURN ENDIF -**** add small increment to circumvent round-off error -c DATE = DATE + 0.0004D0 -**** - IYEAR = DATE + + IYEAR = IDINT(DATE) !Convert DATE to integer (v3.3.0) CALL IYMDMJ(IYEAR, 1, 1, MJD0) IYEAR1 = IYEAR + 1 CALL IYMDMJ(IYEAR1, 1, 1, MJD1) @@ -7513,7 +7773,7 @@ c DATE = DATE + 0.0004D0 ITOTAL = 31 IF (REMDAY .LT. ITOTAL) then MONTH = 1 - IDAY = REMDAY - IBEGIN + 1 + IDAY = IDINT(REMDAY) - IBEGIN + 1 !Convert REMDAY to integer (v3.3.0) CALL IYMDMJ(IYEAR,MONTH,IDAY,MJD) MINS = MJD * 24 * 60 TEST = .FALSE. @@ -7525,7 +7785,7 @@ c DATE = DATE + 0.0004D0 ITOTAL = ITOTAL + M(I) IF (REMDAY .LT. ITOTAL) then MONTH = I - IDAY = REMDAY - IBEGIN + 1 + IDAY = IDINT(REMDAY) - IBEGIN + 1 !Convert REMDAY to integer (v3.3.0) TEST = .FALSE. CALL IYMDMJ(IYEAR,MONTH,IDAY,MJD) MINS = MJD * 24 * 60 @@ -7636,13 +7896,30 @@ C END IF C C........ 1.0 CALCULATION -C - A= IYRP*0.01D0 - B= 2 - A + DINT( A*0.25D0 ) - C= 365.25D0*IYRP - D= 30.6001D0*(IMOP + 1) - MJD = (B + C + D + IDAY - 679006) +C +*** Following 4 lines edited in v3.3.0 to address compiler issue with +*** real to integer conversion. +C A= IYRP*0.01D0 +C B= 2 - A + DINT( A*0.25D0 ) +C C= 365.25D0*IYRP +C D= 30.6001D0*(IMOP + 1) +C MJD = (B + C + D + IDAY - 679006) C +C WRITE(*,*)'Old values:' +C WRITE(*,*)'A, B, C, D =', A, B, C, D +C WRITE(*,*)'MJD =', MJD + + A = IDINT(DBLE(IYRP)*0.01D0) + B = 2 - A + IDINT(DBLE(A)*0.25D0) + C = IDINT(365.25D0*DBLE(IYRP)) + D = IDINT(30.6001D0*(DBLE(IMOP + 1))) + MJD = B + C + D + IDAY - 679006 + +C WRITE(*,*)'New values:' +C WRITE(*,*)'A, B, C, D =', A, B, C, D +C WRITE(*,*)'MJD =', MJD +C WRITE(*,*) + RETURN END ***************************************************** @@ -8511,7 +8788,7 @@ C Done return end -C----------------------------------------------------------------------------------- +C----------------------------------------------------------------------- C On 5/10/2020 removed all code and comments entered on 7/22/2015 that limited use C of HTDP to US and its territories. That functionality was disabled on 5/11/2017. C The removed code still exists in HTDP v3.2.8 and 3.2.7, commented out to disable. -- cgit v1.2.3