C-------------------------------------------------------------------- C PROGRAM FillGrid 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 DrawGrid.exe, C and (if needed) grids of elevation and heat-flow with C regular spacing of points in (longitude, latitude) space, C and either: C * a crustal-thickness grid (for FAULTS .feg files), C with regular spacing of points in (lon,lat); OR C * a parameter file (for PLATES .feg files) C with 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; and (for PLATES .feg files only): C -mantle lithosphere thickness (not including crust). C All units are SI: m, W m**2 (NOT mW m**2), m, m. 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 WHEN PROCESSING .FEG FILES FOR USE WITH -PLATES-, C HEAT-FLOW (FROM ANY SOURCE) IS SUBJECT TO A MINIMUM VALUE C DESIGNED TO GUARANTEE EXISTENCE OF A SOLUTION; AND C ISOSTASY WITH RESPECT TO A SPREADING RISE C IS ACHIEVED BY ADJUSTING THE LAYER THICKNESSES, C LEAVING ELEVATION AND HEAT-FLOW ALONE. 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") INTEGER, PARAMETER :: MAXNOD = 5000 C C MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES). INTEGER, PARAMETER :: MAXEL = 10000 C C MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS); INTEGER, PARAMETER :: MAXFEL = 2000 C C MAXDAT = MAXIMUM NUMBER OF ROWS AND COLUMNS IN 2-D ARRAYS OF C ELEVATION, HEAT-FLOW DATA, AND CRUSTAL-THICKNESS DATA INTEGER, PARAMETER :: MAXDAT = 900 C C MAXBN = MAXIMUM NUMBER OF BOUNDARY NODES (BOTH "REAL" AND "FAKE"). INTEGER, PARAMETER :: MAXBN = 162 C C MAXATP = MAXIMUM NUMBER OF NODES WHICH MAY OVERLAP AT A FAULT- C INTERSECTION POINT. INTEGER, PARAMETER :: MAXATP = 20 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*80 RECORD,TITLE1,TITLE3 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 INTEGER I,IMODE,IROW,IUNITL,IUNITT, + JCOL,IZONE,LIST,MAPTYP, + MAXITR,MXBN,MXDATA,MXEL,MXFEL,MXNODE,MXSTAR, + NCOND,NCX,NCY,NEX,NEY,NODCON,NODTYP,NQX,NQY, + NEMAX,NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,N1000 C LOGICAL BRIEF,CHECKE,CHECKF,CHECKN, + DIMERR,DOTHIS,EDGEFS,EDGETS,EVERYP, + FAULTS,INPUT,LLIMIT,MODIFY, + more_feg,more_grd, + NEEDC,NEEDE,NEEDQ, + OUTPUT,PLATES,USEALI,WORK C REAL ACREEP,ALPHAT,AREA, + BCREEP,BIOT,BY,BYERLY, + CARRAY,CDX,CDY,CX1,CX2,CY1,CY2,CENLON, + CCREEP,CFRIC,CLIMIT,CONDUC,CPNLAT, + DCREEP,DETJ,DIPMAX,DQDTDA,DXS,DYS, + EARRAY,EDX,EDY,EX1,EX2,EY1,EY2, + EASTOF,ECREEP,ELEV,ELEVAT, + FAZ,FDIP,FFRIC,FLEN,FTAN, + GMEAN, + HCMAX,HEATFL,HLMAX, + LRANGE,LX, + NLAT,NLAT0,NLAT1, + OFFMAX,OFFSET,OKDELV,OKTOQT,ONEKM, + PLAT,PLON, + QARRAY,QDX,QDY,QX1,QX2,QY1,QY2, + QLIMIT, + RADIO,RADIUS,REFSTR,RHOAST,RHOBAR,RHOH2O,RX, + TASTHC,TASTHK,TAUMAX,TEMLIM,THICKC,THICKM, + TLNODE,TRHMAX,TSURF,TY, + WEDGE, + XNODE,XOF500,X0ELON, + YNODE,Y0NLAT, + ZMNODE C C--------------------------------------------------------------------- C DIMENSION STATMENTS C C DIMENSIONS USING PARAMETER MAXNOD: DIMENSION CHECKN(MAXNOD), DQDTDA(MAXNOD), + ELEV (MAXNOD), NODTYP(MAXNOD), + TLNODE(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), + EDGETS(3,MAXEL), NODES(6,MAXEL) C C DIMENSIONS USING PARAMETER MAXFEL: DIMENSION CHECKF (MAXFEL), EDGEFS(2,MAXFEL), + FAZ (2,MAXFEL), FDIP (3,MAXFEL), + FLEN (MAXFEL), FTAN (7,MAXFEL), + NODEF(6,MAXFEL), + OFFSET (MAXFEL) C C DIMENSIONS USING PARAMETER MAXDAT: DIMENSION CARRAY (MAXDAT,MAXDAT), + EARRAY (MAXDAT,MAXDAT), + QARRAY (MAXDAT,MAXDAT) C C DIMENSIONS USING PARAMETER MAXBN: DIMENSION NODCON (MAXBN) C C DIMENSIONS USING PARAMETER MAXATP: DIMENSION LIST (MAXATP) 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 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 (DATA ABOUT TRIANGULAR ELEMENTS): COMMON /S1S2S3/ POINTS COMMON /PHITAB/ PHI COMMON /WGTVEC/ WEIGHT C C ENTRIES CORRESPONDING TO BD2 (DATA ABOUT FAULT ELEMENTS): COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS C C------------------------------------------------------------------- C DATA STATEMENTS 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. A CAUTIONARY NOTICE WILL BE OUTPUT. DATA DIPMAX /75./ C C "QLIMIT" IS A LOWER LIMIT ON HEAT-FLOW DESIGNED TO GUARUNTEE C EXISTANCE OF A SOLUTION FOR LAYER THICKNESSES: DATA QLIMIT /0.043/ C C "CLIMIT" IS A LOWER LIMIT ON CRUSTAL THICKNESS: DATA CLIMIT /5000./ C C "HCMAX" IS AN UPPER LIMIT ON CRUSTAL THICKNESS: DATA HCMAX /75000./ C C "HLMAX" IS AN UPPER LIMIT ON TOTAL LITHOSPHERE THICKNESS: DATA HLMAX /200000./ 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 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 MXEL =MAXEL MXFEL =MAXFEL MXBN =MAXBN MXSTAR=MAXATP MXDATA=MAXDAT C WEDGE=ABS(90.-ABS(DIPMAX))*0.017453293 C C LITHOSPHERE/ASTHENOSPHERE BOUNDARY TEMPERATURE, IN C OR K, IS: TASTHC = 1250.0 TASTHK = TASTHC + 273.0 C BUT NOTE THAT THIS IS ONLY USED IN FILLING PLATES .FEG FILES. C C--------------------------------------------------------------------- C C INTRODUCTION / HEADER SECTION: C WRITE (IUNITT,1) 1 FORMAT ( + ' -------------------------------------------------------------' + /' Program FillGrid (version of 13 March 2000)' +//' by Peter Bird, Dept. of Earth & Space Sciences,' + /' University of California, Los Angeles, CA 90095-1567.' +//' Reads a finite element grid produced by DrawGrid.exe,' + /' and (if needed) grids of elevation and heat-flow with' + /' regular spacing of points in (longitude, latitude) space,' + /' and either:' + /' * a crustal-thickness grid (for FAULTS .feg files), OR' + /' * a parameter file (for PLATES .feg files)' + /' 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 (for PLATES .feg only).' + /' All units must be SI: m, W/m**2 (NOT mW/m**2), m, m.') WRITE (IUNITT,*) CALL Prompt_for_Logical( + "Do you want information about .feg files?",.TRUE.,more_feg) IF (more_feg) THEN WRITE (*,"( +//' ------------------------------------------------------------' +,'----------' + /' About Finite Element Grid (.feg) Files' +//' Finite element grids may be either spherical-Earth or flat-E' +,'arth.' + /' Spherical-Earth grids used with SHELLS are produced by OrbWe' +,'ave' + /' and have 3-node spherical triangles and 4-node great-circl' +,'e faults.' + /' Flat-Earth grids used with PLATES or FAULTS are produced by ' +,'DrawGrid' + /' and have 6-node isoparametric triangles and 6-node curved ' +,'faults.' + /' Flat-Earth grids used by LARAMY are created inside that prog' +,'ram' + /' and have only 6-node isoparametric triangles (no faults).' + /' Any .feg file contains nodal locations and nodal data' + /' (elevation, heat-flow, crustal thickness, lithosphere thi' +,'ckness).' + /' They also contain topological connections of the nodes to' +,' form' + /' elements. If there are faults, their dips are specified.' +//' ------------------------------------------------------------' +,'----------')") END IF CALL Prompt_for_Logical( + "Do you want information about .grd files?",.TRUE.,more_grd) IF (more_grd) THEN WRITE (*,"( +//' -------------------------------------------------------------' +,'---------' + /' About Gridded Data (.grd) Files' +//' These files contain scalar data values on a regular rectangul' +,'ar grid,' + /' defined in (longitude,latitude) space.' +//' The first line has 3 numbers: lon_min, d_lon, lon_max;' + /' the 2nd also has 3 numbers: lat_min, d_lat, lat_max.')") WRITE (*,"( + ' All of these numbers are in decimal degrees, with N = +,' + /' S = -, E = +, W = -.' + /' Following lines give the gridded data in text order, i.e., ' +,'beginning' + /' with the NorthWest corner, going W->E along the top row, th' +,'en W->E ' + /' along the 2nd row, etc. The number and position of line br' +,'eaks' + /' is not important in this part of the file.' +//' -------------------------------------------------------------' +,'---------')") END IF WRITE(IUNITT,19) 19 FORMAT (/' ENTER 1 for FAULTS mode or 2 for PLATES mode: '\) READ(*,*) IMODE FAULTS=(IMODE.EQ.1) PLATES=(IMODE.EQ.2) IF (.NOT.(FAULTS.OR.PLATES)) STOP 'Illegal choice.' WRITE (IUNITT,20) 20 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.') IF (FAULTS) WRITE (IUNITT,26) 26 FORMAT ( + ' If a nodal crustal thickness is non-zero, it will be left' + /' unchanged, to preserve effects of hand-editing.' + /' If any crustal-thickness is zero, a new values will be' + /' interpolated from a gridded dataset.') IF (PLATES) WRITE (IUNITT,28) TASTHC,QLIMIT 28 FORMAT ( + /' 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' + /' and is currently set to TASTHC = ',F6.1,' Centigrade.' +//' Heat-flow (from any source) is subject to a minimum value' + /' (',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.') WRITE (IUNITT,30) 30 FORMAT ( + ' ------------------------------------------------------------' +//' PRESS [Enter] to continue...'\) READ(*,"(A)") C1 WRITE(*,31) 31 FORMAT(/ /' It is possible to limit the update of nodal values' + /' to a rectangular region of (x, y)' + /' Do you wish to use this option? (T/F): '\) READ(*,*) LLIMIT IF (LLIMIT) THEN 36 WRITE(*,"(' Enter left-side (minimum) x: '\)") READ(*,*) LX WRITE(*,"(' Enter right-side (maximum) x: '\)") READ(*,*) RX IF (RX.LE.LX) THEN WRITE (*,"(' ERROR: Maximum x must ', + 'be >= minimum.')") GO TO 36 END IF 37 WRITE(*,"(' Enter bottom (minimum) y: '\)") READ(*,*) BY WRITE(*,"(' Enter top (maximum) y: '\)") READ(*,*) TY IF (TY.LE.BY) THEN WRITE (*,"(' ERROR: Maximum y must ', + 'be >= minimum.')") GO TO 37 END IF END IF C IF (PLATES) WRITE (IUNITT,40) QLIMIT, CLIMIT, HCMAX, HLMAX 40 FORMAT(/' THE FOLLOWING LIMITS APPLY IN THIS RUN:' +/' LOWER 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 IF (PLATES) THEN C C READ PARAMETERS ON UNIT 1 C CALL READPM (INPUT, 1, 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) END IF C C READ FINITE ELEMENT GRID ON UNIT 2 C CALL GETNET (INPUT, 2,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) WRITE (IUNITT,50) 50 FORMAT (/' SUCCESSFULLY READ FE GRID.') C C CHECK GRID TOPOLOGY AND COMPUTE GEOMETRIC PROPERTIES C WRITE (IUNITT,51) 51 FORMAT (/' CHECKING TOPOLOGY OF FE GRID...') CALL SQUARE (INPUT,BRIEF,FDIP,IUNITT, + MXBN,MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,WEDGE, + MODIFY,FAZ,XNODE,YNODE, + OUTPUT,AREA,DETJ,DXS,DYS,EDGEFS,EDGETS, + FLEN,FTAN,NCOND,NODCON, + WORK,CHECKN,LIST,NODTYP) WRITE (IUNITT,52) 52 FORMAT (' TOPOLOGY CHECKING COMPLETED.') C C READ IN ELEVATION ARRAY ON UNIT 3, IF NEEDED C NEEDE=.FALSE. DO 60 I=1,NUMNOD IF (LLIMIT) THEN IF ((XNODE(I).GE.LX).AND. + (XNODE(I).LE.RX).AND. + (YNODE(I).GE.BY).AND. + (YNODE(I).LE.TY)) 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 IF ((NEX.GT.MXDATA).OR.(NEY.GT.MXDATA)) THEN NEMAX=MAX(NEX,NEY) WRITE (IUNITT,65) NEMAX 65 FORMAT (/' INCREASE PARAMETER MAXDAT TO ',I5, + /' AND RECOMPILE.') STOP ENDIF READ (3,*) ((EARRAY(IROW,JCOL),JCOL=1,NEX),IROW=1,NEY) ENDIF C C READ IN HEAT-FLOW ARRAY ON UNIT 4, IF NEEDED C NEEDQ=.FALSE. DO 70 I=1,NUMNOD IF (LLIMIT) THEN IF ((XNODE(I).GE.LX).AND. + (XNODE(I).LE.RX).AND. + (YNODE(I).GE.BY).AND. + (YNODE(I).LE.TY)) 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 IF ((NQX.GT.MXDATA).OR.(NQY.GT.MXDATA)) THEN NEMAX=MAX(NQX,NQY) WRITE (IUNITT,65) NEMAX STOP ENDIF READ (4,*) ((QARRAY(IROW,JCOL),JCOL=1,NQX),IROW=1,NQY) CLOSE(4) ELSE C C IMPOSE MINIMUM HEAT-FLOW AT NODES, TO PREVENT UNREASONABLY C THICK AND STIFF LITHOSPHERE ANYWHERE C DO 80 I=1,NUMNOD IF (DQDTDA(I).NE.0.0) DQDTDA(I)=MAX(DQDTDA(I),QLIMIT) 80 CONTINUE END IF C C READ IN CRUSTAL-THICKNESS ARRAY ON UNIT 4, IF NEEDED C NEEDC=.FALSE. IF (FAULTS) THEN DO 81 I=1,NUMNOD IF (LLIMIT) THEN IF ((XNODE(I).GE.LX).AND. + (XNODE(I).LE.RX).AND. + (YNODE(I).GE.BY).AND. + (YNODE(I).LE.TY)) THEN IF (ZMNODE(I).EQ.0.0) NEEDC=.TRUE. END IF ELSE IF (ZMNODE(I).EQ.0.0) NEEDC=.TRUE. END IF 81 CONTINUE END IF IF (NEEDC) THEN WRITE(*,82) 82 FORMAT(/ /' ATTEMPTING TO READ GRIDDED CRUSTAL-THICKNESS:'/) READ (4,*) CX1,CDX,CX2 READ (4,*) CY1,CDY,CY2 NCX=(CX2-CX1)/CDX + 1.5 NCY=(CY2-CY1)/CDY + 1.5 IF ((NCX.GT.MXDATA).OR.(NCY.GT.MXDATA)) THEN NEMAX=MAX(NCX,NCY) WRITE (IUNITT,65) NEMAX STOP ENDIF READ (4,*) ((CARRAY(IROW,JCOL),JCOL=1,NCX),IROW=1,NCY) END IF C C DEFINE THE MAP PROJECTION: C 89 WRITE (IUNITT,90) 90 FORMAT (/' Choose the map projection type' + /' (for conversions between (x,y) and (lon,lat))' + /' from the following list:' + /' 1: conic projection' + /' 2: Universal Transverse Mercator projection' + /' Please enter the integer for your choice: '\) READ (*,*) MAPTYP IF ((MAPTYP.LT.1).OR.(MAPTYP.GT.2)) THEN WRITE (IUNITT,"(' ERROR: Illegal choice.')") GO TO 89 END IF IF (MAPTYP.EQ.1) THEN C conic projection: WRITE (IUNITT,91) 91 FORMAT (/' Please enter the following parameters to define' + /' the conic projection that relates (x,y) to' + /' (longitude, latitude):' + /' Radius of the planet, in m: '\) READ (*,*) RADIUS WRITE (IUNITT,92) 92 FORMAT ( ' Latitude of the tangent parallel, in degrees' + /' (N = +, S = -): '\) READ (*,*) CPNLAT WRITE (IUNITT,93) 93 FORMAT ( ' Latitude of the (x,y) origin, in degrees' + /' (N = +, S = -): '\) READ (*,*) Y0NLAT WRITE (IUNITT,94) 94 FORMAT ( ' Longitude of the (x,y) origin, in degrees' + /' (E = +, W = -): '\) READ (*,*) X0ELON ELSE IF (MAPTYP.EQ.2) THEN C Universal Transverse Mercator WRITE (IUNITT,95) 95 FORMAT (/' UTM zones are 6 degrees wide, and the first zone' + /' extends from -180E to -174E, centered at -177E.' + /' Please enter the number of your zone: '\) READ (*,*) IZONE CENLON=-183.0+6.0*IZONE WRITE (IUNITT,96) CENLON 96 FORMAT (' Central meridian will be at ',F7.2,' degrees E.') WRITE (IUNITT,97) 97 FORMAT(/' Please enter the radius of the Earth, ' + /' in the same units (m?) used in your .feg file: '\) READ (*,*) RADIUS WRITE (IUNITT,98) 98 FORMAT (/' The UTM formula assumes that X = 500 km at the' + /' central meridian for this zone. Please convert' + /' the distance of 500 km to the units (m?) that' + /' are used in your .feg file: '\) READ (*,*) XOF500 END IF C C INITIATE LOG FILE FOR -ASSIGN- OUTPUT: C WRITE (IUNITT,110) IUNITL 110 FORMAT (/ /' ATTEMPTING TO CREATE DETAILED LOG FILE ON UNIT', + I3/) WRITE (IUNITL,120) 120 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...') DO 680 I=1,NUMNOD CALL XYTOLL (INPUT,XNODE(I),YNODE(I), + MAPTYP, + RADIUS,CPNLAT,Y0NLAT,X0ELON, + IZONE,XOF500, + OUTPUT,PLAT,PLON) ELEVAT=ELEV(I) HEATFL=DQDTDA(I) THICKC=ZMNODE(I) THICKM=TLNODE(I) IF (LLIMIT) THEN DOTHIS=((XNODE(I).GE.LX).AND. + (XNODE(I).LE.RX).AND. + (YNODE(I).GE.BY).AND. + (YNODE(I).LE.TY)) ELSE DOTHIS=.TRUE. END IF IF (DOTHIS) THEN CALL ASSIGN (INPUT,FAULTS,PLATES, + ALPHAT, + CARRAY,CX1,CDX,CX2,NCX,CDY,CY2,NCY, + CLIMIT,CONDUC, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT,MXDATA, + ONEKM, + PLON,PLAT,QLIMIT, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + TASTHK,TEMLIM,TSURF, + MODIFY,ELEVAT,HEATFL,THICKC,THICKM) ELEV(I)=ELEVAT DQDTDA(I)=HEATFL ZMNODE(I)=MAX(THICKC,CLIMIT) TLNODE(I)=MAX(THICKM,0.) IF (FAULTS) THEN WRITE (IUNITL,671) I,ELEV(I),DQDTDA(I),ZMNODE(I) 671 FORMAT (' ',I10,1P,3E10.2) ELSE IF (PLATES) THEN WRITE (IUNITL,672) I,ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 672 FORMAT (' ',I10,1P,4E10.2) END IF 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,FAULTS,PLATES, + BRIEF,DQDTDA,ELEV,FAZ,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE) WRITE (IUNITT,800) 800 FORMAT (/' Remember to edit the title line of ', + 'this new .feg file!' + /' Press [Enter] to finish...'\) READ (*,"(A)") C1 STOP ' ' END C C C SUBROUTINE AREAS (INPUT,NODES,NUMEL,NUMNOD,XNODE,YNODE, + 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 INTEGER I,I1,I2,I3,NODES,NUMEL,NUMNOD LOGICAL INPUT,OUTPUT REAL AREA,XNODE,YNODE 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 SUBROUTINE AREAS C C C SUBROUTINE ASSIGN (INPUT,FAULTS,PLATES, + ALPHAT, + CARRAY,CX1,CDX,CX2,NCX,CDY,CY2,NCY, + CLIMIT,CONDUC, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT,MXDATA, + ONEKM, + PLON,PLAT,QLIMIT, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + TASTHK,TEMLIM,TSURF, + MODIFY,ELEVAT,HEATFL,THICKC,THICKM) C C DETERMINES ELEVATION (IF 0.0), HEAT FLOW (IF 0.0), AND EITHER: C IF (FAULTS) THEN C DETERMINES CRUSTAL THICKNESS (IF 0.0) FROM A GRID; C ELSE IF (PLATES) THEN C DETERMINES THE ISOSTATIC C CRUSTAL AND MANTLE LITHOSPHERE THICKNESSES C END IF C FOR ONE POINT: C (PLON = EAST LONGITUDE, PLAT = NORTH LATITUDE), BOTH IN DEGREES. C 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 C NUMBER OF ITERATIONS PERMITTED: INTEGER, PARAMETER :: MAXITR = 25 C DAMPING FACTOR (<1.) FOR STABILITY OF ITERATIONS: REAL, PARAMETER :: DAMP = 0.300 C INTEGER I,IC1,IC2,IR1,IR2,ITER,IUNITL,IUNITT, + MXDATA, + NCX,NCY,NEX,NEY,NQX,NQY LOGICAL BADP,BADT,DIVERG,FAULTS,INPUT,MODIFY, + NEEDC,NEEDE,NEEDQ,OUTPUT,OUTSID,PLATES, + WARNC1, WARNC2, WARNM1, WARNL2, WAYOUT REAL ALPHAT,BOT,CARRAY,CDX,CDY,CLIMIT,CONDUC,COSDEG,CX1,CX2,CY2, + DC,DC24,DC25,DM24,DM25,DTL, + EARRAY,EDX,EDY,ELEVAT,EX1,EX2,EY2,FC,FR, + GEOTH1,GEOTH2,GEOTH3,GEOTH4,GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN,HC,HCMAX,HEATFL,HLMAX,HM, + QARRAY,QDX,QDY,QLIMIT,QRED,QX1,QX2,QY2, + ONEKM,PLAT,PLON,RADIO,RHOAST,RHOBAR,RHOH2O, + SIGZZB,SINDEG, + TA,TANDEG,TASTHK,TAUZZ,TEMLIM,TERR0R,TEST, + THICKC,THICKM,TM,TMOHO,TOP,TSURF C DIMENSION ALPHAT(2),CONDUC(2),RADIO(2),RHOBAR(2),TEMLIM(2) DIMENSION CARRAY(MXDATA,MXDATA),EARRAY(MXDATA,MXDATA), + QARRAY(MXDATA,MXDATA) C LOCAL STORAGE FOR HISTORY OF ITERATION: DIMENSION HC(MAXITR),HM(MAXITR),TM(MAXITR),TA(MAXITR) C--------------------------------------------------------------------- C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) 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 IC1=((PLON-EX1)/EDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NEX-1) IC2=IC1+1 FC=(PLON-(EX1+EDX*(IC1-1)))/EDX 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 IC1=((PLON-QX1)/QDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NQX-1) IC2=IC1+1 FC=(PLON-(QX1+QDX*(IC1-1)))/QDX 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) ENDIF C C DETERMINE CRUSTAL THICKNESS, IF NEEDED C NEEDC=(FAULTS.AND.(THICKC.EQ.0.0)) IF (NEEDC) THEN IR1=((CY2-PLAT)/CDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NCY-1) IR2=IR1+1 FR=((CY2-CDY*(IR1-1))-PLAT)/CDY IC1=((PLON-CX1)/CDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NCX-1) IC2=IC1+1 FC=(PLON-(CX1+CDX*(IC1-1)))/CDX 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=CARRAY(IR1,IC1)+FC*(CARRAY(IR1,IC2)-CARRAY(IR1,IC1)) BOT=CARRAY(IR2,IC1)+FC*(CARRAY(IR2,IC2)-CARRAY(IR2,IC1)) THICKC=TOP+FR*(BOT-TOP) ENDIF 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).') STOP END IF C IF (FAULTS) RETURN C C PROVIDE MINIMUM REASONABLE Q TO GUARANTEE A STEADY-STATE SOLUTION: C HEATFL=MAX(HEATFL,QLIMIT) 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.100.) 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 (' ', I5, 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 RETURN END SUBROUTINE ASSIGN C C C REAL FUNCTION ATAN2F (Y,X) C C CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS C REAL X,Y IF ((X.NE.0.).OR.(Y.NE.0.)) THEN ATAN2F=ATAN2(Y,X) ELSE ATAN2F=0. END IF RETURN END FUNCTION ATAN2F 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 BLOCK DATA BD1 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 BLOCK DATA BD2 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 INTEGER I,J,M,NODE,NODES,NUMEL,NUMNOD LOGICAL INPUT,OUTPUT REAL AI2,AJ11,AJ11S,AJ12,AJ21,AJ22,AREA, + B,C,DETJ,DETJAC,DN,DXS,DYS,X,XNODE,Y,YNODE 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 SUBROUTINE DERIV 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 INTEGER I,IFN,INDEX,IOS,IUNIT7,IUNIT8,IUNITP, + J,K,L, + MXEL,MXFEL,MXNODE, + N,NBASE,NFAKEN,NFL, + NODEF,NODES,NREALN,NTOP,NUMEL,NUMNOD,N1000 LOGICAL ALLOK,BRIEF,CHECKE,CHECKF,CHECKN,INPUT,OUTPUT,WORK REAL AZ1,AZ3, + DIPS,DQDTDA,ELEV,ELEVI,FAZ,FDIP, + OFF,QI,OFFMAX,OFFSET, + TLI,TLNODE,VECTOR, + XI,XNODE,YI,YNODE,ZMI,ZMNODE 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/) TITLE1=' '// + ' ' READ (IUNIT7,2,IOSTAT=IOS) 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 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 SUBROUTINE GETNET 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 INTEGER I,J,K,KELE,KFAULT,L,M1,M2,M3,M4,M5,M6,MXEL,MXFEL, + N1,N2,N3,NFL,NODEF,NODES,NUMEL LOGICAL FOUNDE,FOUNDF,INPUT,OUTPUT DIMENSION NODEF(6,MXFEL),NODES(6,MXEL) C C THREE NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE: N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J, 3)+4,I) N3=NODES(MOD(J+1,3)+1,I) C CHECK FOR ADJACENT FAULT ELEMENT FIRST: FOUNDF=.FALSE. DO 10 K=1,NFL M1=NODEF(1,K) M2=NODEF(2,K) M3=NODEF(3,K) M4=NODEF(4,K) M5=NODEF(5,K) M6=NODEF(6,K) IF (((M1.EQ.N3).AND.(M2.EQ.N2).AND.(M3.EQ.N1)).OR. + ((M4.EQ.N3).AND.(M5.EQ.N2).AND.(M6.EQ.N1))) THEN FOUNDF=.TRUE. KFAULT=K GO TO 11 ENDIF 10 CONTINUE 11 IF (.NOT.FOUNDF) KFAULT=0 C IF THERE WAS A FAULT, REPLACE 3 NODE NUMBERS THAT WE SEARCH FOR: IF (FOUNDF) THEN IF (M2.EQ.N2) THEN N1=M4 N2=M5 N3=M6 ELSE N1=M1 N2=M2 N3=M3 ENDIF ENDIF C SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT: FOUNDE=.FALSE. DO 20 K=1,NUMEL IF (K.NE.I) THEN DO 15 L=1,3 M1=NODES(MOD(L, 3)+1,K) M2=NODES(MOD(L, 3)+4,K) M3=NODES(MOD(L+1,3)+1,K) IF ((M3.EQ.N1).AND.(M2.EQ.N2).AND.(M1.EQ.N3)) THEN FOUNDE=.TRUE. KELE=K GO TO 21 ENDIF 15 CONTINUE ENDIF 20 CONTINUE 21 IF (.NOT.FOUNDE) KELE=0 RETURN END SUBROUTINE NEXT C C C SUBROUTINE Prompt_for_Logical (prompt_text, default, answer) C Writes a line to the default (*) unit, with: C "prompt_text" ["default"]: C and accepts an answer with a logical value. C If [Enter] is pressed with nothing preceding it, C then "answer" takes the value "default". C Note that prompt_text should usually end with '?'. C It can be more than 70 bytes long if necessary, but cannot C contain line breaks, tabs, or other formatting. C Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") + prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE(*, +"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE(*, +"(' (Only the first letter of your answer is used.)')") WRITE(*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE(*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE(*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE Prompt_for_Logical C C C SUBROUTINE PUTNET (INPUT,IUNITO,FAULTS,PLATES, + BRIEF,DQDTDA,ELEV,FAZ,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 INTEGER I,IPRINT,IUNITO,J,K,KM,KP,MXEL,MXFEL,MXNODE, + N1000,NFAKEN,NFL,NODEF,NODES,NP,NREALN,NUMEL,NUMNOD LOGICAL BRIEF,FAULTS,INPUT,PLATES REAL ARGS,AZ1,AZ3,DIPS,DQDTDA,ELEV,FAZ,FDIP, + OFF2,OFFSET,SIDE2,TLNODE, + XC,XM,XNODE,XP,YC,YM,YNODE,YP,ZMNODE C DIMENSION DQDTDA(MXNODE),ELEV(MXNODE), + FAZ(2,MXFEL),FDIP(3,MXFEL),NP(6), + NODEF(6,MXFEL),NODES(6,MXEL), + OFFSET(MXFEL),TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(3),ARGS(3) C C RETURN COORDINATES OF STRAIGHT MID-POINT NODES TO ZERO C DO 100 I=1,NUMNOD DO 90 J=1,NUMEL DO 80 K=4,6 IF (I.EQ.NODES(K,J)) THEN KM=K-3 XM=XNODE(NODES(KM,J)) YM=YNODE(NODES(KM,J)) KP=MOD(K-3,3)+1 XP=XNODE(NODES(KP,J)) YP=YNODE(NODES(KP,J)) SIDE2=(XP-XM)**2+(YP-YM)**2 XC=(XM+XP)*0.5 YC=(YM+YP)*0.5 OFF2=(XNODE(I)-XC)**2+(YNODE(I)-YC)**2 IF (OFF2/SIDE2.LE.0.0001) THEN XNODE(I)=0.0 YNODE(I)=0.0 END IF GO TO 99 END IF 80 CONTINUE 90 CONTINUE 99 CONTINUE 100 CONTINUE C WRITE(*,"(/ /' READY TO CREATE OUTPUT .FEG FILE ON UNIT ',I3)") WRITE(*,"(' Please enter a (new) filename with .feg ending: ')") WRITE (IUNITO,101) TITLE1 101 FORMAT (A80) C WRITE (IUNITO,102) NUMNOD, NREALN,NFAKEN,N1000,BRIEF 102 FORMAT(4I8,L8,' (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)') C DO 108 I=1,NUMNOD IF (I.LE.NREALN) THEN IPRINT=I ELSE IPRINT=N1000+(I-NREALN) ENDIF IF (FAULTS) THEN WRITE (IUNITO,105) I,XNODE(I),YNODE(I), + ELEV(I),DQDTDA(I),ZMNODE(I) 105 FORMAT (I5,1P,2E13.5,3E10.2) END IF IF (PLATES) THEN WRITE (IUNITO,106) I,XNODE(I),YNODE(I), + ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 106 FORMAT (I5,1P,2E13.5,4E10.2) END IF 108 CONTINUE C WRITE (IUNITO,110) NUMEL 110 FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') C DO 200 I=1,NUMEL DO 150 K=1,6 IF (NODES(K,I).LE.NREALN) THEN NP(K)=NODES(K,I) ELSE NP(K)=N1000+(NODES(K,I)-NREALN) ENDIF 150 CONTINUE WRITE (IUNITO,160) I,(NP(K),K=1,6) 160 FORMAT (I5,6I5) 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,6 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,3 DIPS(K)=FDIP(K,I) DIPS(K)=DIPS(K)*57.2958 IF (DIPS(K).GT.90.01) DIPS(K)=DIPS(K)-180. 230 CONTINUE AZ1=FAZ(1,I)*57.2958 AZ3=FAZ(2,I)*57.2958 WRITE (IUNITO,250) I,(NP(K),K=1,6), + (DIPS(K),K=1,3),AZ1,AZ3,OFFSET(I) 250 FORMAT (I5,6I6,3F6.1,2F7.1,1P,E10.2) 300 CONTINUE C RETURN END SUBROUTINE PUTNET 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 INTEGER I,IUNITP,IUNITT,N,NUMBER LOGICAL ANYIN,DOTTED,EXPON,INPUT,OUTPUT,SIGNED REAL VECTOR 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 SUBROUTINE READN 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*80 TITLE3 INTEGER IOS,IUNIT7,IUNIT8,MAXITR REAL ACREEP,ALPHAT,BCREEP,BIOT,BYERLY,CCREEP,CFRIC,CONDUC, + DCREEP,ECREEP,FFRIC,GMEAN,OFFMAX,OKDELV,OKTOQT,ONEKM, + RADIO,REFSTR,RHOAST,RHOBAR,RHOH2O,TEMLIM,TRHMAX,TSURF LOGICAL EVERYP,INPUT,OUTPUT C 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/) TITLE3=' '// + ' ' READ (IUNIT7,2,IOSTAT=IOS) 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 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,198) 198 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 SUBROUTINE READPM C C C SUBROUTINE SQUARE (INPUT,BRIEF,FDIP,IUNIT8, + MXBN,MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,WEDGE, + MODIFY,FAZ,XNODE,YNODE, + OUTPUT,AREA,DETJ,DXS,DYS,EDGEFS,EDGETS, + FLEN,FTAN,NCOND,NODCON, + WORK,CHECKN,LIST,NODTYP) C C CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID C CHARACTER*21 OBLIQU,TAG1,TAG2,TAG3,VERTIC INTEGER I,I1,I2,I3,I4,I5,I6,IO,ITYPE,IUNIT8, + J,J1,J2,JM,JP, + K,KELE,KFAULT, + L,LIST,LM,LP, + M,MXBN,MXEL,MXFEL,MXNODE,MXSTAR, + N,N1,N2,N3,N6,NAZI,NAZL,NCOND,NDIP,NDONE,NFL, + NINSUM,NJ1,NJ2,NL1,NL3,NL4,NL6,NLEFT, + NODCON,NODE,NODE1,NODE6,NODEF,NODES,NODTYP,NP1,NP6, + NREALN,NTOFIX,NUMBER,NUMEL,NUMNOD,NVPART,N1000 LOGICAL AGREED,ALLOK,BRIEF,CHECKE,CHECKF,CHECKN, + EDGEFS,EDGETS,FOUND,INPUT,MATCH,MODIFY, + OUTPUT,SWITCH,VERT1,VERT2,VERT3,WORK REAL AREA,AZ,AZI,AZIMUT,AZL, + DAZ,DAZI,DAZL,DELD,DELD1,DELD2,DELD3,DETJ, + DF1DS,DF2DS,DF3DS, + DIP1,DIP2,DIP3,DX,DXDS,DXS,DY,DYDS,DYS, + F1,F2,F3,FACTOR,FAZ,FDIP,FLEN,FTAN, + OLDX,OLDY, + PARRAL,PERPEN,PHI1,PHI2, + R,R2,RMAX, + S,SHORT, + T1,T2,TEST,TOLER,WEDGE, + X,X1,X2,X3,XMEAN,XNODE,XSUM, + Y,Y1,Y2,Y3,YMEAN,YNODE,YSUM C 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), + EDGEFS(2,MXFEL),EDGETS(3,MXEL),FDIP(3,MXFEL), + FAZ(2,MXFEL),FLEN(MXFEL),FTAN(7,MXFEL), + LIST(MXSTAR),NODCON(MXBN), + 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 (9) 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 WRITE (IUNIT8,800) 800 FORMAT (/' ----------------------------------------' + / /' COMPILING AN ORDERED LIST OF BOUNDARY NODES...'/) 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, 3)+4,I) N3=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 IF (.NOT.CHECKN(N3)) THEN NCOND=NCOND+1 CHECKN(N3)=.TRUE. ENDIF ELSE C (TRIANGULAR ELEMENT HAS AN EXTERIOR FAULT ELEMENT C ADJACENT TO IT) EDGETS(J,I)=.FALSE. N2=NODES(MOD(J, 3)+4,I) IF (NODEF(2,KFAULT).EQ.N2) THEN EDGEFS(2,KFAULT)=.TRUE. DO 806 K=4,6 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,3 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 (NUMNOD.GT.NREALN) THEN DO 824 I=NREALN+1,NUMNOD IF (.NOT.CHECKN(I)) THEN IO=N1000+I-NREALN WRITE(IUNIT8,822) IO 822 FORMAT(' BAD GRID TOPOLOGY; FAKE NODES ARE NOT', + ' PERMITTED IN THE INTERIOR.'/' CHECK NODE ',I6) STOP ENDIF 824 CONTINUE ENDIF 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 FIRST, BE SURE THAT WE ARE NOT STARTING ON A MIDPOINT: IF (NODTYP(I).EQ.2) THEN DO 833 K=1,NUMEL DO 832 L=1,3 IF (EDGETS(L,K)) THEN N2=NODES(MOD(L, 3)+4,K) IF (N2.EQ.I) THEN J=NODES(MOD(L+1,3)+1,K) GO TO 839 ENDIF ENDIF 832 CONTINUE 833 CONTINUE DO 835 K=1,NFL IF (EDGEFS(1,K)) THEN IF (NODEF(2,K).EQ.I) THEN J=NODEF(3,K) GO TO 839 ENDIF ELSE IF (EDGEFS(2,K)) THEN IF (NODEF(5,K).EQ.I) THEN J=NODEF(6,K) GO TO 839 ENDIF ENDIF 835 CONTINUE 839 NDONE=2 NODCON(2)=J NLEFT=NCOND-1 ENDIF 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,3)+4,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. N3=NODES(MOD(J+1,3)+1,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N3 CHECKN(N3)=.FALSE. NLEFT=NLEFT-2 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) N3=NODEF(3,I) GO TO 856 ENDIF ELSE IF (EDGEFS(2,I)) THEN IF (NODEF(4,I).EQ.NODE) THEN N2=NODEF(5,I) N3=NODEF(6,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. NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N3 CHECKN(N3)=.FALSE. NLEFT=NLEFT-2 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 ((NODTYP(I).EQ.1).AND. + ((XNODE(I).EQ.X).AND.(YNODE(I).EQ.Y)))GO TO 867 ENDIF 865 CONTINUE WRITE(IUNIT8,866) NDONE, NODE, X, Y 866 FORMAT(/' AFTER CONNECTING ',I6,' NODES AROUND THE', + ' PERIMETER, '/ + ' PROCESS WAS STOPPED BY BAD GRID TOPOLOGY;'/ + ' COULD NOT FIND ANY WAY TO CONTINUE FROM NODE ',I6/ + ' AT (X=',1P,E10.3,',Y=',E10.3,')'/ + ' EITHER THROUGH SHARED BOUNDARY ELEMENTS, OR'/ + ' THROUGH OTHER BOUNDARY NODES SHARING THE SAME ', + 'POSITION.') 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 X Y') DO 890 I=1,NCOND N=NODCON(I) IF (N.GT.NREALN) N=N1000+N-NREALN WRITE(IUNIT8,882) I, N, XNODE(NODCON(I)), + YNODE(NODCON(I)) 882 FORMAT(' ',2I6,1P,2E11.3) 890 CONTINUE N=NODCON(1) IF (N.GT.NREALN) N=N1000+N-NREALN WRITE (IUNIT8,892) N 892 FORMAT(' (NOTE: NODE ',I6,' COMPLETES THE LOOP, BUT WILL', + ' NOT BE LISTED TWICE.)') ENDIF 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 SUBROUTINE SQUARE 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 INTEGER, PARAMETER :: NDREF = 300 INTEGER I,J,IUNITT,LASTDR,LAYER,LAYER1,LAYER2,N1,N2,NSTEP LOGICAL CALLED,INPUT,OUTPUT REAL ALPHAT, + DENSE,DENSE1,DENSE2,DREF, + ELEVAT, + FRAC,FRAC1,FRAC2, + GEOTH1,GEOTH2,GEOTH3,GEOTH4,GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN,H, + OLDSZZ,OLDPR,ONEKM, + PR,PREF, + RESID,RHOAST,RHOBAR,RHOH2O,RHOTOP, + SIGZZ,SIGZZB, + T,TAUZZ,TEMLIM,TEMPC,TEMPM, + Z,ZBASE,ZM,ZSTOP,ZTOP 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) 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 SUBROUTINE SQUEEZ C C C SUBROUTINE XYTOLL (INPUT,X,Y, + MAPTYP, + RADIUS,CPNLAT,Y0NLAT,X0ELON, + IZONE,XOF500, + OUTPUT,PLAT,PLON) C C Convert points expressed as (x,y) to (longitude, latitude) C in degrees East and North (considered positive), C using either: C MAPTYP = 1: Conic projection, with C tangent latitude CPNLAT and origin at (Y0NLAT,X0ELON) C Necessary constants are global data (RADIUS,CPNLAT, C Y0NLAT,X0ELON). C MAPTYP = 2: Universal Transverse Mercator, with C UTM zone IZONE. C The only global constants required are RADIUS C of the Earth and the length of 500 km, both converted C to the same units as the .feg file: XOF500. C LOGICAL INPUT,IZONE,MAPTYP,OUTPUT REAL ANGLE,ARADS,ATAN2F,CENLON,COSDEG,CPNLAT, + DOT1,DOT2,DOT3,EQUAT, + PLAT,PLON,R,RADIUS,RTAN, + S,SINDEG,TANDEG, + UVEC,UVEC1,UVEC2,UVEC3, + X,XOF500,XRADS,XT,X0ELON,Y,Y0NLAT,YPOLE,YRP DIMENSION UVEC(3),UVEC1(3),UVEC2(3),UVEC3(3) C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C IF (MAPTYP.EQ.1) THEN C Conic projection: RTAN=RADIUS*TANDEG(90.-CPNLAT) YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) 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. ELSE IF (MAPTYP.EQ.2) THEN C Universal Transverse Mercator CENLON=-183.0+6.0*IZONE C C Find 3 orthogonal axes: C Uvec1 points to the equator at the CENLON longitude. UVEC1(1)=COSDEG(CENLON) UVEC1(2)=SINDEG(CENLON) UVEC1(3)=0.0 C Uvec2 points to the South pole. UVEC2(1)=0.0 UVEC2(2)=0.0 UVEC2(3)=-1.0 C Uvec3 points to equator 90 d. E of CENLON. UVEC3(1)=COSDEG(CENLON+90.0) UVEC3(2)=SINDEG(CENLON+90.0) UVEC3(3)=0.0 C C Convert (x,y) to (lon,lat): ARADS=Y/RADIUS XT=X-XOF500 XRADS=1.570796327D0-2.0D0*ATAN(EXP(-(XT*1.0D0)/RADIUS)) C DOT1= COS(XRADS)*COS(ARADS) DOT2= -COS(XRADS)*SIN(ARADS) DOT3= SIN(XRADS) C UVEC(1)=DOT1*UVEC1(1)+DOT2*UVEC2(1)+DOT3*UVEC3(1) UVEC(2)=DOT1*UVEC1(2)+DOT2*UVEC2(2)+DOT3*UVEC3(2) UVEC(3)=DOT1*UVEC1(3)+DOT2*UVEC2(3)+DOT3*UVEC3(3) EQUAT=SQRT(UVEC(1)**2+UVEC(2)**2) PLAT=57.29578*ATAN2(UVEC(3),EQUAT) PLON=57.29578*ATAN2(UVEC(2),UVEC(1)) IF (PLON.LT.-180.) PLON=PLON+360.0 IF (PLON.GT.+180.) PLON=PLON-360.0 END IF RETURN END SUBROUTINE XYTOLL