C PROGRAM SCORER.GPB.ALASKA C C COMPUTES NUMERICAL SCORES OF FINITE ELEMENT VELOCITY MODELS FROM C PROGRAM -PLATES- BY USING SLIP-RATE, STRESS, AND GEODETIC DATA. C C THIS VERSION CUSTOMIZED FOR THE ALASKAN PROBLEM IN THE PLACES C MARKED BY "C****** ALASKA (BEGIN SECTION) ******************" C "C****** ALASKA (END SECTION) ******************" C C *READS FINITE ELEMENT GRID GEOMETRY FROM FORTRAN DEVICE "IUNITF". C *READS DATA ON FAULT SLIP RATES FROM FORTRAN DEVICE "IUNITD". C *READS DATA ON STRESS DIRECTIONS FROM FORTRAN DEVICE "IUNITS". C *READS DATA ON GEODETIC BASELINE LENGTH RATES FROM FORTRAN DEVICE C "IUNITG". C *THEN, IN AN INDEFINATE LOOP: C -READS PARAMETER SET(S) OF FINITE ELEMENT MODELS FROM FORTRAN C DEVICE "IUNITP". C -READS CONVERGED SOLUTION(S) (3 TITLES AND NODE VELOCITIES ONLY) C FROM FORTRAN DEVICE "IUNITV"; C -COMPUTES FAULT SLIP RATES, AND COMPARES THEM TO GEOLOGIC DATA. C -COMPUTES PREDICTED STRESS DIRECTIONS, AND COMPARES TO DATA. C -COMPUTES PREDICTED BASELINE LENGTH RATES (WITH FRICTIONAL PARTS C OF EACH FAULT TEMPORARILY LOCKED), AND COMPARES TO DATA. C -OUTPUTS SCORES TO FORTRAN DEVICE "IUNITT". C C C BY C PETER BIRD, C DEPARTMENT OF EARTH AND SPACE SCIENCES, C UNIVERSITY OF CALIFORNIA, LOS ANGELES, CALIFORNIA 90024 C (WITH ASSISTANCE FROM XIANGHONG KONG). C *** SET THE VERSION DATE IN FORMAT STATEMENT 1 BELOW *** C C C THIS PROGRAM WAS DEVELOPED WITH SUPPORT FROM THE UNIVERSITY OF C CALIFORNIA, THE UNITED STATES GEOLOGIC SURVEY, AND THE NATIONAL C SCIENCE FOUNGSDATION. IT IS IN THE PUBLIC DOMAIN, AND MAY BE COPIED C AND USED WITHOUT RESTRICTION. HOWEVER, SCIENTIFIC ETHICS AND C COURTESY REQUIRE THE SOURCE OF THE PROGRAM TO BE STATED IN ANY C RESULTING PUBLICATIONS. (THE AUTHOR WOULD ALSO LIKE TO BE INFORMED C OF THESE PROJECTS.) 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=1935) C C MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES). PARAMETER (MAXEL=722) C C MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS); PARAMETER (MAXFEL=196) C C MAXATP = MAXIMUM NUMBER OF NODES WHICH MAY OVERLAP AT A FAULT- C INTERSECTION POINT. PARAMETER (MAXATP=20) C C MAXDAT = MAXIMUM NUMBER OF DATA OF ANY PARTICULAR TYPE (SET LARGE): PARAMETER (MAXDAT=1000) C C MAXBEN = MAXIMUM NUMBER OF GEODETIC BENCHMARKS: PARAMETER (MAXBEN=200) C C MAXBAS = MAXIMUM NUMBER OF GEODETIC BASELINES: PARAMETER (MAXBAS=1000) C---------------------------------------------------------------------- C C TYPE STATEMENTS C (NOTE: IMPLICIT CONVENTION OF I-N = INTEGER IS FOLLOWED) C CHARACTER*2 QUALIT CHARACTER*4 REGIME,REGSTR CHARACTER*6 SITE,TAGSTR,TYPE CHARACTER*7 DEPTH CHARACTER*28 REFERE,MORERE CHARACTER*31 BENCH,BENCHN,SITE1,SITE2,SITE1N,SITE2N CHARACTER*59 TAGGSR,TEST,BLANKS CHARACTER*80 TITLE1,TITLE2,TITLE3,TITL1G,TITL2G,TITL3G,TITL4G C DOUBLE PRECISION V,VM 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 ALDONE,BRIEF,EVERYP,OK1,OK2 C C NOTE: THE FOLLOWING ARRAYS COULD BE COMPRESSED WITH "LOGICAL*1" C IN VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN,FSLIPS,PULLED C C---------------------------------------------------------------------- C C DIMENSION STATEMENTS C C DIMENSIONS USING PARAMETER MAXNOD: DIMENSION CHECKN(MAXNOD), DQDTDA(MAXNOD), + ELEV (MAXNOD), NFTRE (MAXNOD), + NODTYP(MAXNOD), TLNODE (MAXNOD), + V (2,MAXNOD), VM (2,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), + ERATE(3,7,MAXEL), + NODES (6,MAXEL), + PULLED (7,MAXEL), + TLINT (7,MAXEL), ZMOHO(7,MAXEL) C C DIMENSIONS USING PARAMETER MAXFEL: DIMENSION CHECKF (MAXFEL), + FAZ (2,MAXFEL), FC(2,2,7,MAXFEL), + FDIP (3,MAXFEL), + FIMUDZ (7,MAXFEL), FLEN (MAXFEL), + FPEAKS (MAXFEL), FSLIPS (MAXFEL), + FTSTAR(2,7,MAXFEL), FTAN (7,MAXFEL), NODEF (6,MAXFEL), + OFFSET (MAXFEL), ZTRANF (7,MAXFEL) C C DIMENSIONS USING PARAMETER MAXATP: DIMENSION LIST (MAXATP) C C DIMENSIONS USING PARAMETER MAXDAT: C GEOLOGIC SLIP-RATE-RELATED: DIMENSION IFDATM(MAXDAT),DATUMS(MAXDAT),IENDND(MAXDAT), + DELVZ(2,MAXDAT),RLAT(2,MAXDAT),TAGGSR(MAXDAT) C STRESS-DIRECTION-RELATED: DIMENSION IESTR (MAXDAT),REGSTR(MAXDAT),WGTSTR (MAXDAT), + ARGSTR(MAXDAT),TAGSTR(MAXDAT),IASIGN(3,MAXDAT) C C DIMENSIONS USING PARAMETER MAXBEN (MAXIMUM NUMBER OF BENCHMARKS): DIMENSION BENCH (MAXBEN), BENCHX(MAXBEN), BENCHY(MAXBEN), + BENS1 (MAXBEN), BENS2 (MAXBEN), + IBENCH (MAXBEN), + UX (MAXBEN), UY (MAXBEN), + VGCRE(2,MAXBEN), VG (2,MAXBEN) C C DIMENSIONS USING PARAMETER MAXBAS (MAXIMUM NUMBER OF BASELINES): DIMENSION BERR (MAXBAS), BRATE (MAXBAS), + INDFST(MAXBAS), INDLST(MAXBAS), + SITE1 (MAXBAS), SITE2 (MAXBAS) C C DIMENSIONS OF FIXED SIZE: DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),RADIO (2),RHOBAR(2),TEMLIM(2) DIMENSION IMSIGN(3) 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 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 DATA STATEMENTS C C THE FOLLOWING ARE THE FORTRAN INPUT AND OUTPUT DEVICE NUMBERS: C C "IUNITF"= DEVICE NUMBER ASSOCIATED WITH THE FINITE ELEMENT GRID C INPUT FILE. DATA IUNITF /1/ C C "IUNITD"= DEVICE NUMBER ASSOCIATED WITH THE GEOLOGIC SLIP RATE DATA. DATA IUNITD /2/ C C "IUNITS"= DEVICE NUMBER ASSOCIATED WITH THE STRESS DIRECTION DATA. DATA IUNITS /3/ C C "IUNITG"= DEVICE NUMBER ASSOCIATED WITH GEODETIC DATA INPUT DATA IUNITG /4/ C C "IUNITT"= DEVICE NUMBER ASSOCIATED WITH TEXT OUTPUT, INCLUDING C STATUS AND WARNING MESSAGES. C (NOTE: ON SOME SYSTEMS, SYSTEM ERR0R MESSAGES ARE ALWAYS C OUTPUT ON DEVICE 6. IF SO, THEN IUNITT SHOULD BE 6.) DATA IUNITT /6/ C C "IUNITP"= DEVICE NUMBER ASSOCIATED WITH THE PARAMETER INPUT FILE(S) C CORRESPONDING TO THE SOLUTION(S) TO BE SCORED. DATA IUNITP /11/ C C "IUNITV"= DEVICE NUMBER ASSOCIATED WITH THE VELOCITY C SOLUTION(S), TO BE SCORED. DATA IUNITV /12/ 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 "UNITL" IS THE NUMBER OF MILLIMETERS PER UNIT LENGTH IN THE C SYSTEM OF UNITS IN USE ON FORTRAN DEVICES IUNITF (FINITE C ELEMENT GRID), AND DEVICE IUNITP (PARAMETERS), C AND IUNITV (VELOCITY SOLUTIONS). FOR EXAMPLE, IF THESE ARE IN C SI UNITS (METERS), THEN UNITL = 1000. DATA UNITL /1000./ C C "UNITT" IS THE NUMBER OF YEARS PER UNIT OF TIME IN THE C SYSTEM OF UNITS IN USE ON FORTRAN DEVICES IUNITF (FINITE C ELEMENT GRID), AND DEVICE IUNITP (PARAMETERS), C AND IUNITV (VELOCITY SOLUTIONS). FOR EXAMPLE, IF THESE ARE IN C SI UNITS (SECONDS), THEN UNITT = 3.1688087E-08 DATA UNITT /3.1688087E-08/ C C FOLLOWING NUMBER IS "VERY BIG", AND PRINTS AS "****" WHEN OUTPUT C WITH AN F FORMAT: DATA STARS /1.E50/ C C---------------------------------------------------------------------- C C STATEMENT FUNCTIONS: C C INTERPOLATION WITHIN A TRIANGULAR ELEMENT: 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 SUM OF BASIS FUNCTIONS #1, #2, AND #3 (POSITIVE ONES) WITHIN A C FAULT ELEMENT: FHIVAL(S,F1,F2,F3)= + F1*(1.-3.*S+2.*S**2)+ + F2*(4.*S*(1.-S))+ + F3*(2.*S**2-S) C C SUM OF DERIVITIVES (D/DS) OF BASIS FUNCTIONS #1, #2, AND #3 C (POSITIVE ONES) WITHIN A FAULT ELEMENT: DERIF(S,F1,F2,F3)= + F1*(4.*S-3.)+ + F2*(4.-8.*S)+ + F3*(4.*S-1.) C---------------------------------------------------------------------- C C BEGIN EXECUTABLE CODE: 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 MXDATA=MAXDAT MXBEN=MAXBEN MXBAS=MAXBAS C ****************************************************************** C C WRITE HEADER ON OUTPUT: WRITE (IUNITT,1) 1 FORMAT ( +/' PROGRAM SCORER.GPB.ALASKA' +/' ' +/' COMPUTES SCORES OF FINITE ELEMENT VELOCITY MODELS FROM' +/' PROGRAM -PLATES- BY USING SLIP-RATE, STRESS, AND GEODETIC DATA' +/' ' +/' THIS VERSION CUSTOMIZED FOR THE ALASKAN PROBLEM IN THE PLACES' +/' MARKED BY C****** ALASKA (BEGIN SECTION) ******************' +/' C****** ALASKA (END SECTION) ******************' +/' ' +/' BY' +/' PETER BIRD,' +/' DEPARTMENT OF EARTH AND SPACE SCIENCES,' +/' UNIVERSITY OF CALIFORNIA, LOS ANGELES, CALIFORNIA 90024' +/' (WITH ASSISTANCE FROM XIANGHONG KONG).' +/' VERSION OF MAY 6, 1993' +/' --------------------------------------------------------' +/) C WEDGE=ABS(90.-ABS(DIPMAX))*0.017453293 C BLANKS=' ' DO 2 I=2,59 BLANKS=BLANKS//' ' 2 CONTINUE C C C INPUT FINITE ELEMENT GRID AND DATA VALUES AT NODE POINTS C CALL GETNET (INPUT,IUNITF,IUNITT, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP, + NFAKEN,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,OFFMAX,OFFSET, + TITLE1,TLNODE,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 READ TEST DATASET OF GEOLOGIC SLIP-RATE LIMITS FOR FAULTS C NFDATA=0 DO 106 I=1,MXDATA READ (IUNITD,105,END=111) + IFDATM(I),DATUMS(I),IENDND(I),RLAT(1,I),RLAT(2,I), + DELVZ(1,I),DELVZ(2,I),TAGGSR(I) 105 FORMAT (I10,F10.2,I10,4F10.2,2X,A59) NFDATA=NFDATA+1 106 CONTINUE WRITE (IUNITT,107) MXDATA 107 FORMAT (/' ERR0R: PARAMETER MAXDAT = ',I6,' IS NOT LARGE'/ + ' ENOUGH FOR THE NUMBER OF GEOLOGIC SLIP-RATE DATA.') STOP 111 IF (NFDATA.EQ.0) THEN WRITE (IUNITT,112) IUNITD 112 FORMAT (/' NO TEST DATA COULD BE READ FROM UNIT ',I2) STOP ELSE WRITE (IUNITT,113) NFDATA,IUNITD 113 FORMAT (/' ',I4,' DATA ON FAULT SLIP RATES', + ' WERE READ FROM UNIT ',I2,':'/) WRITE (IUNITT,114) 114 FORMAT (30X,'-IN MILLIMETERS/YEAR-'/ + ' FAULT POSITION TOWARD MIN. MAX. MIN. MAX.'/ + ' ELEMNT FRACTION NODE R.LAT. R.LAT. DELVZ DELVZ') DO 150 I=1,NFDATA IF (RLAT(1,I).NE.0.0) THEN RLT1=RLAT(1,I) ELSE RLT1=STARS ENDIF IF (RLAT(2,I).NE.0.0) THEN RLT2=RLAT(2,I) ELSE RLT2=STARS ENDIF IF (DELVZ(1,I).NE.0.0) THEN DZT1=DELVZ(1,I) ELSE DZT1=STARS ENDIF IF (DELVZ(2,I).NE.0.0) THEN DZT2=DELVZ(2,I) ELSE DZT2=STARS ENDIF WRITE(IUNITT,145) + IFDATM(I),DATUMS(I),IENDND(I), + RLT1,RLT2,DZT1,DZT2, + TAGGSR(I) 145 FORMAT(' ',I6,F10.3,I8,1X,4F7.3,2X,A59) 150 CONTINUE IF ((DATUMS(I).LT.0.).OR.(DATUMS(I).GT.1.)) THEN WRITE (IUNITT,160) I,DATUMS(I) 160 FORMAT (/' ERR0R IN TEST DATUM ',I4,':'/ + ' S VALUE OF ',F10.6,' IS OUTSIDE ALLOWED ', + 'RANGE OF 0.0 TO 1.0') STOP ENDIF WRITE (IUNITT,170) 170 FORMAT (/ ' NOTE: DELVZ VALUES ARE ALWAYS POSITIVE IN ', + 'CASES OF CONVERGENCE (THRUSTING),'/ + ' AND NEGATIVE IN CASES OF DIVERGENCE ', + '(NORMAL FAULTING).') ENDIF C NCHOP=N1000-(NREALN+1) DO 190 I=1,NFDATA IF (IENDND(I).GT.NREALN) IENDND(I)=IENDND(I)-NCHOP 190 CONTINUE C C READ IN DATA ON ARGUMENT (ANGLE COUNTERCLOCKWISE FROM +X, IN C DEGREES) OF DIRECTION OF GREATEST HORIZONTAL COMPRESSION; C ALSO ELEMENT NUMBER, S1, S2, DATA QUALITY, TAG-NAME, C AND TYPE OF STRESS FIELD (I.E., THRUST). C NSDATA=0 WRITE (IUNITT,201) 201 FORMAT (/ /'-----------------------------------', + '-----------------------------------'/ + / ' DATA ON STRESS DIRECTIONS (IN DEGREES COUNTER', + 'CLOCKWISE FROM +X)'/ / + ' SITE ELE# S1 S2 A-E TYPE ARG REGIME', + ' DEPTH'/) 210 READ (IUNITS,211,ERR=212,END=299) SITE, IE, S1, S2, QUALIT, + TYPE, ARG, REGIME, DEPTH, REFERE 211 FORMAT (A6,1X,I4,2F6.3,1X,A2,3X,A6,2X,F4.0,2X,A4,3X, + A7,2X,A28) 212 WRITE (IUNITT,213) SITE, IE, S1, S2, QUALIT, TYPE, ARG, + REGIME, DEPTH, REFERE 213 FORMAT (' ',A6,1X,I4,2F6.3,1X,A2,3X,A6,2X,F4.0,2X,A4,3X, + A7,2X,A28) NSDATA=NSDATA+1 IF (NSDATA.GT.MXDATA) THEN WRITE (IUNITT,220) MXDATA 220 FORMAT (/' ERR0R: PARAMETER MAXDAT =',I6,' IS NOT LARGE'/ + ' ENOUGH FOR THE NUMBER OF STRESS DATA.') STOP ENDIF IESTR(NSDATA)=IE ARGSTR(NSDATA)=ARG TAGSTR(NSDATA)=SITE REGSTR(NSDATA)=REGIME IF (QUALIT(1:1).EQ.'A') THEN WGTSTR(NSDATA)=1.00 ELSE IF (QUALIT(1:1).EQ.'B') THEN WGTSTR(NSDATA)=0.75 ELSE IF (QUALIT(1:1).EQ.'C') THEN WGTSTR(NSDATA)=0.50 ELSE IF (QUALIT(1:1).EQ.'D') THEN WGTSTR(NSDATA)=0.25 ELSE WGTSTR(NSDATA)=0.00 ENDIF IF (REGIME.EQ.'N ') THEN IASIGN(1,NSDATA)= 0 IASIGN(2,NSDATA)= +1 IASIGN(3,NSDATA)= -1 ELSE IF (REGIME.EQ.'SS ') THEN IASIGN(1,NSDATA)= -1 IASIGN(2,NSDATA)= +1 IASIGN(3,NSDATA)= 0 ELSE IF (REGIME.EQ.'T ') THEN IASIGN(1,NSDATA)= -1 IASIGN(2,NSDATA)= 0 IASIGN(3,NSDATA)= +1 ELSE IF (REGIME.EQ.'SS/T') THEN IASIGN(1,NSDATA)= -1 IASIGN(2,NSDATA)= +1 IASIGN(3,NSDATA)= +1 ELSE IF (REGIME.EQ.'T/SS') THEN IASIGN(1,NSDATA)= -1 IASIGN(2,NSDATA)= +1 IASIGN(3,NSDATA)= +1 ELSE IF (REGIME.EQ.'N/SS') THEN IASIGN(1,NSDATA)= -1 IASIGN(2,NSDATA)= +1 IASIGN(3,NSDATA)= -1 ELSE IF (REGIME.EQ.'SS/N') THEN IASIGN(1,NSDATA)= -1 IASIGN(2,NSDATA)= +1 IASIGN(3,NSDATA)= -1 ELSE IASIGN(1,NSDATA)= 0 IASIGN(2,NSDATA)= 0 IASIGN(3,NSDATA)= 0 ENDIF C C READ ANY EXTENSIONS OF THE REFERENCE ONTO ANOTHER RECORD C 240 READ (IUNITS,242,END=299) TEST,MORERE 242 FORMAT (A59,A28) IF (TEST.EQ.BLANKS) THEN WRITE (IUNITT,243) BLANKS, MORERE 243 FORMAT(' ',A59,A28) GO TO 240 ELSE BACKSPACE (IUNITS) ENDIF C GO TO 210 C 299 CONTINUE C C READ IN LOCATIONS OF GEODETIC BENCHMARKS C WRITE (IUNITT,300) IUNITG 300 FORMAT (/' GEODETIC BENCHMARK LOCATIONS FROM INPUT UNIT ',I3/) C (FIRST, TWO HEADER LINES ARE COPIED) READ (IUNITG,301) TITL1G 301 FORMAT (A80) WRITE (IUNITT,302) TITL1G 302 FORMAT (' ',A80) READ (IUNITG,301) TITL2G WRITE (IUNITT,302) TITL2G C (NOW, WE READ UNTIL '=====' IS ENCOUNTERED IN BENCHMARK NAME FIELD) NBENCH=0 310 READ (IUNITG,311,END=320,ERR=312) BENCHN,IE,S1,S2 311 FORMAT ( A31,1X,I6,1X,F6.3,1X,F6.3) 312 IF (BENCHN(1:5).EQ.'=====') GO TO 320 NBENCH=NBENCH+1 IF (NBENCH.GT.MXBEN) THEN WRITE (IUNITT,313) MXBEN 313 FORMAT (' ERR0R: NUMBER OF BENCHMARKS IN INPUT FILE'/ + ' EXCEEDS PARAMETER MAXBEN =',I5) STOP ENDIF BENCH(NBENCH)=BENCHN WRITE (IUNITT,315) BENCHN,IE,S1,S2 315 FORMAT (' ',A31,1X,I6,1X,F6.3,1X,F6.3) IBENCH(NBENCH)=IE BENS1(NBENCH)=S1 BENS2(NBENCH)=S2 S3=1.00-S1-S2 N1=NODES(1,IE) N2=NODES(2,IE) N3=NODES(3,IE) N4=NODES(4,IE) N5=NODES(5,IE) N6=NODES(6,IE) X1=XNODE(N1) X2=XNODE(N2) X3=XNODE(N3) X4=XNODE(N4) X5=XNODE(N5) X6=XNODE(N6) Y1=YNODE(N1) Y2=YNODE(N2) Y3=YNODE(N3) Y4=YNODE(N4) Y5=YNODE(N5) Y6=YNODE(N6) BENCHX(NBENCH)=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) BENCHY(NBENCH)=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) GO TO 310 C C READ IN DEFINITIONS OF BASELINES, THEIR RATES, AND ERR0RS C 320 WRITE (IUNITT,321) IUNITG 321 FORMAT (/' GEODETIC BASLINE RATES FROM INPUT UNIT ',I3/) C (FIRST, TWO HEADER LINES ARE COPIED) READ (IUNITG,301) TITL3G WRITE (IUNITT,302) TITL3G READ (IUNITG,301) TITL4G WRITE (IUNITT,302) TITL4G C (NOW, WE READ UNTIL '=====' IS ENCOUNTERED) NBASE=0 330 READ (IUNITG,331,END=350,ERR=332) SITE1N,SITE2N,RATE,ERR 331 FORMAT (A31,1X,A31,1X,F6.0,1X,F6.0) 332 IF (SITE1N(1:5).EQ.'=====') GO TO 350 NBASE=NBASE+1 IF (NBASE.GT.MXBAS) THEN WRITE (IUNITT,333) MXBAS 333 FORMAT (' ERR0R: NUMBER OF BASELINES IN INPUT FILE'/ + ' EXCEEDS PARAMETER MAXBAS =',I5) STOP ENDIF SITE1(NBASE)=SITE1N SITE2(NBASE)=SITE2N BRATE(NBASE)=RATE BERR(NBASE)=ERR WRITE (IUNITT,334) SITE1N,SITE2N,RATE,ERR 334 FORMAT (' ',A31,1X,A31,1X,F6.2,1X,F6.2) BRATE(NBASE)=BRATE(NBASE)*UNITT/UNITL BERR(NBASE)=BERR(NBASE)*UNITT/UNITL OK1=.FALSE. DO 340 I=1,NBENCH IF (SITE1(NBASE).EQ.BENCH(I)) THEN OK1=.TRUE. INDFST(NBASE)=I ENDIF 340 CONTINUE OK2=.FALSE. DO 345 I=1,NBENCH IF (SITE2(NBASE).EQ.BENCH(I)) THEN OK2=.TRUE. INDLST(NBASE)=I ENDIF 345 CONTINUE IF ((.NOT.OK1).OR.(.NOT.OK2)) THEN NBASE=NBASE-1 WRITE (IUNITT,346) 346 FORMAT (' BASELINE ABOVE WAS NOT USED BECAUSE', + ' BENCHMARK(S) NOT RECOGNIZABLE.') ENDIF GO TO 330 350 CONTINUE C C ==============BEGINNING OF INDEFINATE LOOP IN INPUT MODELS ========= C C READ SOLUTION (NODE VELOCITIES, PRECEDED BY THREE TITLES) C 1001 CALL OLDVEL (INPUT,IUNITV,MXNODE,NUMNOD, + OUTPUT,ALDONE,TITLE1,TITLE2,TITLE3,V) IF (ALDONE) STOP CALL EDOT (INPUT,DXS,DYS,MXEL,MXNODE,NODES,NUMEL,V, + OUTPUT,ERATE) C WRITE (IUNITT,1002) TITLE1,TITLE2,TITLE3 1002 FORMAT (/ /' ===================================', + '==================================='/ / + ' ',A80/' ',A80/' ',A80/) C C READ PARAMETERS C CALL READPM (INPUT,IUNITP, IUNITT, OFFMAX, + OUTPUT,ACREEP, ALPHAT, BCREEP, BIOT, + BYERLY, CCREEP, CFRIC , CONDUC, + DCREEP, ECREEP, EVERYP, FFRIC , GMEAN , + MAXITR, OKDELV, OKTOQT, ONEKM, RADIO , + REFSTR, RHOAST, RHOBAR, RHOH2O, + TEMLIM, TITLE3, TRHMAX, TSURF) C C COMPUTE SLIP RATES (IN MM/YEAR) AT SPECIFIED POINTS C AND ASSIGN SCORES AND PROBABILITIES C WRITE (IUNITT,1100) 1100 FORMAT (/' GEOLOGIC SLIPRATE PREDICTIONS AND THEIR ERR0RS:'/) WRITE (IUNITT,114) WRITE (IUNITT,1110) 1110 FORMAT('+',55X,'MODEL-R.LAT ERR0R MODEL-DELVZ ERR0R', + ' CHANCES-OF-CONSISTENCY-W/-GEOLOGY'/) FACTOR=1. NGSDAT=0 NB=0 SUMB=0. SUMB2=0. BIGB=0. DO 1190 J=1,NFDATA I=IFDATM(J) IF ((IENDND(J).EQ.NODEF(3,I)) .OR. + (IENDND(J).EQ.NODEF(4,I))) THEN S=DATUMS(J) ELSE IF ((IENDND(J).EQ.NODEF(1,I)) .OR. + (IENDND(J).EQ.NODEF(6,I))) THEN S=1.00-DATUMS(J) ELSE WRITE (IUNITT,1120) J,IENDND(J),I 1120 FORMAT (/' ERR0R IN GEOLOGIC SLIPRATE DATUM ',I4,':'/ + ' NODE ',I6,' IS NOT AN END NODE OF FAULT ', + I4) STOP ENDIF C N1=NODEF(1,I) N2=NODEF(2,I) N3=NODEF(3,I) N4=NODEF(4,I) N5=NODEF(5,I) N6=NODEF(6,I) C X1=XNODE(N1) X2=XNODE(N2) X3=XNODE(N3) X4=XNODE(N4) X5=XNODE(N5) X6=XNODE(N6) C Y1=YNODE(N1) Y2=YNODE(N2) Y3=YNODE(N3) Y4=YNODE(N4) Y5=YNODE(N5) Y6=YNODE(N6) C VX1=V(1,N1) VX2=V(1,N2) VX3=V(1,N3) VX4=V(1,N4) VX5=V(1,N5) VX6=V(1,N6) C VY1=V(2,N1) VY2=V(2,N2) VY3=V(2,N3) VY4=V(2,N4) VY5=V(2,N5) VY6=V(2,N6) C DELVX=FHIVAL(S,VX1,VX2,VX3)-FHIVAL(S,VX6,VX5,VX4) DELVY=FHIVAL(S,VY1,VY2,VY3)-FHIVAL(S,VY6,VY5,VY4) DXDS=DERIF(S,X1,X2,X3) DYDS=DERIF(S,Y1,Y2,Y3) ARG=ATAN2(DYDS,DXDS) UNITX=COS(ARG) UNITY=SIN(ARG) CROSSX= -UNITY CROSSY= UNITX SINIST=DELVX*UNITX+DELVY*UNITY DEXTRA= -SINIST C C NOTE CONVERSION OF UNITS TO MM/YEAR IN NEXT LINE RLT=DEXTRA*UNITL/UNITT C CLOSES=DELVX*CROSSX+DELVY*CROSSY OPENS= -CLOSES DIP=FHIVAL(S,FDIP(1,I),FDIP(2,I),FDIP(3,I)) IF (ABS(DIP-1.570796).LT.WEDGE) THEN RELVZ=0. ELSE C C NOTE CONVERSION OF UNITS TO MM/YEAR IN NEXT LINE RELVZ=CLOSES*ABS(TAN(DIP))*UNITL/UNITT C ENDIF C************ C WRITE(6,777) I,N1,N2,N3,N4,N5,N6, C + X1,X2,X3,X4,X5,X6, C + Y1,Y2,Y3,Y4,Y5,Y6, C + VX1,VX2,VX3,VX4,VX5,VX6, C + VY1,VY2,VY3,VY4,VY5,VY6, C + S,DIP,ARG, C + DELVX,DELVY,RLT,RELVZ C 777 FORMAT(/' I=',I3/' N=',6I10/' X=',1P,6E10.2/' Y=',6E10.2/ C + ' VX=',6E10.2/' VY=',6E10.2/ C + ' S=',0P,F10.5,' DIP=',F10.5,' ARG=',F10.5/ C + ' DELVX=',1P,E10.2,' DELVY=',E10.2,' RLT=',E10.2, C + ' RELVZ=',E10.2) C************ C CALL LIKELY (INPUT,RLAT(1,J),RLAT(2,J),RLT, + OUTPUT,PROBRL) CALL LIKELY (INPUT,DELVZ(1,J),DELVZ(2,J),RELVZ, + OUTPUT,PROBVZ) C FACTOR=FACTOR*PROBRL*PROBVZ C RERR=0. IF ((RLAT(1,J)+RLAT(2,J)).NE.0.) THEN NGSDAT=NGSDAT+1 IF ((RLT.LT.RLAT(1,J)).AND.(RLAT(1,J).NE.0.)) THEN NB=NB+1 RERR=RLAT(1,J)-RLT SUMB=SUMB+RERR SUMB2=SUMB2+RERR**2 BIGB=AMAX1(BIGB,RERR) ENDIF IF ((RLT.GT.RLAT(2,J)).AND.(RLAT(2,J).NE.0.)) THEN NB=NB+1 RERR=RLT-RLAT(2,J) SUMB=SUMB+RERR SUMB2=SUMB2+RERR**2 BIGB=AMAX1(BIGB,RERR) ENDIF ENDIF C ZERR=0. IF ((DELVZ(1,J)+DELVZ(2,J)).NE.0.) THEN NGSDAT=NGSDAT+1 IF ((RELVZ.LT.DELVZ(1,J)).AND.(DELVZ(1,J).NE.0.)) THEN NB=NB+1 ZERR=DELVZ(1,J)-RELVZ SUMB=SUMB+ZERR SUMB2=SUMB2+ZERR**2 BIGB=AMAX1(BIGB,ZERR) ENDIF IF ((RELVZ.GT.DELVZ(2,J)).AND.(DELVZ(2,J).NE.0.)) THEN NB=NB+1 ZERR=RELVZ-DELVZ(2,J) SUMB=SUMB+ZERR SUMB2=SUMB2+ZERR**2 BIGB=AMAX1(BIGB,ZERR) ENDIF ENDIF C C PREPARE ONE LINE IN THE TABLE C IF(IENDND(J).GT.NREALN) THEN IEND=IENDND(J)+NCHOP ELSE IEND=IENDND(J) ENDIF C IF (RLAT(1,J).NE.0.) THEN PRLAT1=RLAT(1,J) ELSE PRLAT1=STARS ENDIF C IF (RLAT(2,J).NE.0.) THEN PRLAT2=RLAT(2,J) ELSE PRLAT2=STARS ENDIF C IF (DELVZ(1,J).NE.0.) THEN PDELV1=DELVZ(1,J) ELSE PDELV1=STARS ENDIF C IF (DELVZ(2,J).NE.0.) THEN PDELV2=DELVZ(2,J) ELSE PDELV2=STARS ENDIF C IF ((RLAT(1,J)+RLAT(2,J)).EQ.0.) THEN XPRLT=STARS ELSE XPRLT=RLT ENDIF C IF (XPRLT.GE.STARS) THEN XPRERR=STARS ELSE XPRERR=RERR ENDIF C IF ((DELVZ(1,J)+DELVZ(2,J)).EQ.0.) THEN XPRELV=STARS ELSE XPRELV=RELVZ ENDIF C IF (XPRELV.GE.STARS) THEN XPZERR=STARS ELSE XPZERR=ZERR ENDIF C IF ((PRLAT1+PRLAT2).GT.STARS) THEN PPRL=STARS ELSE PPRL=PROBRL ENDIF C IF ((PDELV1+PDELV2).GT.STARS) THEN PPVZ=STARS ELSE PPVZ=PROBVZ ENDIF C WRITE(IUNITT,1180)IFDATM(J),DATUMS(J),IEND,PRLAT1,PRLAT2, + PDELV1,PDELV2,XPRLT,XPRERR,XPRELV,XPZERR,PPRL,PPVZ 1180 FORMAT(' ',I6,F10.3,I8,1X,4F7.3, + 5X,F7.3,2X,F7.3,4X,F7.3,2X,F7.3,10X,F7.3,10X,F7.3) 1190 CONTINUE WRITE (IUNITT,170) C C SUMMARY STATISTICS FOR THIS MODEL C WRITE (IUNITT,1192) TITLE1,TITLE2,TITLE3 1192 FORMAT (/ / + ' ',A80/' ',A80/' ',A80/) ENB=FLOAT(NB)/FLOAT(NGSDAT) SUMB=SUMB/FLOAT(NGSDAT) SUMB2=SQRT(SUMB2/FLOAT(NGSDAT)) WRITE (IUNITT,1198) ENB,SUMB,SUMB2,BIGB,FACTOR 1198 FORMAT(/ + ' FRACTION OF PREDICTIONS IN ERR0R =',F6.3/ + ' AVERAGE ERR0R OF ALL PREDICTIONS =',F10.3,' MM/YEAR'/ + ' ROOT-MEAN-SQUARE ERR0R OF ALL PREDICTIONS =',F10.3,' MM/YEAR'/ + ' LARGEST ERR0R OF ANY PREDICTION =',F10.3,' MM/YEAR'/ + ' PROBABILITY THAT MODEL INPUT PARAMETERS AND PHYSICAL EQUATIO', + 'NS ARE CONSISTENT'/' WITH THE TRUE GEOLOGY REPRESENTED BY', + ' THIS DATASET IS ',1PE11.3) C RESULT=FACTOR**(1./NGSDAT) WRITE (IUNITT,1199) RESULT 1199 FORMAT( + ' ',12('####')/' ### ESTIMATE, BASED ON THIS DATASET,OF', + ' ###'/' ### PROBABILITY THAT, IN A GIVEN CASE, ###'/ + ' ### THESE PHYSICAL EQUATIONS AND PARAMETER ###'/ + ' ### VALUES LEAD TO A PREDICTED SLIP RATE ###'/ + ' ### CONSISTENT WITH GEOLOGIC DATA (WITHIN ###'/ + ' ### ITS UNCERTAINTIES) IS ',F7.4,' ###'/ + ' ',12('####')) C C COMPUTE STRESS-DIRECTION SCORES AND PRINT TABLE C WRITE (IUNITT,1201) 1201 FORMAT (/ /' -----------------------------------', + '-----------------------------------'/ + /' SCORES WITH REGARD TO STRESS DIRECTION DATA:'/ + /' DATA DATA DATA DATA DATA', + ' MODEL MODEL MODEL MODEL'/ + ' SITE ELEMENT ARG(S1H) REGIME E-SIGNS', + ' E-SIGNS #WRONG ARG(S1H) ERR0R WEIGHT'/) 1202 FORMAT (' ',A6,I8,F9.0,3X,A4,3I3,3I3,I7,F9.0,F6.0,F7.2) SUMA=0.0 SUMS=0.0 DENOMA=0.0 DENOMS=0.0 M=1 DO 1290 I=1,NSDATA IE=IESTR(I) EXX=ERATE(1,M,IE) EYY=ERATE(2,M,IE) EXY=ERATE(3,M,IE) CALL PRINCE (INPUT,EXX,EYY,EXY, + OUTPUT,E1,E2,U1X,U1Y,U2X,U2Y) ARGMOD=ATAN2F(U1Y,U1X)*57.2957795 EZ= -(EXX+EYY) IF (E1.LT.0.) THEN IMSIGN(1)= -1 ELSE IF (E1.GT.0.) THEN IMSIGN(1)= +1 ELSE IMSIGN(1)= 0 ENDIF IF (E2.LT.0.) THEN IMSIGN(2)= -1 ELSE IF (E2.GT.0.) THEN IMSIGN(2)= +1 ELSE IMSIGN(2)= 0 ENDIF IF (EZ.LT.0.) THEN IMSIGN(3)= -1 ELSE IF (EZ.GT.0.) THEN IMSIGN(3)= +1 ELSE IMSIGN(3)= 0 ENDIF DEGERR=ABS(ARGMOD-ARGSTR(I)) IF (DEGERR.GT.180.) DEGERR=DEGERR-180. IF (DEGERR.GT.180.) DEGERR=DEGERR-180. IF (DEGERR.GT.180.) DEGERR=DEGERR-180. IF (DEGERR.GT.180.) DEGERR=DEGERR-180. IF (DEGERR.GT.90.) DEGERR=180.-DEGERR NWRONG=0 DO 1270 K=1,3 IF (IASIGN(K,I).NE.0) THEN IF (IMSIGN(K).NE.IASIGN(K,I)) THEN NWRONG=NWRONG+1 ENDIF ENDIF 1270 CONTINUE WRITE (IUNITT,1202) TAGSTR(I),IESTR(I),ARGSTR(I), + REGSTR(I),(IASIGN(K,I),K=1,3), + (IMSIGN(K),K=1,3),NWRONG,ARGMOD, + DEGERR,WGTSTR(I) SUMS=SUMS+NWRONG*WGTSTR(I) SUMA=SUMA+DEGERR*WGTSTR(I) DENOMA=DENOMA+WGTSTR(I) DO 1280 K=1,3 IF (IASIGN(K,I).NE.0) DENOMS=DENOMS+WGTSTR(I) 1280 CONTINUE 1290 CONTINUE IF (DENOMA.GT.0.) THEN AVEDEG=SUMA/DENOMA ELSE AVEDEG=STARS ENDIF IF (DENOMS.GT.0.) THEN AVEBAD=SUMS/DENOMS ELSE AVEBAD=STARS ENDIF WRITE (IUNITT,1192) TITLE1,TITLE2,TITLE3 WRITE (IUNITT,1299) AVEDEG, AVEBAD 1299 FORMAT (/ /' WEIGHTED MEAN ERR0R IN STRESS ORIENTATION = ', + F6.2,' DEGREES.'/ / + ' WEIGHTED MEAN NUMBER OF INCORRECT AXIS SIGNS = ', + F7.5/ /) C C COMPUTE GEODETIC LENGTH CHANGE RATE SCORES C CALL GEOSCO(INPUT,ACREEP,ALPHAT,AREA,BCREEP,BIOT, + BYERLY,CCREEP,CFRIC,CONDUC, + DCREEP,DQDTDA,DETJ,DXS,DYS,ECREEP,FAZ,FDIP, + FFRIC,FTAN,BENS1,BENS2,BENCHX,BENCHY, + BERR,BRATE,GMEAN,INDFST,INDLST,IUNITT, + MXEL,MXFEL,MXNODE,IBENCH,NFL, + NBENCH,NBASE, + NODEF,NODES,NUMEL, + OFFMAX,OFFSET,OKDELV,ONEKM, + RADIO,REFSTR,RHOH2O,RHOBAR,SITE1,SITE2, + TITLE1,TITLE2,TITLE3, + TLNODE,TRHMAX,TSURF,UNITL,UNITT, + UX,UY,V,VG,VGCRE,WEDGE, + XNODE,YNODE,ZMNODE,ZMOHO, + WORK,NFTRE, + OUTPUT,FC,FIMUDZ,FPEAKS,FSLIPS,FTSTAR, + SMSCOR,TLINT,ZTRANF) C C---------------------------------------------- C C PRINT SUMMARY SCORE FOR MODEL: C EPSILO=( SUMB2+ + SMSCOR+ + (AVEDEG-17.)/2. )/3. WRITE (IUNITT,2001) TITLE1,TITLE2,TITLE3, + SUMB2,AVEDEG,SMSCOR,EPSILO 2001 FORMAT(/' ================================================' + /' BOTTOM LINE:' + /' ' + /' ',A80 + /' ',A80 + /' ',A80 + /' ' + /' RMS GEOLOGIC SLIP-RATE ERR0R= ',F10.3,' MM/YEAR' + /' WEIGHTED MEAN STRESS ERR0R = ',F10.3,' DEGREES' + /' WEIGHTED RMS GEODETIC ERR0R = ',F10.3,' MM/YEAR' + /' ' + /' EPSILON (COMBINED ERR0R) = ',F10.3,' MM/YEAR' + /' ================================================') C C================= END OF INDEFINATE LOOP ON DATASETS ================= GO TO 1001 C 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 REAL FUNCTION ATAN2F (Y,X) C C CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS C IF ((Y.NE.0.).OR.(X.NE.0.)) THEN ATAN2F=ATAN2(Y,X) ELSE ATAN2F=0. ENDIF RETURN END C C 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 SUBROUTINE EDOT (INPUT,DXS,DYS,MXEL,MXNODE,NODES,NUMEL,V, + OUTPUT,ERATE) C C COMPUTE STRAIN-RATE COMPONENTS EDOTXX, EDOTYY, AND C EDOTXY (TENSOR FORM; EQUAL TO C (1/2) * ((DVX/DY)+(DVY/DX)) C AT THE INTEGRATION POINTS OF TRIANGULAR CONTINUUM ELEMENTS. C DOUBLE PRECISION V DIMENSION DXS(6,7,MXEL),DYS(6,7,MXEL), + ERATE(3,7,MXEL),NODES(6,MXEL), + V(2,MXNODE) C DO 1000 M=1,7 DO 900 I=1,NUMEL EXX=0. EYY=0. EXY=0. DO 800 J=1,6 NODE=NODES(J,I) VX=V(1,NODE) VY=V(2,NODE) DX=DXS(J,M,I) DY=DYS(J,M,I) EXX=EXX+VX*DX EYY=EYY+VY*DY EXY=EXY+(VX*DY+VY*DX)*0.5 800 CONTINUE ERATE(1,M,I)=EXX ERATE(2,M,I)=EYY ERATE(3,M,I)=EXY 900 CONTINUE 1000 CONTINUE RETURN END C C C SUBROUTINE GDTRAN(INPUT,FDIP,BENCHX,BENCHY,MXNODE,NBENCH, + NFL,NODEF,V,WEDGE, + XNODE,YNODE,ZTRANF, + WORK, NFTRE, + OUTPUT, UX,UY) PARAMETER (NSEGME=3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL BENCHX,BENCHY,FDIP,OUTPUT, + UX,UY,WEDGE,WORK,XNODE,YNODE,ZTRANF LOGICAL PROBLE,DIS05 DIMENSION FDIP(3,NFL), BENCHX(NBENCH), BENCHY(NBENCH), + NODEF(6,NFL), NOUTL(10), NOUTR(10), NFTRE(MXNODE), + SLIP(3), UX(NBENCH), UY(NBENCH) DIMENSION XNODE(MXNODE),YNODE(MXNODE),V(2,MXNODE), + ZTRANF(7,NFL) NSEG=NSEGME P18=0.017453292 DP=0.0D0 PI2=3.14159 UN=1.0 NMID=NSEG/2+1 DO 100 I=1,NBENCH UX(I)=0.0 UY(I)=0.0 100 CONTINUE DO 550 I=1,NFL I1=NODEF(1,I) I2=NODEF(2,I) I3=NODEF(3,I) I4=NODEF(4,I) I5=NODEF(5,I) I6=NODEF(6,I) XX1=XNODE(I1) XX2=XNODE(I2) XX3=XNODE(I3) YY1=YNODE(I1) YY2=YNODE(I2) YY3=YNODE(I3) XBEGIN=XX1 YBEGIN=YY1 DO 540 M=1,NSEG SS=REAL(M)/REAL(NSEG) SMID=(REAL(M)-0.5)/REAL(NSEG) IF(M.EQ.NSEG) THEN XEND=XX3 YEND=YY3 ELSE PHI1=1-3.0*SS+2.0*SS**2 PHI2=4.0*SS*(1.0-SS) PHI3=2.0*SS**2-SS XEND=XX1*PHI1+XX2*PHI2+XX3*PHI3 YEND=YY1*PHI1+YY2*PHI2+YY3*PHI3 ENDIF DPHI1S=4.0*SMID-3.0 DPHI2S=4.0-8.0*SMID DPHI3S=4.0*SMID-1.0 DXDS=XX1*DPHI1S+XX2*DPHI2S+XX3*DPHI3S DYDS=YY1*DPHI1S+YY2*DPHI2S+YY3*DPHI3S ARG=ATAN2(DYDS,DXDS) AZIMHS=ARG PHI1M=1-3.0*SMID+2.0*SMID**2 PHI2M=4.0*SMID*(1.0-SMID) PHI3M=2.0*SMID**2-SMID DIP=PHI1M*FDIP(1,I)+PHI2M*FDIP(2,I)+PHI3M*FDIP(3,I) FDEPTH=PHI1M*ZTRANF(1,I)+PHI2M*ZTRANF(4,I)+ + PHI3M*ZTRANF(7,I) XMID=PHI1M*XX1+PHI2M*XX2+PHI3M*XX3 YMID=PHI1M*YY1+PHI2M*YY2+PHI3M*YY3 DUX=(V(1,I1)-V(1,I6))*PHI1M+(V(1,I2)-V(1,I5))*PHI2M + +(V(1,I3)-V(1,I4))*PHI3M DVY=(V(2,I1)-V(2,I6))*PHI1M+(V(2,I2)-V(2,I5))*PHI2M + +(V(2,I3)-V(2,I4))*PHI3M UNITX=COS(ARG) UNITY=SIN(ARG) CROSSX= -UNITY CROSSY= +UNITX IF (ABS(DIP-1.570796).LT.WEDGE) THEN CLOSE=0. VUPDIP=0. ELSE CLOSE=DUX*CROSSX+DVY*CROSSY RELV=CLOSE/TAN(DIP) ENDIF C C PREPARING PARAMETERS C XXOO=XEND YYOO=YEND ABX = XEND - XBEGIN ABY = YEND - YBEGIN AL = SQRT(ABX*ABX+ABY*ABY) C (AL IS THE LENGTH OF THE FAULT) AZIMF=AZIMHS ZTOP=DP DO 400 K=1,10 NOUTL(K)=0 NOUTR(K)=0 400 CONTINUE C DO 500 J=1,NBENCH C C TO DECIDE IF BENCHMARK IS ON THE FAULT NODE WHICH IS END OF C OF FAULT SEGMENT AND OF WHICH SLIP IS USING C IN THIS FAULT LOOP TO AVOID SINGULAR C IF(M.EQ.1) THEN DX=ABS(XX1-BENCHX(J)) DY=ABS(YY1-BENCHY(J)) ELSE IF(M.EQ.NMID) THEN DX=ABS(XX2-BENCHX(J)) DY=ABS(YY2-BENCHY(J)) ELSE IF(M.EQ.NSEG) THEN DX=ABS(XX3-BENCHX(J)) DY=ABS(YY3-BENCHY(J)) ENDIF PROBLE= (M.EQ.1.OR.M.EQ.NMID.OR.M.EQ.NSEG). + AND.(DX.LT.100.0.AND.DY.LT.100.0) IF(PROBLE) THEN IF ((DX*CROSSX+DY*CROSSY).GT.0.0) THEN BX=BX-100.0*CROSSX BY=BY-100.0*CROSSX ELSE BX=BX+100.0*CROSSX BY=BY+100.0*CROSSX ENDIF IF(M.EQ.1) THEN FL=AL SLIP(1)=0.5*(V(1,I1)-V(1,I6)) SLIP(2)=0.5*(V(2,I1)-V(2,I6)) XF=XX1 YF=YY1 ZBOT=ZTRANF(1,I) DIPF=FDIP(1,I) ELSE IF(M.EQ.NMID) THEN FL=0.5*AL SLIP(1)=DUX SLIP(2)=DVY XF=XMID YF=YMID ZBOT=FDEPTH DIPF=DIP ELSE IF(M.EQ.NSEG) THEN FL=AL SLIP(1)=0.5*(V(1,I3)-V(1,I4)) SLIP(2)=0.5*(V(2,I3)-V(2,I4)) XF=XX3 YF=YY3 ZBOT=ZTRANF(7,I) DIPF=FDIP(3,I) ENDIF IF (ABS(DIPF-1.570796).LT.WEDGE) THEN VUPDP1=0. ELSE CLOSE1=SLIP(1)*CROSSX+SLIP(2)*CROSSY VUPDP1=CLOSE1/TAN(DIPF) ENDIF SLIP(3)=VUPDP1 ELSE ZBOT=FDEPTH DIPF=DIP XF=XMID YF=YMID FL=0.5*AL SLIP(1)=DUX SLIP(2)=DVY SLIP(3)=RELV BX=BENCHX(J) BY=BENCHY(J) ENDIF C C TO DECIDE DIRECTION OF SLIP(3) BY THE DIP C C C COMPUTE THEORETICAL DEFORMATIONS AND PARTIALS FOR EACH STATION C C CALL CHANGE (INPUT,AZIMF,BX,BY,DIPF, + FL,SLIP,XF,YF,ZBOT,ZTOP, + OUTPUT,VXB,VYB,DIS05) UX(J) = UX(J) + VXB UY(J) = UY(J) + VYB 500 CONTINUE XBEGIN=XXOO YBEGIN=YYOO 540 CONTINUE 550 CONTINUE RETURN END C C C SUBROUTINE CHANGE (INPUT,AZIMF,BX,BY,DIPF, + LF,SLIP,XF,YF,ZBOT,ZTOP, + OUTPUT,VXB,VYB,DIS05) C C CONVERTS TO MANSINHA-AND-SMYLIE'S COORDINATE SYSTEM C AND USES THE APPROPRIATE DISLOCATION ROUTINE FOR DISPLACEMENTS. C C CONVENTIONS: XF AND YF ARE COORDINATES OF THE MIDPOINT OF THE C TRACE OF THE FAULT SEGMENT. BX AND BY ARE COORDINATES C OF THE BENCHMARK WHERE THE DISPLACEMENT IS NEEDED. C BOTH ARE IN THE GRID COORDINATE C SYSTEM (Y-DIRECTION 90 DEG COUNTERCLOCKWISE FROM X-DIRECTION) C LF IS THE HALF-LENGTH (CENTER TO END) OF THE TRACE OF THE C FAULT SEGMENT. C AZIMF IS THE ARGUMENT OF THE FAULT TRACE, MEASURED IN RADIANS C COUNTERCLOCKWISE FROM X. C DIPF IS FAULT DIP IN RADIANS MEASURED CLOCKWISE C (INITIALLY,DOWN) FROM HORIZONTAL ON THE RIGHT SIDE OF THE C FAULT (WHEN LOOKING IN DIRECTION AZIMF). C ZTOP AND ZBOT ARE DEPTHS TO TOP AND BOTTOM OF THE DISLOCATION, C SLIP IS THE VECTOR OF MOTION OF THE RIGHT SIDE OF THE FAULT C (SEE ABOVE) RELATIVE TO THE LEFT SIDE (") WITH THE VERTICAL C (THIRD) COMPONENT POSITIVE UPWARD; UNITS SAME AS VXB AND VYB. C VXB AND VYB ARE THE BENCHMARK DISPLACEMENTS OR VELOCITIES. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL OUTPUT INTEGER INPUT DOUBLE PRECISION LF DIMENSION SLIP(3),TSLIP(3),UCAP(3),Z(3),U(3) LOGICAL VERT,DIS05 C IF (DIPF.LE.1.5708) THEN TAZIMF=AZIMF THETA=DIPF DO 10 I=1,3 TSLIP(I)=SLIP(I) 10 CONTINUE ELSE TAZIMF=AZIMF+3.141593 THETA=3.141593-DIPF DO 20 I=1,3 TSLIP(I)= -SLIP(I) 20 CONTINUE ENDIF SINAZ=SIN(TAZIMF) COSAZ=COS(TAZIMF) Z(1)=((BX-XF)*COSAZ+(BY-YF)*SINAZ) Z(2)=((BX-XF)*SINAZ-(BY-YF)*COSAZ) Z(3)=0. DTOP=ZTOP/SIN(THETA) C**********SPECIAL STATEMENT TO APPROXIMATE SURFACE FAULT-BREAKOUT*** DTOP=MAX(DTOP,(ABS(Z(2))/20.)) C******************************************************************** DBOT=ZBOT/SIN(THETA) VERT=ABS(THETA-1.5708).LE.0.16 DEP05=0.05*ZBOT DIS05=(ABS(Z(2)).LT.DEP05) IF(DIS05) THEN Z(2)=DEP05 ENDIF UCAP(1)=TSLIP(1)*COSAZ+TSLIP(2)*SINAZ UCAP(2)=TSLIP(1)*SINAZ-TSLIP(2)*COSAZ C IF (VERT) THEN UCAP(2)=0. UCAP(3)=0. CALL HALO (INPUT,UCAP,Z,LF,DBOT,DTOP, + OUTPUT,U) ELSE UCAP(3)= -TSLIP(3) CALL AURA (INPUT,UCAP,THETA,Z,LF,DBOT,DTOP, + OUTPUT,U) ENDIF C VXB= U(1)*COSAZ+U(2)*SINAZ VYB= U(1)*SINAZ-U(2)*COSAZ C RETURN END C C C SUBROUTINE HALO (INPUT,UCAP,Z,L,DBOT,DTOP, + OUTPUT,U) C C COMPUTES DISPLACEMENTS IN A UNIFORM ELASTIC HALF-SPACE C OF POISSON'S RATION 0.25 DUE TO UNIFORM SLIP ON A C RECTANGULAR DISLOCATION IN THE VERTICAL PLANE. C C EQUATIONS BASED ON MANSINHA AND SMYLIE (1967), C J. GEOPHYS. RES., 72, 4731-4743. C IMPLICIT DOUBLE PRECISION (A-Z) REAL OUTPUT INTEGER I,INPUT DIMENSION UCAP(3),Z(3),U(3) C C THE COORDINATE SYSTEM IS CARTESIAN AND RIGHT-HANDED, C WITH Z3 DOWN, Z1 PARALLEL TO THE FAULT TRACE, AND Z2 C PERPENDICULAR. THE DISLOCATION SURFACE IS AT Z2=0, C DTOP.LE.Z3.LE.DBOT (NOTE THAT DTOP.GT.0.), AND C -L.LE.Z1.LE.L. C UCAP CONTAINS THE THREE COMPONENTS OF THE SLIP OF THE C POSITIVE-Z2 SIDE RELATIVE TO THE OTHER SIDE; UCAP(2) C SHOULD THEREFORE BE ZERO (BUT IT IS NOT USED HERE). C U CONTAINS THE THREE COMPONENTS OF THE DISPLACEMENT C AT THE OBSERVATION POINT Z. C DO 10 I=1,3 10 U(I)=0. DO 100 I=1,4 FACTOR=1. IF(I.EQ.2.OR.I.EQ.3) FACTOR= -1. ZETA1=L IF(I.EQ.3.OR.I.EQ.4) ZETA1= -L ZETA3=DBOT IF(I.EQ.2.OR.I.EQ.4) ZETA3= DTOP R=SQRT((Z(1)-ZETA1)**2+Z(2)**2+(Z(3)-ZETA3)**2) Q=SQRT((Z(1)-ZETA1)**2+Z(2)**2+(Z(3)+ZETA3)**2) TERM1=Z(2)*(Z(1)-ZETA1)*(2./(R*(R+Z(3)-ZETA3)) 1 -(5.*Q+8.*ZETA3)/(2.*Q*(Q+Z(3)+ZETA3)**2) 2 +4.*Z(3)*ZETA3*(2.*Q+Z(3)+ZETA3)/ 3 (Q**3*(Q+Z(3)+ZETA3)**2)) 4 +3.*ATANPV((Z(1)-ZETA1)*(Z(3)-ZETA3),(R*Z(2))) 5 -3.*ATANPV((Z(1)-ZETA1)*(Z(3)+ZETA3),(Q*Z(2))) C TERM2= -LOG(R+Z(3)-ZETA3) 1 +0.5*LOG(Q+Z(3)+ZETA3) 2 -(5.*Z(3)-3.*ZETA3)/(2.*(Q+Z(3)+ZETA3)) 3 -4.*Z(3)*ZETA3/(Q*(Q+Z(3)+ZETA3)) 4 +Z(2)**2*(2./(R*(R+Z(3)-ZETA3)) 5 -(5.*Q+8.*ZETA3)/(2.*Q*(Q+Z(3)+ZETA3)**2) 6 +4.*Z(3)*ZETA3*(2.*Q+Z(3)+ZETA3)/ 7 (Q**3*(Q+Z(3)+ZETA3)**2)) C C TERM3=Z(2)*(2./R-2./Q+4.*Z(3)*ZETA3/Q**3 C 1 +3./(Q+Z(3)+ZETA3)+2.*(Z(3)+3.*ZETA3)/ C 2 (Q*(Q+Z(3)+ZETA3))) C C TERM4=Z(2)*(2./R+4./Q-4.*Z(3)*ZETA3/Q**3 C 1 -6.*Z(3)/(Q*(Q+Z(3)+ZETA3))) C C TERM5= -LOG(R+Z(1)-ZETA1)+LOG(Q+Z(1)-ZETA1) C 1 +(6.*Z(3)**2+10.*Z(3)*ZETA3)/(Q*(Q+Z(1)-ZETA1)) C 2 +(6.*Z(3)*(Z(1)-ZETA1))/(Q*(Q+Z(3)+ZETA3)) C 3 +Z(2)**2*(2./(R*(R+Z(1)-ZETA1)) C 4 +4./(Q*(Q+Z(1)-ZETA1)) C 5 -4.*Z(3)*ZETA3*(2.*Q+Z(1)-ZETA1)/(Q**3*(Q+Z(1)-ZETA1)**2)) C C TERM6=Z(2)*(2.*(Z(3)-ZETA3)/(R*(R+Z(1)-ZETA1)) C 1 -2.*(Z(3)+2.*ZETA3)/(Q*(Q+Z(1)-ZETA1)) C 2 -4.*Z(3)*ZETA3*(Z(3)+ZETA3)*(2.*Q+Z(1)-ZETA1)/ C 3 (Q**3*(Q+Z(1)-ZETA1)**2)) C 4 +3.*ATANPV((Z(1)-ZETA1)*(Z(3)-ZETA3),(R*Z(2))) C 5 -3.*ATANPV((Z(1)-ZETA1)*(Z(3)+ZETA3),(Q*Z(2))) C U(1)=U(1)+FACTOR*TERM1*UCAP(1)/37.7 U(2)=U(2)+FACTOR*TERM2*UCAP(1)/37.7 C U(1)=U(1)+FACTOR*TERM4*UCAP(3)/37.7 C U(2)=U(2)+FACTOR*TERM5*UCAP(3)/37.7 C U(3)=U(3)+FACTOR*(TERM3*UCAP(1)+TERM6*UCAP(3))/37.7 100 CONTINUE RETURN END C C C C SUBROUTINE AURA (INPUT,U4,THETA4,X4,L4,DBOT4,DTOP4, + OUTPUT,V4) C C COMPUTES DISPLACEMENTS IN AN INFINITE UNIFORM ELASTIC HALFSPACE C OF POISSON'S RATIO 0.25 CAUSED BY UNIFORM SLIP OVER THE SURFACE C OF A BURIED RECTANGULAR DISLOCATION WITH HORIZONTAL AND PLUNGING C SIDES. TO CALCULATE DISPLACEMENTS CAUSED BY MORE GENERAL SLIP, C USE THIS ROUTINE MANY TIMES IN THE KERNEL OF AN INTEGRAL, C ASSIGNING THE LOCAL AVERAGE SLIP OVER EACH VERY SMALL RECTANGLE. C C EQUATIONS BASED ON MANSINHA AND SMYLIE (1971), C BULL.SEISM.SOC.AMER., 61, 1433-1440. C IMPLICIT DOUBLE PRECISION (A-Z) REAL OUTPUT INTEGER I,INPUT DIMENSION U(3),U4(3),V(3),V4(3),X(3),X4(3) DO 1 I=1,3 U(I)=U4(I) X(I)=X4(I) 1 CONTINUE THETA=THETA4 L=L4 DBOT=DBOT4 DTOP=DTOP4 C C U(3) IS THE DISLOCATION VECTOR(I.E. THE TOTAL SLIP VECTOR C OF THE HANGING WALL RELATIVE TO THE FOOTWALL), C IN A CARTESIAN COORDINATE SYSTEM WHERE X1 IS PARALLEL TO THE C FAULT STRIKE, X2 IS PERPENDICULAR TO IT AND POINTING TO HANGING C WALL, AND Z IS DOWN (TAKE CARE TO MAKE THESE RIGHT-HANDED). C THETA IS THE FAULT DIP ANGLE (IN RADIANS) FROM X2 AXIS. C X(3) GIVES THE OBSERVATION POINT IN THESE COORDINATES. C L IS THE HALF-LENGTH OF THE RECTANGLE, ALONG STRIKE. C DBOT AND DTOP ARE OBLIQUE DISTANCES IN THE FAULT PLANE (AND C DOWN THE DIP) FROM THE SURFACE TO THE BOTTOM AND TOP OF THE C RECTANGULAR DISLOCATION, RESPECTIVELY. C NOTE: FAULT DOES NOT BREAK THE SURFACE, SO DTOP > 0., C ALTHOUGH IT MAY BE SMALL. C V(3) IS THE ANSWER: THE 3-COMPONENT DISPLACEMENT AT THE C OBSERVATION POINT IN SAME UNITS AS U AND IN X1,X2,X3 COORDINATES. C SINTH=SIN(THETA) COSTH=COS(THETA) TANTH=SINTH/COSTH C NOTE: FAULT MAY NOT BE VERTICAL IN THIS ROUTINE. SECTH=1./COSTH R2=X(2)*SINTH-X(3)*COSTH R3=X(2)*COSTH+X(3)*SINTH Q2=X(2)*SINTH+X(3)*COSTH Q3= -X(2)*COSTH+X(3)*SINTH DO 10 I=1,3 10 V(I)=0.D0 DO 100 I=1,4 FACTOR = 1.00D0 IF (I.EQ.2.OR.I.EQ.3) FACTOR = -1.00D0 XI1=L IF(I.GE.3) XI1 = -L XI=DBOT IF(I.EQ.2.OR.I.EQ.4) XI = DTOP XI2=XI*COSTH XI3=XI*SINTH RSQUAR=(X(1)-XI1)**2+(X(2)-XI2)**2+(X(3)-XI3)**2 QSQUAR=(X(1)-XI1)**2+(X(2)-XI2)**2+(X(3)+XI3)**2 K=SQRT(QSQUAR-(Q3+XI)**2) H=SQRT(QSQUAR-(X(1)-XI1)**2) UCAP1=U(1) UCAP=SQRT(U(2)**2+U(3)**2) IF(U(2).LT.0.)UCAP= -UCAP Q=SQRT(QSQUAR) R=SQRT(RSQUAR) Q=MAX(Q,1.00001D0*ABS(X(1)-XI1)) R=MAX(R,1.00001D0*ABS(X(1)-XI1)) TERM1=(X(1)-XI1)*(2.*R2/(R*(R+R3-XI)) 1 -(4.*Q2-2.*X(3)*COSTH)/(Q*(Q+Q3+XI)) 2 -3.*TANTH/(Q+X(3)+XI3) 3 +4.*Q2*X(3)*SINTH/Q**3 4 -4.*Q2*Q3*X(3)*SINTH*(2.*Q+Q3+XI)/ 5 (Q**3*(Q+Q3+XI)**2)) 6 -6.*TANTH**2*ATANPV(((K-Q2*COSTH)*(Q-K)+ 7 (Q3+XI)*K*SINTH),((X(1)-XI1)*(Q3+XI)*COSTH)) 8 +3.*ATANPV((X(1)-XI1)*(R3-XI),(R2*R)) 9 -3.*ATANPV((X(1)-XI1)*(Q3+XI),(Q2*Q)) C TERM2=SINTH*(3.*TANTH*SECTH*LOG(Q+X(3)+XI3)- 1 LOG(R+R3-XI)-(1.+3.*TANTH**2)*LOG 2 (Q+Q3+XI))+2.*R2**2*SINTH/(R*(R+R3-XI)) 3 +2.*R2*COSTH/R-2.*SINTH* 4 (2.*X(3)*(Q2*COSTH-Q3*SINTH)+Q2*(Q2+X(2)*SINTH)) 5 /(Q*(Q+Q3+XI))-3.*TANTH*(X(2)-XI2)/ 6 (Q+X(3)+XI3)+2.*(Q2*COSTH-Q3*SINTH- 7 X(3)*SINTH**2)/Q+4.*Q2*X(3)*SINTH* 8 ((X(2)-XI2)+Q3*COSTH)/Q**3-4.*Q2**2*Q3*X(3) 9 *SINTH**2*(2.*Q+Q3+XI)/(Q**3*(Q+Q3+XI)**2) C C TERM3=COSTH*(LOG(R+R3-XI)+(1.+3.*TANTH**2)*LOG( C 1 Q+Q3+XI)-3.*TANTH*SECTH*LOG(Q+X(3)+XI3)) C 2 +2.*R2*SINTH/R+2.*SINTH*(Q2+X(2)*SINTH)/Q C 3 -2.*R2**2*COSTH/(R*(R+R3-XI))+ C 4 (4.*Q2*X(3)*SINTH**2-2.*(Q2+X(2)*SINTH)*(X(3)+Q3*SINTH))/ C 5 (Q*(Q+Q3+XI))+ C 6 4.*Q2*X(3)*SINTH*(X(3)+XI3-Q3*SINTH)/Q**3 C 7 -4.*Q2**2*Q3*X(3)*COSTH*SINTH*(2.*Q+Q3+XI)/ C 8 (Q**3*(Q+Q3+XI)**2) C TERM4=(X(2)-XI2)*SINTH*(2./R+4./Q-4.*XI3*X(3)/Q**3 1 -3./(Q+X(3)+XI3))-COSTH*(3.*LOG(Q+X(3)+XI3) 2 +2.*(X(3)-XI3)/R+4.*(X(3)-XI3)/Q+4.*XI3*X(3)*(X(3)+XI3) 3 /Q**3)+3.*(LOG(Q+X(3)+XI3)-SINTH*LOG(Q+ 4 Q3+XI))/COSTH+6.*X(3)*(COSTH/Q-Q2*SINTH/ 5 (Q*(Q+Q3+XI))) C TERM5=SINTH*(-LOG(R+X(1)-XI1)+LOG(Q+X(1)-XI1) 1 +4.*XI3*X(3)/(Q*(Q+X(1)-XI1))+3.*(X(1)-XI1)/ 2 (Q+X(3)+XI3)+(X(2)-XI2)**2*(2./(R*(R+X(1)-XI1)) 3 +4./(Q*(Q+X(1)-XI1))-4.*XI3*X(3)*(2.*Q+X(1)-XI1)/ 4 (Q**3*(Q+X(1)-XI1)**2))) 5 -COSTH*((X(2)-XI2)*(2.*(X(3)-XI3)/(R*(R+X(1)-XI1)) 6 +4.*(X(3)-XI3)/(Q*(Q+X(1)-XI1))+4.*XI3*X(3)* 7 (X(3)+XI3)*(2.*Q+X(1)-XI1)/(Q**3*(Q+X(1)-XI1)**2)) 8 +6.*ATANPV((X(1)-XI1)*(X(2)-XI2),((H+X(3)+XI3)*(Q+H))) 9 -3.*ATANPV((X(1)-XI1)*(R3-XI),(R2*R))+ A 6.*ATANPV((X(1)-XI1)*(Q3+XI),(Q2*Q)))+ B 6.*(SECTH*ATANPV(((K-Q2*COSTH)*(Q-K)+(Q3+XI)*K*SINTH) C ,((X(1)-XI1)*(Q3+XI)*COSTH)) D +X(3)*(((SINTH**2-COSTH**2)*(Q3+XI)+2.*Q2* E COSTH*SINTH)/(Q*(Q+X(1)-XI1))+(X(1)-XI1)*SINTH**2 F /(Q*(Q+Q3+XI)))) C C TERM6=SINTH*((X(2)-XI2)*(2.*(X(3)-XI3)/(R*(R+X(1)-XI1)) C 1 +4.*(X(3)-XI3)/(Q*(Q+X(1)-XI1))-4.*XI3*X(3)* C 2 (X(3)+XI3)*(2.*Q+X(1)-XI1)/(Q**3*(Q+X(1)-XI1)**2)) C 3 -6.*ATANPV((X(1)-XI1)*(X(2)-XI2),((H+X(3)+XI3)*(Q+H))) C 4 +3.*ATANPV((X(1)-XI1)*(R3-XI),(R2*R))- C 5 6.*ATANPV((X(1)-XI1)*(Q3+XI),(Q2*Q))) C 6 +COSTH*(LOG(R+X(1)-XI1)-LOG(Q+X(1)-XI1)- C 7 2.*(X(3)-XI3)**2/(R*(R+X(1)-XI1)) C 8 -4.*((X(3)+XI3)**2-XI3*X(3))/(Q*(Q+X(1)-XI1)) C 9 -4.*XI3*X(3)*(X(3)+XI3)**2*(2.*Q+X(1)-XI1)/ C A (Q**3*(Q+X(1)-XI1)**2)) C B +6.*X(3)*(COSTH*SINTH*(2.*(Q3+XI)/(Q*(Q+X(1)-XI1)) C C +(X(1)-XI1)/(Q*(Q+Q3+XI))) C D -Q2*(SINTH**2-COSTH**2)/(Q*(Q+X(1)-XI1))) C DELD=DBOT-DTOP V(1)=V(1)+FACTOR*(TERM1*UCAP1+TERM4*UCAP)/37.69911184D0 V(2)=V(2)+FACTOR*(TERM2*UCAP1+TERM5*UCAP)/37.69911184D0 C V(3)=V(3)+FACTOR*(TERM3*UCAP1+TERM6*UCAP)/37.69911184D0 100 CONTINUE V4(1)=V(1) V4(2)=V(2) C V4(3)=V(3) RETURN END C C C C DOUBLE PRECISION FUNCTION ATANPV (Y,X) C C RETURNS PRINCIPAL VALUE OF INVERSE TANGENT OF (Y,X) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) IF (X.GE.0.) THEN XX=X YY=Y ELSE XX= -X YY= -Y ENDIF IF(YY.EQ.0.0.OR.XX.EQ.0.0) THEN IF(YY.EQ.0.0) ATANPV = 0.0 IF(XX.EQ.0.0) ATANPV = 1.570796327D0 ELSE ATANPV = ATAN2(YY,XX) ENDIF RETURN END C C C SUBROUTINE GEOSCO(INPUT,ACREEP,ALPHAT,AREA,BCREEP,BIOT, + BYERLY,CCREEP,CFRIC,CONDUC, + DCREEP,DQDTDA,DETJ,DXS,DYS, + ECREEP,FAZ,FDIP,FFRIC,FTAN, + BENS1,BENS2,BENCHX,BENCHY, + BERR,BRATE,GMEAN,INDFST,INDLST,IUNITT, + MXEL,MXFEL,MXNODE,IBENCH,NFL, + NBENCH,NBASE, + NODEF,NODES,NUMEL, + OFFMAX,OFFSET,OKDELV,ONEKM, + RADIO,REFSTR,RHOH2O,RHOBAR,SITE1,SITE2, + TITLE1,TITLE2,TITLE3, + TLNODE,TRHMAX,TSURF,UNITL,UNITT, + UX,UY,V,VG,VGCRE,WEDGE, + XNODE,YNODE,ZMNODE,ZMOHO, + WORK,NFTRE, + OUTPUT,FC,FIMUDZ,FPEAKS,FSLIPS,FTSTAR,SMSCOR, + TLINT,ZTRANF) CHARACTER*31 SITE1,SITE2 CHARACTER*80 TITLE1,TITLE2,TITLE3 DOUBLE PRECISION V LOGICAL FSLIPS DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),RADIO (2),RHOBAR(2) DIMENSION AREA(MXEL),DETJ(7,MXEL),DQDTDA(MXNODE), + DXS(6,7,MXEL),DYS(6,7,MXEL), + FAZ(2,MXFEL), FC(2,2,7,MXFEL), FDIP(3,MXFEL), + FIMUDZ(7,MXFEL), FPEAKS(MXFEL), FSLIPS(MXFEL), + FTAN(7,MXFEL), FTSTAR(2,7,MXFEL), + BENS1(NBENCH), BENS2(NBENCH), + BENCHX(NBENCH),BENCHY(NBENCH),BERR(NBASE),BRATE(NBASE), + INDFST(NBASE),INDLST(NBASE),IBENCH(NBENCH), + NFTRE(MXNODE), + NODEF(6,MXFEL),NODES(6,MXEL),OFFSET(MXFEL), + SITE1(NBASE),SITE2(NBASE), + TLINT(7,MXEL), + TLNODE(MXNODE),UX(MXNODE),UY(MXNODE), + V(2,MXNODE),VG(2,NBENCH),VGCRE(2,NBENCH), + XNODE(MXNODE),YNODE(MXNODE), + ZMNODE(MXNODE),ZMOHO(7,MXEL),ZTRANF(7,MXFEL) C C COMPUTE LONG-TERM-AVERAGE BENCHMARK VELOCITIES FROM FE MODEL: C CALL INTEPG (INPUT,BENS1,BENS2,MXEL,MXNODE,IBENCH,NBENCH,NODES, + V, + OUTPUT,VG) C C THICKNESS OF CRUST AND MANTLE LITHOSPHERE: C CALL INTERP (INPUT,ZMNODE,MXEL,MXNODE,NODES,NUMEL, + OUTPUT,ZMOHO) CALL INTERP (INPUT,TLNODE,MXEL,MXNODE,NODES,NUMEL, + OUTPUT,TLINT) C C COMPUTE TACTICAL VALUES OF LIMITS ON VISCOSITY, AND WEIGHTS FOR C IMPOSITION OF CONSTRAINTS IN LINEAR SYSTEMS: C CALL LIMITS (INPUT,AREA,DETJ,IUNITT,MXEL, + NUMEL,OKDELV,REFSTR,TLINT,ZMOHO, + OUTPUT,CONSTR,ETAMAX,FMUMAX,VISMAX) C C LOCATE THE BRITTLE/DUCTILE TRANSITION IN EACH FAULT C DO 100 M=1,7 DO 90 I=1,NFL ZTRANF(M,I)=MAX(ZMNODE(NODEF(2,I))/2.,ONEKM) 90 CONTINUE 100 CONTINUE DO 200 IMOHR=1,3 CALL MOHR (INPUT,ACREEP,ALPHAT,BCREEP,BIOT,BYERLY, + CCREEP,CFRIC,CONDUC,CONSTR,DCREEP,DQDTDA, + ECREEP,FDIP,FFRIC,FMUMAX,FTAN,GMEAN,MXFEL, + MXNODE,NFL,NODEF, + OFFMAX,OFFSET, + ONEKM,RADIO,RHOH2O, + RHOBAR,TLNODE,TSURF,V,WEDGE, + ZMNODE, + MODIFY,ZTRANF, + OUTPUT,FC,FIMUDZ,FPEAKS,FSLIPS,FTSTAR) 200 CONTINUE C C COMPUTE BENCHMARK VELOCITIES DUE TO ELASTIC DISLOCATIONS C ACCUMULATING IN THE BRITTLE PART OF EACH FAULT AT THE RATE PREDICTED C BY THE FINITE ELEMENT MODEL. C CALL GDTRAN(INPUT,FDIP,BENCHX,BENCHY,MXNODE,NBENCH, + NFL,NODEF,V,WEDGE, + XNODE,YNODE,ZTRANF, + WORK,NFTRE, + OUTPUT,UX,UY) C C****** ALASKA (BEGIN SECTION) ****************** C C AUGMENT (UX,UY) VELOCITIES AT BENCHMARKS C (DUE TO POSITIVE ELASTIC DISLOCATIONS ACCUMULATING C IN THE ALEUTIAN TRENCH AT THE RATE IMPLIED BY C THE RELATIVE MOTION BETWEEN THE PACIFIC PLATE AND C THE NORTH AMERICAN PLATE: C CALL BIGDIS (INPUT,AREA,BENCHX,BENCHY, + DXS,DYS, + MXNODE,MXEL, + NBENCH,NODES,NUMEL, + TLNODE,TRHMAX,V, + XNODE,YNODE,ZMNODE, + MODIFY,UX,UY) C C****** ALASKA (END SECTION) ****************** C C COMBINE THE TWO COMPONENTS OF BENCHMARK VELOCITIES, ASSUMING AN C ASEISMIC YEAR IN WHICH THE BRITTLE PARTS OF FAULTS DONT SLIP: C DO 300 I=1,NBENCH VGCRE(1,I)=VG(1,I)-UX(I) VGCRE(2,I)=VG(2,I)-UY(I) 300 CONTINUE C C PRINT TABLE OF OBSERVED AND PREDICTED BASELINE RATES AND ERR0RS: C WRITE(IUNITT,350) 350 FORMAT (' --------------------------------' +,'---------------------------------------' +/' ' +/' GEODETIC BASELINE RATE PREDICTIO' +,'NS AND THEIR ERR0RS:' +/' ' +/' ' +,' OBSERVED PREDICTED MODEL MODEL UNCORRECTED' +/' ' +,' RATE RATE ERR0R SIGMA ERR0R MODEL RATE' +/' BASELINE ' +,' (MM/YR) (MM/YR) (MM/YR) (MM/YR) (SIGMAS) (MM/YR)' +/' ------------------------------- ' +,' -------- -------- ------- ------- ------- --------') 351 FORMAT (' ',A31,'-TO-'/ + ' ',A31,F12.3,F10.3,F9.3,F9.3,2X,F7.2,F12.3) DRSGMA=0.0 REVSG2=0.0 REVSG1=0.0 DRSGM2=0.0 DO 500 I=1,NBASE NEL1=INDFST(I) NEL2=INDLST(I) DX=BENCHX(NEL1)-BENCHX(NEL2) DY=BENCHY(NEL1)-BENCHY(NEL2) DU=VGCRE(1,NEL1)-VGCRE(1,NEL2) DV=VGCRE(2,NEL1)-VGCRE(2,NEL2) DS1=VG(1,NEL1)-VG(1,NEL2) DS2=VG(2,NEL1)-VG(2,NEL2) DXY2=SQRT(DX**2+DY**2) UVSSL=(DS1*DX+DS2*DY)/DXY2 UVLENG=(DU*DX+DV*DY)/DXY2 DRSIG=ABS(UVLENG-BRATE(I))/BERR(I) ERRSI=ABS(UVLENG-BRATE(I)) DRSGMA=DRSGMA+DRSIG DRSGM2=DRSGM2+DRSIG**2 REVSG2=REVSG2+1/BERR(I)**2 REVSG1=REVSG1+1/BERR(I) GRRA=BRATE(I)*UNITL/UNITT UVLN=UVLENG*UNITL/UNITT UVSSL=UVSSL*UNITL/UNITT ERRMPY=ERRSI*UNITL/UNITT BERMPY=BERR(I)*UNITL/UNITT WRITE(IUNITT,351) SITE1(I),SITE2(I),GRRA,UVLN,ERRMPY, + BERMPY,DRSIG,UVSSL 500 CONTINUE AVSCOR=DRSGMA/REVSG1*UNITL/UNITT SMSCOR=SQRT(DRSGM2/REVSG2)*UNITL/UNITT WRITE (IUNITT,599) TITLE1,TITLE2,TITLE3 599 FORMAT (/ / + ' ',A80/' ',A80/' ',A80/) WRITE(IUNITT,600) AVSCOR,SMSCOR 600 FORMAT(/// + ' ',' WEIGHTED MEAN PREDICTION ERR0R RELATIVE TO'/, + ' ',' GEODETIC BASELINE RATES IS ',F8.3,' (MM/YR).'/, + ' ',' WEIGHTED RMS ERR0R IS ',F8.3,' (MM/YR).') RETURN END C C C SUBROUTINE GETNET (INPUT,IUNIT7,IUNIT8, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP, + NFAKEN,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,OFFMAX,OFFSET, + TITLE1,TLNODE,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 CHARACTER*80 TITLE1 LOGICAL ALLOK,BRIEF C C NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN C C EXTERNAL ARRAYS: DIMENSION CHECKE(MXEL),CHECKF(MXFEL),CHECKN(MXNODE), + DQDTDA(MXNODE),ELEV(MXNODE), + FAZ(2,MXFEL),FDIP(3,MXFEL), + NODEF(6,MXFEL),NODES(6,MXEL),OFFSET(MXFEL), + TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) C INTERNAL ARRAYS: DIMENSION DIPS(3),IFN(6),VECTOR(7) 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 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 (OPTION "BRIEF" SUPPRESSES MOST OUTPUT) C READ (IUNIT7,*) NUMNOD,NREALN,NFAKEN,N1000,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 N2=2*NUMNOD 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 (/' ', + ' MANTLE'/ + ' ', + ' CRUSTAL LITHOSPHERE'/ + ' NODE X Y ELEVATION ', + 'HEAT-FLOW THICKNESS THICKNESS'/) 55 FORMAT (' ',I10,1P,2E11.3,4E10.2) ENDIF DO 90 K=1,NUMNOD CHECKN(K)=.FALSE. 90 CONTINUE DO 100 K=1,NUMNOD CALL READN (INPUT,IUNIT7,IUNIT8,7, + OUTPUT,VECTOR) INDEX=VECTOR(1)+0.5 XI=VECTOR(2) YI=VECTOR(3) ELEVI=VECTOR(4) QI=VECTOR(5) ZMI=VECTOR(6) TLI=VECTOR(7) 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,96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') STOP ENDIF XNODE(I)=XI YNODE(I)=YI IF (ZMI.LT.0.) THEN WRITE (IUNIT8,97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') STOP ENDIF ZMNODE(I)=ZMI IF (TLI.LT.0.) THEN WRITE (IUNIT8,98) 98 FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS IS', + ' NON-PHYSICAL.') STOP ENDIF TLNODE(I)=TLI IF (.NOT.BRIEF) THEN WRITE (IUNIT8,55) INDEX,XI,YI,ELEVI,QI,ZMI,TLI 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 OFFMAX=0. 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.)'/ + ' OFFSET IS THE TOTAL PAST SLIP OF THE FAULT.'/ / + ' ELEMENT N1 N2 N3 N4 N5 N6 DIP1 DIP2 DIP3', + ' ARG1 ARG3 OFFSET'/) 240 FORMAT (' ',I6,':',6I5,1X,3F6.1,1X,2F5.0,F9.0) DO 300 K=1,NFL OFF=0. READ (IUNIT7,*,ERR=242) I,(IFN(J),J=1,6),(DIPS(L),L=1,3), + AZ1,AZ3,OFF 242 CHECKF(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,240) I,(IFN(J),J=1,6), + (DIPS(L),L=1,3),AZ1,AZ3,OFF 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 IF (OFF.LT.0.) THEN WRITE (IUNIT8,280) OFF 280 FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2, + ' FOR FAULT ELEMENT',I6/ + ' OFFSETS MAY NOT BE NEGATIVE.') STOP ENDIF OFFSET(I)=OFF OFFMAX=MAX(OFFMAX,OFF) 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 ELSE IF (OFFMAX.GT.0.) THEN WRITE (IUNIT8,400) OFFMAX 400 FORMAT (/' GREATEST FAULT OFFSET READ WAS ',1P,E10.2) ELSE WRITE (IUNIT8,401) 401 FORMAT (/' SINCE FAULT OFFSETS ARE ALL ZERO,', + ' INPUT PARAMETER BYERLY WILL HAVE NO EFFECT.') ENDIF ENDIF IF (.NOT. BRIEF) WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C C C SUBROUTINE INTERP (INPUT,FATNOD,MXEL,MXNODE,NODES,NUMEL, + OUTPUT,FATIP) C C INTERPOLATE A SCALAR FUNCTION KNOWN AT THE NODES ("FATNOD") C TO VALUES AT THE 7 INTEGRATION POINTS IN EACH TRIANGULAR C CONTINUUM ELEMENT. C DOUBLE PRECISION PHI COMMON /PHITAB/ PHI DIMENSION PHI(6,7) DIMENSION FATNOD(MXNODE),FATIP(7,MXEL),NODES(6,MXEL) C DO 100 M=1,7 DO 90 I=1,NUMEL FATIP(M,I)=PHI(1,M)*FATNOD(NODES(1,I))+ + PHI(2,M)*FATNOD(NODES(2,I))+ + PHI(3,M)*FATNOD(NODES(3,I))+ + PHI(4,M)*FATNOD(NODES(4,I))+ + PHI(5,M)*FATNOD(NODES(5,I))+ + PHI(6,M)*FATNOD(NODES(6,I)) 90 CONTINUE 100 CONTINUE RETURN END C C C SUBROUTINE INTEPG (INPUT,BENS1,BENS2,MXEL,MXNODE,IBENCH,NBENCH, + NODES,V, + OUTPUT,VG) C C INTERPOLATE A VECTOR FUNCTION KNOWN AT THE NODES ("V") C TO VALUES AT BENCHMARKS. C DOUBLE PRECISION V DIMENSION BENS1(NBENCH),BENS2(NBENCH),IBENCH(NBENCH) DIMENSION NODES(6,MXEL),VG(2,NBENCH),V(2,MXNODE) C 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 DO 90 I=1,NBENCH T1=BENS1(I) T2=BENS2(I) T3=1.00-T1-T2 IE=IBENCH(I) U1=V(1,NODES(1,IE)) U2=V(1,NODES(2,IE)) U3=V(1,NODES(3,IE)) U4=V(1,NODES(4,IE)) U5=V(1,NODES(5,IE)) U6=V(1,NODES(6,IE)) V1=V(2,NODES(1,IE)) V2=V(2,NODES(2,IE)) V3=V(2,NODES(3,IE)) V4=V(2,NODES(4,IE)) V5=V(2,NODES(5,IE)) V6=V(2,NODES(6,IE)) VG(1,I)=PHIVAL(T1,T2,T3,U1,U2,U3,U4,U5,U6) VG(2,I)=PHIVAL(T1,T2,T3,V1,V2,V3,V4,V5,V6) 90 CONTINUE RETURN END C C C SUBROUTINE LIKELY (INPUT,VMIN,VMAX,V, + OUTPUT,PROB) C C COMPUTE THE PROBABILITY THAT VELOCITY V IS CONSISTENT WITH LIMITS C VMIN AND VMAX (BOTH OF WHICH CONTAIN RANDOM ERR0R). C LOGICAL BOXCAR,POST,POINT DATA BADATA/0.20/,SDMODV/0.10/,SDAGE/0.15/,E/2.71828/, + TPIM1H/0.398942/ C PROB=1. IF ((VMIN.EQ.0.) .AND. (VMAX.EQ.0.)) RETURN V1=VMIN IF (VMIN.EQ.0.) V1= -1.E50 V2=VMAX IF (VMAX.EQ.0.) V2= +1.E50 VCENT=(V1+V2)/2. SPREAD=ABS((V2-V1)/(V1+V2)) BOXCAR=(SPREAD.GE.1.3*SDAGE).OR.(V1.EQ.-1.E50).OR.(V2.EQ.1.E50) POINT=SPREAD.EQ.0. POST=(.NOT.BOXCAR).AND.(.NOT.POINT) VU=V*(1.+2.*SDMODV) VL=V*(1.-2.*SDMODV) IF(V.GE.0.)GO TO 4 TEMP=VU VU=VL VL=TEMP 4 V1U=V1*(1.+2.*SDAGE) V1L=V1*(1.-2.*SDAGE) IF(V1.GE.0.)GO TO 6 TEMP=V1U V1U=V1L V1L=TEMP 6 V2U=V2*(1.+2.*SDAGE) V2L=V2*(1.-2.*SDAGE) IF(V2.GE.0.)GO TO 8 TEMP=V2U V2U=V2L V2L=TEMP 8 IF(VU.GT.V1L)GO TO 10 PROB=BADATA RETURN 10 IF(VL.LT.V2U) GO TO 20 PROB=BADATA RETURN 20 IF((VL.LT.V1U).OR.(VU.GT.V2L)) GO TO 30 PROB=1. RETURN 30 PROB=0. B= -3.1 DB=0.10 D1=ABS(V1)*SDAGE*1.41421 D2=ABS(V2)*SDAGE*1.41421 DW=DB*SDMODV*ABS(V) W=V-31.*DW DO 40 I=1,61 B=B+DB PV=TPIM1H*E**(-0.5*B**2) W=W+DW IF(POINT) PGF=E**(-MIN(99.,0.5*((W-V1)/(SDAGE*V1))**2)) IF(POST)PGF=E**(-MIN(99.,(0.832*(W-VCENT)/(SPREAD*VCENT))**2)) IF(BOXCAR)PGF=MIN(0.5*(1.+ERF((W-V1)/D1)), 1 1.,0.5*(1.+ERF((V2-W)/D2))) PGF=AMAX1(PGF,1.E-50) PG=BADATA+(1.-BADATA)*PGF 40 PROB=PROB+PV*PG*DB RETURN END C C C SUBROUTINE LIMITS (INPUT,AREA,DETJ,IUNITT,MXEL, + NUMEL,OKDELV,REFSTR,TLINT,ZMOHO, + OUTPUT,CONSTR,ETAMAX,FMUMAX,VISMAX) C C COMPUTE AREA, MEAN THICKNESS, AND OTHER DIMENSIONAL PARAMETERS C OF THE PLATE, THEN DETERMINE VALUES OF STIFFNESS LIMITS NEEDED C TO KEEP VELOCITY ERR0RS DOWN TO ORDER "OKDELV" AT SHEAR STRESS C LEVEL "REFSTR". C PARAMETER (PART=0.1) C "PART" IS A DIMENSIONLESS NUMBER, BETWEEN 0.0 AND 1.0, WHICH C EXPRESSES WHAT FRACTION OF THE BOTTOM BOUNDARY OF THE MODEL IS C LIKELY TO BE EXPOSED TO VELOCITY BOUNDARY CONDITIONS (E.G., THE C TOP OF A SUBDUCTING SLAB). PRECISION IS NOT NEEDED. HOWEVER, C IF "PART" IS 1.0 WHILE THE ACTUAL AREA OF BASAL BOUNDARY CONDITIONS C IS MUCH LESS, THEN "ETAMAX" WILL BE SET TOO LOW, AND MAY C ARTIFICIALLY LIMIT THE BASAL TRACTION THAT YOU CAN APPLY. C CONVERSELY, IF BOTTOM VELOCITY BOUNDARY CONDITIONS CONVER THE C WHOLE AREA, BUT "PART" IS SET TO 0.1, THEN "ETAMAX" WILL BE TOO C HIGH, AND THERE IS A POSSIBILITY THAT THE SOLUTION WILL JUST C OSCILLATE WITHOUT CONVERGING. C DOUBLE PRECISION WEIGHT COMMON /WGTVEC/ WEIGHT DIMENSION WEIGHT(7) DIMENSION AREA(MXEL),DETJ(7,MXEL),TLINT(7,MXEL),ZMOHO(7,MXEL) C C DATA ITEM "NFAULT" GIVES THE TYPICAL NUMBER OF FAULTS WHICH ARE C CROSSED BY ANY STRAIGHT LINE RUNNING ACROSS THE MODEL. IT DOES C NOT NEED TO BE VERY ACCURATE] DATA NFAULT /5/ C TOTALA=0. TOTALV=0. DO 20 M=1,7 DO 10 I=1,NUMEL DA=AREA(I)*DETJ(M,I)*WEIGHT(M) TOTALA=TOTALA+DA TOTALV=TOTALV+DA*(ZMOHO(M,I)+TLINT(M,I)) 10 CONTINUE 20 CONTINUE THICK=TOTALV/TOTALA SIDE=SQRT(TOTALA) CONSTR=NFAULT*REFSTR*THICK/OKDELV ETAMAX=0.5*REFSTR*THICK/(PART*SIDE*OKDELV) FMUMAX=NFAULT*REFSTR/OKDELV VISMAX=REFSTR*SIDE/OKDELV WRITE (IUNITT,50) TOTALA,TOTALV,THICK,SIDE,CONSTR,ETAMAX, + FMUMAX,VISMAX 50 FORMAT (/ /' SUBPROGRAM -LIMITS- PERFORMS DIMENSIONAL ANALYSIS'/ + ' AND ESTIMATES NECESSARY STIFFNESS LIMITS TO BALANCE'/1P, + ' THE CONFLICTING OBJECTIVES OF ACCURACY AND PRECISION:'/ / + ' AREA OF MODEL = ',E10.3,' LENGTH**2'/ + ' VOLUME OF CRUST = ',E10.3,' LENGTH**3'/ + ' TYPICAL THICKNESS = ',E10.3,' LENGTH'/ + ' TYPICAL WIDTH = ',E10.3,' LENGTH'/ + ' CONSTR (CONSTRAINT WEIGHT) = ',E10.3,' FORCE-SEC/LENGTH**2'/ + ' ETAMAX (MAX. BASAL COUPLING) = ',E10.3,' FORCE-SEC/LENGTH**3'/ + ' FMUMAX (MAX. FAULT STIFFNESS) = ',E10.3,' FORCE-SEC/LENGTH**3'/ + ' VISMAX (MAX. BLOCK VISCOSITY) = ',E10.3,' FORCE-SEC/LENGTH**2') RETURN END C C C SUBROUTINE MOHR (INPUT,ACREEP,ALPHAT,BCREEP,BIOT,BYERLY, + CCREEP,CFRIC,CONDUC,CONSTR,DCREEP,DQDTDA, + ECREEP,FDIP,FFRIC,FMUMAX,FTAN,GMEAN, + MXFEL,MXNODE,NFL,NODEF, + OFFMAX,OFFSET, + ONEKM,RADIO,RHOH2O, + RHOBAR,TLNODE,TSURF,V,WEDGE, + ZMNODE, + MODIFY,ZTRANF, + OUTPUT,FC,FIMUDZ,FPEAKS,FSLIPS,FTSTAR) C C THIS SUBPROGRAM CONTAINS THE NONLINEAR RHEOLOGY OF THE FAULTS. C FOR EACH OF 7 INTEGRATION POINTS ALONG THE LENGTH OF EACH FAULT C ELEMENT, IT: C C (1) COMPUTES THE SLIP-RATE VECTOR ON THE FAULT SURFACE, C (2) DETERMINES THE SHEAR STRESS ON THE FAULT SURFACE BY MOHR/ C COULOMB/NAVIER THEORY (THIS STRESS IS PROPORTIONAL TO DEPTH, C SO THE CALCULATION IS ACTUALLY DONE AT UNIT DEPTH AND THEN C SCALED), C (3) PROCEEDS DOWN THE DIP OF THE FAULT, CHECKING TEMPERATURE, C STRAIN-RATE, AND PRESSURE TO SEE IF FRICTIONAL OR CREEP C SHEAR STRESS IS LOWER, C (4) REPORTS THE VERTICAL INTEGRAL OF "MU" (THE RATIO OF SHEAR C STRESS TO SLIP RATE) DOWN THE FAULT AS "FIMUDZ". C (NOTE THAT THE INTEGRAL IS VERTICAL, NOT ON A SLANT, EVEN THOUGH C CONDITIONS ARE EVALUATED ALONG A SLANT PATH.) C (5) FOR DIPPING, OBLIQUE-SLIP FAULT ONLY, ALSO REPORTS RECOMMENDED C TACTICAL VALUES FOR THE MATRIX "FC" AND THE VECTOR "FTSTAR" C WHICH JOINTLY DESCRIBE A LINEARIZED RHEOLOGY STIFFER THAN C THE ACTUAL NONLINEAR RHEOLOGY. C (6) "ZTRANF" IS THE LATEST ESTIMATE OF THE DEPTH C TO THE BRITTLE/DUCTILE TRANSITION, AT THE FAULT MIDPOINT. C (7) LOGICAL VARIABLE "FSLIPS" INDICATES WHETHER THE FAULT IS C SLIPPING AT ITS MIDPOINT. OTHERWISE, IT IS IN THE ARTIFICIAL C LINEARIZED REGIME, WITH STIFFNESS "FMUMAX". C (8) "FPEAKS" GIVES THE PEAK SHEAR STRESS AT THE MIDPOINT OF EACH C FAULT, EVALUATED AT THE BRITTLE/DUCTILE TRANSITION. C C NOTE THAT PORE PRESSURES ARE CONSIDERED IN THE CALCULATION OF C FRICTIONAL STRENGTH: C *NORMAL PORE PRESSURES REDUCE THE EFFECTIVE NORMAL STRESS ON THE C FAULT SURFACE BY THE AMOUNT C -BIOT*GMEAN*RHOH20*Z C *IF (OFFMAX.GT.0.), THEN THE REMAINING EFFECTIVE FRICTIONAL STRENGTH C OF THE FAULT IS MULTIPLIED BY THE REDUCING FACTOR C *(1.-BYERLY*OFFSET(I)/OFFMAX). C THIS IS ALSO A PORE PRESSURE EFFECT, BECAUSE BYERLY'S MODEL IS C THAT GOUGE LAYERS HAVE THICKNESS IN PROPORTION TO OFFSET, AND C THAT THEY SUPPORT NON-DARCY STATIC PORE PRESSURE GRADIENTS WHICH C REDUCE THE EFFECTIVE FRICTION OF THE FAULT. C C FOLLOWING PARAMETER GIVES NUMBER OF STEPS IN VERTICAL INTEGRAL C OF CREEP SHEAR STRESS ON DUCTILE PARTS OF FAULTS: PARAMETER (NSTEP=20) C HIGHER VALUES OBVIOUSLY COST MORE. ON THE OTHER HAND, SMALL VALUES C DO NOT MERELY APPROXIMATE THE CREEP LAW; THEY ALSO INTRODUCE C SOME RANDOM ERR0R WHICH CAN PUT A FLOOR ON CONVERGENCE. C C NOTE: IN VS-FORTRAN, FOLLOWING TYPE COULD BE LOGICAL*1: LOGICAL FSLIPS C LOGICAL LOCKED,PURESS,SLOPED DOUBLE PRECISION V DOUBLE PRECISION FPHI COMMON /FPHIS/ FPHI REAL MANTLE,NORMAL C DIMENSIONS PER COMMON BLOCK: DIMENSION FPHI(6,7) C DIMENSIONS OF INTERNAL CONVENIENCE ARRAYS: DIMENSION DLEPDZ(2),DSFDZ(2),RHO(2),SHEART(2),TMEAN(2),ZTRANS(2) C DIMENSIONS OF EXTERNAL ARGUMENTS ARRAYS: DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),DQDTDA(MXNODE), + FC(2,2,7,MXFEL),FDIP(3,MXFEL), + FIMUDZ(7,MXFEL),FPEAKS(2,MXFEL),FSLIPS(MXFEL), + FTAN(7,MXFEL),FTSTAR(2,7,MXFEL),NODEF(6,MXFEL), + OFFSET(MXFEL),RADIO(2),RHOBAR(2),TLNODE(MXNODE), + V(2,MXNODE),ZMNODE(MXNODE),ZTRANF(2,MXFEL) C C FOLLOWING TWO NUMBERS ARE "VERY SMALL" AND "VERY LARGE", BUT NOT C SO EXTREME AS TO CAUSE UNDERFLOW OR OVERFLOW. THEY MAY NEED TO C BE ADJUSTED, DEPENDING ON THE COMPUTER AND COMPILER BEING USED. DATA TINY /1.E-50/ DATA HUGE /1.E+50/ C CGAMMA=(1.+SIN(ATAN(CFRIC)))/(1.-SIN(ATAN(CFRIC))) DO 100 I=1,NFL IF (OFFMAX.LE.0.) THEN FRIC=FFRIC ELSE FRIC=FFRIC*(1.-BYERLY*OFFSET(I)/OFFMAX) ENDIF N1=NODEF(1,I) N2=NODEF(2,I) N3=NODEF(3,I) N4=NODEF(4,I) N5=NODEF(5,I) N6=NODEF(6,I) C C IS THIS A PURELY STRIKE-SLIP FAULT ELEMENT? PURESS=(ABS(FDIP(1,I)-1.570796).LE.WEDGE).AND. + (ABS(FDIP(2,I)-1.570796).LE.WEDGE).AND. + (ABS(FDIP(3,I)-1.570796).LE.WEDGE) C C IF SO, COMPUTE ESTIMATE OF RELATIVE NORMAL STRESS C (RELATIVE TO VERTICAL STRESS) BY USING AMOUNT OF DIVERGENCE C BETWEEN NODES N2 AND N5 (IN SPITE OF CONSTRAINT EQUATION): IF (PURESS) THEN ANGLE=FTAN(4,I) UNITBX=SIN(ANGLE) UNITBY= -COS(ANGLE) DELVX=V(1,N2)-V(1,N5) DELVY=V(2,N2)-V(2,N5) SPREAD=DELVX*UNITBX+DELVY*UNITBY DELTAU=CONSTR*SPREAD IF ((TLNODE(N2).LE.0.).OR.(ZTRANF(2,I).LE.0.)) THEN C CRUST ALONE RESISTS CONVERGENCE: DPMAX= -2.*DELTAU/ZTRANF(1,I) DDPNDZ=DPMAX/ZTRANF(1,I) ELSE C MANTLE LITHOSPHERE HELPS TO RESIST CONVERGENCE: DDPNDZ= -DELTAU/ + (0.5*ZTRANF(1,I)**2+ZTRANF(2,I)*ZMNODE(N2)+ + 0.5*ZTRANF(2,I)**2) ENDIF C DDPNDZ IS THE GRADIENT OF EXCESS NORMAL PRESSURE (IN C EXCESS OF VERTICAL PRESSURE) WITH DEPTH ON THIS FAULT; C CHECK THAT IT LIES WITHIN FRICTIONAL LIMITS OF BLOCKS: Q=0.5*(DQDTDA(N2)+DQDTDA(N5)) TTRANS=TSURF+ZTRANF(1,I)*Q/CONDUC(1)- + ZTRANF(1,I)**2*RADIO(1)/(2.*CONDUC(1)) TMEANC=(TSURF+TTRANS)/2. RHOC=RHOBAR(1)*(1.-ALPHAT(1)*TMEANC) DLEPDC=GMEAN*(RHOC-RHOH2O*BIOT) THRUST=DLEPDC*CGAMMA NORMAL=DLEPDC/CGAMMA DDPNDZ=MAX(DDPNDZ,NORMAL-DLEPDC) DDPNDZ=MIN(DDPNDZ,THRUST-DLEPDC) C ELSE C DIFFERENT LOGIC WILL BE USED; THIS PARAMETER IS NOT C REALLY NEEDED. ZERO IT JUST TO BE CAREFUL. DDPNDZ=0. ENDIF C DO 90 M=1,7 C HEAT-FLOW: Q=DQDTDA(N1)*FPHI(1,M)+DQDTDA(N2)*FPHI(2,M)+ + DQDTDA(N3)*FPHI(3,M) C C CRUSTAL THICKNESS: CRUST=ZMNODE(N1)*FPHI(1,M)+ZMNODE(N2)*FPHI(2,M)+ + ZMNODE(N3)*FPHI(3,M) C C MANTLE LITHOSPHERE THICKNESS: MANTLE=TLNODE(N1)*FPHI(1,M)+TLNODE(N2)*FPHI(2,M)+ + TLNODE(N3)*FPHI(3,M) MANTLE=MAX(MANTLE,0.) C C MOHO TEMPERATURE: TMOHO=TSURF+CRUST*Q/CONDUC(1)- + CRUST**2*RADIO(1)/(2.*CONDUC(1)) C C TEMPERATURE AT BASE OF PLATE: TASTH=TMOHO+MANTLE*(Q-CRUST*RADIO(1))/CONDUC(2)- + MANTLE**2*RADIO(2)/(2.*CONDUC(2)) C C MEAN TEMPERATURES: TMEAN(1)=(TSURF+TMOHO)/2. TMEAN(2)=(TMOHO+TASTH)/2. C C MEAN DENSITIES: RHO(1)=RHOBAR(1)*(1.-ALPHAT(1)*TMEAN(1)) RHO(2)=RHOBAR(2)*(1.-ALPHAT(2)*TMEAN(2)) C C DERIVITIVES OF LITHOSTATIC EFFECTIVE PRESSURE WRT DEPTH DLEPDZ(1)=GMEAN*(RHO(1)-RHOH2O*BIOT) EPMOHO=DLEPDZ(1)*CRUST DLEPDZ(2)=GMEAN*(RHO(2)-RHOH2O*BIOT) C C ANGLE IS THE FAULT STRIKE, IN RADIANS CCLKWS FROM +X. ANGLE=FTAN(M,I) C C UNITA IS A UNIT VECTOR ALONG THE FAULT, FROM N1 TO N3. UNITAX=COS(ANGLE) UNITAY=SIN(ANGLE) C C UNITB IS A PERPENDICULAR UNIT VECTOR, POINTING OUT C TOWARD THE N6-N4 SIDE. UNITBX= -UNITAY UNITBY= +UNITAX C C RELATIVE VELOCITIES ARE FOR N1-3 SIDE RELATIVE TO C THE N6-4 SIDE: DELVX=V(1,N1)*FPHI(1,M)+V(1,N2)*FPHI(2,M)+ + V(1,N3)*FPHI(3,M)+V(1,N4)*FPHI(4,M)+ + V(1,N5)*FPHI(5,M)+V(1,N6)*FPHI(6,M) DELVY=V(2,N1)*FPHI(1,M)+V(2,N2)*FPHI(2,M)+ + V(2,N3)*FPHI(3,M)+V(2,N4)*FPHI(4,M)+ + V(2,N5)*FPHI(5,M)+V(2,N6)*FPHI(6,M) C C SINISTRAL STRIKE-SLIP RATE COMPONENT: SINIST=DELVX*UNITAX+DELVY*UNITAY C C CONVERGENCE RATE COMPONENT (IN HORIZONTAL PLANE): CLOSE =DELVX*UNITBX+DELVY*UNITBY C C DIP OF THE FAULT (FROM HORIZONTAL ON THE N1-3 SIDE): DIP=FDIP(1,I)*FPHI(1,M)+FDIP(2,I)*FPHI(2,M)+ + FDIP(3,I)*FPHI(3,M) SLOPED=ABS(DIP-1.570796).GT.WEDGE C IF (.NOT.SLOPED) THEN C CASE OF A NEAR-VERTICAL FAULT: DSFDZ(1)=(DLEPDZ(1)+DDPNDZ)*FRIC SFMOHO=DSFDZ(2)*CRUST DSFDZ(2)=(DLEPDZ(2)+DDPNDZ)*FRIC SLIP=ABS(SINIST) LOCKED=.FALSE. ELSE C CASE OF A SHALLOW-DIPPING FAULT: C C VUPDIP IS THE UP-DIP VELOCITY COMPONENT, IN THE C FAULT PLANE, OF THE BLOCK ON THE N1-N3 SIDE. VUPDIP=CLOSE/COS(DIP) C C RAKE ANGLE IS MEASURED COUNTERCLOCKWISE IN C FAULT PLANE FROM HORIZONTAL & PARALLEL TO ANGLE. RAKE=ATAN2F(VUPDIP,SINIST) C C DERIVITIVE OF EFFECTIVE NORMAL PRESSURE C WITH RESPECT TO SHEAR TRACTION ON FAULT: DEPDST=TAN(DIP)*SIN(RAKE) C (NOTICE THAT WHEN SENSE OF DIP REVERSES, SIGN C CHANGE CAUSED BY TAN(DIP) IS CANCELLED BY SIGN C CHANGE CAUSED BY SIN(RAKE).) C C ACCORDING TO THEORY, THE EQUATION TO SOLVE IS: C D(SHEAR_TRACTION)/DZ = C "FRIC"*("DLEPDZ"+"DEPDST"*D(SHEAR_TRACTION)/DZ) C THIS MAY HAVE A PHYSICAL SOLUTION (ONE WITH C POSITIVE SHEAR_TRACTION). IF NOT, THE C FAULT IS LOCKED. LOCKED=(FRIC*DEPDST).GE.1.00 IF (LOCKED) THEN DSFDZ(1)=HUGE DSFDZ(2)=HUGE ELSE DSFDZ(1)=FRIC*DLEPDZ(1)/(1.00-FRIC*DEPDST) SFMOHO=DSFDZ(2)*CRUST DSFDZ(2)=FRIC*DLEPDZ(2)/(1.00-FRIC*DEPDST) ENDIF C SLIP=SQRT(SINIST**2+VUPDIP**2) ENDIF SLIP=MAX(SLIP,TINY) C C LOCATE PLASTIC/CREEP TRANSITION(S) C BY ITERATED HALVING OF DOMAIN: C IF (MANTLE.GT.0.) THEN LIMIT=2 ELSE LIMIT=1 ZTRANS(2)=0. SHEART(2)=0. ENDIF DO 60 LAYER=1,LIMIT TOPZ=0. IF (LAYER.EQ.1) THEN BASEZ=CRUST SF0=0. T0=TSURF Q0=Q Z0=0. ELSE BASEZ=MANTLE SF0=SFMOHO T0=TMOHO Q0=Q-CRUST*RADIO(1) Z0=CRUST ENDIF DO 50 KITER=1,15 Z=0.5*(TOPZ+BASEZ) ZABS=Z+Z0 SHEARF=Z*DSFDZ(LAYER)+SF0 SHEARP=MIN(SHEARF,DCREEP(LAYER)) T=T0+Q0*Z/CONDUC(LAYER)-(RADIO(LAYER)/ + (2.*CONDUC(LAYER)))*Z**2 IF (ZABS.LE.17.*ONEKM) THEN T90PC=0.5*ZABS ELSE T90PC=25.*ONEKM-2.83*ZABS+ + 0.11111*ONEKM*(ZABS/ONEKM)**2 ENDIF C SEE TURCOTTE ET AL (1980) JGR, 85, B11, 6224-6230 STRAIN=SLIP/T90PC SHEARC=ACREEP(LAYER)*(STRAIN**ECREEP)* + EXP((BCREEP(LAYER)+CCREEP(LAYER)*Z)/T) IF (SHEARC.LT.SHEARP) THEN BASEZ=Z ELSE TOPZ=Z ENDIF 50 CONTINUE ZTRANS(LAYER)=0.5*(TOPZ+BASEZ) SHEART(LAYER)=ZTRANS(LAYER)*DSFDZ(LAYER)+SF0 60 CONTINUE C C PLASTIC PART OF VERTICAL INTEGRAL(S) OF TRACTION: C (A) CRUST: IF (SHEART(1).LE.DCREEP(1)) THEN VITDZ=0.5*SHEART(1)*ZTRANS(1) ELSE ZP=ZTRANS(1)*DCREEP(1)/SHEART(1) VITDZ=DCREEP(1)*(ZTRANS(1)-0.5*ZP) ENDIF C (B) MANTLE LITHOSPHERE: IF ((MANTLE.GT.0.).AND.(SHEART(2).GT.SFMOHO)) THEN IF (SHEART(2).LE.DCREEP(2)) THEN VITDZ=VITDZ+0.5*(SFMOHO+SHEART(2))*ZTRANS(2) ELSE ZP=ZTRANS(2)*(DCREEP(2)-SFMOHO)/ + (SHEART(2)-SFMOHO) ZP=MAX(ZP,0.) VITDZ=VITDZ+0.5*(SFMOHO+SHEART(2))*ZP+ + DCREEP(2)*(ZTRANS(2)-ZP) ENDIF ENDIF C C ADD CREEP PART(S) OF INTEGRAL, USING PARABOLIC RULE C SUM=0. DO 80 LAYER=1,LIMIT IF (LAYER.EQ.1) THEN THICK=CRUST T0=TSURF Q0=Q ZABS=0. ELSE THICK=MANTLE T0=TMOHO Q0=Q-CRUST*RADIO(1) ZABS=CRUST ENDIF DZ=(THICK-ZTRANS(LAYER))/NSTEP OLDSC=SHEART(LAYER) OLDSC=MIN(OLDSC,DCREEP(LAYER)) Z0=ZTRANS(LAYER) DO 70 J=1,NSTEP ZHALF=Z0+0.5*DZ ZFULL=Z0+DZ AZHALF=ZHALF+ZABS AZFULL=ZFULL+ZABS THALF=T0+Q0*ZHALF/CONDUC(LAYER)- + (RADIO(LAYER)/ + (2.*CONDUC(LAYER)))*ZHALF**2 TFULL=T0+Q0*ZFULL/CONDUC(LAYER)- + (RADIO(LAYER)/ + (2.*CONDUC(LAYER)))*ZFULL**2 IF (AZHALF.LE.17.*ONEKM) THEN WHALF=0.5*AZHALF ELSE WHALF=25.*ONEKM-2.83*AZHALF+ + 0.11111*ONEKM*(AZHALF/ONEKM)**2 ENDIF IF (AZFULL.LE.17.*ONEKM) THEN WFULL=0.5*AZFULL ELSE WFULL=25.*ONEKM-2.83*AZFULL+ + 0.11111*ONEKM*(AZFULL/ONEKM)**2 ENDIF C SEE TURCOTTE ET AL (1980) JGR, 85, B11, 6224-6230 EHALF=SLIP/WHALF EFULL=SLIP/WFULL SCHALF=ACREEP(LAYER)*(EHALF**ECREEP)* + EXP((BCREEP(LAYER)+CCREEP(LAYER)*ZHALF) + /THALF) SCHALF=MIN(SCHALF,DCREEP(LAYER)) SCFULL=ACREEP(LAYER)*(EFULL**ECREEP)* + EXP((BCREEP(LAYER)+CCREEP(LAYER)*ZFULL) + /TFULL) SCFULL=MIN(SCFULL,DCREEP(LAYER)) SUM=SUM+DZ*(0.1666667*OLDSC+ + 0.6666667*SCHALF+ + 0.1666666*SCFULL) Z0=ZFULL OLDSC=SCFULL 70 CONTINUE 80 CONTINUE C VITDZ=VITDZ+SUM C VIMUDZ=VITDZ/SLIP C FIMUDZ(M,I)=MIN(VIMUDZ,FMUMAX*(CRUST+MANTLE)) C C DIPPING, OBLIQUE-SLIP INTEGRATION C POINTS ARE ALSO CHARACTERIZED C BY "FC" AND "FTSTAR": C IF (SLOPED) THEN TS=SINIST*FIMUDZ(M,I) TU=VUPDIP*FIMUDZ(M,I) IF (LOCKED) THEN FC(1,1,M,I)=FIMUDZ(M,I) FC(1,2,M,I)=0. FC(2,1,M,I)=0. FC(2,2,M,I)=FIMUDZ(M,I) ELSE SINR=SIN(RAKE) COSR=COS(RAKE) TAND=TAN(DIP) C C *** IMPORTANT NOTE: *** C THE FOLLOWING 7 STATEMENTS ARE -NOT- THE C RESULT OF THEORY, BUT A TACTICAL CHOICE C WHICH ATTEMPTS TO COMPROMISE BETWEEN C STABILITY OF THE LINEAR SYSTEM, STABILITY C OF THE ITERATION, AND EFFICIENCY. C THEY MAY BE CHANGED IF THE PROGRAM DOES C NOT CONVERGE SATISFACTORILY] C TUNE=2. FC(1,1,M,I)=FIMUDZ(M,I)* + (1.-TUNE*SINR*COSR**2*TAND) FC(1,2,M,I)=FIMUDZ(M,I)* + (TUNE*COSR**3*TAND) FC(2,1,M,I)=FIMUDZ(M,I)* + (-TUNE*SINR**2*COSR*TAND) FC(2,2,M,I)=FIMUDZ(M,I)* + (1.+TUNE*SINR*COSR**2*TAND) C (OFTEN, FC(1,2) IS THE BIGGEST TERM. C IN SOME CASES, DIAGONALS BECOME NEGATIVE. C FOR STABILITY, BE SURE THAT THE FC C MATRIX REMAINS POSITIVE DEFINITE: FC(1,1,M,I)=MAX(FC(1,1,M,I),ABS(FC(1,2,M,I))) FC(2,2,M,I)=MAX(FC(2,2,M,I),ABS(FC(1,2,M,I))) ENDIF FTSTAR(1,M,I)=TS-FC(1,1,M,I)*SINIST- + FC(1,2,M,I)*VUPDIP FTSTAR(2,M,I)=TU-FC(2,1,M,I)*SINIST- + FC(2,2,M,I)*VUPDIP ENDIF C C PROVIDE INTERESTING DIAGNOSTIC DATA AT MIDPOINTS ONLY: C IF (M.EQ.4) THEN FSLIPS(I)=(.NOT.LOCKED).AND. + (FIMUDZ(M,I).LT.(0.99*FMUMAX*(CRUST+MANTLE))) ZTRANF(1,I)=ZTRANS(1) FPEAKS(1,I)=SHEART(1) ZTRANF(2,I)=ZTRANS(2) FPEAKS(2,I)=SHEART(2) ENDIF C 90 CONTINUE 100 CONTINUE RETURN END C C C SUBROUTINE OLDVEL (INPUT,IUNITV,MXNODE,NUMNOD, + OUTPUT,ALDONE,TITLE1,TITLE2,TITLE3,V) C C READ OLD VELOCITY SOLUTION FROM UNIT IUNITV, OR ELSE SET FLAG C "ALDONE". C LOGICAL ALDONE CHARACTER*80 TITLE1,TITLE2,TITLE3 DOUBLE PRECISION V DIMENSION V(2,MXNODE) C READ (IUNITV,'(A80)',END=100,ERR=100) TITLE1 READ (IUNITV,'(A80)',END=100,ERR=100) TITLE2 READ (IUNITV,'(A80)',END=100,ERR=100) TITLE3 READ (IUNITV,*,END=100,ERR=100) ((V(J,I),J=1,2),I=1,NUMNOD) ALDONE=.FALSE. RETURN C ------------------(THIS SECTION EXECUTED ONLY IF READ FAILS)--------- 100 ALDONE=.TRUE. RETURN END C C C SUBROUTINE PRINCE (INPUT,E11,E22,E12, + OUTPUT,E1,E2,U1X,U1Y,U2X,U2Y) C C FIND PRINCIPAL VALUES (E1,E2) OF THE SYMMETRIC 2X2 TENSOR E11 E12 C E12 E22 C AND ALSO THE ASSOCIATED EIGENVECTORS #1=(U1X,U1Y),#2=(U2X,U2Y). C THE CONVENTION IS THAT E1 <= E2. C R=SQRT(((E11-E22)/2.)**2+E12**2) C=(E11+E22)/2. E1=C-R E2=C+R THETA=ATAN2F(E11-E1, -E12) U1X=COS(THETA) U1Y=SIN(THETA) U2X=U1Y U2Y= -U1X RETURN END C C C SUBROUTINE READPM (INPUT,IUNIT7, IUNIT8, OFFMAX, + OUTPUT,ACREEP, ALPHAT, BCREEP, BIOT , + BYERLY, CCREEP, CFRIC , CONDUC, + DCREEP, ECREEP, EVERYP, FFRIC , GMEAN , + MAXITR, OKDELV, OKTOQT, ONEKM, RADIO , + REFSTR, RHOAST, RHOBAR, RHOH2O, + TEMLIM, TITLE3, TRHMAX, TSURF) C C READS STRATEGIC AND TACTICAL INPUT PARAMETERS FROM DEVICE IUNIT7, C AND ECHOES THEM ON DEVICE IUNIT8 WITH ANNOTATIONS. C CHARACTER*5 COLOR CHARACTER*80 TITLE3 LOGICAL EVERYP DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),RADIO(2), + RHOBAR(2),TEMLIM(2) C WRITE (IUNIT8,1) IUNIT7 1 FORMAT (/ /' ATTEMPTING TO READ PARAMETERS FROM UNIT',I3) READ (IUNIT7,2) TITLE3 2 FORMAT (A80) WRITE (IUNIT8,3) TITLE3 3 FORMAT (/' TITLE OF PARAMETER SET ='/' ',A80) WRITE (IUNIT8,4) 4 FORMAT (/' **************************************************'/ + ' IT IS THE USERS RESPONSIBILITY TO INPUT ALL OF THE'/ + ' FOLLOWING NUMERICAL QUANTITIES IN CONSISTENT UNITS,'/ + ' SUCH AS SYSTEM-INTERNATIONAL (SI) OR CM-G-S (CGS).'/ + ' NOTE THAT TIME UNIT MUST BE THE SECOND (HARD-CODED).'/ + ' **************************************************'/ + /' ========== STRATEGIC PARAMETERS (DEFINE THE REAL', + '-EARTH PROBLEM) ======'/) READ (IUNIT7,*) FFRIC WRITE (IUNIT8,20) FFRIC 20 FORMAT (' ', F10.3,' COEFFICIENT OF FRICTION ON FAULTS') READ (IUNIT7,*) CFRIC WRITE (IUNIT8,30) CFRIC 30 FORMAT (' ', F10.3,' COEFFICIENT OF FRICTION WITHIN BLOCKS') READ (IUNIT7,*) BIOT BIOT = MAX(0.0,MIN(1.0,BIOT)) WRITE (IUNIT8,40) BIOT 40 FORMAT (' ',F10.4,' EFFECTIVE-PRESSURE (BIOT) COEFFICIENT,', + ' 0.0 TO 1.0') READ (IUNIT7,*) BYERLY BYERLY = MAX(0.0,MIN(0.99,BYERLY)) IF (OFFMAX.GT.0.) THEN WRITE (IUNIT8,41) BYERLY 41 FORMAT (' ',F10.4,' BYERLY COEFFICIENT (0. TO 0.99) ='/ + 11X,' FRACTIONAL FRICTION REDUCTION ON MASTER', + ' FAULT'/ + 11X,' (OTHER FAULTS HAVE LESS REDUCTION, IN', + ' PROPORTION TO'/ + 11X,' THEIR TOTAL PAST OFFSETS)') ELSE WRITE (IUNIT8,42) BYERLY 42 FORMAT (' ',F10.4,' BYERLY COEFFICIENT (NOT USED IN', + ' THIS RUN,'/ + 11X,' AS ALL OFFSETS ARE ZERO)') ENDIF CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,ACREEP) IF (ACREEP(2).EQ.0.) ACREEP(2)=ACREEP(1) WRITE (IUNIT8,50) ACREEP(1),ACREEP(2) 50 FORMAT (' ',1P, E10.2,'/',E10.2,' A FOR CREEP = ', + 'PRE-EXPONENTIAL SHEAR', + ' STRESS CONSTANT FOR CREEP. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,BCREEP) IF (BCREEP(2).EQ.0.) BCREEP(2)=BCREEP(1) WRITE (IUNIT8,60) BCREEP(1),BCREEP(2) 60 FORMAT (' ', F10.0,'/',F10.0,' B FOR CREEP =(ACTIVATION ', + 'ENERGY)/R/N', + ' IN K. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,CCREEP) IF (CCREEP(2).EQ.0.) CCREEP(2)=CCREEP(1) WRITE (IUNIT8,70) CCREEP(1),CCREEP(2) 70 FORMAT (' ',1P, E10.2,'/',E10.2,' C FOR CREEP = DERIVATIVE OF B', + ' WITH RESPECT TO DEPTH. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,DCREEP) IF (DCREEP(2).EQ.0.) DCREEP(2)=DCREEP(1) WRITE (IUNIT8,80) DCREEP(1),DCREEP(2) 80 FORMAT (' ',1P, E10.2,'/',E10.2,' D FOR CREEP = MAXIMUM SHEAR ', + 'STRESS UNDER ANY CONDITIONS. (CRUST/MANTLE)') READ (IUNIT7,*) ECREEP WRITE (IUNIT8,90) ECREEP 90 FORMAT (' ', F10.6,' E FOR CREEP = STRAIN-RATE EXPONENT FOR', + ' CREEP (1/N). (SAME FOR CRUST AND MANTLE])') READ (IUNIT7,*) TRHMAX IF (TRHMAX.EQ.1.) TRHMAX=9.9E39 WRITE (IUNIT8,101) TRHMAX 101 FORMAT (' ',1P,E10.2,' LIMIT ON HORIZONTAL TRACTIONS', + ' APPLIED TO BASE OF PLATE') READ (IUNIT7,*) RHOH2O WRITE (IUNIT8,110) RHOH2O 110 FORMAT (' ',1P,E10.3,' DENSITY OF GROUNDWATER, LAKES, & OCEANS') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,RHOBAR) IF (RHOBAR(2).EQ.0.) RHOBAR(2)=RHOBAR(1) WRITE (IUNIT8,120) RHOBAR(1),RHOBAR(2) 120 FORMAT (' ',1P,E10.3,'/',E10.3,' MEAN DENSITY,', + ' CORRECTED TO 0 DEGREES KELVIN. (CRUST/MANTLE)') READ (IUNIT7,*) RHOAST WRITE (IUNIT8,130) RHOAST 130 FORMAT (' ',1P,E10.3,' DENSITY OF ASTHENOSPHERE') READ (IUNIT7,*) GMEAN WRITE (IUNIT8,140) GMEAN 140 FORMAT (' ',1P,E10.3,' MEAN GRAVITATIONAL ACCELERATION', + ' (LENGTH/SEC**2)') READ (IUNIT7,*) ONEKM WRITE (IUNIT8,150) ONEKM 150 FORMAT (' ',1P,E10.3,' NUMBER OF LENGTH UNITS NEEDED TO', + ' MAKE 1 KILOMETER'/11X, + ' (E.G., 1000. IN SI, 1.E5 IN CGS)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,ALPHAT) IF (ALPHAT(2).EQ.0.) ALPHAT(2)=ALPHAT(1) WRITE (IUNIT8,160) ALPHAT(1),ALPHAT(2) 160 FORMAT (' ',1P,E10.2,'/',E10.2,' VOLUMETERIC THERMAL ', + 'EXPANSION OF CRUST', + ' (1/VOL)*(D.VOL/D.T). (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,CONDUC) IF (CONDUC(2).EQ.0.) CONDUC(2)=CONDUC(1) WRITE (IUNIT8,170) CONDUC(1),CONDUC(2) 170 FORMAT (' ',1P,E10.2,'/',E10.2,' THERMAL CONDUCTIVITY, ENERGY/', + 'LENGTH/SEC/DEG. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,RADIO) IF (RADIO(2).EQ.0.) RADIO(2)=RADIO(1) WRITE (IUNIT8,180) RADIO(1),RADIO(2) 180 FORMAT (' ',1P,E10.2,'/',E10.2,' RADIOACTIVE HEAT PRODUCTION', + ' ENERGY/VOLUME/SEC. (CRUST/MANTLE)') READ (IUNIT7,*) TSURF WRITE (IUNIT8,185) TSURF 185 FORMAT (' ', F10.0,' SURFACE TEMPERATURE, ON', + ' ABSOLUTE SCALE') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,TEMLIM) IF (TEMLIM(2).EQ.0.) TEMLIM(2)=TEMLIM(1) WRITE (IUNIT8,190) TEMLIM(1),TEMLIM(2) 190 FORMAT (' ', F10.0,'/',F10.0,' CONVECTING TEMPERATURE (TMAX), ON' + ,' ABSOLUTE SCALE. (CRUST/MANTLE)') WRITE (IUNIT8,199) 199 FORMAT (/' ========== TACTICAL PARAMETERS (HOW TO REACH ', + 'THE SOLUTION) =========='/) READ (IUNIT7,*) MAXITR WRITE (IUNIT8,200) MAXITR 200 FORMAT (' ',I10,' MAXIMUM ITERATIONS WITHIN VELOCITY SOLUTION') READ (IUNIT7,*) OKTOQT WRITE (IUNIT8,210) OKTOQT 210 FORMAT (' ',F10.6,' ACCEPTABLE FRACTIONAL CHANGE IN VELOCITY ', + '(STOPS ITERATION EARLY)') READ (IUNIT7,*) REFSTR WRITE (IUNIT8,220) REFSTR 220 FORMAT (' ',1P,E10.2,' EXPECTED MEAN VALUE OF SHEAR STRESS IN', + ' CRUST'/' ',10X, + ' (USED TO INITIALIZE AND SET STIFFNESS LIMITS)') READ (IUNIT7,*) OKDELV WRITE (IUNIT8,230) OKDELV 230 FORMAT (' ',1P,E10.2,' MAGNITUDE OF VELOCITY ERR0RS ALLOWED', + ' DUE TO FINITE STIFFNESS'/11X, + '(SUCH ERR0RS MAY APPEAR IN SUCH FORMS AS:'/11X, + ' 1. FICTICIOUS BASAL SLIP OF CRUST OVER MANTLE'/11X, + ' 2. ERRONEOUS CONVERGENCE/DIVERGENCE AT VERTICAL FAULTS'/ + 11X, + ' 3. VELOCITY EFFECT OF FICTICIOUS VISCOUS COMPLIANCES'/11X, + ' HOWEVER, VALUES WHICH ARE TOO SMALL WILL CAUSE ILL-CONDITIONED' + /11X, + ' LINEAR SYSTEMS AND STRESS ERR0RS, ', + 'AND MAY PREVENT CONVERGENCE])') READ (IUNIT7,*) EVERYP WRITE (IUNIT8,240) EVERYP 240 FORMAT (' ',L10,' SHOULD NODAL VELOCITIES BE OUTPUT EVERY STE', + 'P? (FOR CONVERGENCE STUDIES)') WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') DO 2000 I=1,999 READ (IUNIT7,1001) COLOR 1001 FORMAT (11X,A5) IF (COLOR.EQ.'COLOR') THEN WRITE (IUNIT8,1002) 1002 FORMAT (' SKIPPED OVER PLOT-CONTROL LINES') WRITE (IUNIT8,999) RETURN ENDIF 2000 CONTINUE 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 AGREED,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) ELEMENT OR FAULT ELEMENT; C DO 110 I=1,NREALN CHECKN(I)=.FALSE. 110 CONTINUE DO 130 I=1,NUMEL DO 120 J=1,6 CHECKN(NODES(J,I))=.TRUE. 120 CONTINUE 130 CONTINUE DO 132 I=1,NFL DO 131 J=1,6 CHECKN(NODEF(J,I))=.TRUE. 131 CONTINUE 132 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 CONTINUUM ELEMENT OR FAULT:') 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. C (MEANS "NOT YET INVOLVED IN AVERAGING') 410 CONTINUE DO 490 I=1,NFL DO 480 J1=1,3,2 NJ1=NODEF(J1,I) C (FAULT ENDS ARE THE ONLY PLACES THAT CAN HAVE PROBLEMS) IF (.NOT.CHECKN(NJ1)) THEN LIST(1)=NJ1 CHECKN(NJ1)=.TRUE. C BEGIN LIST OF NEIGHBORS WITH PAIRED NODE J2=7-J1 NJ2=NODEF(J2,I) LIST(2)=NJ2 CHECKN(NJ2)=.TRUE. NINSUM=2 C FIND SHORTEST FAULT CONNECTED TO EITHER ONE SHORT=SQRT( + (XNODE(NODEF(1,I))-XNODE(NODEF(3,I)))**2+ + (YNODE(NODEF(1,I))-YNODE(NODEF(3,I)))**2) DO 470 K=1,NFL NL1=NODEF(1,K) NL3=NODEF(3,K) NL4=NODEF(4,K) NL6=NODEF(6,K) IF ((NJ1.EQ.NL1).OR.(NJ2.EQ.NL1).OR. + (NJ1.EQ.NL3).OR.(NJ2.EQ.NL3).OR. + (NJ1.EQ.NL4).OR.(NJ2.EQ.NL4).OR. + (NJ1.EQ.NL6).OR.(NJ2.EQ.NL6)) THEN TEST=SQRT( + (XNODE(NL1)-XNODE(NL3))**2+ + (YNODE(NL1)-YNODE(NL3))**2) SHORT=MIN(SHORT,TEST) ENDIF 470 CONTINUE C COLLECT ALL CORNER NODES WITHIN 10% OF THIS TOLER=SHORT/10. T2=TOLER**2 DO 471 K=1,NUMNOD IF (.NOT.CHECKN(K)) THEN IF (NODTYP(K).EQ.1) THEN R2=(XNODE(NJ1)-XNODE(K))**2+ + (YNODE(NJ1)-YNODE(K))**2 IF (R2.LT.T2) THEN NINSUM=NINSUM+1 LIST(NINSUM)=K CHECKN(K)=.TRUE. ENDIF ENDIF ENDIF 471 CONTINUE C (QUICK EXIT IF ALL NODES IN SAME PLACE) AGREED=.TRUE. DO 472 K=2,NINSUM AGREED=AGREED.AND. + (XNODE(K).EQ.XNODE(1)).AND. + (YNODE(K).EQ.YNODE(1)) 472 CONTINUE IF (AGREED) GO TO 480 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,476) NINSUM, + (LIST(N),N=1,NINSUM) 476 FORMAT(/ + ' AVERAGING TOGETHER THE POSITIONS OF', + ' THESE ',I6,' NODES:',(/' ',12I6)) WRITE (IUNIT8,477) RMAX 477 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: NDIP=1+J/2 IF (ABS(FDIP(NDIP,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(3,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 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 C C C SUBROUTINE READN (INPUT,IUNITP,IUNITT,N, + OUTPUT,VECTOR) C C A UTILITY ROUTINE DESIGNED TO PERMIT -FAULTS- INPUT FILES C TO ALSO BE USED BY -PLATES-, WHICH EXPECTS MORE NUMBERS C IN SOME RECORDS. C THIS ROUTINE ATTEMPTS TO READ 'N' FLOATING-POINT VALUES C (USING * FORMAT) FROM THE NEXT RECORD ON DEVICE 'IUNITP'. C IF ANYTHING GOES WRONG, THE MISSING VALUES ARE SET TO ZERO. C CHARACTER*1 C CHARACTER*80 LINE LOGICAL ANYIN,DOTTED,EXPON,SIGNED DIMENSION VECTOR(N) C READ (IUNITP,1) LINE 1 FORMAT (A80) C NUMBER=0 ANYIN=.FALSE. EXPON=.FALSE. SIGNED=.FALSE. DOTTED=.FALSE. DO 10 I=1,80 C=LINE(I:I) IF ((C.EQ.' ').OR.(C.EQ.',').OR.(C.EQ.'/')) THEN SIGNED=.FALSE. EXPON=.FALSE. DOTTED=.FALSE. IF (ANYIN) THEN NUMBER=NUMBER+1 ANYIN=.FALSE. ENDIF ELSE IF ((C.EQ.'+').OR.(C.EQ.'-')) THEN IF (SIGNED) THEN GO TO 50 ELSE SIGNED=.TRUE. ENDIF ELSE IF ((C.EQ.'E').OR.(C.EQ.'D').OR. + (C.EQ.'e').OR.(C.EQ.'d')) THEN IF (EXPON) THEN GO TO 50 ELSE EXPON=.TRUE. SIGNED=.FALSE. DOTTED=.TRUE. ENDIF ELSE IF (C.EQ.'.') THEN IF (DOTTED) THEN GO TO 50 ELSE DOTTED=.TRUE. ENDIF ELSE IF ((C.EQ.'0').OR.(C.EQ.'1').OR.(C.EQ.'2').OR. + (C.EQ.'3').OR.(C.EQ.'4').OR.(C.EQ.'5').OR. + (C.EQ.'6').OR.(C.EQ.'7').OR.(C.EQ.'8').OR. + (C.EQ.'9')) THEN SIGNED=.TRUE. ANYIN=.TRUE. ELSE GO TO 50 ENDIF 10 CONTINUE IF (ANYIN) NUMBER=NUMBER+1 C 50 IF (NUMBER.EQ.0) THEN WRITE (IUNITT,91) N,LINE 91 FORMAT (/' ERR0R: A LINE OF PARAMETER INPUT WHICH', + ' WAS SUPPOSED TO CONTAIN 1-',I2,' NUMBERS'/ + ' COULD NOT BE INTERPRETED. LINE FOLLOWS:'/ + ' ',A80) STOP ELSE NUMBER=MIN(NUMBER,N) BACKSPACE IUNITP READ (IUNITP,*) (VECTOR(I),I=1,NUMBER) IF (NUMBER.LT.N) THEN DO 99 I=NUMBER+1,N VECTOR(I)=0. 99 CONTINUE ENDIF ENDIF RETURN END C C C SUBROUTINE LLTOXY (INPUT,CPNLAT, + PLAT,PLON, + RADIUS,RTAN,X0ELON,YPOLE, + OUTPUT,X,Y) C C CONVERT A (NORTH LATITUDE=PLAT, EAST LONGITUDE=PLON) POSITION C INTO AN (X,Y) POSITION ON A CONIC PROJECTION WITH TANGENT C LATITUDE CPNLAT, 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 TANGENT LATITUDE (CPNLAT), ALONG A MERIDIAN WHICH C IS ON THE OPPOSITE SIDE OF THE EARTH FROM X0ELON. C IF PLAT IS MORE THAN 90 DEGREES DIFFERENT FROM CPNLAT, THE C POINT DOES NOT FALL ONTO THE PROJECTION AT ALL. TO PREVENT C CRASHES, WE MERELY PLACE IT VERY FAR OUT ON THE PROJECTION. C C NOTE: FOLLOWING TWO LINES ARE PRECOMPUTED AND PASSED TO SAVE TIME: C RTAN=RADIUS*TANDEG(90.-CPNLAT) C YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C IF (ABS(PLAT-CPNLAT).GE.90.) PLAT=CPNLAT+89.*(PLAT-CPNLAT)/ + ABS(PLAT-CPNLAT) R=RTAN-RADIUS*TANDEG(PLAT-CPNLAT) DPLON=PLON-X0ELON IF (DPLON.LT.-180.) DPLON=DPLON+360. IF (DPLON.GT. 180.) DPLON=DPLON-360. ANGLE=DPLON*SINDEG(CPNLAT) X=R*SINDEG(ANGLE) Y=YPOLE-R*COSDEG(ANGLE) RETURN END C C C SUBROUTINE XYTOLL (INPUT,CPNLAT,RADIUS,RTAN,X0ELON,YPOLE, + X,Y, + OUTPUT,PLAT,PLON) C C CONVERT POINTS EXPRESSED AS (X,Y) ON A CONIC PROJECTION PLANE C WITH TANGENT LATITUDE CPNLAT AND ORIGIN AT (Y0NLAT,X0ELON) C TO (PLAT = NORTH_LATITUDE, PLON = EAST_LONGITUDE) C IN DEGREES C C NOTE: FOLLOWING TWO VARIABLES ARE PRECOMPUTED TO SAVE TIME: C RTAN=RADIUS*TANDEG(90.-CPNLAT) C YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C YRP=Y-YPOLE R=SQRT(X**2+YRP**2) ANGLE=57.29578*ATAN2F(X,-YRP) PLON=X0ELON+ANGLE/SINDEG(CPNLAT) PLAT=CPNLAT+57.29578*ATAN((RTAN-R)/RADIUS) PLAT=MIN(90.,MAX(PLAT,-90.)) IF ((PLON-X0ELON).GT. 180.) PLON=PLON-360. IF ((PLON-X0ELON).GT. 180.) PLON=PLON-360. IF ((PLON-X0ELON).LT.-180.) PLON=PLON+360. IF ((PLON-X0ELON).LT.-180.) PLON=PLON+360. RETURN END C****** ALASKA (BEGIN SECTION) ****************** C C C SUBROUTINE BIGDIS (INPUT,AREA,BENCHX,BENCHY, + DXS,DYS, + MXNODE,MXEL, + NBENCH,NODES,NUMEL, + TLNODE,TRHMAX,V, + XNODE,YNODE,ZMNODE, + MODIFY,UX,UY) C C AUGMENT (UX,UY) VELOCITIES AT BENCHMARKS C (DUE TO POSITIVE ELASTIC DISLOCATIONS ACCUMULATING C IN THE ALEUTIAN TRENCH AT THE RATE IMPLIED BY C THE RELATIVE MOTION BETWEEN THE PACIFIC PLATE AND C THE NORTH AMERICAN PLATE: C DOUBLE PRECISION AZIMF,BX,BY,DIPF,LF,SLIP,XF,YF,ZBOT,ZTOP, + V,VXB,VYB LOGICAL DIS05,RESIST DIMENSION AREA(MXEL),BENCHX(NBENCH),BENCHY(NBENCH), + DXS(6,7,MXEL),DYS(6,7,MXEL),NODES(6,MXEL), + SLIP(3),TLNODE(MXNODE),ZMNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE), + UX(NBENCH),UY(NBENCH),V(2,MXNODE) 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 DO 1000 I=1,NUMEL DO 900 M=1,7 N1=NODES(1,I) N2=NODES(2,I) N3=NODES(3,I) N4=NODES(4,I) N5=NODES(5,I) N6=NODES(6,I) X=PHI(1,M)*XNODE(N1)+ + PHI(2,M)*XNODE(N2)+ + PHI(3,M)*XNODE(N3)+ + PHI(4,M)*XNODE(N4)+ + PHI(5,M)*XNODE(N5)+ + PHI(6,M)*XNODE(N6) Y=PHI(1,M)*YNODE(N1)+ + PHI(2,M)*YNODE(N2)+ + PHI(3,M)*YNODE(N3)+ + PHI(4,M)*YNODE(N4)+ + PHI(5,M)*YNODE(N5)+ + PHI(6,M)*YNODE(N6) VCX=PHI(1,M)*V(1,N1)+ + PHI(2,M)*V(1,N2)+ + PHI(3,M)*V(1,N3)+ + PHI(4,M)*V(1,N4)+ + PHI(5,M)*V(1,N5)+ + PHI(6,M)*V(1,N6) VCY=PHI(1,M)*V(2,N1)+ + PHI(2,M)*V(2,N2)+ + PHI(3,M)*V(2,N3)+ + PHI(4,M)*V(2,N4)+ + PHI(5,M)*V(2,N5)+ + PHI(6,M)*V(2,N6) CALL BOTTOM (INPUT,TRHMAX,VCX,VCY,X,Y, + OUTPUT,RESIST,VMX,VMY) IF (RESIST) THEN DEEP=PHI(1,M)*(ZMNODE(N1)+TLNODE(N1))+ + PHI(2,M)*(ZMNODE(N2)+TLNODE(N2))+ + PHI(3,M)*(ZMNODE(N3)+TLNODE(N3))+ + PHI(4,M)*(ZMNODE(N4)+TLNODE(N4))+ + PHI(5,M)*(ZMNODE(N5)+TLNODE(N5))+ + PHI(6,M)*(ZMNODE(N6)+TLNODE(N6)) PATCH=AREA(I)*WEIGHT(M) SIDE=SQRT(PATCH) DTHDX=DXS(1,M,I)*(ZMNODE(N1)+TLNODE(N1))+ + DXS(2,M,I)*(ZMNODE(N2)+TLNODE(N2))+ + DXS(3,M,I)*(ZMNODE(N3)+TLNODE(N3))+ + DXS(4,M,I)*(ZMNODE(N4)+TLNODE(N4))+ + DXS(5,M,I)*(ZMNODE(N5)+TLNODE(N5))+ + DXS(6,M,I)*(ZMNODE(N6)+TLNODE(N6)) DTHDY=DYS(1,M,I)*(ZMNODE(N1)+TLNODE(N1))+ + DYS(2,M,I)*(ZMNODE(N2)+TLNODE(N2))+ + DYS(3,M,I)*(ZMNODE(N3)+TLNODE(N3))+ + DYS(4,M,I)*(ZMNODE(N4)+TLNODE(N4))+ + DYS(5,M,I)*(ZMNODE(N5)+TLNODE(N5))+ + DYS(6,M,I)*(ZMNODE(N6)+TLNODE(N6)) DIP=ATAN(SQRT(DTHDX**2+DTHDY**2)) TOWARD=ATAN2F(DTHDY,DTHDX) AZIMF=TOWARD-1.570796327D0 DIPF=3.141592654D0-DIP LF=0.5D0*SIDE SLIP(1)=VMX-VCX SLIP(2)=VMY-VCY SLIP(3)=(-SLIP(1)*COS(TOWARD) + -SLIP(2)*SIN(TOWARD))*TAN(DIP) AWAY=DEEP/TAN(DIP) XF=X-AWAY*COS(TOWARD) YF=Y-AWAY*SIN(TOWARD) ZBOT=DEEP+0.5*SIDE*SIN(DIP) ZTOP=DEEP-0.5*SIDE*SIN(DIP) DO 800 N=1,NBENCH BX=BENCHX(N) BY=BENCHY(N) CALL CHANGE (INPUT,AZIMF,BX,BY,DIPF, + LF,SLIP,XF,YF,ZBOT,ZTOP, + OUTPUT,VXB,VYB,DIS05) UX(N)=UX(N)+VXB UY(N)=UY(N)+VYB IF (N.EQ.76) THEN CCCC WRITE (6,799)I,M,AZIMF,BX,BY,DIPF,LF, CCCC + SLIP(1),SLIP(2),SLIP(3),XF,YF, CCCC + ZBOT,ZTOP,VXB,VYB,UX(N),UY(N) C 799 FORMAT (/' I=',I5,', M=',I1 C + /' AZIMF=',1P,D10.2 C + /' BX=',D10.2 C + /' BY=',D10.2 C + /' DIPF=',D10.2 C + /' LF=',D10.2 C + /' SLIP(1)=',D10.2 C + /' SLIP(2)=',D10.2 C + /' SLIP(3)=',D10.2 C + /' XF=',D10.2 C + /' YF=',D10.2 C + /' ZBOT=',D10.2 C + /' ZTOP=',D10.2 C + /' VXB=',D10.2 C + /' VYB=',D10.2 C + /' UXB=',E10.2 CCCCC+ /' UYB=',E10.2) ENDIF 800 CONTINUE ENDIF 900 CONTINUE 1000 CONTINUE C RETURN END C C C SUBROUTINE BOTTOM (INPUT,TRHMAX,VCX,VCY,X,Y, + OUTPUT,RESIST,VMX,VMY) C C COMPUTES HORIZONTAL COMPONENTS OF FLOW AT TOP OF ANY SUBDUCTING C SLABS OR OTHER STRONG ELEMENTS WHICH MAY APPLY TRACTIONS TO THE C BASE OF THE PLATE. C C (BUT, ALL SUCH TRACTIONS ARE TURNED OFF IF TRHMAX=0.) C C**************************************************************** C CAVEAT HACKER ]]] C UNLIKE OTHER SUBPROGRAMS IN THIS PACKAGE, "BOTTOM" IS VERY C SPECIFIC TO A PARTICULAR PROBLEM: C -IT ONLY DESCRIBES THE PACIFIC/NORTH AMERICAN BOUNDARY IN THE C REGION OF THE ALASKAN-ALEUTIAN ARC. C -IT ASSUMES A PARTICULAR ORIGIN AND ORIENTATION OF THE X-AXIS. C (ORIGIN AT 61 N, 147 W, WITH +X POINTING E) C -IT ASSUMES THAT INPUT COORDINATES ARE IN METERS. C -IT IS BASED ON A PARTICULAR PLATE MODEL (ROTATION POLE): C THAT OF DEMETS ET. AL. (1990): NUVEL-1. C C YOU WILL PROBABLY NEED TO REPLACE THE CODE GIVEN HERE WITH C NEW CODE OF YOUR OWN ]]] C**************************************************************** C PARAMETER (LAST=17) LOGICAL CONTAC,RESIST C XYARC HOLDS PAIRS OF (X, Y) OF VOLCANIC ARC, WITH X INCREASING: DIMENSION XYARC(2,LAST) DATA ((XYARC(I,J),I=1,2),J=1,17) / + -2.688E+06,+2.382E+05, + -2.564E+06,-9.064E+04, + -2.343E+06,-3.474E+05, + -2.088E+06,-5.215E+05, + -1.682E+06,-6.322E+05, + -1.420E+06,-6.290E+05, + -1.064E+06,-5.516E+05, + -8.054E+05,-4.764E+05, + -6.004E+05,-3.647E+05, + -3.868E+05,-2.185E+05, + -2.812E+05,+5.439E+04, + -2.080E+05,+2.306E+05, + -1.141E+05,+3.241E+05, + +7.915E+04,+1.160E+05, + +2.325E+05,+4.720E+04, + +4.203E+05,-7.314E+04, + +6.269E+05,-2.666E+05/ C IF (TRHMAX.LE.0.) THEN C NO-DRAG OPTION: RESIST=.FALSE. VMX=VCX VMY=VCY ELSE IF ((X.LT.XYARC(1,1)).OR.(X.GT.XYARC(1,LAST)).OR. + (X.GT.(-90000.-Y*1.28)).OR. + (X.GT.(287000.+Y*0.860))) THEN CONTAC=.FALSE. ELSE DO 10 I=2,LAST IF ((X.GE.XYARC(1,I-1)).AND. + (X.LE.XYARC(1,I))) THEN I1=I-1 I2=I FRAC=(X-XYARC(1,I1))/ + (XYARC(1,I2)-XYARC(1,I1)) GO TO 11 ENDIF 10 CONTINUE 11 YARC=XYARC(2,I1)+FRAC*(XYARC(2,I2)-XYARC(2,I1)) CONTAC=Y.LT.YARC ENDIF IF (CONTAC) THEN RESIST=.TRUE. CALL SIDES (INPUT,X,Y, + OUTPUT,VMX,VMY) ELSE RESIST=.FALSE. VMX=VCX VMY=VCY ENDIF ENDIF RETURN END C C C SUBROUTINE SIDES (INPUT,X,Y, + OUTPUT,VX,VY) C C COMPUTES HORIZONTAL COMPONENTS OF FLOW AT SIDES OF MODEL, IF C NEEDED FOR IMPOSITION OF TYPE-3 VELOCITY BOUNDARY CONDITIONS. C C**************************************************************** C CAVEAT HACKER ]]] C UNLIKE OTHER SUBPROGRAMS IN THIS PACKAGE, "SIDES" IS VERY C SPECIFIC TO A PARTICULAR PROBLEM: C -IT ONLY DESCRIBES THE PLATE BOUNDARIES IN THE ALASKA/BERING SEA C REGION. C -IT ASSUMES A PARTICULAR ORIGIN AND ORIENTATION OF THE X-AXIS. C (ORIGIN AT 61 N, 147 W, WITH +X POINTING E) C -IT IS BASED ON A PARTICULAR PLATE MODELS (ROTATION POLES): C THOSE OF DEMETS ET. AL. (1990): NUVEL-1, AND C COOK ET AL. (1986): OKHOTSK/NORTH AMERICA. C -NOTE THAT ALL VELOCITIES ARE RELATIVE TO NORTH AMERICA. C C YOU WILL PROBABLY NEED TO REPLACE THE CODE GIVEN HERE WITH C NEW CODE OF YOUR OWN ]]] C**************************************************************** C LOGICAL PACIFI,EURASI,OKHOTS REAL N1,N2,N3 C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C IF (Y.LT. +5.617E5) THEN PACIFI=.TRUE. EURASI=.FALSE. OKHOTS=.FALSE. ELSE DX=X-(-1.200E6) DY=Y-(+5.617E5) ARG=ATAN2F(DY,DX) IF (ARG.GT.2.3999) THEN OKHOTS=.TRUE. EURASI=.FALSE. PACIFI=.FALSE. ELSE EURASI=.TRUE. OKHOTS=.FALSE. PACIFI=.FALSE. ENDIF ENDIF C RADIUS=6371000. CPNLAT=61. Y0NLAT=61. X0ELON=-147. C RTAN=RADIUS*TANDEG(90.-CPNLAT) YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) CALL XYTOLL (INPUT,CPNLAT,RADIUS,RTAN,X0ELON,YPOLE, + X,Y, + OUTPUT,PLAT,PLON) R1=RADIUS*COSDEG(PLON)*COSDEG(PLAT) R2=RADIUS*SINDEG(PLON)*COSDEG(PLAT) R3=RADIUS*SINDEG(PLAT) N1= -COSDEG(PLON)*SINDEG(PLAT) N2= -SINDEG(PLON)*SINDEG(PLAT) N3=COSDEG(PLAT) E1= -SINDEG(PLON) E2=COSDEG(PLON) E3=0. C IF (PACIFI) THEN C (DEMETS ET AL.: C (48.7S, 101.8E, 0.78 DEG./MILLION-YEARS) RATE=SINDEG(0.78)/(1.E6*3.15576E7) P1=COSDEG(101.8)*COSDEG(-48.7)*RATE P2=SINDEG(101.8)*COSDEG(-48.7)*RATE P3=SINDEG(-48.7)*RATE ELSE IF (EURASI) THEN C (DEMETS ET AL.: C (62.4N, 135.8E, 0.22 DEG./MILLION-YEARS) RATE=SINDEG(0.22)/(1.E6*3.15576E7) P1=COSDEG(135.8)*COSDEG(+62.4)*RATE P2=SINDEG(135.8)*COSDEG(+62.4)*RATE P3=SINDEG(+62.4)*RATE ELSE IF (OKHOTS) THEN C (COOK ET AL., 1986; RATE IS REALLY A GUESS]) C (72.4N, 169.8E, 0.4 DEG./MILLION-YEARS) RATE=SINDEG(0.4 )/(1.E6*3.15576E7) P1=COSDEG(169.8)*COSDEG(+72.4)*RATE P2=SINDEG(169.8)*COSDEG(+72.4)*RATE P3=SINDEG(+72.4)*RATE ELSE VX=0. VY=0. RETURN ENDIF V1=P2*R3-P3*R2 V2=P3*R1-P1*R3 V3=P1*R2-P2*R1 VN=N1*V1+N2*V2+N3*V3 VE=E1*V1+E2*V2+E3*V3 VSQ=VN*VN+VE*VE IF (VSQ.GT.0.) THEN TIME=RADIUS*SINDEG(1.)/SQRT(VSQ) ELSE VX=0. VY=0. RETURN ENDIF QLAT=PLAT+57.296*VN*TIME/RADIUS QLON=PLON+57.296*VE*TIME/(RADIUS*COSDEG(PLAT)) CALL LLTOXY (INPUT,CPNLAT, + QLAT,QLON, + RADIUS,RTAN,X0ELON,YPOLE, + OUTPUT,XP,YP) VX=(XP-X)/TIME VY=(YP-Y)/TIME RETURN C****** ALASKA (END SECTION) ****************** END