C-------------------------------------------------------------------- C PROGRAM ORBDATA C written in Fortran 77 (fixed-column) syntax, C but the need for dynamic memory allocation (ALLOCATE) C for the (up to) 3 huge arrays of elevation, C heat-flow and (possibly) seafloor ages C has now forced this code into Fortran 90. C Therefore, it now has ALLOCATE in main, C REAL, DIMENSION(:,:) in both main and ASSIGN, C an interface block for ASSIGN, and matching C Fortran 90 parameter declarations in ASSIGN. C C BY PETER BIRD, DEPT. OF EARTH & SPACE SCIENCES, C UNIVERSITY OF CALIFORNIA, LOS ANGELES, CA 90095-1567. C (FOR VERSION DATE, SEE 1st WRITE BELOW.) C C READS A FINITE-ELEMENT GRID PRODUCED BY ORBWEAVE C (IN EITHER SHELLS OR RESTORE3 MODE) C AND A PARAMETER FILE (FORMATTED FOR SHELLS) WITH C THERMAL CONDUCTIVITIES, DENSITIES, ET CETERA; C THEN, FILLS IN OR ADJUSTS NODAL DATA: C -ELEVATION (+ ABOVE SEA LEVEL, - BELOW); C -HEAT-FLOW (ASSUMING STEADY-STATE GEOTHERM); C -CRUSTAL THICKNESS; C -MANTLE LITHOSPHERE THICKNESS (NOT INCLUDING CRUST). C C ALL UNITS ARE SI: m, W/m**2 (not mW/m**2), m, m. C C IN THE OUTPUT DATA, THE BASE OF THE MANTLE LITHOSPHERE C WILL BE AN ISOTHERMAL SURFACE. THE ISOTHERM VALUE C IS CHOSEN TO LIE ON THE ASTHENOSPHERE ADIABAT C (EVALUATED AT AN ARBITRARY DEPTH OF 100 km). C C IF A NODAL ELEVATION IS NON-ZERO, IT WILL ALWAYS BE LEFT C UNCHANGED, TO PRESERVE EFFECTS OF HAND-EDITING. C IF ANY ELEVATION(S) IS/ARE ZERO, NEW VALUE(S) WILL BE C INTERPOLATED FROM A GRIDDED DATASET. C C IF A NODAL HEAT-FLOW IS NON-ZERO, IT WILL ALWAYS BE LEFT C UNCHANGED, TO PRESERVE EFFECTS OF HAND-EDITING. C IF ANY HEAT-FLOW(S) IS/ARE ZERO, NEW VALUE(S) WILL BE C INTERPOLATED FROM A GRIDDED DATASET. C C THERE IS AN OPTION TO BASE HEAT-FLOW ON AGE OF SEAFLOOR, C WHERE KNOWN, FROM A GRIDDED-AGE DATASET. C C HEAT-FLOW (FROM ANY SOURCE) IS SUBJECT TO A MINIMUM VALUE C DESIGNED TO GUARANTEE EXISTENCE OF A SOLUTION. C C FINALLY, ISOSTASY WITH RESPECT TO A SPREADING RISE C IS ACHIEVED BY ADJUSTING THE LAYER THICKNESSES, C LEAVING ELEVATION AND HEAT-FLOW ALONE. C--------------------------------------------------------------------- INTERFACE SUBROUTINE ASSIGN (INPUT,AARRAY,AX1,ADX,AX2,NAX,ADY,AY2,NAY, + ALPHAT,CLIMIT,CONDUC, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT, + ONEKM, + PLON,PLAT, + QLIM0,DQL_DE,QLIM1, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + TASTHK,TEMLIM,TSURF,USEAGE, + MODIFY,ELEVAT,HEATFL,THICKC,THICKM) INTEGER, INTENT(IN) :: INPUT,IUNITL,IUNITT,MODIFY,NAX,NAY, & NEX,NEY,NQX,NQY LOGICAL, INTENT(IN) :: USEAGE REAL, INTENT(IN) :: AX1,ADX,ADY,AX2,AY2,CLIMIT,EX1,EDX,EX2, & EDY,EY2,GMEAN,HCMAX,HLMAX,ONEKM,PLON,PLAT, & QLIM0,DQL_DE,QLIM1,QX1,QDX,QX2,QDY,QY2, & RHOAST,RHOH2O,TASTHK,TSURF REAL, DIMENSION(2), INTENT(IN) :: ALPHAT,CONDUC,RADIO,RHOBAR, & TEMLIM REAL, DIMENSION(:,:), INTENT(IN) :: AARRAY,EARRAY,QARRAY REAL, INTENT(INOUT) :: ELEVAT,HEATFL,THICKC,THICKM END SUBROUTINE ASSIGN END INTERFACE C--------------------------------------------------------------------- C PARAMETER (ARRAY-SIZE) STATEMENTS C C SET THE FOLLOWING PARAMETERS AT LEAST AS LARGE AS YOUR PROBLEM: C C MAXNOD = MAXIMUM NUMBER OF NODES (INCLUDES BOTH "REAL"AND & "FAKE") PARAMETER (MAXNOD=20000) C C MAXBN = MAXIMUM NUMBER OF BOUNDARY NODES (BOTH "REAL" AND "FAKE"). PARAMETER (MAXBN=6000) C C MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES). PARAMETER (MAXEL=30000) C C MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS); PARAMETER (MAXFEL=6000) C C MAXATP = MAXIMUM NUMBER OF NODES WHICH MAY OVERLAP AT A FAULT- C INTERSECTION POINT. PARAMETER (MAXATP=10) C C PARAMETER GIVING THE (EXACT) NUMBER OF PB2002 PLATES: PARAMETER (NUMPLT=52) C C--------------------------------------------------------------------- C TYPE STATEMENTS C C (NOTE: THE IMPLICIT TYPING OF I-N = INTEGER, AND A-H,O-Z = REAL C IS OBSERVED IN THIS PROGRAM. NO TYPES ARE STATED FOR SUCH NAMES.) C CHARACTER*1 C1 CHARACTER*2 NAMES CHARACTER*80 TITLE1,TITLE3 C LOGICAL BRIEF,DOTHIS,EVERYP,LLIMIT,NEEDE,NEEDQ C REAL LRANGE,NLAT,NLAT0,NLAT1 C C NOTE: THE FOLLOWING ARRAYS COULD BE COMPRESSED WITH "LOGICAL*1" C IN VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN,EDGETS,EDGEFS,USEAGE C C THE FOLLOWING TO AGREE WITH BLOCK DATA BD1: DOUBLE PRECISION POINTS,WEIGHT C THE FOLLOWING TO AGREE WITH BLOCK DATA BD2: DOUBLE PRECISION FPHI,FPOINT,FGAUSS C C--------------------------------------------------------------------- C DIMENSION STATMENTS C C DIMENSIONs using dynamic memory allocation: REAL, DIMENSION(:,:), ALLOCATABLE :: EARRAY,QARRAY,AARRAY C C DIMENSIONS USING PARAMETER MAXNOD: DIMENSION CHECKN(MAXNOD), DQDTDA(MAXNOD), + ELEV (MAXNOD), + TLNODE(MAXNOD), + XNODE (MAXNOD), YNODE (MAXNOD), ZMNODE(MAXNOD) C C DIMENSIONS USING PARAMETER MAXBN: DIMENSION NODCON (MAXBN) C C DIMENSIONS USING PARAMETER MAXEL: DIMENSION AREA (MAXEL), CHECKE (MAXEL), + DETJ (7,MAXEL), + DXS(2,2,3,7,MAXEL), DYS (2,2,3,7,MAXEL), + DXSP (3,7,MAXEL), DYSP (3,7,MAXEL), + EDGETS (3,MAXEL), ELEDAT (3,MAXEL), + FPSFER(2,2,3,7,MAXEL), + NODES (3,MAXEL), + SITA (7,MAXEL) C C DIMENSIONS USING PARAMETER MAXFEL: DIMENSION CHECKF (MAXFEL), EDGEFS (2,MAXFEL), + FDIP (2,MAXFEL), + FLEN (MAXFEL), + FPFLT(2,2,2,7,MAXFEL), + FARG (2,MAXFEL), NODEF (4,MAXFEL), + OFFSET (MAXFEL) C C DIMENSIONS USING PARAMETER MAXATP: DIMENSION LIST (MAXATP) C C DIMENSIONS USING PARAMETER NUMPLT: DIMENSION NAMES (NUMPLT) C C FIXED DIMENSIONS: DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),RADIO(2), + RHOBAR(2),TAUMAX(2),TEMLIM(2) C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: DIMENSION POINTS(3,7),WEIGHT(7) C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: DIMENSION FPHI(4,7),FPOINT(7),FGAUSS(7) C C--------------------------------------------------------------------- C COMMON STATEMENTS C C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: COMMON /S1S2S3/ POINTS COMMON /WGTVEC/ WEIGHT C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS C------------------------------------------------------------------- C DATA STATEMENTS C C "QLIM0" is the lower limit on heat-flow for points C with an elevation of zero. C (A minimum heat-flow limit is necessary to find C a converged solution.) DATA QLIM0 /0.045/ C C "DQL_DE" is the derivitive d(QLIMIT)/d(elevation), C which adjusts the minimum heat-flow for elevation. C (A minimum heat-flow limit is necessary to find C a converged solution.) DATA DQL_DE /1.43E-06/ C C "QLIM1 is the upper limit on heat-flow for all points. DATA QLIM1 /0.300/ C C "CLIMIT" IS A LOWER LIMIT ON CRUSTAL THICKNESS: DATA CLIMIT /6570./ C Changed from 5000. to 6570. on 2005.05.24, to agree with CRUST2.grd. C C "HCMAX" IS AN UPPER LIMIT ON CRUSTAL THICKNESS: DATA HCMAX /75000./ C Already agreed with CRUST2.grd; no need to change! C C "HLMAX" IS AN UPPER LIMIT ON TOTAL LITHOSPHERE THICKNESS: DATA HLMAX /300000./ C C (NOTE THAT ALL OF THE ABOVE ARE IN SI UNITS (W/m**2, m, m, m). C C "IUNITT" = FORTRAN DEVICE NUMBER FOR TEXT MESSAGES: DATA IUNITT /6/ C C "IUNITL" = FORTRAN DEVICE NUMBER FOR -ASSIGN- LOG FILE: DATA IUNITL /11/ C plate names (in alphabetical order): DATA names / 'AF','AM','AN', + 'AP','AR','AS', + 'AT','AU','BH', + 'BR','BS','BU', + 'CA','CL','CO', + 'CR','EA','EU', + 'FT','GP','IN', + 'JF','JZ','KE', + 'MA','MN','MO', + 'MS','NA','NB', + 'ND','NH','NI', + 'NZ','OK','ON', + 'PA','PM','PS', + 'RI','SA','SB', + 'SC','SL','SO', + 'SS','SU','SW', + 'TI','TO','WL', + 'YA' / 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 C STATEMENT FUNCTIONS: EASTOF(X)=MOD((X + 720.),360.) C C--------------------------------------------------------------------- C C BEGINNING OF 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 MXBN =MAXBN MXEL =MAXEL MXFEL =MAXFEL MXSTAR=MAXATP C WEDGE=ABS(90.-ABS(DIPMAX))*0.017453293 C--------------------------------------------------------------------- C C INTRODUCTION / HEADER SECTION: C WRITE (IUNITT,1) 1 FORMAT ( + ' -------------------------------------------------------------' + /' Program OrbData (version of 23 July 2001,' + /' except CLIMIT increased to 6570 m on 24 May 2005)' +//' by Peter Bird, Dept. of Earth & Space Sciences,' + /' University of California, Los Angeles, CA 90095-1567.' +//' Reads a finite element grid produced by OrbWeave' + /' (in either Shells or Restore3 mode)' + /' and a parameter file (formatted for Shells) with' + /' thermal conductivities, densities, et cetera;' + /' Then, fills in or adjusts nodal data:' + /' -elevation (+ above sea level, - below);' + /' -heat-flow (assuming steady-state geotherm);' + /' -crustal thickness;' + /' -mantle lithosphere thickness (NOT including crust).' +//' All units are SI: m, W/m**2 (NOT mW/m**2), m, m.' +//' In the output data, the base of the mantle lithosphere' + /' will be an isothermal surface. The isotherm value' + /' is chosen to lie on the asthenosphere adiabat' + /' (evaluated at an arbitrary depth of 100 km).' +//' PRESS [Enter] to continue...'\) READ(*,"(A)") C1 WRITE (IUNITT,2) QLIM0, DQL_DE, QLIM1 2 FORMAT ( + /' If an elevation is non-zero, this elevation will be left' + /' unchanged, to preserve effects of hand-editing.' + /' If any elevation is zero, a new values will be' + /' interpolated from a gridded dataset.' + /' If a nodal heat-flow is non-zero, this heat-flow will be left' + /' unchanged, to preserve effects of hand-editing.' + /' If any heat-flow is zero, a new values will be' + /' interpolated from a gridded dataset.' + /' There is an option to base heat-flow on age of seafloor,' + /' where known (< 200 Ma), from a gridded-age dataset.' +//' Heat-flow (from any source) is subject to a minimum value of' + /' (',F5.3,' + ',ES10.3,' * elevation)' + /' and a maximum value of ',F5.3 + /' chosen to guaruntee existence of a solution.' +//' Finally, isostasy with respect to a spreading rise' + /' is achieved by adjusting the layer thicknesses,' + /' leaving elevation and heat-flow alone.' + /' ------------------------------------------------------------' +//' PRESS [Enter] to continue...'\) READ(*,"(A)") C1 WRITE(*,3) 3 FORMAT(/ /' It is possible to limit the update of nodal values' + /' to a rectangular region of (longitude,latitude)' + /' Do you wish to use this option? (T/F): '\) READ(*,*) LLIMIT IF (LLIMIT) THEN WRITE(*,"(' Enter West-side longitude (East = + ): '\)") READ(*,*) WELON WRITE(*,"(' Enter East-side longitude (East = + ): '\)") READ(*,*) EELON STELON=EASTOF(WELON) LRANGE=EASTOF(EELON - WELON) 37 WRITE(*,"(' Enter Southernmost latitude (North = +): '\)") READ(*,*) NLAT0 WRITE(*,"(' Enter Northernmost latitude (North = +): '\)") READ(*,*) NLAT1 IF (NLAT1.LE.NLAT0) THEN WRITE (*,"(' ERROR: Maximum latitude must ', + 'be >= minimum.')") GO TO 37 END IF END IF C WRITE (IUNITT,5) QLIM0,DQL_DE, QLIM1, CLIMIT, HCMAX, HLMAX 5 FORMAT(/' THE FOLLOWING LIMITS APPLY IN THIS RUN:' +/' LOWER LIMIT ON HEAT-FLOW = ',F5.3,'+',ES10.3,' * elevation' +/' UPPER LIMIT ON HEAT-FLOW = ',F5.3 +/' LOWER LIMIT ON CRUSTAL THICKNESS = ',F7.0 +/' UPPER LIMIT ON CRUSTAL THICKNESS = ',F7.0 +/' UPPER LIMIT ON TOTAL LITHOSPHERE THICKNESS = ',F7.0) C C READ PARAMETERS ON UNIT 1 C CALL READPM (INPUT, 1, 6, NAMES, NUMPLT, OFFMAX, + OUTPUT,ACREEP, ALPHAT, BCREEP, BIOT , + BYERLY, CCREEP, CFRIC , CONDUC, + DCREEP, ECREEP, EVERYP, FFRIC , GMEAN , + GRADIE, ICONVE, IPVREF, + MAXITR, OKDELV, OKTOQT, ONEKM, RADIO , + RDETH, REFSTR, RHOAST, RHOBAR, RHOH2O, + TADIAB, + TAUMAX, TEMLIM, TITLE3, TRHMAX, TSURF, + ZBASTH) C C LITHOSPHERE/ASTHENOSPHERE TEMPERATURE, IN KELVIN, IS DETERMINED C FROM THE ASTHENOSPHERE ADIABAT IN THE PARAMETER FILE, C EVALUATED (ARBITRARILY) AT 100 KM: C TASTHK = TADIAB + GRADIE * 100.E3 C C READ FINITE ELEMENT GRID ON UNIT 2 C CALL GETNET (INPUT, 2,IUNITT, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + NFAKEN,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,OFFMAX,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE, + WORK,CHECKE,CHECKF,CHECKN) WRITE (IUNITT,40) 40 FORMAT (/' Successfully read FE grid, now verifying topology...') CALL SQUARE (INPUT,BRIEF,FDIP,IUNITT, + MXBN,MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES, + NUMEL,NUMNOD,RDETH,WEDGE, + MODIFY,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,EDGEFS, + EDGETS,FLEN,FPFLT,FPSFER, + FARG,NCOND,NODCON,SITA, + WORK,CHECKN,LIST) WRITE (IUNITT,50) 50 FORMAT (/' Topology of finite element grid verified.'/) C C READ IN ELEVATION ARRAY ON UNIT 3, IF NEEDED C NEEDE=.FALSE. DO 60 I=1,NUMNOD IF (LLIMIT) THEN ELON=YNODE(I)*57.2957795 NLAT=90.-XNODE(I)*57.2957795 IF ((EASTOF(ELON-WELON).LE.LRANGE).AND. + (NLAT.GE.NLAT0).AND. + (NLAT.LE.NLAT1)) THEN IF (ELEV(I).EQ.0.0) NEEDE=.TRUE. END IF ELSE IF (ELEV(I).EQ.0.0) NEEDE=.TRUE. END IF 60 CONTINUE IF (NEEDE) THEN WRITE(*,61) 61 FORMAT(/ /' ATTEMPTING TO READ GRIDDED ELEVATIONS:'/) READ (3,*) EX1,EDX,EX2 READ (3,*) EY1,EDY,EY2 NEX=(EX2-EX1)/EDX + 1.5 NEY=(EY2-EY1)/EDY + 1.5 ALLOCATE ( EARRAY(NEY, NEX) ) READ (3,*) ((EARRAY(IROW,JCOL),JCOL=1,NEX),IROW=1,NEY) ELSE WRITE (*,69) 69 FORMAT (/' ALL NODES HAVE NON-ZERO ELEVATION.' & /' NO INPUT GRID IS NEEDED.') ENDIF C C READ IN HEAT-FLOW ARRAY ON UNIT 4, IF NEEDED C NEEDQ=.FALSE. DO 70 I=1,NUMNOD IF (LLIMIT) THEN ELON=YNODE(I)*57.2957795 NLAT=90.-XNODE(I)*57.2957795 IF ((EASTOF(ELON-WELON).LE.LRANGE).AND. + (NLAT.GE.NLAT0).AND. + (NLAT.LE.NLAT1)) THEN IF (DQDTDA(I).EQ.0.0) NEEDQ=.TRUE. END IF ELSE IF (DQDTDA(I).EQ.0.0) NEEDQ=.TRUE. END IF 70 CONTINUE IF (NEEDQ) THEN WRITE(*,71) 71 FORMAT(/ /' ATTEMPTING TO READ GRIDDED HEAT-FLOW:'/) READ (4,*) QX1,QDX,QX2 READ (4,*) QY1,QDY,QY2 NQX=(QX2-QX1)/QDX + 1.5 NQY=(QY2-QY1)/QDY + 1.5 ALLOCATE ( QARRAY(NQY, NQX) ) READ (4,*) ((QARRAY(IROW,JCOL),JCOL=1,NQX),IROW=1,NQY) C C CHECK IF GRID OF SEA-FLOOR AGES WILL BE USED? WRITE (*,72) 72 FORMAT(/' DO YOU WISH TO USE ' & ,'A GRID OF SEA-FLOOR AGES (T/F)?') READ (*,*) USEAGE IF (USEAGE) THEN C C READ DATASET OF GRIDDED AGES, ON UNIT 7: C (NOTE: AGE > 200 Ma MEANS "UNKNOWN" OR "CONTINENTAL") WRITE(*,73) 73 FORMAT(/ /' ATTEMPTING TO READ GRIDDED AGES:'/) READ (7,*) AX1,ADX,AX2 READ (7,*) AY1,ADY,AY2 NAX=(AX2-AX1)/ADX + 1.5 NAY=(AY2-AY1)/ADY + 1.5 ALLOCATE ( AARRAY(NAY, NAX) ) READ (7,*) ((AARRAY(IROW,JCOL),JCOL=1,NAX),IROW=1,NAY) ELSE C define array size to keep compiler happy: ALLOCATE ( AARRAY(1, 1) ) C ENDIF C ELSE C WRITE (*,79) 79 FORMAT (/' ALL NODES HAVE NON-ZERO HEAT-FLOW.' & /' NO INPUT DATA SETS ARE NEEDED.') C C IMPOSE HEAT-FLOW LIMITS AT NODES, TO PREVENT UNREASONABLY C THICK AND STIFF LITHOSPHERE ANYWHERE C DO 80 I=1,NUMNOD QLIMIT=QLIM0+DQL_DE*ELEV(I) IF (DQDTDA(I).NE.0.0) DQDTDA(I)=MAX(DQDTDA(I),QLIMIT) DQDTDA(I)=MIN(DQDTDA(I),QLIM1) 80 CONTINUE C define array size to keep compiler happy: ALLOCATE ( AARRAY(1, 1) ) END IF C C INITIATE LOG FILE FOR -ASSIGN- OUTPUT: C WRITE (IUNITT,91) IUNITL 91 FORMAT (/ /' ATTEMPTING TO CREATE DETAILED LOG FILE ON UNIT', + I3/) WRITE (IUNITL,92) 92 FORMAT ('DETAILED MESSAGES FROM SUBPROGRAM -ASSIGN-:') C C PROCESS ALL NODES EQUALLY: C WRITE (IUNITT,600) 600 FORMAT (/' COMPUTING LAYER THICKNESSES AT ALL NODES...' + /' 0 nodes completed...') WRITE (IUNITL,601) 601 FORMAT (' NODE ELEV DQDTDA ZMNODE TLNODE') DO 680 I=1,NUMNOD PLON=YNODE(I)*57.2957795 PLAT=90.-XNODE(I)*57.2957795 ELEVAT=ELEV(I) HEATFL=DQDTDA(I) IF (LLIMIT) THEN DOTHIS=((EASTOF(PLON-WELON).LE.LRANGE).AND. + (PLAT.GE.NLAT0).AND. + (PLAT.LE.NLAT1)) ELSE DOTHIS=.TRUE. END IF IF (DOTHIS) THEN CALL ASSIGN (INPUT,AARRAY,AX1,ADX,AX2,NAX,ADY,AY2,NAY, + ALPHAT,CLIMIT,CONDUC, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT, + ONEKM, + PLON,PLAT, + QLIM0,DQL_DE,QLIM1, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + TASTHK,TEMLIM,TSURF,USEAGE, + MODIFY,ELEVAT,HEATFL,THICKC,THICKM) ELEV(I)=ELEVAT DQDTDA(I)=HEATFL ZMNODE(I)=MAX(THICKC,CLIMIT) TLNODE(I)=MAX(THICKM,0.) WRITE (IUNITL,678) I,ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 678 FORMAT (' ',I5,4ES10.2) WRITE (IUNITT,679) I 679 FORMAT ('+',I6,' nodes completed...') END IF 680 CONTINUE WRITE (IUNITT,700) 700 FORMAT (/ /' CHECK THE LOG FILE CAREFULLY FOR PROBLEMS!') C C OUTPUT THE MODIFIED .FEG FILE: C CALL PUTNET (INPUT, 12, + BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE) STOP 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 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) CALL PAUSE() 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 SQUEEZ (INPUT,ALPHAT,ELEVAT, + GEOTH1,GEOTH2,GEOTH3,GEOTH4, + GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN, + IUNITT,ONEKM,RHOAST,RHOBAR,RHOH2O, + TEMLIM,ZM,ZSTOP, + OUTPUT,TAUZZ,SIGZZB) C C CALCULATES "TAUZZ", THE VERTICAL INTEGRAL THROUGH THE PLATE C OF THE VERTICAL STANDARDIZED STRESS ANOMALY, WHICH IS C RELATIVE TO A COLUMN OF MANTLE AT ASTHENOSPHERE TEMPERATURE C WITH A 5 KM CRUST AND A 2.7 KM OCEAN ON TOP, LIKE A MID-OCEAN C RISE. THE INTEGRAL IS FROM EITHER THE LAND SURFACE OR THE C SEA SURFACE, DOWN TO A DEPTH OF "ZSTOP" BELOW THE TOP OF C THE CRUST. C IF "ZSTOP" EXCEEDS MOHO DEPTH "ZM", THEN PROPERTIES OF THE MANTLE C WILL BE USED IN THE LOWER PART OF THE INTEGRAL. C ALSO RETURNS "SIGZZB", THE STANDARDIZED VERTICAL STRESS ANOMALY C AT DEPTH "ZSTOP" BELOW THE SOLID ROCK SURFACE. C NOTE: THIS VERSION IS DIFFERENT FROM THE VERSION FOUND IN THE LARAMY C PROGRAM PACKAGE. FIRST, IT ACTS ON ONLY A SINGLE POINT. C SECOND, IT INFERS SUB-PLATE NORMAL-STRESS ANOMALIES FROM C THE GIVEN TOPOGRAPHY, INSTEAD OF FROM MODEL STRUCTURE. C PARAMETER (NDREF=700) LOGICAL CALLED C INTERNAL ARRAYS: DIMENSION DREF(NDREF),PREF(0:NDREF) C ARGUMENT ARRAYS: DIMENSION ALPHAT(2),RHOBAR(2),TEMLIM(2) SAVE CALLED,DREF,PREF DATA CALLED /.FALSE./ C C STATEMENT FUNCTIONS: TEMPC(H)=MIN(TEMLIM(1),GEOTH1+GEOTH2*H+GEOTH3*H**2 + +GEOTH4*H**3) TEMPM(H)=MIN(TEMLIM(2),GEOTH5+GEOTH6*H+GEOTH7*H**2 + +GEOTH8*H**3) C C CREATE REFERENCE TEMPERATURE & DENSITY PROFILES TO DEPTH OF NDREF KM C IF (.NOT.CALLED) THEN RHOTOP=RHOBAR(1)*(1.-ALPHAT(1)*GEOTH1) DREF(1)=RHOH2O DREF(2)=RHOH2O DREF(3)=0.7*RHOH2O+0.3*RHOTOP DREF(4)=RHOTOP DREF(5)=RHOTOP DREF(6)=RHOTOP DREF(7)=RHOTOP DREF(8)=0.7*RHOTOP+0.3*RHOAST DO 50 J=9,NDREF DREF(J)=RHOAST 50 CONTINUE PREF(0)=0. DO 100 I=1,NDREF PREF(I)=PREF(I-1)+DREF(I)*GMEAN*ONEKM 100 CONTINUE ENDIF C C ROUTINE PROCESSING (ON EVERY CALL): C IF (ELEVAT.GT.0.) THEN C LAND ZTOP= -ELEVAT ZBASE=ZSTOP-ELEVAT DENSE1=RHOBAR(1)*(1.-GEOTH1*ALPHAT(1)) H=0. LAYER1=1 ELSE C OCEAN ZTOP=0. ZBASE=ZSTOP+(-ELEVAT) DENSE1=RHOH2O H=ELEVAT LAYER1=0 ENDIF LASTDR=ZBASE/ONEKM IF (ZBASE.GT.ONEKM*LASTDR) LASTDR=LASTDR+1 IF (LASTDR.GT.NDREF) THEN WRITE(IUNITT,110) LASTDR 110 FORMAT(' IN SUBPROGRAM SQUEEZ, PARAMETER NDREF '/ + ' MUST BE INCREASED TO AT LEAST ',I10) CALL PAUSE() STOP ENDIF NSTEP=(ZBASE-ZTOP)/ONEKM OLDSZZ=0. OLDPR=0. SIGZZ=0. TAUZZ=0. Z=ZTOP DO 200 I=1,NSTEP Z=Z+ONEKM H=H+ONEKM IF (H.GT.0.) THEN IF (H.LE.ZM) THEN T=TEMPC(H) DENSE2=RHOBAR(1)*(1.-T*ALPHAT(1)) LAYER2=1 ELSE T=TEMPM(H-ZM) DENSE2=RHOBAR(2)*(1.-T*ALPHAT(2)) LAYER2=2 ENDIF ELSE DENSE2=RHOH2O LAYER2=0 ENDIF IF ((LAYER1.EQ.0).AND.(LAYER2.EQ.1)) THEN FRAC2=H/ONEKM FRAC1=1.-FRAC2 ELSE IF ((LAYER1.EQ.1).AND.(LAYER2.EQ.2)) THEN FRAC2=(H-ZM)/ONEKM FRAC1=1.-FRAC2 ELSE FRAC1=0.5 FRAC2=0.5 ENDIF DENSE=FRAC1*DENSE1+FRAC2*DENSE2 IF (Z.GT.0.) THEN N1=Z/ONEKM N2=N1+1 FRAC=Z/ONEKM-N1 PR=PREF(N1)+FRAC*(PREF(N2)-PREF(N1)) ELSE PR=0. ENDIF SIGZZ=SIGZZ-DENSE*GMEAN*ONEKM+(PR-OLDPR) TAUZZ=TAUZZ+0.5*(SIGZZ+OLDSZZ)*ONEKM DENSE1=DENSE2 OLDSZZ=SIGZZ OLDPR=PR LAYER1=LAYER2 200 CONTINUE RESID=ZBASE-Z H=ZSTOP Z=ZBASE IF (ZSTOP.LE.ZM) THEN T=TEMPC(H) DENSE2=RHOBAR(1)*(1.-T*ALPHAT(1)) ELSE T=TEMPM(H-ZM) DENSE2=RHOBAR(2)*(1.-T*ALPHAT(2)) ENDIF DENSE=0.5*(DENSE1+DENSE2) IF (Z.GT.0.) THEN N1=Z/ONEKM N2=N1+1 FRAC=Z/ONEKM-N1 PR=PREF(N1)+FRAC*(PREF(N2)-PREF(N1)) ELSE PR=0. ENDIF SIGZZB=SIGZZ-DENSE*GMEAN*RESID+(PR-OLDPR) TAUZZ=TAUZZ+0.5*(SIGZZB+OLDSZZ)*RESID CALLED=.TRUE. RETURN END C C C SUBROUTINE ASSIGN (INPUT,AARRAY,AX1,ADX,AX2,NAX,ADY,AY2,NAY, + ALPHAT,CLIMIT,CONDUC, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT, + ONEKM, + PLON,PLAT, + QLIM0,DQL_DE,QLIM1, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + TASTHK,TEMLIM,TSURF,USEAGE, + MODIFY,ELEVAT,HEATFL,THICKC,THICKM) C C DETERMINES ELEVATION (IF 0.0), HEAT FLOW (IF 0.0), AND ISOSTATIC C VALUES OF CRUSTAL AND MANTLE LITHOSPHERE THICKNESS, AT ONE POINT: C (PLON = EAST LONGITUDE, PLAT = NORTH LATITUDE), BOTH IN DEGREES. C (UNLIKE SUBPROGRAM -ASSIGN- IN ALGRIDDR, THIS VERSION NEVER C DOES REGIONAL AVERAGING. NEITHER DOES IT INTERPOLATE BY C KRIGING; A SIMPLE INTERPOLATION OFF A RECTANGULAR GRID IS C USED INSTEAD.) C INTEGER, INTENT(IN) :: INPUT,IUNITL,IUNITT,MODIFY,NAX,NAY, & NEX,NEY,NQX,NQY LOGICAL, INTENT(IN) :: USEAGE REAL, INTENT(IN) :: AX1,ADX,ADY,AX2,AY2,CLIMIT,EX1,EDX,EX2, & EDY,EY2,GMEAN,HCMAX,HLMAX,ONEKM,PLON,PLAT, & QLIM0,DQL_DE,QLIM1,QX1,QDX,QX2,QDY,QY2, & RHOAST,RHOH2O,TASTHK,TSURF REAL, DIMENSION(2), INTENT(IN) :: ALPHAT,CONDUC,RADIO,RHOBAR, & TEMLIM REAL, DIMENSION(:,:), INTENT(IN) :: AARRAY,EARRAY,QARRAY REAL, INTENT(INOUT) :: ELEVAT,HEATFL,THICKC,THICKM C C NUMBER OF ITERATIONS PERMITTED: PARAMETER (MAXITR = 50) C DAMPING FACTOR (<1.) FOR STABILITY OF ITERATIONS: PARAMETER (DAMP = 0.200) C LOGICAL BADP,BADT,DIVERG,OUTSID,NEEDE,NEEDQ, + WARNC1, WARNC2, WARNM1, WARNL2, WAYOUT C LOCAL STORAGE FOR HISTORY OF ITERATION: DIMENSION HC(MAXITR),HM(MAXITR),TM(MAXITR),TA(MAXITR) C--------------------------------------------------------------------- OUTSID=.FALSE. WAYOUT=.FALSE. WARNC1=.FALSE. WARNC2=.FALSE. WARNM1=.FALSE. WARNL2=.FALSE. C C DETERMINE THE ELEVATIONS, IF NEEDED C NEEDE=(ELEVAT.EQ.0.0) IF (NEEDE) THEN IR1=((EY2-PLAT)/EDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NEY-1) IR2=IR1+1 FR=((EY2-EDY*(IR1-1))-PLAT)/EDY C IF ((PLON.LT.EX1).OR.(PLON.GT.EX2)) THEN IF (((PLON+360.).GE.EX1).AND. + ((PLON+360.).LE.EX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.EX1).AND. + ((PLON-360.).LE.EX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-EX1)/EDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NEX-1) IC2=IC1+1 FC=(USELON-(EX1+EDX*(IC1-1)))/EDX C OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=EARRAY(IR1,IC1)+FC*(EARRAY(IR1,IC2)-EARRAY(IR1,IC1)) BOT=EARRAY(IR2,IC1)+FC*(EARRAY(IR2,IC2)-EARRAY(IR2,IC1)) ELEVAT=TOP+FR*(BOT-TOP) ENDIF C C DETERMINE HEAT FLOW, IF NEEDED C NEEDQ=(HEATFL.EQ.0.0) IF (NEEDQ) THEN IR1=((QY2-PLAT)/QDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NQY-1) IR2=IR1+1 FR=((QY2-QDY*(IR1-1))-PLAT)/QDY C IF ((PLON.LT.QX1).OR.(PLON.GT.QX2)) THEN IF (((PLON+360.).GE.QX1).AND. + ((PLON+360.).LE.QX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.QX1).AND. + ((PLON-360.).LE.QX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-QX1)/QDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NQX-1) IC2=IC1+1 FC=(USELON-(QX1+QDX*(IC1-1)))/QDX C OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=QARRAY(IR1,IC1)+FC*(QARRAY(IR1,IC2)-QARRAY(IR1,IC1)) BOT=QARRAY(IR2,IC1)+FC*(QARRAY(IR2,IC2)-QARRAY(IR2,IC1)) HEATFL=TOP+FR*(BOT-TOP) C C CHECK FOR SEAFLOOR AGE OVERRIDING HEAT-FLOW GRID VALUE: IF (USEAGE) THEN C IR1=((AY2-PLAT)/ADY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NAY-1) IR2=IR1+1 FR=((AY2-ADY*(IR1-1))-PLAT)/ADY C IF ((PLON.LT.AX1).OR.(PLON.GT.AX2)) THEN IF (((PLON+360.).GE.AX1).AND. + ((PLON+360.).LE.AX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.AX1).AND. + ((PLON-360.).LE.AX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-AX1)/ADX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NAX-1) IC2=IC1+1 FC=(USELON-(AX1+ADX*(IC1-1)))/ADX C OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=AARRAY(IR1,IC1)+FC*(AARRAY(IR1,IC2)-AARRAY(IR1,IC1)) BOT=AARRAY(IR2,IC1)+FC*(AARRAY(IR2,IC2)-AARRAY(IR2,IC1)) AGEMA=TOP+FR*(BOT-TOP) IF (AGEMA.LT.200.0) THEN C SEAFLOOR AGE IS VALID (NOT "UNKNOWN" OR "CONTINENTAL") IF (AGEMA.LE.0.0) THEN HEATFL=QLIM1 ELSE C Carol A. Stein & Seth Stein [1992] C A model for the global variation in oceanic C depth and heat flow with lithospheric age, C Nature, v. 359, 10 September, p. 123-129. C According to their preferred GDH1 model: IF (AGEMA.LE.55.0) THEN HEATFL=0.510/SQRT(AGEMA) ELSE HEATFL=0.048+0.096*EXP(-0.0278*AGEMA) END IF END IF END IF C END IF ENDIF C C GIVE WARNING IF ANY NECESSARY DATA COULD NOT BE READ C IF (OUTSID) THEN WRITE (IUNITL,1) PLON,PLAT 1 FORMAT(' WARNING: (',F8.3,'E, ',F7.3,'N) IS ', + 'SLIGHTLY OUTSIDE DATA GRID(S).') END IF IF (WAYOUT) THEN WRITE (IUNITL,2) PLON,PLAT WRITE (IUNITT,2) PLON,PLAT 2 FORMAT(' ERROR: (',F8.3,'E, ',F7.3,'N) IS ', + 'FAR OUTSIDE DATA GRID(S).') CALL PAUSE() CCCCC STOP END IF C C APPLY LIMITS ON Q TO GUARANTEE A STEADY-STATE SOLUTION: C QLIMIT=QLIM0+DQL_DE*ELEVAT HEATFL=MAX(HEATFL,QLIMIT) HEATFL=MIN(HEATFL,QLIM1) C C This is the re-entry point when heat-flow must be reduced, C and the entire iteration process restarted, for this point: 10 CONTINUE C C INFER CRUSTAL THICKNESS FROM ISOSTASY, AND MANTLE LITHOSPHERE C THICKNESS FROM ASTHENOSPHERE ADIABAT TEMPERATURE (@ 100 km). C THICKC=35.E3 THICKM=65.E3 GEOTH1=TSURF GEOTH2=HEATFL/CONDUC(1) GEOTH3= -RADIO(1)/(2.*CONDUC(1)) GEOTH4=0. GEOTH7= -RADIO(2)/(2.*CONDUC(2)) GEOTH8=0. C DO 150 ITER=1,MAXITR QRED=HEATFL-THICKC*RADIO(1) TMOHO=GEOTH1+GEOTH2*THICKC+GEOTH3*THICKC**2 C C REMEMBER VALUE FOR COVERGENCE REPORT: TM(ITER)=TMOHO C GEOTH5=TMOHO GEOTH6=QRED/CONDUC(2) TEST=GEOTH5+GEOTH6*THICKM+GEOTH7*THICKM**2 C C REMEMBER VALUE FOR COVERGENCE REPORT: TA(ITER)=TEST C TERR0R=TEST-TASTHK DTL= -TERR0R/GEOTH6 C C NOTE: NEXT LINE USES DAMPING FACTOR: DTL=DTL*DAMP C DTL=MIN(DTL,+10000.) DTL=MAX(DTL,-10000.) C THICKM=THICKM+DTL C IF ((THICKM.LT.0.).AND.(.NOT.WARNM1)) THEN WRITE (IUNITL,11) PLON,PLAT 11 FORMAT(' WARNING: (',F8.3,'E, ',F7.3,'N) ', + 'HIT MINIMUM MANTLE LITHOSPHERE THICKNESS OF 0' + /' HEAT-FLOW MAY BE UNREASONABLY ' + ,'HIGH.') WARNM1=.TRUE. END IF THICKM=MAX(THICKM,0.) IF (((THICKC+THICKM).GT.HLMAX).AND.(.NOT.WARNL2)) THEN WRITE (IUNITL,12) PLON,PLAT,HLMAX 12 FORMAT(' WARNING: (',F8.3,'E, ',F7.3,'N) ', + 'HIT MAXIMUM TOTAL LITHOSPHERE THICKNESS OF ' + ,F7.0 + /' SO TEMPERATURE AT THIS DEPTH ' + ,'WILL BE TOO LOW') WARNL2=.TRUE. END IF THICKM=MIN(THICKM,HLMAX-THICKC) C C REMEMBER VALUE FOR COVERGENCE REPORT: HM(ITER)=THICKM C CALL SQUEEZ (INPUT,ALPHAT,ELEVAT, + GEOTH1,GEOTH2,GEOTH3,GEOTH4, + GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN,IUNITT, + ONEKM,RHOAST,RHOBAR,RHOH2O, + TEMLIM,THICKC,THICKC+THICKM, + OUTPUT,TAUZZ,SIGZZB) C DC= -SIGZZB/(GMEAN*(RHOBAR(2)-RHOBAR(1))) C C NOTE: NEXT LINE USES DAMPING FACTOR: THICKC=THICKC+DAMP*DC C IF ((THICKC.LT.CLIMIT).AND.(.NOT.WARNC1)) THEN WRITE (IUNITL,13) PLON,PLAT,CLIMIT 13 FORMAT(' WARNING: (',F8.3,'E, ',F7.3,'N) ', + 'HIT MINIMUM CRUSTAL THICKNESS OF ',F7.0 + /' CHECK FOR LOW-PRESSURE ANOMALY AT ' + ,'BASE OF PLATE') WARNC1=.TRUE. END IF THICKC=MAX(THICKC,CLIMIT) IF ((THICKC.GT.HCMAX).AND.(.NOT.WARNC2)) THEN WRITE (IUNITL,14) PLON,PLAT,HCMAX 14 FORMAT(' WARNING: (',F8.3,'E, ',F7.3,'N) ', + 'HIT MAXIMUM CRUSTAL THICKNESS OF ',F7.0 + /' CHECK FOR HIGH-PRESSURE ANOMALY AT ' + ,'BASE OF PLATE') WARNC2=.TRUE. END IF THICKC=MIN(THICKC,HCMAX) C C REMEMBER VALUE FOR COVERGENCE REPORT: HC(ITER)=THICKC C 150 CONTINUE C C TEST FOR SUCCESSFUL CONVERGENCE: C C ISOSTASY IS BAD IF BASAL PRESSURE ANOMALY IS LARGE: BADP=(ABS(SIGZZB).GT.4.E7) C C GEOTHERM MIGHT NOT CONNECT TO ASTHENOSPHERE ADIABAT: BADT=(ABS(TERR0R).GT.20.) C C COMPUTATION MIGHT BE SUBJECT TO DIVERGENCE: DC25=ABS(HC(MAXITR)-HC(MAXITR-1)) DM25=ABS(HM(MAXITR)-HM(MAXITR-1)) DC24=ABS(HC(MAXITR-1)-HC(MAXITR-2)) DM24=ABS(HM(MAXITR-1)-HM(MAXITR-2)) DIVERG=((DC25.GT.100.).AND.(DC25.GT.DC24)).OR. + ((DM25.GT.100.).AND.(DM25.GT.DM24)) C IF (DIVERG) THEN WRITE (IUNITL,101) PLON,PLAT 101 FORMAT(/' ERROR: (',F8.3,'E, ',F7.3,'N) DID NOT ', + 'CONVERGE:') END IF IF (BADP) THEN WRITE (IUNITL,102) PLON,PLAT 102 FORMAT(/' ERROR: (',F8.3,'E, ',F7.3,'N) IS NOT ', + 'ISOSTATIC WITH RISE:') END IF IF (BADT) THEN WRITE (IUNITL,103) PLON,PLAT 103 FORMAT(/' ERROR: (',F8.3,'E, ',F7.3,'N) DOES NOT ' + ,'CONNECT TO ASTHENOSPHERE ADIABAT') END IF IF (BADP.OR.BADT.OR.DIVERG) THEN WRITE (IUNITL,201) 201 FORMAT(/' TRIAL CRUST T.MOHO MANTLE-' + ,'LITHOSPHERE T.ASTH') 202 FORMAT (' ',I15,F8.0,F7.0,F19.0,F7.0) DO 299 ITER=1,MAXITR WRITE(IUNITL,202)ITER,HC(ITER),TM(ITER),HM(ITER), + TA(ITER) 299 CONTINUE END IF C C It may be necessary to repeat the whole computation with C reduced heat-flow if results were hotter than the asthenosphere C temperature: IF (BADT.AND.(TERR0R.GT.0).AND.(HEATFL.GT.(1.01*QLIMIT))) THEN QNEW=0.99*HEATFL WRITE (IUNITL,301) HEATFL,QNEW 301 FORMAT(/' REDUCING HEATFLOW FROM ',ES10.3, + ' TO ',ES10.3,' AND REPEATING COMPUTATION:') HEATFL=QNEW GO TO 10 END IF C RETURN END C C C BLOCK DATA BD1 C C DEFINE "WEIGHT" (GAUSSIAN INTEGRATION WEIGHTS) OF THE C SEVEN INTEGRATION POINTS IN EACH ELEMENT, DEFINED BY INTERNAL C COORDINATES "POINTS(3,7)", WHERE POINTS(1-3,M)=S1-S3 OF C INTEGRATION POINT NUMBER M. 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 POINTS,WEIGHT COMMON /S1S2S3/ POINTS COMMON /WGTVEC/ WEIGHT DIMENSION POINTS(3,7),WEIGHT(7) C 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. "POINTS" IS ALSO THE SET OF NODAL FUNCTIONS C FOR UNPROJECTED SCALAR QUANTITIES WITHIN AN ELEMENT. DATA POINTS / + 0.3333333333333333D0, 0.3333333333333333D0, 0.3333333333333333D0, + 0.0597158733333333D0, 0.4701420633333333D0, 0.4701420633333333D0, + 0.4701420633333333D0, 0.0597158733333333D0, 0.4701420633333333D0, + 0.4701420633333333D0, 0.4701420633333333D0, 0.0597158733333333D0, + 0.7974269866666667D0, 0.1012865066666667D0, 0.1012865066666667D0, + 0.1012865066666667D0, 0.7974269866666667D0, 0.1012865066666667D0, + 0.1012865066666667D0, 0.1012865066666667D0, 0.7974269866666667D0/ 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) AT THE 7 INTEGRATION POINTS IN EACH FAULT ELEMENT, C DEFINED BY INTERNAL COORDINATE "FPOINT(M=1,...,7)" C 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(4,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 NODE2 (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 4 NODAL FUNCTIONS (ONE PER NODE) C AT EACH OF THESE 7 INTEGRATION POINTS IN THE FAULT ELEMENT. C A SPECIAL CONVENTION IS THAT THE NODAL FUNCTION OF NODE 3 C IS THE NEGATIVE OF THAT FOR NODE 2, WHILE THE NODAL FUNCTION C FOR NODE 4 IS THE NEGATIVE OF THAT FOR NODE 1. THIS SIMPLIFIES C MANY EXPRESSIONS IN WHICH WE WOULD OTHERWISE HAVE TO HAVE C A SEPARATE FACTOR OR +1 OR -1 FOR THE TWO SIDES OF THE FAULT. DATA FPHI/ + 0.9745539D0, 0.0254461D0, -0.0254461D0, -0.9745539D0, + 0.8707656D0, 0.1292344D0, -0.1292344D0, -0.8707656D0, + 0.7029226D0, 0.2970774D0, -0.2970774D0, -0.7029226D0, + 0.5000000D0, 0.5000000D0, -0.5000000D0, -0.5000000D0, + 0.2970774D0, 0.7029226D0, -0.7029226D0, -0.2970774D0, + 0.1292344D0, 0.8707656D0, -0.8707656D0, -0.1292344D0, + 0.0254461D0, 0.9745539D0, -0.9745539D0, -0.0254461D0/ C END C C C SUBROUTINE SQUARE (INPUT,BRIEF,FDIP,IUNIT8, + MXBN,MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES, + NUMEL,NUMNOD,RDETH,WEDGE, + MODIFY,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,EDGEFS, + EDGETS,FLEN,FPFLT,FPSFER, + FARG,NCOND,NODCON,SITA, + WORK,CHECKN,LIST) C C CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID C LOGICAL AGREED,ALLOK,BRIEF,FOUND,SWITCH,VERT1,VERT2 C C NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL CHECKN,EDGEFS,EDGETS C CHARACTER*21 OBLIQU,TAG1,TAG2,VERTIC DOUBLE PRECISION FPHI,FPOINT,FGAUSS COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS DIMENSION FANGLE(2),FPHI(4,7),FPOINT(7),FGAUSS(7),PHI(2),THETA(2) DIMENSION AREA(MXEL),CHECKN(MXNODE), + DETJ(7,MXEL), + DXS(2,2,3,7,MXEL),DYS(2,2,3,7,MXEL), + DXSP(3,7,MXEL),DYSP(3,7,MXEL), + EDGEFS(2,MXFEL),EDGETS(3,MXEL),FDIP(2,MXFEL), + FLEN(MXFEL), + FPFLT(2,2,2,7,MXFEL), + FPSFER(2,2,3,7,MXEL),FARG(2,MXFEL), + LIST(MXSTAR),NODCON(MXBN), + NODEF(4,MXFEL),NODES(3,MXEL), + SITA(7,MXEL),XNODE(MXNODE),YNODE(MXNODE) DATA OBLIQU /'(DIP SLIP IS ALLOWED)'/ DATA VERTIC /'(STRIKE-SLIP ONLY) '/ C C (1) CHECK THAT ALL NODES ARE CONNECTED TO AT LEAST ONE C CONTINUUM (TRIANGULAR) ELEMENT OR FAULT ELEMENT; C DO 110 I=1,NUMNOD CHECKN(I)=.FALSE. 110 CONTINUE DO 130 I=1,NUMEL DO 120 J=1,3 CHECKN(NODES(J,I))=.TRUE. 120 CONTINUE 130 CONTINUE DO 136 I=1,NFL DO 134 J=1,4 CHECKN(NODEF(J,I))=.TRUE. 134 CONTINUE 136 CONTINUE ALLOK=.TRUE. DO 140 I=1,NUMNOD ALLOK=ALLOK.AND.CHECKN(I) 140 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,150) 150 FORMAT(' BAD GRID TOPOLOGY: FOLLOWING REAL NODES DO NOT'/ 1 ' BELONG TO ANY TRIANGULAR CONTINUUM ELEMENT'/ 2 ' OR FAULT ELEMENT:') DO 160 I=1,NUMNOD IF (.NOT.CHECKN(I)) WRITE (IUNIT8,155) I 155 FORMAT (' ',43X,I6) 160 CONTINUE CALL PAUSE() STOP ENDIF C C (2) 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,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=5-J1 NJ2=NODEF(J2,I) LIST(2)=NJ2 CHECKN(NJ2)=.TRUE. NINSUM=2 C FIND SHORTEST FAULT CONNECTED TO EITHER ONE DX=XNODE(NJ1)-XNODE(NJ2) DY=YNODE(NJ1)-YNODE(NJ2) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NJ1)) SHORT=SQRT(DX**2+DY**2) DO 470 K=1,NFL NL1=NODEF(1,K) NL2=NODEF(2,K) NL3=NODEF(3,K) NL4=NODEF(4,K) IF ((NJ1.EQ.NL1).OR.(NJ2.EQ.NL1).OR. + (NJ1.EQ.NL2).OR.(NJ2.EQ.NL2).OR. + (NJ1.EQ.NL3).OR.(NJ2.EQ.NL3).OR. + (NJ1.EQ.NL4).OR.(NJ2.EQ.NL4)) THEN DX=XNODE(NL1)-XNODE(NL2) DY=YNODE(NL1)-YNODE(NL2) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NL1)) TEST=SQRT(DX**2+DY**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 DX=XNODE(NJ1)-XNODE(K) DY=YNODE(NJ1)-YNODE(K) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NJ1)) R2=DX**2+DY**2 IF (R2.LT.T2) THEN NINSUM=NINSUM+1 IF (NINSUM.GT.MXSTAR) THEN WRITE(IUNIT8,421) 421 FORMAT(/' INCREASE VALUE' + ,' OF PARAMETER MAXATP.') CALL PAUSE() STOP ENDIF LIST(NINSUM)=K CHECKN(K)=.TRUE. ENDIF ENDIF 471 CONTINUE C (QUICK EXIT IF ALL NODES IN SAME PLACE) AGREED=.TRUE. DO 472 K=2,NINSUM AGREED=AGREED.AND. + (XNODE(LIST(K)).EQ.XNODE(LIST(1))).AND. + (YNODE(LIST(K)).EQ.YNODE(LIST(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 (3) COMPUTE DERIVITIVES OF NODAL C FUNCTIONS AT INTEGRATION POINTS; C THEN CHECK FOR NEGATIVE AREAS C CALL DERIV (INPUT,IUNIT8,MXEL,MXNODE,NODES,NUMEL, + RDETH,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,FPSFER,SITA) ALLOK=.TRUE. DO 620 I=1,NUMEL DO 610 M=1,7 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) WRITE(IUNIT8,606) AREA(I),DETJ(M,I) 606 FORMAT('AREA = ',1P,E14.4,' DETJ: ',0P,F14.6) ALLOK=.FALSE. ENDIF 610 CONTINUE 620 CONTINUE IF (.NOT.ALLOK) THEN CALL PAUSE() STOP ENDIF C C (4) COMPUTE LENGTHS OF FAULT ELEMENTS. C DO 750 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) THETA1=XNODE(N1) THETA2=XNODE(N2) PHI1 =YNODE(N1) PHI2 =YNODE(N2) FLEN(I)=FLTLEN (PHI1,PHI2,RDETH,THETA1,THETA2) 750 CONTINUE C C (5) MAKE A LIST OF NODES THAT ARE ON THE BOUNDARY AND REQUIRE C BOUNDARY CONDITIONS (NODCON); THESE ARE IN COUNTERCLOCKWISE C ORDER. ALSO MAKE A LISTS OF ELEMENT SIDES WHICH CONTAIN THESE C NODES: EDGETS AND EDGEFS. C NCOND=0 DO 801 I=1,NUMNOD CHECKN(I)=.FALSE. 801 CONTINUE DO 802 I=1,NFL EDGEFS(1,I)=.FALSE. EDGEFS(2,I)=.FALSE. 802 CONTINUE DO 810 I=1,NUMEL DO 809 J=1,3 CALL NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) IF (KELE.GT.0) THEN C (ORDINARY INTERIOR SIDE) EDGETS(J,I)=.FALSE. ELSE IF (KFAULT.EQ.0) THEN C (EXTERIOR SIDE) EDGETS(J,I)=.TRUE. N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J+1,3)+1,I) IF (.NOT.CHECKN(N1)) THEN NCOND=NCOND+1 CHECKN(N1)=.TRUE. ENDIF IF (.NOT.CHECKN(N2)) THEN NCOND=NCOND+1 CHECKN(N2)=.TRUE. ENDIF ELSE C (TRIANGULAR ELEMENT HAS AN EXTERIOR FAULT ELEMENT C ADJACENT TO IT) EDGETS(J,I)=.FALSE. N1=NODES(MOD(J, 3)+1,I) IF (NODEF(2,KFAULT).EQ.N1) THEN EDGEFS(2,KFAULT)=.TRUE. DO 806 K=3,4 N=NODEF(K,KFAULT) IF (.NOT.CHECKN(N)) THEN NCOND=NCOND+1 CHECKN(N)=.TRUE. ENDIF 806 CONTINUE ELSE EDGEFS(1,KFAULT)=.TRUE. DO 808 K=1,2 N=NODEF(K,KFAULT) IF (.NOT.CHECKN(N)) THEN NCOND=NCOND+1 CHECKN(N)=.TRUE. ENDIF 808 CONTINUE ENDIF ENDIF 809 CONTINUE 810 CONTINUE IF (NCOND.GT.MXBN) THEN WRITE(IUNIT8,820) NCOND 820 FORMAT(/' INCREASE PARAMETER MAXBN TO',I6,' AND RECOMPILE.') CALL PAUSE() STOP ENDIF C C STOP WORK IF NO BOUNDARY NODES FOUND (GLOBAL GRID) C IF (NCOND.EQ.0) GO TO 899 C C BEGIN CIRCUIT WITH LOWEST-NUMBERED BOUNDARY NODE DO 830 I=1,NUMNOD IF (CHECKN(I)) GO TO 831 830 CONTINUE 831 NODCON(1)=I NDONE=1 NLEFT=NCOND C BEGINNING OF INDEFINATE LOOP WHICH TRACES AROUND THE PERIMETER. C EACH TIME, IT PROGRESSES BY ONE OF 3 STEPS: C -2 NODES AT A TIME ALONG A TRIANGLE SIDE, OR C -2 NODES AT A TIME ALONG A FAULT ELEMENT SIDE, OR C -BY FINDING ANOTHER (CORNER) NODE WHICH SHARES THE SAME LOCATION. C BEGINNING OF MAIN INDEFINATE LOOP: 840 NODE=NODCON(NDONE) X=XNODE(NODE) Y=YNODE(NODE) C LOOK FOR AN ADJACENT TRIANGULAR ELEMENT USING THIS NODE. DO 844 I=1,NUMEL DO 842 J=1,3 IF (EDGETS(J,I)) THEN N1=NODES(MOD(J,3)+1,I) IF (N1.EQ.NODE) GO TO 846 ENDIF 842 CONTINUE 844 CONTINUE GO TO 850 846 N2=NODES(MOD(J+1,3)+1,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR AN ADJACENT FAULT ELEMENT USING THIS NODE. 850 DO 854 I=1,NFL IF (EDGEFS(1,I)) THEN IF (NODEF(1,I).EQ.NODE) THEN N2=NODEF(2,I) GO TO 856 ENDIF ELSE IF (EDGEFS(2,I)) THEN IF (NODEF(3,I).EQ.NODE) THEN N2=NODEF(4,I) GO TO 856 ENDIF ENDIF 854 CONTINUE GO TO 860 856 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR ANOTHER EXTERIOR CORNER NODE AT SAME LOCATION 860 DO 865 I=1,NUMNOD IF ((I.NE.NODE).AND.CHECKN(I)) THEN IF ( (ABS(XNODE(I)-X).LT.1.E-6) .AND. + (ABS(YNODE(I)-Y).LT.1.E-6) ) GO TO 867 ENDIF 865 CONTINUE WRITE(IUNIT8,866) NODE 866 FORMAT(' BAD GRID TOPOLOGY: WHILE TRACING PERIMETER,'/ + ' COULD NOT FIND ANY WAY TO CONTINUE FROM NODE ',I6/ + ' EITHER THROUGH SHARED BOUNDARY ELEMENTS, OR'/ + ' THROUGH OTHER BOUNDARY NODES SHARING THE SAME ', + 'POSITION.') CALL PAUSE() STOP 867 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=I CHECKN(I)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) GO TO 840 C END OF INDEFINATE LOOP WHICH TRACES AROUND PERIMETER. 870 IF (.NOT.BRIEF) THEN WRITE(IUNIT8,880) 880 FORMAT(/ /' HERE FOLLOWS A LIST, IN CONSECUTIVE ORDER,'/ + ' OF THE NODES WHICH DEFINE THE PERIMETER'/ + ' OF THE MODEL; THESE NODES REQUIRE BOUNDARY', + ' CONDITIONS:'/' BC# NODE') DO 890 I=1,NCOND N=NODCON(I) WRITE(IUNIT8,882) I, N 882 FORMAT(' ',2I6) 890 CONTINUE N=NODCON(1) WRITE (IUNIT8,892) N 892 FORMAT(' (NOTE: NODE ',I6,' COMPLETES THE LOOP, BUT WILL', + ' NOT BE LISTED TWICE.)') ENDIF 899 CONTINUE C C (6) 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 VERT1=ABS(DELD1).LE.WEDGE VERT2=ABS(DELD2).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 SWITCH=((NVPART.GT.0).AND.(NVPART.LT.2)) 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. WRITE (IUNIT8,905) I,DIP1,TAG1,DIP2,TAG2 905 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',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) 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 912 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',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 914 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',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 (7) CALCULATE FAULT ARGUMENT (IN RADIANS, MEASURED COUNTERCLOCKWISE C FROM +THETA) FOR EACH INTEGRATION POINT IN EACH FAULT ELEMENT. C C HERE USE CROSS PRODUCT "D" OF POSITION VECTOR "A" AND "B" C OF FAULT ELEMENT NODE "N1" AND "N2" RESPECTIVELY, TO CALCULATE C VECTOR "D" IN PLANE OF "A" AND "B" AND PERPENDICULAR TO "A". C THEN USE "D" AND "A" TO DECIDE THE POSITION VECTOR "GOLD" OF C INTERPOLATION POINT "M" AND "G" OF "M+1". FINALLY USE THE C DOT PRODUCT OF OF "GOLD-POL" AND "G-GOLD" TO DECIDE THE ANGLE C BETWEEN THESE TWO VECTORS, I.E., FARG. HERE POL IS POSITION C VECTOR OF NORTH POLE. DO 1000 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) THETA(1)=XNODE(N1) THETA(2)=XNODE(N2) PHI(1) =YNODE(N1) PHI(2) =YNODE(N2) CALL FANGLS (INPUT,PHI,THETA, + OUTPUT,FANGLE) DO 900 J=1,2 FARG(J,I)=FANGLE(J) FANGLE(J)=FANGLE(J)*57.29577951 900 CONTINUE 1000 CONTINUE C C (8) 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-4, 2-3 (J = 1 OR 2): DO 1900 J=1,2 C DIP MUST BE WITHIN "WEDGE" OF VERTICAL FOR CONSTRAINT: IF (ABS(FDIP(J,I)-1.570796).LE.WEDGE) THEN NAZI=J N1=J IF(J.EQ.1) THEN N4=4 ELSE N4=3 ENDIF NODE1=NODEF(N1,I) NODE4=NODEF(N4,I) C NO CONSTRAINT APPLIED WHERE A FAULT ENDS: IF (NODE1.NE.NODE4) 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. + (NODE4.EQ.NODEF(4,L))).OR. + ((NODE1.EQ.NODEF(4,L)).AND. + (NODE4.EQ.NODEF(1,L)))) THEN FOUND=.TRUE. NUMBER=L NAZL=1 GO TO 1601 ENDIF ENDIF IF (ABS(FDIP(2,L)-1.5708).LE.WEDGE) THEN IF (((NODE1.EQ.NODEF(2,L)).AND. + (NODE4.EQ.NODEF(3,L))).OR. + ((NODE1.EQ.NODEF(3,L)).AND. + (NODE4.EQ.NODEF(2,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): IF(NAZI.EQ.NAZL) THEN AZL=FARG(NAZL,NUMBER)+3.141592654 ELSE AZL=FARG(NAZL,NUMBER) ENDIF AZI=FARG(NAZI,I) COSZ=0.5*(COS(AZI)+COS(AZL)) SINZ=0.5*(SIN(AZI)+SIN(AZL)) AZIMUT=ATAN2F(SINZ,COSZ) FARG(NAZI,I)=AZIMUT IF(NAZL.EQ.NAZI) THEN FARG(NAZL,NUMBER)=AZIMUT-3.141592654 ELSE FARG(NAZL,NUMBER)=AZIMUT 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 C (9) CALCULATE NODAL FUNCTION AT INTERPOLATION POINTS C ON SPHERICAL GREAT CIRCLE FAULT. C CALL FNODAL (INPUT,IUNITT,MXFEL, + MXNODE,NFL,NODEF,XNODE,YNODE, + OUTPUT,FPFLT) C IF (.NOT. BRIEF) WRITE (IUNIT8,9999) 9999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE DERIV (INPUT,IUNITT,MXEL,MXNODE, + NODES,NUMEL, + RDETH,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,FPSFER,SITA) C C SETS UP 6 VECTOR NODAL FUNCTIONS (FPSFER) OF EACH SPHERICAL C TRIANGLE FINITE ELEMENT, AT EACH OF ITS 7 INTEGRATION POINTS. C CALCULATES DXS AND DYS, THE THETA-DERIVITIVE AND PHI-DERIVITIVE C OF EACH OF THESE 6 VECTOR NODAL FUNCTIONS. C ALSO COMPUTES AREA, THE AREAS OF THE PLANE TRIANGLES. C ALSO COMPUTES DETJ, THE LOCAL RATIO OF AREAS ON THE SPHERE C TO AREAS ON THE PLANE TRIANGLES. C DOUBLE PRECISION POINTS DOUBLE PRECISION FFF,SKKC,SKKE,SNCSNE,SNCCSE,CSCCSE,CSCSNE DOUBLE PRECISION XA,XB,XC,YA,YB,YC,ZA,ZB,ZC,XYZP DIMENSION XNODE(MXNODE),YNODE(MXNODE),NODES(3,MXEL),AREA(MXEL) DIMENSION DETJ(7,MXEL) DIMENSION DXS(2,2,3,7,MXEL),DYS(2,2,3,7,MXEL) DIMENSION DXSP(3,7,MXEL),DYSP(3,7,MXEL),POINTS(3,7) DIMENSION PHI(3),THETA(3),SKKC(3),SKKE(3),FFF(3), + SITA(7,MXEL),FPSFER(2,2,3,7,MXEL) DIMENSION XLON(3),YLAT(3) COMMON /S1S2S3/ POINTS DATA OEZOPI /57.29578/ C DO 900 I=1,NUMEL DO 100 J=1,3 THETA(J)=XNODE(NODES(J,I)) PHI(J) =YNODE(NODES(J,I)) 100 CONTINUE X21=SIN(THETA(2))*COS(PHI(2))-SIN(THETA(1))*COS(PHI(1)) X31=SIN(THETA(3))*COS(PHI(3))-SIN(THETA(1))*COS(PHI(1)) Y21=SIN(THETA(2))*SIN(PHI(2))-SIN(THETA(1))*SIN(PHI(1)) Y31=SIN(THETA(3))*SIN(PHI(3))-SIN(THETA(1))*SIN(PHI(1)) Z21=COS(THETA(2))-COS(THETA(1)) Z31=COS(THETA(3))-COS(THETA(1)) A=Y21*Z31-Y31*Z21 B=Z21*X31-Z31*X21 C=X21*Y31-X31*Y21 AREAP=SQRT(A*A+B*B+C*C) AREA(I)=0.5*RDETH*RDETH*AREAP IF (AREA(I).LE.0.0) THEN DO 9 K=1,3 XLON(K)=PHI(K)*OEZOPI YLAT(K)=90.-THETA(K)*OEZOPI 9 CONTINUE WRITE (IUNITT,10) AREA(I), I, + (NODES(K,I), XLON(K), YLAT(K), K=1,3) 10 FORMAT (/' ERR0R: NON-POSITIVE AREA OF ',1P,E10.2, + ' FOR ELEMENT ',I5/ + ' NODES LONGITUDE LATITUDE', + 3(/' ',I7,0P,2F10.3)) CALL PAUSE() STOP ENDIF PNX=A/AREAP PNY=B/AREAP PNZ=C/AREAP DD1=SIN(THETA(1))*COS(PHI(1))*PNX DD2=SIN(THETA(1))*SIN(PHI(1))*PNY DD3=COS(THETA(1))*PNZ DD=DD1+DD2+DD3 C C THIS PART IS TO TEST IF KONG'S METHOD AND BIRD'S METHOD ARE SAME FOR C CALCULATING DEERIVATIVE C XA=SIN(THETA(1))*COS(PHI(1)) XB=SIN(THETA(2))*COS(PHI(2)) XC=SIN(THETA(3))*COS(PHI(3)) YA=SIN(THETA(1))*SIN(PHI(1)) YB=SIN(THETA(2))*SIN(PHI(2)) YC=SIN(THETA(3))*SIN(PHI(3)) ZA=COS(THETA(1)) ZB=COS(THETA(2)) ZC=COS(THETA(3)) CKA=(YB*ZC-ZB*YC)*XA+(ZB*XC-XB*ZC)*YA+(XB*YC-YB*XC)*ZA C C C DO 800 M=1,7 SNCCSE=0.0 SNCSNE=0.0 COSM=0.0 DO 200 J=1,3 SNCCSE=SNCCSE+POINTS(J,M)*SIN(THETA(J))*COS(PHI(J)) SNCSNE=SNCSNE+POINTS(J,M)*SIN(THETA(J))*SIN(PHI(J)) COSM=COSM+POINTS(J,M)*COS(THETA(J)) 200 CONTINUE XYZP=SQRT(SNCCSE*SNCCSE+SNCSNE*SNCSNE+COSM*COSM) SNCCSE=SNCCSE/XYZP SNCSNE=SNCSNE/XYZP COSM=COSM/XYZP SITAJ=ACOS(COSM) TY=SNCSNE TX=SNCCSE PHAIJ=ATAN2F(TY,TX) CSCCSE=COS(SITAJ)*COS(PHAIJ) CSCSNE=COS(SITAJ)*SIN(PHAIJ) C C BIRD'S METHOD C FFF(1)=((YB*ZC-ZB*YC)*SNCCSE+(ZB*XC-XB*ZC)*SNCSNE + +(XB*YC-YB*XC)*COSM)/CKA FFF(2)=((YC*ZA-ZC*YA)*SNCCSE+(ZC*XA-XC*ZA)*SNCSNE + +(XC*YA-YC*XA)*COSM)/CKA FFF(3)=((YA*ZB-ZA*YB)*SNCCSE+(ZA*XB-XA*ZB)*SNCSNE + +(XA*YB-YA*XB)*COSM)/CKA SKKC(1)=((YB*ZC-ZB*YC)*CSCCSE + +(ZB*XC-XB*ZC)*CSCSNE + -(XB*YC-YB*XC)*SIN(SITAJ))/CKA SKKC(2)=((YC*ZA-ZC*YA)*CSCCSE + +(ZC*XA-XC*ZA)*CSCSNE + -(XC*YA-YC*XA)*SIN(SITAJ))/CKA SKKC(3)=((YA*ZB-ZA*YB)*CSCCSE + +(ZA*XB-XA*ZB)*CSCSNE + -(XA*YB-YA*XB)*SIN(SITAJ))/CKA SKKE(1)=(-(YB*ZC-ZB*YC)*SNCSNE + +(ZB*XC-XB*ZC)*SNCCSE)/CKA SKKE(2)=(-(YC*ZA-ZC*YA)*SNCSNE + +(ZC*XA-XC*ZA)*SNCCSE)/CKA SKKE(3)=(-(YA*ZB-ZA*YB)*SNCSNE + +(ZA*XB-XA*ZB)*SNCCSE)/CKA C C C SITA(M,I)=SITAJ RR1=SIN(SITAJ)*COS(PHAIJ) RR2=SIN(SITAJ)*SIN(PHAIJ) RR3=COS(SITAJ) RN=RR1*PNX+RR2*PNY+RR3*PNZ PP=DD/RN DPDC=(COS(SITAJ)*COS(PHAIJ)*PNX+COS(SITAJ)*SIN(PHAIJ)*PNY + -SIN(SITAJ)*PNZ) DPDE=(-SIN(SITAJ)*SIN(PHAIJ)*PNX+ + SIN(SITAJ)*COS(PHAIJ)*PNY) DDPN=PP/RN DPDC=-DDPN*DPDC DPDE=-DDPN*DPDE IF(SITA(M,I).LE.0.0.OR.SITA(M,I).GE.3.141592654) THEN SITAMI=SITA(M,I)*57.29577951 WRITE(IUNITT,220) M,I,SITAMI 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5,' OF ELEMENT', + I5,' IS OUT RANGE', + E14.4) CALL PAUSE() STOP ENDIF DO 500 J=1,3 DXSP(J,M,I)= DPDC*FFF(J)+PP*SKKC(J) DYSP(J,M,I)= DPDE*FFF(J)+PP*SKKE(J) CSCS=COS(THETA(J))*COS(PHI(J)) CSSN=COS(THETA(J))*SIN(PHI(J)) SNC=SIN(THETA(J)) SNE=SIN(PHI(J)) CSE=COS(PHI(J)) FPSFER(1,1,J,M,I)= CSCS*CSCCSE+CSSN*CSCSNE + +SNC*SIN(SITAJ) FPSFER(2,1,J,M,I)=-SNE*CSCCSE+CSE*CSCSNE FPSFER(1,2,J,M,I)=-CSCS*SIN(PHAIJ)+CSSN*COS(PHAIJ) FPSFER(2,2,J,M,I)=SNE*SIN(PHAIJ)+CSE*COS(PHAIJ) DXS(1,1,J,M,I)=(-CSCS*SNCCSE-CSSN*SNCSNE + +SNC*COS(SITAJ))*FFF(J) + +FPSFER(1,1,J,M,I)*SKKC(J) DXS(2,1,J,M,I)=(SNE*SNCCSE-CSE*SNCSNE)*FFF(J) + +FPSFER(2,1,J,M,I)*SKKC(J) DYS(1,1,J,M,I)=(-CSCS*CSCSNE+CSSN*CSCCSE)*FFF(J) + +FPSFER(1,1,J,M,I)*SKKE(J) DYS(2,1,J,M,I)=(SNE*CSCSNE+CSE*CSCCSE)*FFF(J) + +FPSFER(2,1,J,M,I)*SKKE(J) DXS(1,2,J,M,I)=FPSFER(1,2,J,M,I)*SKKC(J) DXS(2,2,J,M,I)=FPSFER(2,2,J,M,I)*SKKC(J) DYS(1,2,J,M,I)=(-CSCS*COS(PHAIJ)-CSSN*SIN(PHAIJ)) + *FFF(J) + +FPSFER(1,2,J,M,I)*SKKE(J) DYS(2,2,J,M,I)=(SNE*COS(PHAIJ)-CSE*SIN(PHAIJ)) + *FFF(J) + +FPSFER(2,2,J,M,I)*SKKE(J) FPSFER(1,1,J,M,I)=FPSFER(1,1,J,M,I)*FFF(J) FPSFER(2,1,J,M,I)=FPSFER(2,1,J,M,I)*FFF(J) FPSFER(1,2,J,M,I)=FPSFER(1,2,J,M,I)*FFF(J) FPSFER(2,2,J,M,I)=FPSFER(2,2,J,M,I)*FFF(J) 500 CONTINUE PFQ=FFF(1)+FFF(2)+FFF(3) DETJ(M,I)=RN**3/(DD*DD) 800 CONTINUE 900 CONTINUE RETURN END C C C SUBROUTINE NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) C C DETERMINE WHETHER THERE ARE MORE ELEMENTS ADJACENT TO SIDE J OF C TRIANGULAR CONTINUUM ELEMENT #I. C J = 1 MEANS THE SIDE OPPOSITE NODE # NODES(1,I). C J = 2 MEANS THE SIDE OPPOSITE NODE # NODES(2,I). C J = 3 MEANS THE SIDE OPPOSITE NODE # NODES(3,I). C IF A FAULT ELEMENT IS ADJACENT, ITS NUMBER IS KFAULT; OTHERWISE, C KFAULT IS SET TO ZERO. C IF ANOTHER TRIANGULAR CONTINUUM ELEMENT IS ADJACENT (EVEN ACROSS C FAULT ELEMENT KFAULT!) THEN ITS NUMBER IS KELE; OTHERWISE, KELE = 0. C LOGICAL FOUNDF DIMENSION NODEF(4,MXFEL),NODES(3,MXEL) C C THREE NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE: N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J+1,3)+1,I) C CHECK FOR ADJACENT FAULT ELEMENT FIRST: FOUNDF=.FALSE. KFAULT=0 IF (NFL.GT.0) THEN DO 10 K=1,NFL M1=NODEF(1,K) M2=NODEF(2,K) M3=NODEF(3,K) M4=NODEF(4,K) IF (((M1.EQ.N2).AND.(M2.EQ.N1)).OR. + ((M3.EQ.N2).AND.(M4.EQ.N1))) THEN FOUNDF=.TRUE. KFAULT=K GO TO 11 ENDIF 10 CONTINUE 11 IF (.NOT.FOUNDF) KFAULT=0 C IF THERE WAS A FAULT, REPLACE 2 NODE NUMBERS THAT WE SEARCH FOR: IF (FOUNDF) THEN IF (M2.EQ.N1) THEN N1=M3 N2=M4 ELSE N1=M1 N2=M2 ENDIF ENDIF ENDIF C SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT: KELE=0 KLOW=I KHIGH=I C --- BEGIN IRREGULAR LOOP, TO SEARCH OUT NEAREST ELEMENTS FIRST --- 100 KLOW=KLOW-1 IF (KLOW.GE.1) THEN DO 110 L=1,3 M1=NODES(MOD(L, 3)+1,KLOW) M2=NODES(MOD(L+1,3)+1,KLOW) IF ((M2.EQ.N1).AND.(M1.EQ.N2)) THEN KELE=KLOW RETURN ENDIF 110 CONTINUE ENDIF KHIGH=KHIGH+1 IF (KHIGH.LE.NUMEL) THEN DO 120 L=1,3 M1=NODES(MOD(L, 3)+1,KHIGH) M2=NODES(MOD(L+1,3)+1,KHIGH) IF ((M2.EQ.N1).AND.(M1.EQ.N2)) THEN KELE=KHIGH RETURN ENDIF 120 CONTINUE ENDIF IF ((KLOW.GT.1).OR.(KHIGH.LT.NUMEL)) GO TO 100 RETURN END C C C REAL FUNCTION FLTLEN (PHI1,PHI2,RDETH,THETA1,THETA2) C C CALCULATES LENGTH OF GREAT CIRCLE SEGMENT BETWEEN C POINT (THETA1,PHI1) AND POINT (THETA2,PHI2), C IN PHYSICAL LENGTH UNITS (RADIANS*PLANET RADIUS). C DOUBLE PRECISION AB AB =SIN(THETA1)*SIN(THETA2)*COS(PHI1)*COS(PHI2)+ + SIN(THETA1)*SIN(THETA2)*SIN(PHI1)*SIN(PHI2)+ + COS(THETA1)*COS(THETA2) AB=ACOS(AB) FLTLEN=AB*RDETH RETURN END C C C SUBROUTINE READPM (INPUT,IUNIT7, IUNIT8, NAMES, NUMPLT, OFFMAX, + OUTPUT,ACREEP, ALPHAT, BCREEP, BIOT , + BYERLY, CCREEP, CFRIC , CONDUC, + DCREEP, ECREEP, EVERYP, FFRIC , GMEAN , + GRADIE, ICONVE, IPVREF, + MAXITR, OKDELV, OKTOQT, ONEKM, RADIO , + RDETH, REFSTR, RHOAST, RHOBAR, RHOH2O, + TADIAB, + TAUMAX, TEMLIM, TITLE3, TRHMAX, TSURF, + ZBASTH) C C READS STRATEGIC AND TACTICAL INPUT PARAMETERS FROM DEVICE IUNIT7, C AND ECHOES THEM ON DEVICE IUNIT8 WITH ANNOTATIONS. C CHARACTER*2 NAMES,PLTREF CHARACTER*80 TITLE3 LOGICAL EVERYP DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),NAMES(NUMPLT),RADIO(2), + RHOBAR(2),TAUMAX(2),TEMLIM(2),VECTOR(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) 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) 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) 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) 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,*) TADIAB, GRADIE WRITE (IUNIT8,92) TADIAB, GRADIE 92 FORMAT (' ',F10.0,1P,E10.2,' INTERCEPT AND GRADIENT OF UPPER' + ,' MANTLE ADIABAT (K, K/M)') READ (IUNIT7,*) ZBASTH WRITE (IUNIT8,94) ZBASTH 94 FORMAT (' ',1P,E10.2,' DEPTH OF BASE OF ASTHENOSPHERE') C READ (IUNIT7,952) PLTREF 952 FORMAT(A2) WRITE (IUNIT8,954) PLTREF 954 FORMAT(' ',A2,9X,'PLTREF: PLATE DEFINING VELOCITY ', + 'REFERENCE FRAME (EU, NA, AF, ...)') IPVREF=0 DO 956 I=1,NUMPLT IF (NAMES(I).EQ.PLTREF) IPVREF=I 956 CONTINUE IF (IPVREF.EQ.0) THEN WRITE (IUNIT8,958) (NAMES(I),I=1,NUMPLT) 958 FORMAT (/' *** ERR','OR in parameter input file: ***' + /' In line 13 (after ZBASTH, before ICONVE),' + /' in the first two columns of the line,' + /' define the velocity reference frame by' + /' entering one of the following plate names:' + /' ',26(A2,1X)) CALL PAUSE() STOP END IF C READ (IUNIT7,*) ICONVE WRITE (IUNIT8,96) ICONVE 96 FORMAT (' ',I10,' ICONVE: CODE FOR LOWER MANTLE FLOW ', + '(0=NONE)') READ (IUNIT7,*) TRHMAX WRITE (IUNIT8,101) TRHMAX 101 FORMAT (' ',1P,E10.2,' LIMIT ON HORIZONTAL TRACTIONS', + ' APPLIED TO BASE OF PLATE') C CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,VECTOR) TAUMAX(1)=VECTOR(1) TAUMAX(2)=VECTOR(2) C PROVIDE FOR OLD PARAMETER FILES WITH ONLY ONE TAUMAX: IF (TAUMAX(2).LE.0.0) TAUMAX(2)=TAUMAX(1) WRITE (IUNIT8,106) TAUMAX(1),TAUMAX(2) 106 FORMAT (' ',1P,E7.1,',',E7.1, + ' TAUMAX: SEA/LAND UPPER LIMIT', + ' ON INTEGRATED SUBDUCTION DRAG (N/m)') IF ((TAUMAX(1).LT.0.).OR.(TAUMAX(2).LT.0.)) THEN WRITE (IUNIT8,107) 107 FORMAT (/' *** ERR','OR in parameter input file: ***' + /' TAUMAX may not be negative.') CALL PAUSE() STOP END IF C 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) 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)') READ (IUNIT7,*) RDETH WRITE (IUNIT8,155) RDETH 155 FORMAT (' ',1P,E10.3,' RADIUS OF THE EARTH') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,ALPHAT) 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) 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) 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) 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 (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE FANGLS (INPUT, PHI,THETA, + OUTPUT, FANGLE) C C CALCULATE THE ARGUMENTS (ANGLES COUNTERCLOCKWISE FROM +THETA) C AT BOTH ENDS OF AN ARC OF A GREAT CIRCLE. C DOUBLE PRECISION FPOINT COMMON /SFAULT/ FPOINT DIMENSION FANGLE(2),FPOINT(7),PHI(2),THETA(2) DG180=3.141592654 A1=SIN(THETA(1))*COS(PHI(1)) A2=SIN(THETA(1))*SIN(PHI(1)) A3=COS(THETA(1)) B1=SIN(THETA(2))*COS(PHI(2)) B2=SIN(THETA(2))*SIN(PHI(2)) B3=COS(THETA(2)) S=0.99 XX=S*A1+(1.0-S)*B1 YY=S*A2+(1.0-S)*B2 ZZ=S*A3+(1.0-S)*B3 XVAL=SQRT(XX*XX+YY*YY+ZZ*ZZ) XX=XX/XVAL YY=YY/XVAL ZZ=ZZ/XVAL DX=XX-A1 DY=YY-A2 DZ=ZZ-A3 SITA=THETA(1) PHAI=PHI(1) S1=COS(SITA)*COS(PHAI) S2=COS(SITA)*SIN(PHAI) S3=-SIN(SITA) P1=-SIN(PHAI) P2=COS(PHAI) DXX=DX*S1+DY*S2+DZ*S3 DYY=DX*P1+DY*P2 FANGLE(1)=ATAN2F(DYY,DXX) S=0.01 XX=S*A1+(1.0-S)*B1 YY=S*A2+(1.0-S)*B2 ZZ=S*A3+(1.0-S)*B3 XVAL=SQRT(XX*XX+YY*YY+ZZ*ZZ) XX=XX/XVAL YY=YY/XVAL ZZ=ZZ/XVAL DX=B1-XX DY=B2-YY DZ=B3-ZZ SITA=ACOS(ZZ) PHAI=ATAN2F(YY,XX) IF(PHAI.LT.0.0) PHAI=2.0*DG180+PHAI S1=COS(SITA)*COS(PHAI) S2=COS(SITA)*SIN(PHAI) S3=-SIN(SITA) P1=-SIN(PHAI) P2=COS(PHAI) DXX=DX*S1+DY*S2+DZ*S3 DYY=DX*P1+DY*P2 FANGLE(2)=ATAN2F(DYY,DXX) RETURN END C C C SUBROUTINE FNODAL (INPUT,IUNITT,MXFEL,MXNODE,NFL,NODEF, + XNODE,YNODE, + OUTPUT,FPFLT) C C CALCULATES VECTOR NODAL FUNCTIONS AT INTERPOLATION POINT ON A C GREAT CIRCLE FAULT ELEMENT C DOUBLE PRECISION FPHI DIMENSION FPHI(4,7),FPFLT(2,2,2,7,MXFEL),FPP(2,2,2,7), + NODEF(4,MXFEL),PHI(2),THETA(2), + XNODE(MXNODE),YNODE(MXNODE) COMMON /FPHIS/ FPHI C SNCCOP: SIN(SITA)*COS(PHAI) AT INTERPOLATION POINT C SNCSNP SIN(SITA)*SIN(PHAI) AT INTERPOLATION POINT C DO 900 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) THETA(1)=XNODE(N1) THETA(2)=XNODE(N2) PHI(1)=YNODE(N1) PHI(2)=YNODE(N2) CALL SNODAL (INPUT,PHI,THETA, + OUTPUT,FPP) DO 800 M=1,7 DO 500 J=1,2 DO 400 K=1,2 DO 300 L=1,2 FPFLT(L,K,J,M,I)=FPP(L,K,J,M) 300 CONTINUE 400 CONTINUE 500 CONTINUE SNCCOP=FPHI(1,M)*SIN(THETA(1))*COS(PHI(1))+ + FPHI(2,M)*SIN(THETA(2))*COS(PHI(2)) SNCSNP=FPHI(1,M)*SIN(THETA(1))*SIN(PHI(1))+ + FPHI(2,M)*SIN(THETA(2))*SIN(PHI(2)) COSM =FPHI(1,M)*COS(THETA(1))+FPHI(2,M)*COS(THETA(2)) PP=SQRT(SNCCOP*SNCCOP+SNCSNP*SNCSNP+COSM*COSM) COSM=COSM/PP SITA=ACOS(COSM) IF(SITA.LE.0.0.OR.SITA.GE.3.141592654) THEN SITA=SITA*57.29577951 WRITE(IUNITT,220) M,I,SITA 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5, + ' OF FAULT ELEMENT', + I5,' IS OUT RANGE', + E14.4) ENDIF 800 CONTINUE 900 CONTINUE RETURN END C C C SUBROUTINE SNODAL (INPUT,PHI,THETA, + OUTPUT,FPP) C C CALCULATES VECTOR NODAL FUNCTION AT INTERPOLATION POINT ON A C GREAT CIRCLE OF A SIDE OF A FINITE ELEMENT C DOUBLE PRECISION FPHI,PP DIMENSION FPHI(4,7),FPP(2,2,2,7),PHI(2),THETA(2) COMMON /FPHIS/ FPHI C SNCCOP: SIN(SITA)*COS(PHAI) AT INTERPOLATION POINT C SNCSNP SIN(SITA)*SIN(PHAI) AT INTERPOLATION POINT C X1=SIN(THETA(1))*COS(PHI(1)) Y1=SIN(THETA(1))*SIN(PHI(1)) Z1=COS(THETA(1)) X2=SIN(THETA(2))*COS(PHI(2)) Y2=SIN(THETA(2))*SIN(PHI(2)) Z2=COS(THETA(2)) XN=X1+X2 YN=Y1+Y2 ZN=Z1+Z2 XYZN=SQRT(XN*XN+YN*YN+ZN*ZN) XN=XN/XYZN YN=YN/XYZN ZN=ZN/XYZN DD=X1*XN+Y1*YN+Z1*ZN DO 800 M=1,7 XX=FPHI(1,M)*X1+FPHI(2,M)*X2 YY=FPHI(1,M)*Y1+FPHI(2,M)*Y2 ZZ=FPHI(1,M)*Z1+FPHI(2,M)*Z2 PP=SQRT(XX*XX+YY*YY+ZZ*ZZ) XX=XX/PP YY=YY/PP ZZ=ZZ/PP SITAJ=ACOS(ZZ) PHAIJ=ATAN2F(YY,XX) RN=XX*XN+YY*YN+ZZ*ZN PPM=RN/DD CSCS=COS(SITAJ)*COS(PHAIJ) CSSN=COS(SITAJ)*SIN(PHAIJ) SNSN=SIN(SITAJ)*SIN(PHAIJ) SNC=SIN(SITAJ) SNP=SIN(PHAIJ) CSP=COS(PHAIJ) DO 500 J=1,2 FP=FPHI(J,M)*PPM FPP(1,1,J,M)=( COS(THETA(J))*COS(PHI(J))*CSCS + +COS(THETA(J))*SIN(PHI(J))*CSSN + +SIN(THETA(J))*SNC)*FP FPP(2,1,J,M)=(-SIN(PHI(J))*CSCS+COS(PHI(J))*CSSN)*FP FPP(1,2,J,M)=(-COS(THETA(J))*COS(PHI(J))*SNP + +COS(THETA(J))*SIN(PHI(J))*CSP)*FP FPP(2,2,J,M)=( SIN(PHI(J))*SNP + +COS(PHI(J))*CSP)*FP 500 CONTINUE 800 CONTINUE RETURN END C C C SUBROUTINE GETNET (INPUT,IUNIT7,IUNIT8, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEDAT,ELEV,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 DIMENSION CHECKE(MXEL),CHECKF(MXFEL),CHECKN(MXNODE), + DQDTDA(MXNODE),ELEDAT(3,MXEL),ELEV(MXNODE), + FDIP(2,MXFEL), NODEF(4,MXFEL), + NODES(3,MXEL),OFFSET(MXFEL), TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(3),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. 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 BRIEF=.TRUE. C IF (NUMNOD.NE.(NREALN+NFAKEN)) THEN WRITE (IUNIT8,4) NUMNOD,NREALN,NFAKEN 4 FORMAT (/' ERR0R: NUMNOD (',I6,') IS NOT EQUAL TO SUM' + /' OF NREALN (',I6,') AND NFAKEN (',I6,').') CALL PAUSE() STOP ENDIF C IF (NREALN.GT.N1000) THEN WRITE (IUNIT8,5) NREALN,N1000 5 FORMAT (/' ERR0R: NREALN (',I6,') IS GREATER THAN' + /' N1000 (',I6,').') CALL PAUSE() STOP ENDIF C 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.') CALL PAUSE() STOP ENDIF C 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 40 FORMAT (/' THERE ARE',I5,' NODES IN THE GRID') WRITE (IUNIT8,50) 50 FORMAT (/ + 77X,' MANTLE'/ + 77X,' CRUSTAL LITHOSPHERE'/ + ' NODE E-LONGITUDE N-LATITUDE', + ' THETA PHI ELEVATION', + ' HEAT-FLOW THICKNESS THICKNESS'/) 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 IF (INDEX.GT.NREALN) THEN IF ((INDEX.LE.N1000).OR. + (INDEX.GT.(N1000+NFAKEN))) THEN WRITE (IUNIT8,91) INDEX 91 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER: ',I6) CALL PAUSE() STOP ENDIF ENDIF PLON=VECTOR(2) PLAT=VECTOR(3) IF (ABS(PLAT).GT.90.01) THEN WRITE (IUNIT8,92) INDEX 92 FORMAT (/' ERR0R: ABS(LATITUDE) > 90 AT NODE ',I6) CALL PAUSE() STOP ENDIF IF (ABS(PLAT).GT.89.99) THEN WRITE (IUNIT8,93) INDEX 93 FORMAT (/' ERR0R: NODE ',I6,' LIES ON A POLE.' + /' THIS IS A SINGULAR POINT OF THE' + ,' SPHERICAL COORDINATE SYSTEM.' + /' MOVE THIS NODE, AT LEAST SLIGHTLY.') CALL PAUSE() STOP ENDIF XI=(90.0-PLAT)*0.017453292 YI=PLON*0.017453292 ELEVI=VECTOR(4) QI=VECTOR(5) ZMI=VECTOR(6) TLI=VECTOR(7) IF (INDEX.LE.NREALN) THEN I=INDEX ELSE I=NREALN+INDEX-N1000 ENDIF CHECKN(I)=.TRUE. XNODE(I)=XI YNODE(I)=YI ELEV(I)=ELEVI DQDTDA(I)=QI IF (QI.LT.0.) THEN WRITE (IUNIT8,96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL PAUSE() STOP ENDIF IF (ZMI.LT.0.) THEN WRITE (IUNIT8,97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') CALL PAUSE() STOP ENDIF ZMNODE(I)=ZMI IF (TLI.LT.0.) THEN WRITE (IUNIT8,98) 98 FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS IS', + ' NON-PHYSICAL.') CALL PAUSE() STOP ENDIF TLNODE(I)=TLI IF (.NOT.BRIEF) THEN WRITE (IUNIT8,99) INDEX,PLON,PLAT,XI,YI,ELEVI, + QI,ZMI,TLI 99 FORMAT (' ',I10,0P,2F12.3,2F11.5,1P,3E10.2,E12.2) ENDIF 100 CONTINUE ALLOK=.TRUE. DO 101 I=1,NUMNOD ALLOK=ALLOK.AND.CHECKN(I) 101 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,102) 102 FORMAT(' THE FOLLOWING NODES WERE NEVER READ:') DO 104 I=1,NUMNOD IF (I.LE.NREALN) THEN INDEX=I ELSE INDEX=N1000+I-NREALN ENDIF IF (.NOT.CHECKN(I)) WRITE(IUNIT8,103)INDEX 103 FORMAT (' ',36X,I6) 104 CONTINUE CALL PAUSE() 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.') CALL PAUSE() 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, COUNTER', + 'CLOCKWISE'/ / + ' ELEMENT C1 C2 C3') DO 200 K=1,NUMEL C (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) C (READN METHOD IS USED TO ALLOW FOR POSSIBLE ELEMENT DATA.) CALL READN (INPUT,IUNIT7,IUNIT8,7, + OUTPUT,VECTOR) I=NINT(VECTOR(1)) IF ((I.LT.1).OR.(I.GT.NUMEL)) THEN WRITE (IUNIT8,111) I 111 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I6) CALL PAUSE() STOP ENDIF NODES(1,I)=NINT(VECTOR(2)) NODES(2,I)=NINT(VECTOR(3)) NODES(3,I)=NINT(VECTOR(4)) ELEDAT(1,I)=VECTOR(5) ELEDAT(2,I)=VECTOR(6) ELEDAT(3,I)=VECTOR(7) CHECKE(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,120) I,(NODES(J,I),J=1,3), + (ELEDAT(J,I),J=1,3) 120 FORMAT (' ',I6,':',3I10,1P,E8.1,0P,F7.1,1P,E8.1) DO 130 J=1,3 N=NODES(J,I) IF (N.GT.NREALN) N=NREALN+(N-N1000) IF ((N.LE.0).OR.(N.GT.NUMNOD)) THEN WRITE (IUNIT8,125) NODES(J,I) 125 FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.') CALL PAUSE() STOP ENDIF 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 CALL PAUSE() 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.') CALL PAUSE() 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,' GREAT CIRCLE FAULT ELEMENTS.') IF ((.NOT.BRIEF).AND.(NFL.GT.0)) WRITE(IUNIT8,231) 231 FORMAT (/' (THE 4 NODE NUMBERS DEFINING EACH ELEMENT MUST BE', + ' IN A COUNTERCLOCKWISE ORDER:'/ + ' N1, AND N2 ARE IN LEFT-TO-RIGHT SEQUENCE ON THE', + ' NEAR SIDE,'/ + ' THEN N3 IS OPPOSITE N2, N4 IS OPPOSITE N1 '/, + ' (FAULT DIPS ARE GIVEN AT N1, N2, ', + ' IN DEGREES FROM HORIZONTAL;'/ + ' POSITIVE DIPS ARE TOWARD N1, AND N2, RESPECTIVELY, '/ + ' WHILE NEGATIVE DIPS ARE TOWARD N4, AND N3.)'/ + ' (THE ARGUMENT OF THE FAULT TRACE IS GIVEN AT N1 AND N2,'/ + ' IN DEGREES COUNTERCLOCKWISE FROM THE THETA AXIS.)'/ + ' OFFSET IS THE TOTAL PAST SLIP OF THE FAULT.'/ / + ' ELEMENT N1 N2 N3 N4 DIP1 DIP2', + ' OFFSET'/) 240 FORMAT (' ',I6,':',4I5,1X,2F6.1,1X,F9.0) DO 300 K=1,NFL OFF=0. READ(IUNIT7,*) I,(NODEF(J,K),J=1,4),(DIPS(L),L=1,2),OFF IF ((I.LT.1).OR.(I.GT.NFL)) THEN WRITE (IUNIT8,241) I 241 FORMAT (/' ERR0R: ILLEGAL FAULT NUMBER: ',I6) CALL PAUSE() STOP ENDIF CHECKF(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,240) I,(NODEF(J,I),J=1,4), + (DIPS(L),L=1,2),OFF DO 250 J=1,4 N=NODEF(J,I) IF (N.GT.NREALN) N=NREALN+(N-N1000) IF ((N.LE.0).OR.(N.GT.NUMNOD)) THEN WRITE (IUNIT8,243) NODEF(J,I),I 243 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER (',I6, + ') IN FAULT ',I6) CALL PAUSE() STOP ENDIF NODEF(J,I)=N 250 CONTINUE DO 260 L=1,2 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 SIDE;'/ + ' A - PREFIX INDICATES A DIP TOWARD', + ' THE N4-N3 SIDE.)') CALL PAUSE() STOP ENDIF IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L) FDIP(L,I)=DIPS(L)*0.017453293 260 CONTINUE 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.') CALL PAUSE() 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 CALL PAUSE() 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 SUBROUTINE PUTNET (INPUT,IUNITO, + BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE) C C WRITES FINITE ELEMENT GRID TO UNIT "IUNITO". C CHARACTER*80 TITLE1 LOGICAL ANYEMU,BRIEF C DIMENSION DQDTDA(MXNODE),ELEDAT(3,MXEL),ELEV(MXNODE), + FDIP(2,MXFEL),NP(4), + NODEF(4,MXFEL),NODES(3,MXEL), + OFFSET(MXFEL),TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(2) C WRITE(*,"(/ /' READY TO CREATE OUTPUT .FEG FILE ON UNIT ',I3)") WRITE(*,"(' Please enter a (new) filename with .feg ending: ')") WRITE (IUNITO,1) TITLE1 1 FORMAT (A80) C WRITE (IUNITO,2) NUMNOD, NREALN,NFAKEN,N1000,BRIEF 2 FORMAT(4I8,L8,' (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)') C DO 100 I=1,NUMNOD IF (I.LE.NREALN) THEN IPRINT=I ELSE IPRINT=N1000+(I-NREALN) ENDIF PLAT=90.-XNODE(I)*57.29577951 PLON=YNODE(I)*57.29577951 WRITE (IUNITO,91) I,PLON,PLAT,ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 91 FORMAT (I6,0P,2F11.5,1P,4E10.2) 100 CONTINUE C WRITE (IUNITO,110) NUMEL 110 FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') C ANYEMU=.FALSE. DO 140 I=1,NUMEL DO 130 K=1,3 IF (ELEDAT(K,I).NE.0.0) THEN ANYEMU=.TRUE. GO TO 141 END IF 130 CONTINUE 140 CONTINUE 141 DO 200 I=1,NUMEL DO 150 K=1,3 IF (NODES(K,I).LE.NREALN) THEN NP(K)=NODES(K,I) ELSE NP(K)=N1000+(NODES(K,I)-NREALN) ENDIF 150 CONTINUE IF (ANYEMU) THEN WRITE (IUNITO,155) I,(NP(K),K=1,3),(ELEDAT(K,I),K=1,3) 155 FORMAT (I5,3I6,1P,E8.1,0P,F7.1,1P,E8.1) ELSE WRITE (IUNITO,160) I,(NP(K),K=1,3) 160 FORMAT (I5,3I6) END IF 200 CONTINUE C WRITE (IUNITO,210) NFL 210 FORMAT (I10,' (NFL = NUMBER OF CURVILINEAR FAULT ELEMENTS)') C DO 300 I=1,NFL DO 220 K=1,4 IF (NODEF(K,I).LE.NREALN) THEN NP(K)=NODEF(K,I) ELSE NP(K)=N1000+(NODEF(K,I)-NREALN) ENDIF 220 CONTINUE DO 230 K=1,2 DIPS(K)=FDIP(K,I) DIPS(K)=DIPS(K)*57.2958 IF (DIPS(K).GT.90.01) DIPS(K)=DIPS(K)-180. 230 CONTINUE WRITE (IUNITO,250) I,(NP(K),K=1,4), + (DIPS(K),K=1,2),OFFSET(I) 250 FORMAT (I5,4I6,2F6.1,1P,E10.2) 300 CONTINUE C RETURN END C C C SUBROUTINE PAUSE() IMPLICIT NONE WRITE(*,"(' Press [Enter]..'\)") READ(*,*) RETURN END