C C PROGRAM STRXY2IS C C READS NORTH AMERICAN STRESS DATA FILE ORIGINALLY FROM THE C GEOPHYSICS OF NORTH AMERICA COMPACT DISK, BUT PREPROCESSED BY C PROGRAM "STRLL2XY" (WHICH PERFORMS MAP PROJECTION AND WINDOWS). C C THIS PROGRAM CONVERTS (X,Y) TO THE TRIANGULAR CONTINUUM ELEMENT C NUMBER (IE) AND THE 3 INTERNAL COORDINATE VALUES (S1,S2,S3). C POINTS WHICH DO NOT LIE IN ANY ELEMENT ARE MERELY DROPPED FROM C THE OUPUT FILE. C C---------------------------------------------------------------------- C C PARAMETER STATEMENTS C (TO SET MAXIMUM ARRAY SIZES AND DIMENSIONS) C C MAXNOD = MAXIMUM NUMBER OF NODES (INCLUDES BOTH "REAL"AND & "FAKE") PARAMETER (MAXNOD=2000) C C MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES). PARAMETER (MAXEL=800) C C MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS); PARAMETER (MAXFEL=200) C C MAXATP = MAXIMUM NUMBER OF NODES WHICH MAY OVERLAP AT A FAULT- C INTERSECTION POINT. PARAMETER (MAXATP=20) C---------------------------------------------------------------------- C C TYPE STATEMENTS C (NOTE: IMPLICIT CONVENTION OF I-N = INTEGER IS FOLLOWED) C C CHARACTER*80 TITLE1 CHARACTER*6 SITE CHARACTER*2 QUALIT CHARACTER*6 TYPE CHARACTER*4 REGIME CHARACTER*7 DEPTH CHARACTER*28 REFERE,MORERE CHARACTER*59 TEST,BLANKS C C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: DOUBLE PRECISION PHI,POINTS,WEIGHT C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: DOUBLE PRECISION FPHI,FPOINT,FGAUSS C C NOTE: THE FOLLOWING CAN BE MADE "INTEGER*2" IN VS-FORTRAN: INTEGER NODTYP C LOGICAL BRIEF LOGICAL ATSEA,SENDIT C C NOTE: THE FOLLOWING ARRAYS COULD BE COMPRESSED WITH "LOGICAL*1" C IN VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN C C---------------------------------------------------------------------- C C DIMENSION STATEMENTS C C DIMENSIONS USING PARAMETER MAXNOD: DIMENSION CHECKN(MAXNOD), DQDTDA(MAXNOD), + ELEV (MAXNOD), + NODTYP(MAXNOD), + XNODE (MAXNOD), YNODE (MAXNOD), ZMNODE(MAXNOD) C C DIMENSIONS USING PARAMETER MAXEL: DIMENSION AREA (MAXEL), CHECKE (MAXEL), + DETJ (7,MAXEL), + DXS (6,7,MAXEL), DYS (6,7,MAXEL), + NODES (6,MAXEL) C C DIMENSIONS USING PARAMETER MAXFEL: DIMENSION CHECKF (MAXFEL), + FAZ (2,MAXFEL), FDIP (3,MAXFEL), + FLEN (MAXFEL), + FTAN (7,MAXFEL), NODEF (6,MAXFEL) C C DIMENSIONS USING PARAMETER MAXATP: DIMENSION LIST (MAXATP) C C DIMENSIONS OF FIXED SIZE: C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: DIMENSION PHI(6,7),POINTS(5,7),WEIGHT(7) C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: DIMENSION FPHI(6,7),FPOINT(7),FGAUSS(7) C-------------------------------------------------------------- C C COMMON STATEMENTS C C NAMED COMMON BLOCKS HOLD THE FIXED VALUES OF THE POSITIONS, C WEIGHTS, AND NODAL FUNCTION VALUES AT THE INTEGRATION POINTS C IN THE ELEMENTS (TRIANGULAR ELEMENTS IN BLOCK DATA BD1, AND C FAULT ELEMENTS IN BLOCK DATA BD2). C C ENTRIES CORRESPONDING TO BD1: COMMON /S1S2S3/ POINTS COMMON /PHITAB/ PHI COMMON /WGTVEC/ WEIGHT C C ENTRIES CORRESPONDING TO BD2: COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS C C---------------------------------------------------------------------- C C DATA STATEMENTS C C FILES IN USE, BY FORTRAN DEVICE NUMBER: C 1 = FINITE ELEMENT GRID FILE C 2 = INPUT STRESS FILE C 3 = OUTPUT STRESS FILE C 6 = ERR0R MESSAGES, ETC. C DATA IUNITG /1/ DATA IUNITI /2/ DATA IUNITO /3/ DATA IUNITT /6/ C C "DIPMAX" IS THE MAXIMUM DIP (FROM HORIZONTAL, IN DEGREES) FOR A C FAULT ELEMENT TO BE TREATED AS A DIP-SLIP FAULT, WITH TWO DEGREES C OF FREEDOM PER NODE-PAIR. AT STEEPER DIPS, THE DEGREE OF FREEDOM C CORRESPONDING TO OPENING OR CONVERGENCE OF THE OPPOSITE SIDES IS C ELIMINATED BY A CONSTRAINT EQUATION, AND THE FAULT IS TREATED AS C A VERTICAL STRIKE-SLIP FAULT. THIS ARBITRARY LIMIT IS NECESSARY C BECAUSE THE EQUATIONS FOR DIP-SLIP FAULTS BECOME SINGULAR AS THE C DIP APPROACHES 90 DEGREES. IN PRACTICE, IT IS BEST TO SPECIFY DIPS C AS EITHER (1) VERTICAL, OR (2) CLEARLY LESS THAN "DIPMAX", WITHIN C EACH FAULT ELEMENT. IF THE DIP VARIES WITHIN AN ELEMENT IN SUCH A C WAY THAT IT PASSES THROUGH THIS LIMIT WITHIN THE ELEMENT, THEN C THE REPRESENTATION OF THAT FAULT ELEMENT IN THE EQUATIONS MAY C BE INACCURATE. DATA DIPMAX /75./ C C "N1000" IS JUST A CONVENIENT "ROUND-NUMBER" INTEGER THAT IS GREATER C THAN OR EQUAL TO THE NUMBER OF REAL NODES. ON INPUT AND ON OUTPUT, C THE NODE NUMBERS OF FAKE NODES (BOUNDARY NODES WITH BOTH COMPONENTS C OF VELOCITY FIXED, AND NO DEGREES OF FREEDOM LEFT) WILL BEGIN C WITH N1000+1 AND GO UP. (HOWEVER, INTERNALLY, THE GAP IN NODE C NUMBERS IS REMOVED FOR PROGRAMMING CONVENIENCE.) DATA N1000 /30000/ C-------------------------------------------------------------- C C STATEMENT FUNCTIONS: C C C BEGIN EXECUTABLE CODE: C C C *** KLUDGE ALERT ************************************************* C CONVERSION OF PARAMETERS (CONSTANTS) TO VARIABLES SHOULD LOGICALLY C HAVE NO EFFECT, BUT IN FACT HELPS TO SUPPRESS SOME SPURIOUS C MESSAGES FROM THE IBM VS-FORTRAN COMPILER. MXNODE=MAXNOD MXEL =MAXEL MXFEL =MAXFEL MXSTAR=MAXATP C ****************************************************************** WEDGE=ABS(90.-ABS(DIPMAX))*0.017453293 C BLANKS=' ' DO 1 I=2,59 BLANKS=BLANKS//' ' 1 CONTINUE C C INPUT FINITE ELEMENT GRID AND DATA VALUES AT NODE POINTS C CALL GETNET (INPUT,IUNITG,IUNITT, + MXEL,MXFEL,MXNODE,N1000, + OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD, + TITLE1,XNODE,YNODE,ZMNODE, + WORK,CHECKE,CHECKF,CHECKN) C C CHECK GRID TOPOLOGY AND COMPUTE GEOMETRIC PROPERTIES C CALL SQUARE (INPUT,BRIEF,FDIP,IUNITT, + MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,WEDGE, + MODIFY,FAZ,XNODE,YNODE, + OUTPUT,AREA,DETJ,DXS,DYS, + FLEN,FTAN, + WORK,CHECKN,LIST,NODTYP) C C C READ DATUM ACCORDING TO MODIFIED ZOBACK FORMAT C (ARG = ARGUMENT OF HORIZONTAL GREATEST COMPRESSION (T1), IN C DEGREES COUNTERCLOCKWISE FROM +X.) C 10 READ (IUNITI,11,ERR=30,END=101) SITE, X, Y, QUALIT, TYPE, ARG, + REGIME, DEPTH, REFERE 11 FORMAT (A6,1X,2F10.0,1X,A2,3X,A6,2X,F4.0,2X,A4,3X, + A7,2X,A28) C C FIND TRIANGULAR CONTINUUM ELEMENT WHOSE CENTROID IS CLOSEST, AND C USE THIS TO INITIALIZE THE SEARCH: C IE=1 RMIN=9.99E59 DO 20 I=1,NUMEL N1=NODES(1,I) N2=NODES(2,I) N3=NODES(3,I) X1=XNODE(N1) X2=XNODE(N2) X3=XNODE(N2) Y1=YNODE(N1) Y2=YNODE(N2) Y3=YNODE(N2) XMID=(X1+X2+X3)/3. YMID=(Y1+Y2+Y3)/3. R=SQRT((X-XMID)**2+(Y-YMID)**2) IF (R.LT.RMIN) THEN IE=I RMIN=R ENDIF 20 CONTINUE IESAVE=IE S1=0.333 S2=0.333 S3=0.334 C 30 CALL LOOKUP (INPUT,DETJ,IUNITT,MXEL,MXFEL,MXNODE, + NFL,NODEF,NODES,NUMEL, + X,XNODE,Y,YNODE, + MODIFY,IE,S1,S2,S3, + OUTPUT,ATSEA) SENDIT=.NOT.ATSEA IF (SENDIT) THEN WRITE (IUNITO,21) SITE, IE, S1, S2, QUALIT, + TYPE, ARG, REGIME, DEPTH, REFERE 21 FORMAT (A6,1X,I4,2F6.3,1X,A2,3X,A6,2X,F4.0,2X,A4,3X, + A7,2X,A28) ELSE WRITE (IUNITT,25) SITE, X, Y, IESAVE, IE, S1, S2, S3 25 FORMAT (/' DATUM ',A6,' AT X = ',1P,E12.4,', Y = ',E12.4/ + ' WAS SOUGHT NEAR ELEMENT ',I5/ + ' BUT FELL OFF GRID EDGE AT IE =',I5, + 0P,', S1 = ',F6.3,', S2 = ',F6.3,', S3 = ',F6.3) ENDIF C C READ ANY EXTENSIONS OF THE REFERENCE ONTO ANOTHER RECORD C 40 READ (IUNITI,12,END=101) TEST,MORERE 12 FORMAT (A59,A28) IF (TEST.EQ.BLANKS) THEN IF (SENDIT) WRITE (IUNITO,22) BLANKS,MORERE 22 FORMAT (A59,A28) GO TO 40 ELSE BACKSPACE (IUNITI) ENDIF C C LOOP C GO TO 10 C C EXIT C 101 CONTINUE STOP END C C C SUBROUTINE LOOKUP (INPUT,DETJ,IUNITT,MXEL,MXFEL,MXNODE, + NFL,NODEF,NODES,NUMEL, + X,XNODE,Y,YNODE, + MODIFY,IE,S1,S2,S3, + OUTPUT,ATSEA) C C FINDS ELEMENT AND INTERNAL COORDINATES IN GRID MATCHING LOCATION C OF A PARTICULAR POINT (X,Y) AND REPORTS THEM AS IE AND S1,S2,S3. C C THE INITIAL VALUES OF THESE VARIABLES DO NOT NEED TO BE ACCURATE, C BUT IT WILL SAVE TIME IF THEY ARE. IF YOU ARE LOCATING A GROUP C OF POINTS THAT ARE NOT FAR APART, JUST LET THE COORDINATES OF THE C LAST POINT PROCESSED BE THE INITIAL ESTIMATE FOR THE NEXT. C C A RETURNED VALUE OF ATSEA INDICATES THAT POINT FELL OFF EDGE C OF THE GRID. C PARAMETER (NTOTRY=50) LOGICAL ATSEA,ISTRAP,TRUBBL REAL M11,M12,M13,M21,M22,M23 DIMENSION DETJ(7,MXEL), + NODEF(6,MXFEL),NODES(6,MXEL), + XNODE(MXNODE),YNODE(MXNODE) DIMENSION IEHIST(NTOTRY),SHIST(3,NTOTRY) C C STATEMENT FUNCTION: C PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6)= + F1*(-S1+2.*S1**2)+ + F2*(-S2+2.*S2**2)+ + F3*(-S3+2.*S3**2)+ + F4*(4.*S1*S2)+ + F5*(4.*S2*S3)+ + F6*(4.*S3*S1) C C NTRIED=0 C C LOOP AS MANY TIMES AS NEEDED: C 100 NTRIED=NTRIED+1 IEHIST(NTRIED)=IE TRUBBL=(NTRIED.GE.(NTOTRY-10)).AND.(IEHIST(NTRIED).EQ. + IEHIST(NTRIED-2)) IF (TRUBBL) THEN WRITE (IUNITT,101) X, Y, NTRIED, + (IEHIST(K),K=1,NTRIED) 101 FORMAT (/' POINT AT X = ',1P,E12.4,', Y = ',E12.4/ + ' CAUSES TROUBLE AFTER SEARCHING ',I3, + ' ELEMENTS:',99(/' ',I5)) ATSEA=.TRUE. RETURN ENDIF I1=NODES(1,IE) I2=NODES(2,IE) I3=NODES(3,IE) I4=NODES(4,IE) I5=NODES(5,IE) I6=NODES(6,IE) X1=XNODE(I1) X2=XNODE(I2) X3=XNODE(I3) Y1=YNODE(I1) Y2=YNODE(I2) Y3=YNODE(I3) ISTRAP=(DETJ(1,IE).LE.0.2).OR. + (DETJ(2,IE).LE.0.2).OR. + (DETJ(3,IE).LE.0.2).OR. + (DETJ(4,IE).LE.0.2).OR. + (DETJ(5,IE).LE.0.2).OR. + (DETJ(6,IE).LE.0.2).OR. + (DETJ(7,IE).LE.0.2) IF (ISTRAP) THEN X4=0.5*(X1+X2) X5=0.5*(X2+X3) X6=0.5*(X3+X1) Y4=0.5*(Y1+Y2) Y5=0.5*(Y2+Y3) Y6=0.5*(Y3+Y1) ELSE X4=XNODE(I4) X5=XNODE(I5) X6=XNODE(I6) Y4=YNODE(I4) Y5=YNODE(I5) Y6=YNODE(I6) ENDIF S3=1.-S1-S2 LIMIT=3 NREFIN=0 C C LOOP TO REFINE ESTIMATE OF INTERNAL COORDINATES C 150 NREFIN=NREFIN+1 XT=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) YT=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) C C COEF:=MAT((DXDS1,DXDS2,DXDS3), C (DYDS1,DYDS2,DYDS3),(1,1,1)); C COEF11=4.*S2*X4+4.*S1*X1+4.*X6*S3-X1 COEF12=4.*S2*X2+4.*S1*X4+4.*X5*S3-X2 COEF13=4.*S2*X5+4.*S1*X6+4.*X3*S3-X3 COEF21=4.*S2*Y4+4.*S1*Y1+4.*Y6*S3-Y1 COEF22=4.*S2*Y2+4.*S1*Y4+4.*Y5*S3-Y2 COEF23=4.*S2*Y5+4.*S1*Y6+4.*Y3*S3-Y3 M11=COEF22-COEF23 M12=COEF21-COEF23 M13=COEF21-COEF22 M21=COEF12-COEF13 M22=COEF11-COEF13 M23=COEF11-COEF12 CF11=+M11 CF12=-M12 CF13=+M13 CF21=-M21 CF22=+M22 CF23=-M23 DET=COEF11*CF11+COEF12*CF12+COEF13*CF13 IF (DET.EQ.0.0) THEN ATSEA=.TRUE. RETURN ENDIF DETI= 1./DET STEP11=CF11 STEP12=CF21 STEP21=CF12 STEP22=CF22 STEP31=CF13 STEP32=CF23 DELX=X-XT DELY=Y-YT DS1=(STEP11*DELX+STEP12*DELY)*DETI DS2=(STEP21*DELX+STEP22*DELY)*DETI DS3=(STEP31*DELX+STEP32*DELY)*DETI ERR=(DS1+DS2+DS3)/3. DS1=DS1-ERR DS2=DS2-ERR DS3=DS3-ERR DSTEP=MAX(ABS(DS1),ABS(DS2),ABS(DS3)) IF (DSTEP.GT. 0.10) THEN LIMIT=LIMIT+1 DS1=DS1*0.1/DSTEP DS2=DS2*0.1/DSTEP DS3=DS3*0.1/DSTEP ENDIF S1=S1+DS1 S2=S2+DS2 S3=S3+DS3 C C LOOP-BACK (WITH SOME CONDITIONS): C IF ((NREFIN.LT.LIMIT.AND.LIMIT.LE.(NTOTRY-10)).AND. + (S1.GE.-0.1.AND.S1.LE.1.1).AND. + (S2.GE.-0.1.AND.S2.LE.1.1).AND. + (S3.GE.-0.1.AND.S3.LE.1.1)) GO TO 150 C C POINT IS NOW AS WELL-LOCATED AS POSSIBLE "IN" THE CURRENT ELEMENT; C HOWEVER, THE INTERNAL COORDINATES MAY NOT ALL BE POSITIVE, SO C POINT MAY BE OUTSIDE, AND WE MAY NEED TO SHIFT TO A NEW ELEMENT. C SHIST(1,NTRIED)=S1 SHIST(2,NTRIED)=S2 SHIST(3,NTRIED)=S3 IF (TRUBBL.OR.NTRIED.GE.NTOTRY) THEN WRITE(IUNITT,201) X,Y 201 FORMAT(' REQUEST FOR VALUE AT LOCATION', + ' (',1P,E10.2,',',E10.2,') CAUSES ', + 'INFINITE LOOP IN LOOKUP.'/ + ' HISTORY OF SEARCH: ELEMENT S1 S2', + ' S3') DO 203 N=1,NTRIED-1 WRITE(IUNITT,202) IEHIST(N),(SHIST(K,N),K=1,3) 202 FORMAT(22X,I3,2X,3F12.4) 203 CONTINUE WRITE(IUNITT,204) IEHIST(NTRIED-1), + (NODES(J,IEHIST(NTRIED-1)),J=1,6), + (XNODE(NODES(J,IEHIST(NTRIED-1))),J=1,6), + (YNODE(NODES(J,IEHIST(NTRIED-1))),J=1,6) WRITE(IUNITT,204) IEHIST(NTRIED), + (NODES(J,IEHIST(NTRIED)),J=1,6), + (XNODE(NODES(J,IEHIST(NTRIED))),J=1,6), + (YNODE(NODES(J,IEHIST(NTRIED))),J=1,6) 204 FORMAT(' ELEMENT',I3,' NODES:',I3,5I10/ + 9X,'X:',1P,6E10.2/9X,'Y:',6E10.2) RETURN ENDIF IF (S1.GT.-0.03) THEN IF (S2.GT.-0.03) THEN IF (S3.GT.-0.03) THEN C POINT HAS BEEN SUCCESSFULLY FOUND] ATSEA=.FALSE. RETURN ELSE CALL NEXT (INPUT,IE,3,MXEL,MXFEL,NFL, + NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) ENDIF ELSE CALL NEXT (INPUT,IE,2,MXEL,MXFEL,NFL, + NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) ENDIF ELSE CALL NEXT (INPUT,IE,1,MXEL,MXFEL,NFL, + NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) ENDIF IF (KELE.GT.0) THEN IE=KELE S1=0.3333 S2=0.3333 S3=0.3334 GO TO 100 ELSE ATSEA=.TRUE. RETURN ENDIF C C NOTE: INDENTATION REFLECTS INDEFINITE LOOP ON TRIAL ELEMENT IE. C END C C C SUBROUTINE NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) C C DETERMINE WHETHER THERE ARE MORE ELEMENTS ADJACENT TO SIDE J OF C TRIANGULAR CONTINUUM ELEMENT #I. C J = 1 MEANS THE SIDE OPPOSITE NODE # NODES(1,I). C J = 2 MEANS THE SIDE OPPOSITE NODE # NODES(2,I). C J = 3 MEANS THE SIDE OPPOSITE NODE # NODES(3,I). C IF A FAULT ELEMENT IS ADJACENT, ITS NUMBER IS KFAULT; OTHERWISE, C KFAULT IS SET TO ZERO. C IF ANOTHER TRIANGULAR CONTINUUM ELEMENT IS ADJACENT (EVEN ACROSS C FAULT ELEMENT KFAULT]) THEN ITS NUMBER IS KELE; OTHERWISE, KELE = 0. C LOGICAL FOUNDE,FOUNDF DIMENSION NODEF(6,MXFEL),NODES(6,MXEL) C C THREE NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE: N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J, 3)+4,I) N3=NODES(MOD(J+1,3)+1,I) C CHECK FOR ADJACENT FAULT ELEMENT FIRST: FOUNDF=.FALSE. DO 10 K=1,NFL M1=NODEF(1,K) M2=NODEF(2,K) M3=NODEF(3,K) M4=NODEF(4,K) M5=NODEF(5,K) M6=NODEF(6,K) IF (((M1.EQ.N3).AND.(M2.EQ.N2).AND.(M3.EQ.N1)).OR. + ((M4.EQ.N3).AND.(M5.EQ.N2).AND.(M6.EQ.N1))) THEN FOUNDF=.TRUE. KFAULT=K GO TO 11 ENDIF 10 CONTINUE 11 IF (.NOT.FOUNDF) KFAULT=0 C IF THERE WAS A FAULT, REPLACE 3 NODE NUMBERS THAT WE SEARCH FOR: IF (FOUNDF) THEN IF (M2.EQ.N2) THEN N1=M4 N2=M5 N3=M6 ELSE N1=M1 N2=M2 N3=M3 ENDIF ENDIF C SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT: FOUNDE=.FALSE. DO 20 K=1,NUMEL IF (K.NE.I) THEN DO 15 L=1,3 M1=NODES(MOD(L, 3)+1,K) M2=NODES(MOD(L, 3)+4,K) M3=NODES(MOD(L+1,3)+1,K) IF ((M3.EQ.N1).AND.(M2.EQ.N2).AND.(M1.EQ.N3)) THEN FOUNDE=.TRUE. KELE=K GO TO 21 ENDIF 15 CONTINUE ENDIF 20 CONTINUE 21 IF (.NOT.FOUNDE) KELE=0 RETURN END C C C SUBROUTINE GETNET (INPUT,IUNIT7,IUNIT8, + MXEL,MXFEL,MXNODE,N1000, + OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD, + TITLE1,XNODE,YNODE,ZMNODE, + WORK,CHECKE,CHECKF,CHECKN) C C READ FINITE ELEMENT GRID FROM UNIT "IUNIT7". C ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8". C C READ NUMBER OF NODES, AND HOW MANY ARE "REAL" VERSUS "FAKE". C INPUT NODAL LOCATIONS (X,Y), ELEVATIONS, HEAT-FLOW, AND ISOSTATIC C GRAVITY ANOMALIES (ONE RECORD PER NODE). C CHARACTER*80 TITLE1 LOGICAL ALLOK,BRIEF C C NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN C DIMENSION CHECKE(MXEL),CHECKF(MXFEL),CHECKN(MXNODE), + DQDTDA(MXNODE),ELEV(MXNODE), + FAZ(2,MXFEL),FDIP(3,MXFEL), + NODEF(6,MXFEL),NODES(6,MXEL), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(3),IFN(6) C WRITE (IUNIT8,1) IUNIT7 1 FORMAT (' ATTEMPTING TO READ FINITE ELEMENT GRID FROM UNIT',I3) READ (IUNIT7,2) TITLE1 2 FORMAT (A80) WRITE (IUNIT8,3) TITLE1 3 FORMAT(/' TITLE OF FINITE ELEMENT GRID ='/' ',A80) C C READ NODES, AND SUBDIVIDE INTO 2 CLASSES C (OPTION "BRIEF" SUPPRESSES MOST OUTPUT) C READ (IUNIT7,*) NUMNOD,NREALN,NFAKEN,BRIEF C IF (NUMNOD.NE.(NREALN+NFAKEN)) THEN WRITE (IUNIT8,5) 5 FORMAT (/' INCONSISTENT DATA:'/ + ' NUMBER OF NODES SHOULD EQUAL TOTAL OF REAL', + ' NODES AND FAKE NODES.') STOP ENDIF IF (NUMNOD.GT.MXNODE) THEN WRITE (IUNIT8,10) NUMNOD 10 FORMAT(/' INCREASE PARAMETER MAXNOD TO BE AT LEAST' + /' THE NUMBER OF NODES (',I6,') AND RECOMPILE.') STOP ENDIF IF (NREALN.GT.N1000) THEN WRITE (IUNIT8,20) NREALN 20 FORMAT (/' INCREASE THE DATA VALUE N1000 TO BE GREATER' + /' OR EQUAL TO NREALN (',I6,') AND RECOMPILE.') STOP ENDIF C NBASE=N1000+1 NTOP=N1000+NFAKEN IF (BRIEF) THEN WRITE (IUNIT8,35) 35 FORMAT(/' (SINCE OPTION BRIEF=.TRUE., GRID WILL NOT BE ', + 'ECHOED HERE. BE CAREFUL]]])') ELSE WRITE (IUNIT8,40) NUMNOD,NREALN,NREALN 40 FORMAT (/' THERE ARE',I5,' NODES IN THE GRID:'/ ' ',I5, + ' OF THESE ARE NUMBERED 1-',I4,' AND THESE REAL NODES', + ' HAVE TWO VELOCITY VARIABLES UNLESS CONSTRAINED.') IF (NFAKEN.GT.0) WRITE (IUNIT8,42) NFAKEN,NBASE,NTOP 42 FORMAT(' ',I5, + ' OF THESE ARE NUMBERED ',I6,'-',I6,' AND THESE ARE', + ' ARTIFICIAL;'/' THEIR VELOCITIES MUST BE', + ' COMPLETELY SPECIFIED.') WRITE (IUNIT8,49) 49 FORMAT(/ + ' (NOTE: X AND Y COORDINATES MAY BE ZERO FOR NODES WHICH' + ,' WILL BE AT MIDPOINTS OF ELEMENT SIDES AND/OR FAULTS.'/ + ' THE PROGRAM WILL INTERPOLATE VALUES FOR THESE.)') WRITE (IUNIT8,50) 50 FORMAT (/' ', + ' CRUSTAL'/ + ' NODE X Y ELEVATION ', + 'HEAT-FLOW THICKNESS'/) ENDIF DO 90 K=1,NUMNOD CHECKN(K)=.FALSE. 90 CONTINUE DO 100 K=1,NUMNOD READ (IUNIT7,*) INDEX,XI,YI,ELEVI,QI,ZMI C (NODES NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) IF (INDEX.LE.NREALN) THEN I=INDEX ELSE I=INDEX-N1000+NREALN ENDIF CHECKN(I)=.TRUE. ELEV(I)=ELEVI DQDTDA(I)=QI IF (QI.LT.0.) THEN WRITE (IUNIT8,91) 91 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') STOP ENDIF XNODE(I)=XI YNODE(I)=YI IF (ZMI.LT.0.) THEN WRITE (IUNIT8,92) 92 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') STOP ENDIF ZMNODE(I)=ZMI IF (.NOT.BRIEF) THEN WRITE (IUNIT8,95) INDEX,XI,YI,ELEVI,QI,ZMI 95 FORMAT (' ',I10,1P,2E11.3,3E10.2) ENDIF 100 CONTINUE ALLOK=.TRUE. DO 101 I=1,NUMNOD ALLOK=ALLOK.AND.CHECKN(I) 101 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,102) 102 FORMAT(' THE FOLLOWING NODES WERE NEVER READ:') DO 104 I=1,NUMNOD IF (.NOT.CHECKN(I)) WRITE(IUNIT8,103)I 103 FORMAT (' ',36X,I6) 104 CONTINUE STOP ENDIF C C READ TRIANGULAR ELEMENTS C READ (IUNIT7,*) NUMEL IF (NUMEL.GT.MXEL) THEN WRITE (IUNIT8,108) NUMEL 108 FORMAT(/' INCREASE PARAMETER MAXEL TO BE AT LEAST EQUAL' + /' TO THE NUMBER OF ELEMENTS (',I6,') AND RECOMPILE.') STOP ENDIF DO 109 K=1,NUMEL CHECKE(K)=.FALSE. 109 CONTINUE IF (.NOT.BRIEF) WRITE (IUNIT8,110) NUMEL 110 FORMAT(/' THERE ARE ',I6,' TRIANGULAR CONTINUUM ELEMENTS.'/ + ' (NODE NUMBERS FOR EACH ARE GIVEN CORNERS-FIRST, COUNTER', + 'CLOCKWISE; THEN'/' MIDPOINTS, COUNTERCLOCKWISE, BEGINNING' + ,' WITH THE MIDPOINT BETWEEN CORNER #1 AND CORNER #2)'/ / + ' ELEMENT C1 C2 C3 M1 M2', + ' M3') DO 200 K=1,NUMEL C (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) READ (IUNIT7,*) I,(IFN(J),J=1,6) CHECKE(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,120) I,(IFN(J),J=1,6) 120 FORMAT (' ',I6,':',6I10) DO 130 J=1,6 N=IFN(J) IF ((N.LE.0).OR.(N.GT.NTOP).OR. + ((N.GT.NREALN).AND.(N.LE.N1000))) THEN WRITE (IUNIT8,125) N 125 FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.') STOP ENDIF IF (N.GT.NREALN) N=N-N1000+NREALN NODES(J,I)=N 130 CONTINUE 200 CONTINUE ALLOK=.TRUE. DO 201 I=1,NUMEL ALLOK=ALLOK.AND.CHECKE(I) 201 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,202) 202 FORMAT (' THE FOLLOWING ELEMENTS WERE NEVER READ:') DO 204 I=1,NUMEL IF (.NOT.CHECKE(I)) WRITE(IUNIT8,203)I 203 FORMAT (' ',39X,I6) 204 CONTINUE STOP ENDIF C C READ FAULT ELEMENTS C READ(IUNIT7,*) NFL IF (NFL.GT.MXFEL) THEN WRITE (IUNIT8,220)NFL 220 FORMAT (/' INCREASE PARAMETER MAXFEL TO BE AT LEAST EQUAL' + /' TO THE NUMBER OF FAULTS (',I6,') AND RECOMPILE.') STOP ENDIF DO 222 I=1,NFL CHECKF(I)=.FALSE. 222 CONTINUE IF (.NOT.BRIEF) WRITE(IUNIT8,230) NFL 230 FORMAT(/ /' THERE ARE ',I6,' CURVILINEAR FAULT ELEMENTS.') IF ((.NOT.BRIEF).AND.(NFL.GT.0)) WRITE(IUNIT8,231) 231 FORMAT (/' (THE 6 NODE NUMBERS DEFINING EACH ELEMENT MUST BE', + ' IN A COUNTERCLOCKWISE ORDER:'/ + ' N1, N2, AND N3 ARE IN LEFT-TO-RIGHT SEQUENCE ON THE', + ' NEAR SIDE,'/ + ' THEN N4 IS OPPOSITE N3, N5 IS OPPOSITE N2, AND ', + 'N6 IS OPPOSITE N1.)'/' (FAULT DIPS ARE GIVEN AT N1, N2, ', + 'AND N3, IN DEGREES FROM HORIZONTAL;'/ + ' POSITIVE DIPS ARE TOWARD N1, N2, AND N3, RESPECTIVELY, '/ + ' WHILE NEGATIVE DIPS ARE TOWARD N6, N5, AND N4.)'/ + ' (THE ARGUMENT OF THE FAULT TRACE IS GIVEN AT N1 AND N3,'/ + ' IN DEGREES COUNTERCLOCKWISE FROM THE X AXIS.)'/ / + ' ELEMENT N1 N2 N3 N4 N5 N6 DIP1 DIP2 DIP3 ', + 'ARG1 ARG3 '/) 240 FORMAT (' ',I6,':',6I5,3F6.1,2F5.0) DO 300 K=1,NFL READ (IUNIT7,*) I,(IFN(J),J=1,6),(DIPS(L),L=1,3), + AZ1,AZ3 CHECKF(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,240) I,(IFN(J),J=1,6), + (DIPS(L),L=1,3),AZ1,AZ3 DO 250 J=1,6 N=IFN(J) IF ((N.LE.0).OR.(N.GT.NTOP).OR. + ((N.GT.NREALN).AND.(N.LE.N1000))) THEN WRITE (IUNIT8,125) N STOP ENDIF IF (N.GT.NREALN) N=N-N1000+NREALN NODEF(J,I)=N 250 CONTINUE DO 260 L=1,3 IF (ABS(DIPS(L)).GT.90.) THEN WRITE(IUNIT8,252) DIPS(L) 252 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN', + ' RANGE OF -90. TO +90. DEGREES.'/ + ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE', + ' HORIZONAL;'/ + ' A + PREFIX (OR NONE) INDICATES A DIP', + ' TOWARD THE N1-N2-N3 SIDE;'/ + ' A - PREFIX INDICATES A DIP TOWARD', + ' THE N6-N5-N4 SIDE.)') STOP ENDIF IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L) FDIP(L,I)=DIPS(L)*0.017453293 260 CONTINUE IF ((ABS(AZ1).GT.361.).OR.(ABS(AZ3).GT.361.)) THEN WRITE(IUNIT8,272) AZ1,AZ3 272 FORMAT(' ILLEGAL ARGUMENT OF ',F10.4,' OR ',F10.4, + '; SHOULD BE IN RANGE -360. TO +360. DEGREES.') STOP ENDIF FAZ(1,I)=AZ1*0.017453293 FAZ(2,I)=AZ3*0.017453293 300 CONTINUE ALLOK=.TRUE. DO 301 I=1,NFL ALLOK=ALLOK.AND.CHECKF(I) 301 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,302) 302 FORMAT(' THE FOLLOWING FAULTS WERE NEVER READ:') DO 304 I=1,NFL IF (.NOT.CHECKF(I)) WRITE(IUNIT8,303)I 303 FORMAT(' ',36X,I6) 304 CONTINUE STOP ENDIF IF (.NOT. BRIEF) WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE SQUARE (INPUT,BRIEF,FDIP,IUNIT8, + MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,WEDGE, + MODIFY,FAZ,XNODE,YNODE, + OUTPUT,AREA,DETJ,DXS,DYS, + FLEN,FTAN, + WORK,CHECKN,LIST,NODTYP) C C CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID C LOGICAL ALLOK,BRIEF,FOUND,MATCH,SWITCH,VERT1,VERT2,VERT3 C C NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL CHECKN C C NOTE: THE FOLLOWING COULD BE MADE "INTEGER*2" IN VS-FORTRAN: INTEGER NODTYP C CHARACTER*21 OBLIQU,TAG1,TAG2,TAG3,VERTIC DOUBLE PRECISION FPHI,FPOINT,FGAUSS COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS DIMENSION FPHI(6,7),FPOINT(7),FGAUSS(7) DIMENSION AREA(MXEL),CHECKN(MXNODE), + DETJ(7,MXEL),DXS(6,7,MXEL),DYS(6,7,MXEL), + FDIP(3,MXFEL), + FAZ(2,MXFEL),FLEN(MXFEL),FTAN(7,MXFEL), + LIST(MXSTAR), + NODEF(6,MXFEL),NODES(6,MXEL),NODTYP(MXNODE), + XNODE(MXNODE),YNODE(MXNODE) DATA OBLIQU /'(DIP SLIP IS ALLOWED)'/ DATA VERTIC /'(STRIKE-SLIP ONLY) '/ C C (1) CHECK THAT ALL REAL NODES ARE CONNECTED TO AT LEAST ONE C CONTINUUM (TRIANGULAR) OR FAULT (LINEAR) ELEMENT: C DO 110 I=1,NREALN CHECKN(I)=.FALSE. 110 CONTINUE DO 120 I=1,NUMEL DO 115 J=1,6 CHECKN(NODES(J,I))=.TRUE. 115 CONTINUE 120 CONTINUE DO 130 I=1,NFL DO 125 J=1,6 CHECKN(NODEF(J,I))=.TRUE. 125 CONTINUE 130 CONTINUE ALLOK=.TRUE. DO 140 I=1,NREALN ALLOK=ALLOK.AND.CHECKN(I) 140 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,150) 150 FORMAT(' BAD GRID TOPOLOGY: FOLLOWING REAL NODES DO NOT'/ 1 ' BELONG TO ANY TRIANGULAR CONTINUUM ELEMENT:') DO 160 I=1,NREALN IF (.NOT.CHECKN(I)) WRITE (IUNIT8,155) I 155 FORMAT (' ',43X,I6) 160 CONTINUE STOP ENDIF C C (2) CHECK THAT EVERY NODE IS EITHER A CORNER OR A MIDPOINT NODE, C BUT NOT BOTH. C DO 210 I=1,NUMNOD NODTYP(I)=0 210 CONTINUE NTOFIX=0 ALLOK=.TRUE. DO 250 I=1,NUMEL DO 240 J=1,6 IF (J.LE.3) THEN ITYPE=1 ELSE ITYPE=2 ENDIF N=NODES(J,I) IF (NODTYP(N).EQ.0) THEN NODTYP(N)=ITYPE IF (ITYPE.EQ.2) THEN IF ((XNODE(N).EQ.0.).AND.(YNODE(N).EQ.0.)) + NTOFIX=NTOFIX+1 ENDIF ELSE IF (NODTYP(N).NE.ITYPE) THEN ALLOK=.FALSE. WRITE (IUNIT8,220) N 220 FORMAT(' BAD GRID TOPOLOGY: NODE ',I6, + ' CANNOT BE AN ELEMENT CORNER AND AN', + ' ELEMENT SIDE-MIDPOINT AT THE SAME', + ' TIME.') ENDIF ENDIF 240 CONTINUE 250 CONTINUE DO 290 I=1,NFL DO 280 J=1,6 IF ((J.EQ.2).OR.(J.EQ.5)) THEN ITYPE=2 ELSE ITYPE=1 ENDIF N=NODEF(J,I) IF (NODTYP(N).EQ.0) THEN NODTYP(N)=ITYPE ELSE IF (NODTYP(N).NE.ITYPE) THEN ALLOK=.FALSE. WRITE (IUNIT8,220) N ENDIF ENDIF 280 CONTINUE 290 CONTINUE IF (.NOT.ALLOK) STOP C C (3) CHECK THAT EACH FAULT SIDE WITH REAL NODES ALONG IT SHARES C THOSE SAME 3 NODES WITH A TRIANGULAR CONTINUUM ELEMENT. C ALLOK=.TRUE. DO 390 I=1,NFL DO 380 J=2,5,3 N=NODEF(J,I) IF (N.LE.NREALN) THEN DO 320 K=1,NUMEL DO 310 L=4,6 IF (NODES(L,K).EQ.N) THEN LP=L-2 IF (LP.EQ.4) LP=1 LM=L-3 MATCH=((NODEF(J-1,I).EQ.NODES(LP,K)) + .OR.(NODEF(J-1,I).GT.NREALN)) + .AND.((NODEF(J+1,I).EQ.NODES(LM,K)) + .OR.(NODEF(J+1,I).GT.NREALN)) IF (.NOT.MATCH) THEN ALLOK=.FALSE. WRITE(IUNIT8,305) I,K 305 FORMAT(' BAD GRID TOPOLOGY:', + ' FAULT ',I6,' IS NOT PROPERL' + ,'Y CONNECTED TO ELEMENT ',I6) ELSE GO TO 380 ENDIF ENDIF 310 CONTINUE 320 CONTINUE ENDIF 380 CONTINUE 390 CONTINUE IF (.NOT.ALLOK) STOP C C (4) AVERAGE TOGETHER THE COORDINATES OF ALL NODES AT ONE "POINT" C DO 410 I=1,NUMNOD CHECKN(I)=.FALSE. 410 CONTINUE DO 490 I=1,NFL DO 480 J1=1,3,2 NJ1=NODEF(J1,I) IF (.NOT.CHECKN(NJ1)) THEN LIST(1)=NJ1 CHECKN(NJ1)=.TRUE. J2=7-J1 NJ2=NODEF(J2,I) LIST(2)=NJ2 CHECKN(NJ2)=.TRUE. NINSUM=2 DO 470 K=I,NFL DO 460 L1=1,3,2 NL1=NODEF(L1,K) IF (.NOT.CHECKN(NL1)) THEN MATCH=.FALSE. DO 420 M=1,NINSUM MATCH=MATCH.OR.(NL1.EQ.LIST(M)) 420 CONTINUE IF (MATCH) THEN NINSUM=NINSUM+1 IF (NINSUM.GT.MXSTAR) THEN WRITE(IUNIT8,421) 421 FORMAT(/' INCREASE VALUE' + ,' OF PARAMETER MAXATP.') STOP ENDIF LIST(NINSUM)=NL1 CHECKN(NL1)=.TRUE. ENDIF L2=7-L1 NL2=NODEF(L2,K) MATCH=.FALSE. DO 430 M=1,NINSUM MATCH=MATCH.OR.(NL2.EQ.LIST(M)) 430 CONTINUE IF (MATCH) THEN NINSUM=NINSUM+1 IF (NINSUM.GT.MXSTAR) THEN WRITE(IUNIT8,421) STOP ENDIF LIST(NINSUM)=NL2 CHECKN(NL2)=.TRUE. ENDIF ENDIF 460 CONTINUE 470 CONTINUE XSUM=0. YSUM=0. DO 473 K=1,NINSUM XSUM=XSUM+XNODE(LIST(K)) YSUM=YSUM+YNODE(LIST(K)) 473 CONTINUE XMEAN=XSUM/NINSUM YMEAN=YSUM/NINSUM RMAX=0. DO 474 K=1,NINSUM R=SQRT((XNODE(LIST(K))-XMEAN)**2+ + (YNODE(LIST(K))-YMEAN)**2) RMAX=MAX(RMAX,R) 474 CONTINUE DO 475 K=1,NINSUM XNODE(LIST(K))=XMEAN YNODE(LIST(K))=YMEAN 475 CONTINUE IF (.NOT.BRIEF) THEN IF (RMAX.GT.0.) THEN WRITE(IUNIT8,472) NINSUM, + (LIST(N),N=1,NINSUM) 472 FORMAT(/ + ' AVERAGING TOGETHER THE POSITIONS OF', + ' THESE ',I6,' NODES:',(/' ',12I6)) WRITE (IUNIT8,476) RMAX 476 FORMAT (' MAXIMUM CORRECTION TO ', + 'ANY POSITION IS',1P,E10.2/ + ' YOU ARE RESPONSIBLE FOR ', + ' DECIDING WHETHER THIS IS A', + ' SERIOUS ERR0R]') ENDIF ENDIF ENDIF 480 CONTINUE 490 CONTINUE C C (5) SURVEY STRIKE-SLIP (VERTICAL) FAULTS TO CHECK FOR CONFLICTS IN C ARGUMENT THAT WOULD LOCK THE FAULT. C C LOOP ON ALL FAULT ELEMENTS (I): DO 2000 I=1,NFL C LOOP ON 2 TERMINAL NODE PAIRS, 1-6, 4-3 (J = 1 OR 4): DO 1900 J=1,4,3 C DIP MUST BE WITHIN "WEDGE" OF VERTICAL FOR CONSTRAINT: IF (ABS(FDIP(J,I)-1.570796).LE.WEDGE) THEN NAZI=1+J/4 N1=J N6=7-J NODE1=NODEF(N1,I) NODE6=NODEF(N6,I) C NO CONSTRAINT APPLIED WHERE A FAULT ENDS: IF (NODE1.NE.NODE6) THEN C ENDPOINT PAIRS MUST BE CHECKED FOR DUPLICATION: C LOOK FOR OTHER STRIKE-SLIP FAULTS SHARING THIS C PAIR OF NODES, AT EITHER END: FOUND=.FALSE. DO 1600 L=1,NFL IF (L.NE.I) THEN IF (ABS(FDIP(1,L)-1.5708).LE.WEDGE) THEN IF (((NODE1.EQ.NODEF(1,L)).AND. + (NODE6.EQ.NODEF(6,L))).OR. + ((NODE1.EQ.NODEF(6,L)).AND. + (NODE6.EQ.NODEF(1,L)))) THEN FOUND=.TRUE. NUMBER=L NAZL=1 GO TO 1601 ENDIF ENDIF IF (ABS(FDIP(2,L)-1.5708).LE.WEDGE) THEN IF (((NODE1.EQ.NODEF(3,L)).AND. + (NODE6.EQ.NODEF(4,L))).OR. + ((NODE1.EQ.NODEF(4,L)).AND. + (NODE6.EQ.NODEF(3,L)))) THEN FOUND=.TRUE. NUMBER=L NAZL=2 GO TO 1601 ENDIF ENDIF ENDIF 1600 CONTINUE C DON'T WORRY IF THIS PAIR ALREADY CHECKED] 1601 IF (FOUND.AND.(NUMBER.GT.I)) THEN C AVERAGE ARGUMENTS TOGETHER (AVOID CYCLE SHIFTS): AZI=MOD(FAZ(NAZI,I)+1.570796,3.14159265) + -1.570796 AZL=MOD(FAZ(NAZL,NUMBER)+1.570796,3.14159265) + -1.570796 AZIMUT=0.5*(AZI+AZL) FAZ(NAZI,I)=AZIMUT FAZ(NAZL,NUMBER)=AZIMUT IF (ABS(AZI-AZL).GT.0.02) THEN DAZI=AZI*57.2957795 DAZL=AZL*57.2957795 DAZ=AZIMUT*57.2957795 IF (NODE1.LE.NREALN) THEN NP1=NODE1 ELSE NP1=N1000+NODE1-NREALN ENDIF IF (NODE6.LE.NREALN) THEN NP6=NODE6 ELSE NP6=N1000+NODE6-NREALN ENDIF WRITE (IUNIT8,1610) I,NUMBER,NP1,NP6, + DAZI,DAZL,DAZ 1610 FORMAT(/' WARNING: STRIKE-SLIP FAULT ELEMENTS' + ,I7,' AND',I7/' SHARE NODES',I7,' AND', + I7/' BUT THEIR ARGUMENTS OF ',F6.1, + ' AND ',F6.1,' DEGREES DIFFER SUBSTAN', + 'TIALLY.'/' THE ARGUMENTS WILL BE AVERAGED,' + ,' AND A VALUE OF ',F6.1,' WILL BE USED.' + /' THIS IS NECESSARY TO PREVENT FAULT', + ' LOCKING;'/' IF YOU -WANT- THE FAULT LOCKED' + ,', THEN USE A SINGLE NODE AT THIS POINT.') ENDIF ENDIF C ^END BLOCK WHICH LOOKS FOR CONSTRAINTS ON REAL NODES ENDIF C ^END BLOCK WHICH CHECKS FOR DISTINCT NODE NUMBERS ENDIF C ^END BLOCK WHICH CHECKS FOR DIP OF OVER 75 DEGREES 1900 CONTINUE C ^END LOOP ON 2 NODE PAIRS IN FAULT ELEMENT 2000 CONTINUE C ^END LOOP ON FAULT ELEMENTS C C (6) COMPUTE COORDINATES OF MIDPOINT NODES THAT WERE NOT INPUT. C C FIRST, FAULTS: IF ((.NOT.BRIEF).AND.(NFL.GT.0)) WRITE (IUNIT8,540) 540 FORMAT(/ /' FOLLOWING FAULT MID-POINT POSITIONS WERE COMPUTED:'/ + /' FAULT NODE2 NODE5 X Y'/) DO 550 I=1,NFL I1=NODEF(1,I) I2=NODEF(2,I) I3=NODEF(3,I) I5=NODEF(5,I) DX= XNODE(I3)- XNODE(I1) DY= YNODE(I3)- YNODE(I1) AZ=ATAN2(DY,DX) PHI1=FAZ(1,I)-AZ PHI1=MOD(PHI1+1.570796,3.14159265)-1.570796 PHI2=AZ-FAZ(2,I) PHI2=MOD(PHI2+1.570796,3.14159265)-1.570796 IF ((ABS(PHI1).GT.0.001).OR.(ABS(PHI2).GT.0.001)) THEN T1=TAN(PHI1) T2=TAN(PHI2) IF (ABS(T2-T1).GE.ABS(T1+T2)) THEN FACTOR=0.99*ABS(T1+T2)/ABS(T2-T1) IF (ABS(T1).GT.ABS(T2)) THEN T2=T1+FACTOR*(T2-T1) ELSE T1=T2+FACTOR*(T1-T2) ENDIF ENDIF PARRAL=(T2-T1)/(4.*(T1+T2)) PERPEN= T1*T2 /(2.*(T1+T2)) XNODE(I2)=XNODE(I1)+DX/2.+PARRAL*DX-PERPEN*DY YNODE(I2)=YNODE(I1)+DY/2.+PERPEN*DX+PARRAL*DY ELSE XNODE(I2)=(XNODE(I1)+XNODE(I3))/2. YNODE(I2)=(YNODE(I1)+YNODE(I3))/2. ENDIF XNODE(I5)= XNODE(I2) YNODE(I5)= YNODE(I2) NTOFIX=NTOFIX-1 IF (.NOT.BRIEF) WRITE (IUNIT8,549) I,I2,I5,XNODE(I2), 1 YNODE(I2) 549 FORMAT(' ',I6,2I10,1P,2E12.4) 550 CONTINUE C C NEXT, OTHER ELEMENT SIDES, IF NEEDED: IF ((.NOT.BRIEF).AND.(NTOFIX.GT.0)) WRITE (IUNIT8,551) 551 FORMAT(/ /' FOLLOWING MID-POINTS OF CONTINUUM ELEMENT SIDES', + ' THAT WERE 0.0 IN THE' + / ' INPUT DATASET ARE NOW COMPUTED, AS FOLLOWS:' + / / ' ELEMENT NODE X Y'/) DO 590 I=1,NUMEL DO 580 J=4,6 N=NODES(J,I) IF ((XNODE(N).EQ.0.).AND.(YNODE(N).EQ.0.)) THEN JP=J-2 IF (J.EQ.6) JP=1 JM=J-3 XNODE(N)=0.5* + (XNODE(NODES(JP,I))+XNODE(NODES(JM,I))) YNODE(N)=0.5* + (YNODE(NODES(JP,I))+YNODE(NODES(JM,I))) IF (.NOT.BRIEF) + WRITE (IUNIT8,579) I,N,XNODE(N),YNODE(N) 579 FORMAT(' ',I6,I10,1P,2E12.4) ENDIF 580 CONTINUE 590 CONTINUE C C (7) COMPUTE AREAS OF ELEMENTS AND COMPUTE DERIVITIVES OF NODAL C FUNCTIONS AT INTEGRATION POINTS; C THEN CHECK FOR NEGATIVE AREAS C CALL AREAS (INPUT,NODES,NUMEL,NUMNOD,XNODE,YNODE, + OUTPUT,AREA) CALL DERIV (INPUT,AREA,NODES,NUMEL,NUMNOD,XNODE,YNODE, + OUTPUT,DETJ,DXS,DYS) ALLOK=.TRUE. DO 620 M=1,7 DO 610 I=1,NUMEL TEST=AREA(I)*DETJ(M,I) IF (TEST.LE.0.) THEN WRITE(IUNIT8,605) M,I 605 FORMAT(/' EXCESSIVELY DISTORTED ELEMENT LEADS TO ' + ,'NEGATIVE AREA AT POINT ',I1,' IN ELEMENT ', + I5) ALLOK=.FALSE. ENDIF 610 CONTINUE 620 CONTINUE IF (.NOT.ALLOK) STOP C C (8) COMPUTE LENGTHS OF FAULT ELEMENTS. C DO 750 I=1,NFL FLEN(I)=0. X1=XNODE(NODEF(1,I)) X2=XNODE(NODEF(2,I)) X3=XNODE(NODEF(3,I)) Y1=YNODE(NODEF(1,I)) Y2=YNODE(NODEF(2,I)) Y3=YNODE(NODEF(3,I)) OLDX=X1 OLDY=Y1 DO 740 J=1,20 S=J/20. F1=1.-3.*S+2.*S**2 F2=4.*S*(1.-S) F3= -S+2.*S**2 X=X1*F1+X2*F2+X3*F3 Y=Y1*F1+Y2*F2+Y3*F3 FLEN(I)=FLEN(I)+SQRT((X-OLDX)**2+(Y-OLDY)**2) OLDX=X OLDY=Y 740 CONTINUE 750 CONTINUE C C (10) SURVEY FAULT ELEMENTS AND ISSUE WARNING IF ANY ELEMENT IS OF C MIXED TYPE (PART STRIKE-SLIP, AND PART SHALLOW-DIPPING: C DO 920 I=1,NFL DELD1=FDIP(1,I)-1.570796 DELD2=FDIP(2,I)-1.570796 DELD3=FDIP(3,I)-1.570796 VERT1=ABS(DELD1).LE.WEDGE VERT2=ABS(DELD2).LE.WEDGE VERT3=ABS(DELD3).LE.WEDGE NVPART=0 IF (VERT1) THEN NVPART=NVPART+1 TAG1=VERTIC ELSE TAG1=OBLIQU ENDIF IF (VERT2) THEN NVPART=NVPART+1 TAG2=VERTIC ELSE TAG2=OBLIQU ENDIF IF (VERT3) THEN NVPART=NVPART+1 TAG3=VERTIC ELSE TAG3=OBLIQU ENDIF SWITCH=((NVPART.GT.0).AND.(NVPART.LT.3)) IF (SWITCH) THEN DIP1=FDIP(1,I)*57.2957795 IF (DIP1.GT.90.) DIP1=DIP1-180. DIP2=FDIP(2,I)*57.2957795 IF (DIP2.GT.90.) DIP2=DIP2-180. DIP3=FDIP(3,I)*57.2957795 IF (DIP3.GT.90.) DIP3=DIP3-180. WRITE (IUNIT8,905) I,DIP1,TAG1,DIP2,TAG2,DIP3,TAG3 905 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',F7.2,' DEGREES ',A21/ + ' ',F7.2,' DEGREES ',A21/ + ' ',F7.2,' DEGREES ',A21/ + ' WHICH MAKES IT MIXED-MODE.'/ + ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ + ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ + ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', + ' FAULT ELEMENT TO THE NEXT.)') ELSE NVPART=0 DO 910 M=1,7 DELD=DELD1*FPHI(1,M)+DELD2*FPHI(2,M)+ + DELD3*FPHI(3,M) IF (ABS(DELD).LE.WEDGE) NVPART=NVPART+1 910 CONTINUE IF ((NVPART.GT.0).AND.(NVPART.LT.7)) THEN IF (NVPART.GE.4) THEN WRITE (IUNIT8,912) I,DIP1,DIP2,DIP3 912 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',F7.2,' DEGREES,'/ + ' ',F7.2,' DEGREES, AND'/ + ' ',F7.2,' DEGREES'/ + ' WHICH APPEAR TO MAKE IT STRIKE-SLIP.'/ + ' HOWEVER, THESE VALUES ARE SUCH THAT DIP-SLIP'/ + ' IS PERMITTED AT ONE OR MORE INTEGRATION POINTS.'/ + ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ + ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ + ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', + ' FAULT ELEMENT TO THE NEXT.)') ELSE WRITE (IUNIT8,914) I,DIP1,DIP2,DIP3 914 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',F7.2,' DEGREES,'/ + ' ',F7.2,' DEGREES, AND'/ + ' ',F7.2,' DEGREES'/ + ' WHICH APPEAR TO MAKE IT FREE-SLIPPING.'/ + ' HOWEVER, THESE VALUES ARE SUCH THAT DIP-SLIP'/ + ' IS PROHIBITED AT ONE OR MORE INTEGRATION POINTS.'/ + ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ + ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ + ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', + ' FAULT ELEMENT TO THE NEXT.)') ENDIF ENDIF ENDIF 920 CONTINUE C C (11) CALCULATE FAULT ARGUMENT (IN RADIANS, MEASURED COUNTERCLOCKWISE C FROM +X) FOR EACH INTEGRATION POINT IN EACH FAULT ELEMENT. C DO 1000 M=1,7 S=FPOINT(M) DF1DS= -3.+4.*S DF2DS=4.-8.*S DF3DS= -1.+4.*S DO 900 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) N3=NODEF(3,I) X1=XNODE(N1) X2=XNODE(N2) X3=XNODE(N3) Y1=YNODE(N1) Y2=YNODE(N2) Y3=YNODE(N3) DXDS=X1*DF1DS+X2*DF2DS+X3*DF3DS DYDS=Y1*DF1DS+Y2*DF2DS+Y3*DF3DS FTAN(M,I)=ATAN2(DYDS,DXDS) 900 CONTINUE 1000 CONTINUE C IF (.NOT. BRIEF) WRITE (IUNIT8,9999) 9999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE AREAS (INPUT,NODES,NUMEL,NUMNOD,XNODE,YNODE, 1 OUTPUT,AREA) C C COMPUTE AREAS OF ELEMENTS IN GRID AS IF THEY HAD STRAIGHT C SIDES. EFFECT OF SIDE CURVATURE WILL BE HANDLED LATER BY C MULTIPLYING BY DETERMINANT OF JACOBIAN MATRIX FOR THE SIDE- C BENDING MAPPING. NOTE THAT AREA MAY BE NEGATIVE, BUT ELEMENT C IS OK IF DETERMINANT IN DERIV IS ALSO NEGATIVE. C DIMENSION AREA(NUMEL),NODES(6,NUMEL),XNODE(NUMNOD),YNODE(NUMNOD) DO 100 I=1,NUMEL I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) AREA(I)= 0.5*(XNODE(I1)*YNODE(I2)-XNODE(I2)*YNODE(I1) + +XNODE(I2)*YNODE(I3)-XNODE(I3)*YNODE(I2) + +XNODE(I3)*YNODE(I1)-XNODE(I1)*YNODE(I3)) 100 CONTINUE RETURN END C C C SUBROUTINE DERIV (INPUT,AREA,NODES,NUMEL,NUMNOD,XNODE,YNODE, + OUTPUT,DETJ,DXS,DYS) C C CALCULATES DXS AND DYS, THE X-DERIVITIVE AND Y-DERIVITIVE C OF EACH OF THE 6 NODAL FUNCTIONS OF A DEFORMED-TRIANGLE C FINITE ELEMENT, AT EACH OF THE 7 INTEGRATION POINTS IN C THAT ELEMENT. ALSO PROVIDES DETJ, THE DETERMINANT OF THE C JACOBIAN MATRIX FOR THE TRANSFORMATION IN WHICH INTERNAL C POINTS OF A TRIANGLE WITH STRAIGHT SIDES ARE MAPPED INTO C NEW LOCATIONS AS SIDES BEND (BUT CORNERS STAY FIXED). C DOUBLE PRECISION POINTS DIMENSION AREA(NUMEL),DETJ(7,NUMEL), + DXS(6,7,NUMEL),DYS(6,7,NUMEL), + NODES(6,NUMEL), + XNODE(NUMNOD),YNODE(NUMNOD) DIMENSION B(4),C(4),DN(6,2),POINTS(5,7),X(6),Y(6) COMMON /S1S2S3/ POINTS DO 500 I=1,NUMEL DO 100 J=1,6 NODE=NODES(J,I) X(J)=XNODE(NODE) Y(J)=YNODE(NODE) 100 CONTINUE B(1)=Y(2)-Y(3) B(2)=Y(3)-Y(1) B(3)=Y(1)-Y(2) B(4)=B(1) C(1)=X(3)-X(2) C(2)=X(1)-X(3) C(3)=X(2)-X(1) C(4)=C(1) AI2=1./(2.*AREA(I)) DO 400 M=1,7 DO 200 J=1,3 DN(J,1)=AI2*B(J)*(4.*POINTS(J,M)-1.) DN(J+3,1)=AI2*4.*(B(J)*POINTS(J+1,M) + +B(J+1)*POINTS(J,M)) DN(J,2)=AI2*C(J)*(4.*POINTS(J,M)-1.) DN(J+3,2)=AI2*4.*(C(J)*POINTS(J+1,M) + +C(J+1)*POINTS(J,M)) 200 CONTINUE AJ11=0. AJ12=0. AJ21=0. AJ22=0. DO 300 J=1,6 AJ11=AJ11+DN(J,1)*X(J) AJ12=AJ12+DN(J,1)*Y(J) AJ21=AJ21+DN(J,2)*X(J) AJ22=AJ22+DN(J,2)*Y(J) 300 CONTINUE DETJAC=AJ11*AJ22-AJ12*AJ21 DETJ(M,I)=DETJAC AJ11S=AJ11 AJ11=AJ22/DETJAC AJ12=-AJ12/DETJAC AJ21=-AJ21/DETJAC AJ22=AJ11S/DETJAC DO 350 J=1,6 DXS(J,M,I)=AJ11*DN(J,1)+AJ12*DN(J,2) DYS(J,M,I)=AJ21*DN(J,1)+AJ22*DN(J,2) 350 CONTINUE 400 CONTINUE 500 CONTINUE RETURN END C C C BLOCK DATA BD1 C C DEFINE "PHI" (NODAL FUNCTIONS) AND "WEIGHT" (GAUSSIAN INTEGRATION C WEIGHTS) OF THE 6-NODE TRIANGULAR FINITE ELEMENT FOR THE C SEVEN INTEGRATION POINTS IN EACH ELEMENT, DEFINED BY INTERNAL C COORDINATES "POINTS(5,7)", WHERE POINTS(1-3,M)=S1-S3 OF C INTEGRATION POINT NUMBER M. (NOTE: POINTS(4,M)=POINTS(1,M) AND C POINTS(5,M)=POINTS(2,M), FOR PROGRAMMING CONVENIENCE, AS IN C SUBPROGRAM "DERIV".) C BECAUSE ALL OF THESE ARRAYS ARE FUNCTIONS OF INTERNAL C COORDINATES, THEY ARE NOT AFFECTED BY SCALING OR DEFORMATION OF C THE ELEMENTS. C DOUBLE PRECISION PHI,POINTS,WEIGHT COMMON /S1S2S3/ POINTS COMMON /PHITAB/ PHI COMMON /WGTVEC/ WEIGHT DIMENSION PHI(6,7),POINTS(5,7),WEIGHT(7) C C "PHI" CONTAINS THE VALUES OF THE 6 NODAL FUNCTIONS AT THE 7 C GAUSSIAN INTEGRATION POINTS (FOR AREA INTEGRALS) OF THE C TRIANGULAR ELEMENTS. DATA PHI / +-0.1111111111111111D0,-0.1111111111111111D0,-0.1111111111111111D0, + 0.4444444444444444D0, 0.4444444444444444D0, 0.4444444444444444D0, +-0.0525839022774079D0,-0.0280749439026853D0,-0.0280749439026853D0, + 0.1122997756107412D0, 0.8841342388612960D0, 0.1122997756107412D0, +-0.0280749439026853D0,-0.0525839022774079D0,-0.0280749439026853D0, + 0.1122997756107412D0, 0.1122997756107412D0, 0.8841342388612960D0, +-0.0280749439026853D0,-0.0280749439026853D0,-0.0525839022774079D0, + 0.8841342388612960D0, 0.1122997756107412D0, 0.1122997756107412D0, + 0.4743526114618935D0,-0.0807685938011933D0,-0.0807685938011933D0, + 0.3230743752047730D0, 0.0410358257309469D0, 0.3230743752047730D0, +-0.0807685938011933D0, 0.4743526114618935D0,-0.0807685938011933D0, + 0.3230743752047730D0, 0.3230743752047730D0, 0.0410358257309469D0, +-0.0807685938011933D0,-0.0807685938011933D0, 0.4743526114618935D0, + 0.0410358257309469D0, 0.3230743752047730D0, 0.3230743752047730D0/ C C "POINTS" CONTAINS THE INTERNAL COORDINATES (S1,S2,S3) OF THE 7 C GAUSSIAN INTEGRATION POINTS (FOR AREA INTEGRALS) OF THE C TRIANGULAR ELEMENTS. DATA POINTS / + 0.3333333333333333D0, 0.3333333333333333D0, 0.3333333333333333D0, + 0.3333333333333333D0, 0.3333333333333333D0, + 0.0597158733333333D0, 0.4701420633333333D0, 0.4701420633333333D0, + 0.0597158733333333D0, 0.4701420633333333D0, + 0.4701420633333333D0, 0.0597158733333333D0, 0.4701420633333333D0, + 0.4701420633333333D0, 0.0597158733333333D0, + 0.4701420633333333D0, 0.4701420633333333D0, 0.0597158733333333D0, + 0.4701420633333333D0, 0.4701420633333333D0, + 0.7974269866666667D0, 0.1012865066666667D0, 0.1012865066666667D0, + 0.7974269866666667D0, 0.1012865066666667D0, + 0.1012865066666667D0, 0.7974269866666667D0, 0.1012865066666667D0, + 0.1012865066666667D0, 0.7974269866666667D0, + 0.1012865066666667D0, 0.1012865066666667D0, 0.7974269866666667D0, + 0.1012865066666667D0, 0.1012865066666667D0/ C C "WEIGHT" IS THE GAUSSIAN WEIGHT (FOR AREA INTEGRALS) OF THE 7 C INTEGRATION POINTS IN EACH TRIANGULAR ELEMENT. DATA WEIGHT / 0.2250000000000000D0, + 0.1323941500000000D0, 0.1323941500000000D0, 0.1323941500000000D0, + 0.1259391833333333D0, 0.1259391833333333D0, 0.1259391833333333D0/ C END C C C BLOCK DATA BD2 C C DEFINE "FPHI" (NODAL FUNCTIONS) AND "FGAUSS" (GAUSSIAN INTEGRATION C WEIGHTS) OF THE 6-NODE LINEAR FAULT ELEMENT FOR THE SEVEN C INTEGRATION POINTS IN EACH ELEMENT, DEFINED BY INTERNAL C COORDINATE "FPOINT(M=1-7)", WHICH CONTAINS THE RELATIVE POSITION C (FRACTIONAL LENGTH) OF THE INTEGRATION POINTS. C BECAUSE ALL OF THESE ARRAYS ARE FUNCTIONS OF INTERNAL C COORDINATES, THEY ARE NOT AFFECTED BY SCALING OR DEFORMATION OF C THE ELEMENTS. C DOUBLE PRECISION FPHI,FPOINT,FGAUSS COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS DIMENSION FPHI(6,7),FPOINT(7),FGAUSS(7) C C "FPOINT" CONTAINS THE SEVEN INTEGRATION POINT LOCATIONS FOR THE FAULT C ELEMENTS. EACH VALUE GIVES A POSITION AS A FRACTION OF TOTAL LENGTH C MEASURED FROM NODE1 TO NODE3 (OF ARRAY "NODEF"). DATA FPOINT/ 1 0.0254461D0, 2 0.1292344D0, 3 0.2970774D0, 4 0.5000000D0, 5 0.7029226D0, 6 0.8707656D0, 7 0.9745539D0 / C C "FGAUSS" CONTAINS THE SEVEN CORRESPONDING WEIGHT FACTORS FOR USE IN C LINE INTEGRALS. DATA FGAUSS/ 1 0.0647425D0, 2 0.1398527D0, 3 0.1909150D0, 4 0.2089796D0, 5 0.1909150D0, 6 0.1398527D0, 7 0.0647425D0/ C C "FPHI" CONTAINS THE VALUES OF THE 6 NODAL FUNCTIONS (ONE PER NODE) C AT EACH OF THESE 7 INTEGRATION POINTS IN THE FAULT ELEMENT. DATA FPHI/ + .92495670801042D0, .09919438397916D0,-.02415109198958D0, + .02415109198958D0,-.09919438397916D0,-.92495670801042D0, + .64569986028672D0, .45013147942656D0,-.09583133971328D0, + .09583133971328D0,-.45013147942656D0,-.64569986028672D0, + .28527776318152D0, .83528967363696D0,-.12056743681848D0, + .12056743681848D0,-.83528967363696D0,-.28527776318152D0, + 0.0D0, 1.0D0, 0.0D0, + 0.0D0, -1.0D0, 0.0D0, + -.12056743681848D0, .83528967363696D0, .28527776318152D0, + -.28527776318152D0,-.83528967363696D0, .12056743681848D0, + -.09583133971328D0, .45013147942656D0, .64569986028672D0, + -.64569986028672D0,-.45013147942656D0, .09583133971328D0, + -.02415109198958D0, .09919438397916D0, .92495670801042D0, + -.92495670801042D0,-.09919438397916D0, .02415109198958D0/ C END