//EIQ2GPBS JOB TIME=(0,05) CONVERT STRESS DATA FROM LAT,LON TO X,Y /*SCHEDULE PRIORITY=1.5 //CLEAR EXEC PGM=IEFBR14 //TRASH DD DISP=(OLD,DELETE),DSN=EIQ2GPB.STRESSES.ALXY // EXEC FORTCLG,LIBE='APP1.FLIB', // OPTION='FIPS(F),ICA(CVAR,MSGOFF(84)),', // OPTIONS='MAP,NOOPTIMIZE,NOVEC,SDUMP(ISN),TEST,XREF,NOIL', // TYPE=LINKLIB //FORT.SYSIN DD * C PROGRAM STRLL2XY C C READS NORTH AMERICAN STRESS DATA FILE FROM THE C GEOPHYSICS OF NORTH AMERICA COMPACT DISK C (BY MARY LOU ZOBACK) C AND CONVERTS (LATITUDE,LONGITUDE) -> (X,Y) IN THE FINITE ELEMENT C MODEL PLANE (ASSUMING A CONFORMAL OR SIMPLE CONIC PROJECTION), C AND CONVERTS AZIMITH OF COMPRESSION (CLOCKWISE, FROM N) C TO ARGUMENT OF COMPRESSION (COUNTERCLOCKWISE, FROM +X). C C ONLY VALUES WITHIN A REASONABLE DISTANCE (PARAMETER RING) C OF THE NEW (X,Y) ORIGIN ARE OUTPUT. PARAMETER (RING=3000.E3) C C ALASKA VERSION: DIGITIZED TRENCH AND VOLCANIC ARC ARE USED TO C SEPARATE OUT EVENTS FROM NORTH AMERICAN PLATE ONLY. C C NINARC = NUMBER OF POINTS DEFINING TRENCH AND ARC PARAMETER (NINARC=16) C C LOGICAL SENDIT C CHARACTER*6 SITE CHARACTER*2 QUALIT CHARACTER*6 TYPE CHARACTER*3 AZIMUT CHARACTER*4 REGIME CHARACTER*7 DEPTH CHARACTER*28 REFERE,MORERE CHARACTER*59 TEST,BLANKS C DIMENSION XYARC(3,NINARC) DATA ((XYARC(I,J),I=1,3),J=1,NINARC) / + -2.687E+06,-1.034E+05,+2.382E+05, + -2.564E+06,-3.140E+05,-9.064E+04, + -2.342E+06,-5.558E+05,-3.474E+05, + -2.087E+06,-6.976E+05,-5.215E+05, + -1.680E+06,-7.804E+05,-6.322E+05, + -1.419E+06,-7.911E+05,-6.290E+05, + -1.065E+06,-7.600E+05,-5.516E+05, + -8.055E+05,-7.374E+05,-4.764E+05, + -6.005E+05,-6.966E+05,-3.647E+05, + -3.846E+05,-5.816E+05,-2.185E+05, + -2.797E+05,-5.161E+05,+5.439E+04, + -2.046E+05,-4.898E+05,+2.306E+05, + -1.131E+05,-4.350E+05,+3.241E+05, + +1.452E+05,-2.029E+05,+1.285E+05, + +2.954E+05,-2.222E+05,+3.182E+04, + +6.240E+05,-3.286E+05,+3.874E+03 + / C C NOTE: CPNLA2.GE.CPNLA1. C BOTH ARE CONIC PROJECTION PARALLELS IN DEGREES NORTH. C FOR CONFORMAL CONIC, THEIR VALUES DIFFER. C FOR SIMPLE CONIC, THEY ARE THE SAME. DATA CPNLA1/61.000/ DATA CPNLA2/61.000/ DATA RADIUS /6371.E+03/ DATA Y0NLAT/61.0000/ DATA X0ELON/-147.0000/ C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C BLANKS=' ' DO 1 I=2,59 BLANKS=BLANKS//' ' 1 CONTINUE C C 1 = INPUT FILE C 2 = OUTPUT FILE C C C READ ACCORDING TO ZOBACK'S FORMAT C (NOTE: AZIMUTH READ AS TEXT BECAUSE IT IS LEFT-JUSTIFIED INTEGER]) C 10 READ (1,11,ERR=30,END=101) SITE,PLAT,PLON,QUALIT,TYPE,AZIMUT, + REGIME,DEPTH,REFERE 11 FORMAT (A6,1X,F6.3,2X,F8.3,1X,A2,3X,A6,3X,A3,2X, + A4,3X,A7,2X,A28) READ (DEPTH,'(F7.0)') DEEP DEEP=DEEP*1000. C C CONVERT POSITION TO (XT,YT) IN A CONFORMAL CONIC PROJECTION C 30 CALL LLTOXY (INPUT,CPNLA1,CPNLA2, + PLAT,PLON, + RADIUS,X0ELON,Y0NLAT, + OUTPUT,X,Y) C C DECIDE WHETHER TO KEEP THIS POINT: C IF ((X .LT.XYARC(1,1)).OR. + (X .GT.XYARC(1,NINARC))) THEN FARC=99. ELSE IF (X .LE.-1.54E6) THEN ACRUST=20.E3 ELSE IF (X .LE.-1.8E5) THEN ACRUST=20.E3+20.E3*(X -(-1.54E6))/1.36E6 ELSE IF (X .LE.1.7E4) THEN ACRUST=40.E3+15.E3*(X -(-1.8E5))/1.97E5 ELSE IF (X .LE.6.2E5) THEN ACRUST=55.E3-10.E3*(X -1.7E4)/6.03E5 ELSE ACRUST=45.E3 ENDIF DO 190 K=2,NINARC IF ((X .GE.XYARC(1,K-1)).AND. + (X .LE.XYARC(1,K))) THEN K1=K-1 K2=K FRAC=(X -XYARC(1,K1))/ + (XYARC(1,K2)-XYARC(1,K1)) GO TO 191 ENDIF 190 CONTINUE 191 YTRE=XYARC(2,K1)+FRAC*(XYARC(2,K2)-XYARC(2,K1)) YARC=XYARC(3,K1)+FRAC*(XYARC(3,K2)-XYARC(3,K1)) FARC=(Y -YTRE)/(YARC-YTRE) ENDIF IF (FARC.LE.1.5) THEN C C PACIFIC PLATE (REJECT POINT): C IF (FARC.LE.0.0) THEN PLATE= -1000. C C ARC PLATE THICKNESS: C ELSE IF (FARC.LT.0.8) THEN PLATE=10.E3+30.E3*FARC+50.E3*FARC**3 ELSE IF (FARC.LT.1.2) THEN PLATE=ACRUST+14.E3 ELSE PLATE=65.E3 ENDIF ENDIF SENDIT=((ABS(X).LE.RING).AND.(ABS(Y).LE.RING)) .AND. + (AZIMUT(2:2).NE.'-').AND.(DEEP.LT.PLATE) C C CONVERT AZIMUT TO FLOATING-POINT NUMBER, IN DEGREES C IF (SENDIT) THEN READ (AZIMUT,*) IAZ AZ=IAZ C C PERFORM SAME OPERATIONS ON TAIL OF SIGMA1 VECTOR, WITH LENGTH 1 DEGREE C TAILNL=PLAT+COSDEG(AZ) TAILEL=PLON+SINDEG(AZ)/COSDEG(PLAT) CALL LLTOXY (INPUT,CPNLA1,CPNLA2, + TAILNL,TAILEL, + RADIUS,X0ELON,Y0NLAT, + OUTPUT,X2,Y2) ARG=57.298*ATAN2((Y2-Y),(X2-X)) IF (ARG.LT.0.) ARG=ARG+180. WRITE(2,21) SITE, X, Y, QUALIT, TYPE, ARG, REGIME, + DEPTH, REFERE 21 FORMAT (A6,1X,F10.0,F10.0,1X,A2,3X,A6,2X,F4.0,2X,A4,3X, + A7,2X,A28) ENDIF C C READ ANY EXTENSIONS OF THE REFERENCE ONTO ANOTHER RECORD C 40 READ (1,12,END=101) TEST,MORERE 12 FORMAT (A59,A28) IF (TEST.EQ.BLANKS) THEN IF (SENDIT) WRITE (2,22) BLANKS,MORERE 22 FORMAT (A63,A28) GO TO 40 ELSE BACKSPACE (1) ENDIF C C LOOP C GO TO 10 101 CONTINUE C STOP END C C C SUBROUTINE LLTOXY (INPUT,CPNLA1,CPNLA2, + PLAT,PLON, + RADIUS,X0ELON,Y0NLAT, + OUTPUT,X,Y) C C CONVERT A (NORTH LATITUDE=PLAT, EAST LONGITUDE=PLON) POSITION C IN UNITS OF DEGREES C INTO AN (X,Y) POSITION ON A CONFORMAL CONIC PROJECTION WITH C STANDARD LATITUDES CPNLA1 AND CPNLA2, WHEN THE (X,Y) ORIGIN IS AT C (NORTH LATITUDE=Y0NLAT, EAST LONGITUDE=X0ELON). C THE CUT NECESSARY IN THIS PROJECTION IS FROM THE POLE NEAREST C TO THE STANDARD LATITUDE CPNLA2, ALONG A MERIDIAN WHICH C IS ON THE OPPOSITE SIDE OF THE EARTH FROM X0ELON. C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C PHI0=0.5*(CPNLA1+CPNLA2) SINPH0=SINDEG(PHI0) C= ( RADIUS*COSDEG(CPNLA1) ) / + ( SINPH0*(TANDEG(45.-0.5*CPNLA1))**SINPH0 ) RHO0=C*(TANDEG(45.-0.5*Y0NLAT))**SINPH0 RHO=C*(TANDEG(45.-0.5*PLAT))**SINPH0 IF (PLON-X0ELON.GT.180.) THEN PLONP=PLON-360. ELSE IF (X0ELON-PLON.GT.180.) THEN PLONP=PLON+360. ELSE PLONP=PLON ENDIF ANGLE=(PLONP-X0ELON)*SINPH0 X=RHO*SINDEG(ANGLE) Y=RHO0-RHO*COSDEG(ANGLE) RETURN END //GO.FT01F001 DD DISP=SHR,DSN=EIQ2GPB.STRESSES //GO.FT02F001 DD DISP=(NEW,CATLG),UNIT=DATA, // SPACE=(TRK,(1,1),RLSE),DCB=(RECFM=FB,LRECL=91,BLKSIZE=9100), // DSN=EIQ2GPB.STRESSES.ALXY //