!-------------------------------------------------------------------- PROGRAM OrbData ! OrbData.f90 has been updated to Fortran 90 (free-format) syntax. ! By Peter Bird, Dept. of Earth, Planetary, and Space Sciences, ! University of California, Los Angeles, CA 90095-1567. ! (For version data, see 1st WRITE below.) ! Reads a finite-element grid produced by OrbWin ! (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). ! If any nodal elevation is non-zero, it will always be left ! unchanged, to preserve effects of manual-editing. ! If any elevation(s) is/are zero, new value(s) will be ! interpolated from a gridded dataset. ! If any nodal heat-flow is non-zero, it will always be left ! unchanged, to preserve effects of manual-editing. ! If any heat-flow(s) is/are zero, new value(s) will be ! interpolated from a gridded dataset. ! There is an option to base heat-flow on age of seafloor, ! where known, from a gridded-age dataset. ! Heat-flow (from any source) is subject to a minimum value ! designed to guarantee 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. !--------------------------------------------------------------------- ! TYPE STATEMENTS ! (Note: The IMPLICIT typing of i~n = INTEGER, and a~h, o~z = REAL*8 ! is observed in this program. No types are stated for such names.) IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) ! PARAMETER (ARRAY-SIZE) STATEMENTS ! Set the following PARAMETERs at least as large as your problem: ! maxNod = maximim number of nodes (includes both "real" and "fake") INTEGER, PARAMETER :: maxnod = 50000 ! maxBn = maximum number of boundary nodes (both "real" and "fake"). INTEGER, PARAMETER :: maxbn = 10000 ! maxEl = maximum number of continuum elements (triangles). INTEGER, PARAMETER :: maxel = 100000 ! maxFEl = maximum number of fault elements (line segments); INTEGER, PARAMETER :: maxfel = 12000 ! maxAtP = maximum number of nodes which may overlap at a fault- ! intersection point. INTEGER, PARAMETER :: maxatp = 10 ! PARAMETER giving the (exact) number of PB2002 plates: INTEGER, PARAMETER :: numplt = 52 !--------------------------------------------------------------------- CHARACTER*1 :: c1 CHARACTER*2 :: names CHARACTER*80 :: title1, title3 INTEGER :: aArray_columns, aArray_rows, & & eArray_columns, eArray_rows, & & qArray_columns, qArray_rows INTEGER :: continuum_LRi, fault_LRi LOGICAL :: brief, dothis, everyp, llimit, neede, needq LOGICAL*1 :: checke, checkf, checkn, edgets, edgefs, useage REAL*8 :: lrange, nlat, nlat0, nlat1 ! THE FOLLOWING TO AGREE WITH BLOCK DATA BD1: DOUBLE PRECISION :: points, weight ! THE FOLLOWING TO AGREE WITH BLOCK DATA BD2: DOUBLE PRECISION :: fphi, fpoint, fgauss !--------------------------------------------------------------------- ! DIMENSION STATMENTS ! DIMENSIONs using dynamic memory allocation: REAL*8, DIMENSION(:, :), allocatable :: eArray, qArray, aArray ! DIMENSIONS USING PARAMETER MAXNOD: DIMENSION checkn(maxnod), dqdtda(maxnod), & & elev (maxnod), & & tlnode(maxnod), & & xnode (maxnod), ynode (maxnod), zmnode(maxnod) ! DIMENSIONS USING PARAMETER MAXBN: DIMENSION nodcon (maxbn) ! DIMENSIONS USING PARAMETER MAXEL: DIMENSION area (maxel), checke (maxel), & & continuum_LRi(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) ! DIMENSIONS USING PARAMETER MAXFEL: DIMENSION checkf (maxfel), edgefs (2, maxfel), & & fault_LRi (maxfel), & & fdip (2, maxfel), & & flen (maxfel), & & fpflt(2, 2, 2, 7, maxfel), & & farg (2, maxfel), nodef (4, maxfel), & & offset (maxfel) ! DIMENSIONS USING PARAMETER MAXATP: DIMENSION list (maxatp) ! DIMENSIONS USING PARAMETER NUMPLT: DIMENSION names (numplt) ! FIXED DIMENSIONS: DIMENSION acreep(2), alphat(2), bcreep(2), ccreep(2), conduc(2), & & dcreep(2), radio(2), & & rhobar(2), taumax(2), temlim(2) ! FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: DIMENSION points(3, 7), weight(7) ! FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: DIMENSION fphi(4, 7), fpoint(7), fgauss(7) !--------------------------------------------------------------------- ! COMMON STATEMENTS ! FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: COMMON / s1s2s3 / points COMMON / wgtvec / weight ! FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: COMMON / sfault / fpoint COMMON / fphis / fphi COMMON / fglist / fgauss !------------------------------------------------------------------- ! DATA STATEMENTS ! "QLIM0" is the lower limit on heat-flow for points ! with an elevation of zero. ! (A minimum heat-flow limit is necessary to find ! a converged solution.) DATA qlim0 / 0.045D0 / ! "DQL_DE" is the derivitive d(QLIMIT)/d(elevation), ! which adjusts the minimum heat-flow for elevation. ! (A minimum heat-flow limit is necessary to find ! a converged solution.) DATA dql_de / 1.43D-06 / ! "QLIM1 is the upper limit on heat-flow for all points. DATA qlim1 / 0.300D0 / ! "CLIMIT" IS A LOWER LIMIT ON CRUSTAL THICKNESS: DATA climit / 6570.0D0 / ! Changed from 5000. to 6570. on 2005.05.24, to agree with CRUST2.grd ! "HCMAX" IS AN UPPER LIMIT ON CRUSTAL THICKNESS: DATA hcmax / 75000.0D0 / ! Already agreed with CRUST2.grd; no need to change! ! "HLMAX" IS AN UPPER LIMIT ON TOTAL LITHOSPHERE THICKNESS: DATA hlmax / 300000.0D0 / ! (NOTE THAT ALL OF THE ABOVE ARE IN SI UNITS (W/m**2, m, m, m). ! "IUNITT" = FORTRAN DEVICE NUMBER FOR TEXT MESSAGES: DATA iunitt / 6 / ! "IUNITL" = FORTRAN DEVICE NUMBER FOR -Assign- LOG FILE: DATA iunitl / 11 / ! 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' / ! "DIPMAX" IS THE MAXIMUM DIP (FROM HORIZONTAL, IN DEGREES) FOR A ! FAULT ELEMENT TO BE TREATED AS A DIP-SLIP FAULT, WITH TWO DEGREES ! OF FREEDOM PER NODE-PAIR. AT STEEPER DIPS, THE DEGREE OF FREEDOM ! CORRESPONDING TO OPENING OR CONVERGENCE OF THE OPPOSITE SIDES IS ! ELIMINATED BY A CONSTRAINT EQUATION, AND THE FAULT IS TREATED AS ! A VERTICAL STRIKE-SLIP FAULT. THIS ARBITRARY LIMIT IS NECESSARY ! BECAUSE THE EQUATIONS FOR DIP-SLIP FAULTS BECOME SINGULAR AS THE ! DIP APPROACHES 90 DEGREES. IN PRACTICE, IT IS BEST TO SPECIFY DIPS ! AS EITHER (1) VERTICAL, OR (2) CLEARLY LESS THAN "DIPMAX", WITHIN ! EACH FAULT ELEMENT. IF THE DIP VARIES WITHIN AN ELEMENT IN SUCH A ! WAY THAT IT PASSES THROUGH THIS LIMIT WITHIN THE ELEMENT, THEN ! THE REPRESENTATION OF THAT FAULT ELEMENT IN THE EQUATIONS MAY ! BE INACCURATE. DATA dipmax / 75.0D0 / !--------------------------------------------------------------------- ! STATEMENT FUNCTIONS: EastOf(x) = MOD((x + 720.0D0), 360.0D0) !--------------------------------------------------------------------- ! BEGINNING OF EXECUTABLE CODE ! *** KLUDGE ALERT *** ! CONVERSION OF PARAMETERS (CONSTANTS) TO VARIABLES SHOULD LOGICALLY ! HAVE NO EFFECT, BUT IN FACT HELPS TO SUPPRESS SOME SPURIOUS ! MESSAGES FROM THE IBM VS-FORTRAN COMPILER. mxnode = maxnod mxbn = maxbn mxel = maxel mxfel = maxfel mxstar = maxatp wedge = ABS(90.0D0 - ABS(dipmax)) * 0.017453293D0 !--------------------------------------------------------------------- ! INTRODUCTION / HEADER SECTION: WRITE (iunitt, 1) 1 FORMAT ( & & ' -------------------------------------------------------------' & & /' Program OrbData (version of 14 February 2019)' & &//' by Peter Bird, Dept. of Earth, Planetary, & Space Sciences,' & & /' University of California, Los Angeles, CA 90095-1567.' & &//' Reads a finite element grid produced by OrbWin' & & /' (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 <= nlat0) THEN WRITE (*, "(' ERROR: Maximum latitude must ', & & 'be >= minimum.')") GO TO 37 END IF END IF 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) ! Read parameters on UNIT 1 CALL ReadPm ( 1, 6, names, numplt, offmax, & ! inputs & acreep, alphat, bcreep, Biot , & ! outputs & Byerly, ccreep, cfric , conduc, & ! outputs & dcreep, ecreep, everyp, ffric , gmean , & ! outputs & gradie, iconve, ipvref, & ! outputs & maxitr, okdelv, oktoqt, onekm, radio , & ! outputs & rdeth, refstr, rhoast, rhobar, rhoh2o, & ! outputs & tadiab, & ! outputs & taumax, temlim, title3, trhmax, tsurf, & ! outputs & zbasth) ! outputs ! Lithosphere/asthenosphere temperature, in Kelvin, is determined ! from the asthenosphere adiabat in the parameter file, ! evaluated (arbitrarily) at 100 km depth below sea level: tasthk = tadiab + gradie * 100.E3 ! Read finite element grid on UNIT 2: CALL GetNet ( 2, iunitt, & ! inputs & mxel, mxfel, mxnode, & ! inputs & brief, & ! outputs & continuum_LRi, & ! outputs & dqdtda, eledat, elev, & ! outputs & fault_LRi, fdip, & ! outputs & nfaken, nfl, nodef, nodes, nrealn, & ! outputs & numel, numnod, n1000, offmax, offset, & ! outputs & title1, tlnode, xnode, ynode, zmnode, & ! outputs & checke, checkf, checkn) ! work-space arrays WRITE (iunitt, 40) 40 FORMAT (/' Successfully read FE grid, now verifying topology...') CALL Square (brief, fdip, iunitt, & ! inputs & mxbn, mxel, mxfel, mxnode, & ! inputs & mxstar, nfl, nodef, nodes, & ! inputs & numel, numnod, rdeth, wedge, & ! inputs & xnode, ynode, & ! modify & area, detj, & ! outputs & dxs, dys, dxsp, dysp, edgefs, & ! outputs & edgets, flen, fpflt, fpsfer, & ! outputs & farg, ncond, nodcon, sita, & ! outputs & checkn, list) ! work-space arrays WRITE (iunitt, 50) 50 FORMAT (/' Topology of finite element grid verified.'/) ! Read in elevation array on UNIT 3, if needed: neede = .FALSE. DO 60 i = 1, numnod IF (llimit) THEN elon = ynode(i) * 57.2957795D0 nlat = 90.0D0 - xnode(i) * 57.2957795D0 IF ((EastOf(elon - welon) <= lrange).AND. & & (nlat >= nlat0).AND. & & (nlat <= nlat1)) THEN IF (elev(i) == 0.0D0) neede = .TRUE. END IF ELSE IF (elev(i) == 0.0D0) 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 eArray_rows = NINT(1.0D0 + (ey2 - ey1) / edy) eArray_columns = NINT(1.0D0 + (ex2 - ex1) / edx) ALLOCATE ( eArray(eArray_rows, eArray_columns) ) READ (3, * ) ((eArray(irow, jcol), jcol = 1, eArray_columns), irow = 1, eArray_rows) ELSE WRITE (*, 69) 69 FORMAT (/' All nodes have non-zero elevation.' & & /' No input elevation grid is needed.') END IF ! Read in heat-flow array on UNIT 4, if needed: needq = .FALSE. DO 70 i = 1, numnod IF (llimit) THEN elon = ynode(i) * 57.2957795D0 nlat = 90.0D0 - xnode(i) * 57.2957795D0 IF ((EastOf(elon - welon) <= lrange).AND. & & (nlat >= nlat0).AND. & & (nlat <= nlat1)) THEN IF (dqdtda(i) == 0.0D0) needq = .TRUE. END IF ELSE IF (dqdtda(i) == 0.0D0) 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 qArray_rows = NINT(1.0D0 + (qy2 - qy1) / qdy) qArray_columns = NINT(1.0D0 + (qx2 - qx1) / qdx) ALLOCATE ( qArray(qArray_rows, qArray_columns) ) READ (4, * ) ((qArray(irow, jcol), jcol = 1, qArray_columns), irow = 1, qArray_rows) ! Check whether 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 ! Read dataset of gridded ages, on UNIT 7: ! (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 aArray_rows = NINT(1.0D0 + (ay2 - ay1) / ady) aArray_columns = NINT(1.0D0 + (ax2 - ax1) / adx) ALLOCATE ( aArray(aArray_rows, aArray_columns) ) READ (7, * ) ((aArray(irow, jcol), jcol = 1, aArray_columns), irow = 1, aArray_rows) ELSE ! define dummy array size to keep compiler happy: ALLOCATE ( aArray(1, 1) ) END IF ELSE WRITE (*, 79) 79 FORMAT (/' All nodes have non-zero heat-flow.' & & /' No input heat-flow data sets are needed.') ! Impose heat-flow limits at nodes, to prevent unreasonably ! thick and stiff lithosphere anywhere: DO 80 i = 1, numnod qlimit = qlim0 + dql_de * elev(i) IF (dqdtda(i) /= 0.0D0) dqdtda(i) = MAX(dqdtda(i), qlimit) dqdtda(i) = MIN(dqdtda(i), qlim1) 80 CONTINUE ! define dummy array size to keep compiler happy: ALLOCATE ( aArray(1, 1) ) END IF ! INITIATE LOG FILE FOR -Assign- OUTPUT: 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-:') ! PROCESS ALL NODES EQUALLY: 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.2957795D0 plat = 90.0D0 - xnode(i) * 57.2957795D0 elevat = elev(i) heatfl = dqdtda(i) IF (llimit) THEN dothis = ((EastOf(plon - welon) <= lrange).AND. & & (plat >= nlat0).AND. & & (plat <= nlat1)) ELSE dothis = .TRUE. END IF IF (dothis) THEN CALL Assign (aArray, aArray_columns, aArray_rows, & ! inputs & ax1, adx, ax2, ady, ay2, & ! inputs & alphat, climit, conduc, & ! inputs & eArray, eArray_columns, eArray_rows, & ! inputs & ex1, edx, ex2, edy, ey2, & ! inputs & gmean, hcmax, hlmax, & ! inputs & iunitl, iunitt, & ! inputs & onekm, & ! inputs & plon, plat, & ! inputs & qlim0, dql_de, qlim1, & ! inputs & qArray, qArray_columns, qArray_rows, & ! inputs & qx1, qdx, qx2, qdy, qy2, & ! inputs & radio, rhoast, rhobar, rhoh2o, & ! inputs & tasthk, temlim, tsurf, useage, & ! inputs & elevat, heatfl, thickc, thickm) ! modify elev(i) = elevat dqdtda(i) = heatfl zmnode(i) = MAX(thickc, climit) tlnode(i) = MAX(thickm, 0.0D0) WRITE (iunitl, 678) i, elev(i), dqdtda(i), zmnode(i), tlnode(i) 678 FORMAT (' ', I8, 4ES10.2) WRITE (iunitt, 679) i 679 FORMAT ('+', I8, ' nodes completed...') END IF 680 CONTINUE WRITE (iunitt, 700) 700 FORMAT (/ /' CHECK THE LOG FILE CAREFULLY FOR PROBLEMS!') ! Output the modified .FEG file: CALL PutNet ( 12, & ! inputs & brief, & ! inputs & continuum_LRi, & ! inputs & dqdtda, eledat, elev, & ! inputs & fault_LRi, fdip, & ! inputs & mxel, mxfel, mxnode, n1000, & ! inputs & nfaken, nfl, nodef, nodes, & ! inputs & nrealn, numel, numnod, offset, & ! inputs & title1, tlnode, xnode, ynode, zmnode) ! inputs CALL Pause() END PROGRAM OrbData REAL*8 FUNCTION ATan2F (y, x) ! Corrects for problem of two zero arguments REAL*8, INTENT(IN) :: y, x IF ((y /= 0.0D0).OR.(x /= 0.0D0)) THEN ATan2F = ATAN2(y, x) ELSE ATan2F = 0.0D0 END IF END FUNCTION ATan2F SUBROUTINE ReadN (iunitp, iunitt, n, & ! inputs & vector) ! outputs ! A utility routine designed to permit -Faults- input files ! to also be used by -Plates-, which expects more numbers ! in some records. ! This routine attempts to read 'n' floating-point values ! (using * format) from the next record on device 'iUnitP'. ! If anything goes wrong, the missing values are set to zero. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) CHARACTER*1 c CHARACTER*80 line LOGICAL anyin, dotted, expon, signed DIMENSION vector(n) READ (iunitp, 1) line 1 FORMAT (A80) number = 0 anyin = .FALSE. expon = .FALSE. signed = .FALSE. dotted = .FALSE. DO 10 i = 1, 80 c = line(i:i) IF ((c == ' ').OR.(c == ',').OR.(c == '/')) THEN signed = .FALSE. expon = .FALSE. dotted = .FALSE. IF (anyin) THEN number = number + 1 anyin = .FALSE. END IF ELSE IF ((c == '+').OR.(c == '-')) THEN IF (signed) THEN GO TO 50 ELSE signed = .TRUE. END IF ELSE IF ((c == 'E').OR.(c == 'D').OR. & & (c == 'e').OR.(c == 'd')) THEN IF (expon) THEN GO TO 50 ELSE expon = .TRUE. signed = .FALSE. dotted = .TRUE. END IF ELSE IF (c == '.') THEN IF (dotted) THEN GO TO 50 ELSE dotted = .TRUE. END IF ELSE IF ((c == '0').OR.(c == '1').OR.(c == '2').OR. & & (c == '3').OR.(c == '4').OR.(c == '5').OR. & & (c == '6').OR.(c == '7').OR.(c == '8').OR. & & (c == '9')) THEN signed = .TRUE. anyin = .TRUE. ELSE GO TO 50 END IF 10 CONTINUE IF (anyin) number = number + 1 50 IF (number == 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 < n) THEN DO 99 i = (number + 1), n vector(i) = 0.0D0 99 CONTINUE END IF END IF END SUBROUTINE ReadN SUBROUTINE ReadNFromString (passed_line, iUniTt, n, & ! inputs & vector) ! outputs ! Just like subprogram ReadN (above), except that it reads from ! the single-line string "passed_line" argument, instead of ! reading a record from an attached Fortran device number. ! A utility routine designed to permit -Faults- input files ! to also be used by -Plates-, which expects more numbers ! in some records. ! This routine attempts to read 'n' floating-point values ! (using * format) from the next record on device 'iUnitP'. ! If anything goes wrong, the missing values are set to zero. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) CHARACTER*1 c CHARACTER*80 line, passed_line LOGICAL anyin, dotted, expon, signed DIMENSION vector(n) line = passed_line ! copying, to prevent any accidental changes number = 0 anyin = .FALSE. expon = .FALSE. signed = .FALSE. dotted = .FALSE. DO 10 i = 1, 80 c = line(i:i) IF ((c == ' ').OR.(c == ',').OR.(c == '/')) THEN signed = .FALSE. expon = .FALSE. dotted = .FALSE. IF (anyin) THEN number = number + 1 anyin = .FALSE. END IF ELSE IF ((c == '+').OR.(c == '-')) THEN IF (signed) THEN GO TO 50 ELSE signed = .TRUE. END IF ELSE IF ((c == 'E').OR.(c == 'D').OR. & & (c == 'e').OR.(c == 'd')) THEN IF (expon) THEN GO TO 50 ELSE expon = .TRUE. signed = .FALSE. dotted = .TRUE. END IF ELSE IF (c == '.') THEN IF (dotted) THEN GO TO 50 ELSE dotted = .TRUE. END IF ELSE IF ((c == '0').OR.(c == '1').OR.(c == '2').OR. & & (c == '3').OR.(c == '4').OR.(c == '5').OR. & & (c == '6').OR.(c == '7').OR.(c == '8').OR. & & (c == '9')) THEN signed = .TRUE. anyin = .TRUE. ELSE GO TO 50 END IF 10 CONTINUE IF (anyin) number = number + 1 50 IF (number == 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) line = passed_line ! re-load (although perhaps not necessary) READ (line, * ) (vector(i), i = 1, number) IF (number < n) THEN DO 99 i = (number + 1), n vector(i) = 0.0D0 99 CONTINUE END IF END IF END SUBROUTINE ReadNFromString SUBROUTINE Squeez (alphat, elevat, & ! inputs & geoth1, geoth2, geoth3, geoth4, & ! inputs & geoth5, geoth6, geoth7, geoth8, & ! inputs & gmean, & ! inputs & iunitt, onekm, rhoast, rhobar, rhoh2o, & ! inputs & temlim, zm, zstop, & ! inputs & tauzz, sigzzb) ! outputs ! Calculates "tauZZ", the vertical integral through the plate ! of the vertical standardized stress anomaly, which is ! relative to a column of mantle at asthenosphere temperature ! with a 5-km crust and a 2.7-km ocean on top, like a mid-ocean ! rise. The integral is from either the land surface or the ! sea surface, down to a depth of "zStop" below the top of ! the crust. ! If "zStop" exceeds Moho depth "zM", then properties of the mantle ! will be used in the lower part of the integral. ! Also returns "sigZZB", the standardized vertical stress anomaly ! at depth "zStop" below the solid rock surface. ! Note: This version is different from the version found in the Laramy ! program package. First, it acts on only a single point. ! Second, it infers sub-plate normal-stress anomalies from ! the given topography, instead of from model structure. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) INTEGER, PARAMETER :: ndref = 700 LOGICAL called ! Internal arrays: DIMENSION dref(ndref), pref(0:ndref) ! Argument arrays: DIMENSION alphat(2), rhobar(2), temlim(2) SAVE called, dref, pref DATA called / .FALSE. / ! 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) ! Create reference temperature & density profiles to depth of ndref km: IF (.NOT.called) THEN rhotop = rhobar(1) * (1. - alphat(1) * geoth1) dref(1) = rhoh2o dref(2) = rhoh2o dref(3) = 0.7D0 * rhoh2o + 0.3D0 * rhotop dref(4) = rhotop dref(5) = rhotop dref(6) = rhotop dref(7) = rhotop dref(8) = 0.7D0 * rhotop + 0.3D0 * rhoast DO 50 j = 9, ndref dref(j) = rhoast 50 CONTINUE pref(0) = 0.0D0 DO 100 i = 1, ndref pref(i) = pref(i - 1) + dref(i) * gmean * onekm 100 CONTINUE END IF ! ROUTINE PROCESSING (ON EVERY CALL): IF (elevat > 0.0D0) THEN ! LAND ztop = -elevat zbase = zstop - elevat dense1 = rhobar(1) * (1. - geoth1 * alphat(1)) h = 0.0D0 layer1 = 1 ELSE ! OCEAN ztop = 0.0D0 zbase = zstop + (-elevat) dense1 = rhoh2o h = elevat layer1 = 0 END IF lastdr = zbase / onekm IF (zbase > onekm * lastdr) lastdr = lastdr + 1 IF (lastdr > ndref) THEN WRITE(iunitt, 110) lastdr 110 FORMAT(' IN SUBPROGRAM SQUEEZ, PARAMETER NDREF '/ & & ' MUST BE INCREASED TO AT LEAST ',I10) CALL Pause() STOP END IF nstep = (zbase - ztop) / onekm oldszz = 0.0D0 oldpr = 0.0D0 sigzz = 0.0D0 tauzz = 0.0D0 z = ztop DO 200 i = 1, nstep z = z + onekm h = h + onekm IF (h > 0.0D0) THEN IF (h <= zm) THEN t = TempC(h) dense2 = rhobar(1) * (1.0D0 - t * alphat(1)) layer2 = 1 ELSE t = TempM(h - zm) dense2 = rhobar(2) * (1.0D0 - t * alphat(2)) layer2 = 2 END IF ELSE dense2 = rhoh2o layer2 = 0 END IF IF ((layer1 == 0).AND.(layer2 == 1)) THEN frac2 = h / onekm frac1 = 1.0D0 - frac2 ELSE IF ((layer1 == 1).AND.(layer2 == 2)) THEN frac2 = (h - zm) / onekm frac1 = 1.0D0 - frac2 ELSE frac1 = 0.5D0 frac2 = 0.5D0 END IF dense = frac1 * dense1 + frac2 * dense2 IF (z > 0.) THEN n1 = z / onekm n2 = n1 + 1 frac = z / onekm - n1 pr = pref(n1) + frac * (pref(n2) - pref(n1)) ELSE pr = 0.0D0 END IF sigzz = sigzz - dense * gmean * onekm + (pr - oldpr) tauzz = tauzz + 0.5D0 * (sigzz + oldszz) * onekm dense1 = dense2 oldszz = sigzz oldpr = pr layer1 = layer2 200 CONTINUE resid = zbase - z h = zstop z = zbase IF (zstop <= zm) THEN t = TempC(h) dense2 = rhobar(1) * (1.0D0 - t * alphat(1)) ELSE t = TempM(h - zm) dense2 = rhobar(2) * (1.0D0 - t * alphat(2)) END IF dense = 0.5D0 * (dense1 + dense2) IF (z > 0.0D0) THEN n1 = z / onekm n2 = n1 + 1 frac = z / onekm - n1 pr = pref(n1) + frac * (pref(n2) - pref(n1)) ELSE pr = 0.0D0 END IF sigzzb = sigzz - dense * gmean * resid + (pr - oldpr) tauzz = tauzz + 0.5D0 * (sigzzb + oldszz) * resid called = .TRUE. END SUBROUTINE Squeez SUBROUTINE Assign (aArray, aArray_columns, aArray_rows, & ! inputs & ax1, adx, ax2, ady, ay2, & ! inputs & alphat, climit, conduc, & ! inputs & eArray, eArray_columns, eArray_rows, & ! inputs & ex1, edx, ex2, edy, ey2, & ! inputs & gmean, hcmax, hlmax, & ! inputs & iunitl, iunitt, & ! inputs & onekm, & ! inputs & plon, plat, & ! inputs & qlim0, dql_de, qlim1, & ! inputs & qArray, qArray_columns, qArray_rows, & ! inputs & qx1, qdx, qx2, qdy, qy2, & ! inputs & radio, rhoast, rhobar, rhoh2o, & ! inputs & tasthk, temlim, tsurf, useage, & ! inputs & elevat, heatfl, thickc, thickm) ! modify ! Determines elevation (IF 0.0), heat-flow (IF 0.0), and isostatic ! values of crustal and mantle-lithosphere thickness, at one point: ! (pLon = East longitude, pLat = North latitude), both in degrees. ! (Unlike subprogram -Assign- in AlGriddr, this version never ! does regional averaging. Neither does it interpolate by ! Kriging; a simple interpolation off a rectangular grid is ! used instead.) IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) INTEGER, INTENT(IN) :: iunitl, iunitt INTEGER, INTENT(IN) :: aArray_columns, aArray_rows, & & eArray_columns, eArray_rows, & & qArray_columns, qArray_rows LOGICAL, INTENT(IN) :: useage REAL*8, 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*8, DIMENSION(2), INTENT(IN) :: alphat, conduc, radio, rhobar, & & temlim REAL*8, DIMENSION(aArray_rows, aArray_columns), INTENT(IN) :: aArray REAL*8, DIMENSION(eArray_rows, eArray_columns), INTENT(IN) :: eArray REAL*8, DIMENSION(qArray_rows, qArray_columns), INTENT(IN) :: qArray REAL*8, INTENT(INOUT) :: elevat, heatfl, thickc, thickm ! Number of iterations permitted: INTEGER, PARAMETER :: maxitr = 50 ! Damping factor (<1.0D0) for stability of iterations: REAL*8, PARAMETER :: damp = 0.200D0 LOGICAL :: badp, badt, diverg, outsid, neede, needq, & & warnc1, warnc2, warnm1, warnl2, wayout ! Local storage for history of iteration: DIMENSION hc(maxitr), hm(maxitr), tm(maxitr), ta(maxitr) !--------------------------------------------------------------------- outsid = .FALSE. wayout = .FALSE. warnc1 = .FALSE. warnc2 = .FALSE. warnm1 = .FALSE. warnl2 = .FALSE. ! Determine the elevations, if needed: neede = (elevat == 0.0D0) IF (neede) THEN ir1 = ((ey2 - plat) / edy) + 1.00001D0 ir1 = MAX(ir1, 1) ir1 = MIN(ir1, eArray_rows - 1) ir2 = ir1 + 1 fr = ((ey2 - edy * (ir1 - 1)) - plat) / edy IF ((plon < ex1).OR.(plon > ex2)) THEN IF (((plon + 360.0D0) >= ex1).AND. & & ((plon + 360.0D0) <= ex2)) THEN uselon = plon + 360.0D0 ELSE IF (((plon - 360.0D0) >= ex1).AND. & & ((plon - 360.0D0) <= ex2)) THEN uselon = plon - 360.0D0 ELSE uselon = plon END IF ELSE uselon = plon END IF ic1 = ((uselon - ex1) / edx) + 1.00001D0 ic1 = MAX(ic1, 1) ic1 = MIN(ic1, eArray_columns - 1) ic2 = ic1 + 1 fc = (uselon - (ex1 + edx * (ic1 - 1))) / edx outsid = outsid.OR.(fr < -0.01D0).OR.(fr > 1.01D0).OR. & & (fc < -0.01D0).OR.(fc > 1.01D0) wayout = wayout.OR.(fr < -1.01D0).OR.(fr > 2.01D0).OR. & & (fc < -1.01D0).OR.(fc > 2.01D0) fr = MIN(1.0D0, MAX(0.0D0, fr)) fc = MIN(1.0D0, MAX(0.0D0, 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) END IF ! Determine heat-flow, if needed: needq = (heatfl == 0.0D0) IF (needq) THEN ir1 = ((qy2 - plat) / qdy) + 1.00001D0 ir1 = MAX(ir1, 1) ir1 = MIN(ir1, qArray_rows - 1) ir2 = ir1 + 1 fr = ((qy2 - qdy * (ir1 - 1)) - plat) / qdy IF ((plon < qx1).OR.(plon > qx2)) THEN IF (((plon + 360.0D0) >= qx1).AND. & & ((plon + 360.0D0) <= qx2)) THEN uselon = plon + 360.0D0 ELSE IF (((plon - 360.0D0) >= qx1).AND. & & ((plon - 360.0D0) <= qx2)) THEN uselon = plon - 360.0D0 ELSE uselon = plon END IF ELSE uselon = plon END IF ic1 = ((uselon - qx1) / qdx) + 1.00001D0 ic1 = MAX(ic1, 1) ic1 = MIN(ic1, qArray_columns - 1) ic2 = ic1 + 1 fc = (uselon - (qx1 + qdx * (ic1 - 1))) / qdx outsid = outsid.OR.(fr < -0.01D0).OR.(fr > 1.01D0).OR. & & (fc < -0.01D0).OR.(fc > 1.01D0) wayout = wayout.OR.(fr < -1.01D0).OR.(fr > 2.01D0).OR. & & (fc < -1.01D0).OR.(fc > 2.01D0) fr = MIN(1.0D0, MAX(0.0D0, fr)) fc = MIN(1.0D0, MAX(0.0D0, 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) ! Check for seafloow-age overriding heat-flow grid value: IF (useage) THEN ir1 = ((ay2 - plat) / ady) + 1.00001D0 ir1 = MAX(ir1, 1) ir1 = MIN(ir1, aArray_rows - 1) ir2 = ir1 + 1 fr = ((ay2 - ady * (ir1 - 1)) - plat) / ady IF ((plon < ax1).OR.(plon > ax2)) THEN IF (((plon + 360.0D0) >= ax1).AND. & & ((plon + 360.0D0) <= ax2)) THEN uselon = plon + 360.0D0 ELSE IF (((plon - 360.0D0) >= ax1).AND. & & ((plon - 360.0D0) <= ax2)) THEN uselon = plon - 360.0D0 ELSE uselon = plon END IF ELSE uselon = plon END IF ic1 = ((uselon - ax1) / adx) + 1.00001 ic1 = MAX(ic1, 1) ic1 = MIN(ic1, aArray_columns - 1) ic2 = ic1 + 1 fc = (uselon - (ax1 + adx * (ic1 - 1))) / adx outsid = outsid.OR.(fr < -0.01D0).OR.(fr > 1.01D0).OR. & & (fc < -0.01D0).OR.(fc > 1.01D0) wayout = wayout.OR.(fr < -1.01D0).OR.(fr > 2.01D0).OR. & & (fc < -1.01D0).OR.(fc > 2.01D0) fr = MIN(1.0D0, MAX(0.0D0, fr)) fc = MIN(1.0D0, MAX(0.0D0, 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 < 200.0D0) THEN ! SEAFLOOR AGE IS VALID (NOT "UNKNOWN" OR "CONTINENTAL IF (agema <= 0.0D0) THEN heatfl = qlim1 ELSE ! Carol A. Stein & Seth Stein [1992] ! A model for the global variation in oceanic ! depth and heat flow with lithospheric age, ! Nature, v. 359, 10 September, p. 123-129. ! According to their preferred GDH1 model: IF (agema <= 55.0D0) THEN heatfl = 0.510D0 / SQRT(agema) ELSE heatfl = 0.048D0 + 0.096D0 * EXP(-0.0278D0 * agema) END IF END IF END IF END IF END IF ! Give warning if any necessary data could not be read: 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() !CCCC STOP END IF ! APPLY LIMITS ON Q TO GUARANTEE A STEADY-STATE SOLUTION: qlimit = qlim0 + dql_de * elevat heatfl = MAX(heatfl, qlimit) heatfl = MIN(heatfl, qlim1) ! This is the re-entry point when heat-flow must be reduced, ! and the entire iteration process restarted, for this point: 10 CONTINUE ! INFER CRUSTAL THICKNESS FROM ISOSTASY, AND MANTLE LITHOSPHERE ! THICKNESS FROM ASTHENOSPHERE ADIABAT TEMPERATURE (@ 100 km). thickc = 35.0D3 thickm = 65.0D3 geoth1 = tsurf geoth2 = heatfl / conduc(1) geoth3 = -radio(1) / (2.0D0 * conduc(1)) geoth4 = 0.0D0 geoth7 = -radio(2) / (2.0D0 * conduc(2)) geoth8 = 0.0D0 DO 150 iter = 1, maxitr qred = heatfl - thickc * radio(1) tmoho = geoth1 + geoth2 * thickc + geoth3 * thickc**2 ! REMEMBER VALUE FOR COVERGENCE REPORT: tm(iter) = tmoho geoth5 = tmoho geoth6 = qred / conduc(2) test = geoth5 + geoth6 * thickm + geoth7 * thickm**2 ! REMEMBER VALUE FOR COVERGENCE REPORT: ta(iter) = test terr0r = test - tasthk dtl = -terr0r / geoth6 ! NOTE: NEXT LINE USES DAMPING FACTOR: dtl = dtl * damp dtl = MIN(dtl, + 10000.0D0) dtl = MAX(dtl, -10000.0D0) thickm = thickm + dtl IF ((thickm < 0.0D0).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.0D0) IF (((thickc + thickm) > 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) ! Remember value for convergence report: hm(iter) = thickm CALL Squeez (alphat, elevat, & ! inputs & geoth1, geoth2, geoth3, geoth4, & ! inputs & geoth5, geoth6, geoth7, geoth8, & ! inputs & gmean, iunitt, & ! inputs & onekm, rhoast, rhobar, rhoh2o, & ! inputs & temlim, thickc, thickc + thickm, & ! inputs & tauzz, sigzzb) ! outputs dc = -sigzzb / (gmean * (rhobar(2) - rhobar(1))) ! Note: Next line uses damping factor: thickc = thickc + damp * dc IF ((thickc < 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 > 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) ! Remember value for convergence report: hc(iter) = thickc 150 CONTINUE ! Test for successful convergence: ! Isostasy is BAD if basal pressure anomaly is large: badp = (ABS(sigzzb) > 4.0D7) ! in Pascals (SI) ! Geotherm might not connect to asthenosphere adiabat: badt = (ABS(terr0r) > 20.0D0) ! 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 > 100.0D0).AND.(dc25 > dc24)).OR. & & ((dm25 > 100.0D0).AND.(dm25 > dm24)) 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 ! It may be necessary to repeat the whole computation with ! reduced heat-flow if results were hotter than the asthenosphere ! temperature: IF (badt.AND.(terr0r > 0.0D0).AND.(heatfl > (1.01D0 * qlimit))) THEN qnew = 0.99D0 * 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 END SUBROUTINE Assign BLOCK DATA BD1 ! Define "weight" (Gaussian integration weights) of the ! seven integration points in each element, defined by internal ! coordinates "points(3,7)", where points(1-3,m)=s1~s3 of ! integration point number m. ! Because all of these arrays are functions of internal ! coordinates, they are not affected by scaling or deformation of ! the finite elements. DOUBLE PRECISION points, weight COMMON / s1s2s3 / points COMMON / wgtvec / weight DIMENSION points(3, 7), weight(7) ! "points" contains the internal coordinates (s1,s2,s3) of the 7 ! Gaussian integration points (for area-integrals) of the ! triangular elements. "points" is also the set of nodal functions ! 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 / ! "weight" is the Gaussian weight (for area-integrals) of the 7 ! integration points in each triangular element. DATA weight / 0.2250000000000000D0, & & 0.1323941500000000D0, 0.1323941500000000D0, 0.1323941500000000D0, & & 0.1259391833333333D0, 0.1259391833333333D0, 0.1259391833333333D0 / END BLOCK DATA BD1 BLOCK DATA BD2 ! Define "fPhi" (nodal functions) and "fGauss" (Gaussian integration ! weights) at the 7 integration points in each fault element, ! defined by internal coordinate "fPpointm=1,...,7)" ! which contains the relative position ! (fractional length) at the integration points. ! Because all of these arrays are functions of internal ! coordinates, they are not affected by scaling or deformation of ! the elements. DOUBLE PRECISION fphi, fpoint, fgauss COMMON / sfault / fpoint COMMON / fphis / fphi COMMON / fglist / fgauss DIMENSION fphi(4, 7), fpoint(7), fgauss(7) ! "fPoint" contains the seven integration point locations for the fault ! elements. Each value gives a position as a fraction of total length ! measured from node1 to node2 (of array "nodeF"). DATA fpoint / & & 0.0254461D0, & & 0.1292344D0, & & 0.2970774D0, & & 0.5000000D0, & & 0.7029226D0, & & 0.8707656D0, & & 0.9745539D0 / ! "fGauss" contains the sever corresponding weight factors for use in ! line-integrals. DATA fgauss / & & 0.0647425D0, & & 0.1398527D0, & & 0.1909150D0, & & 0.2089796D0, & & 0.1909150D0, & & 0.1398527D0, & & 0.0647425D0 / ! "fPhi" contains the values of the 4 nodal functions (one per node) ! at each of these 7 integration points in the fault element. ! A special convention is that the nodal function of node 3 ! is the negative of that for node 2, while the nodal function ! for node 4 is the negative of that for node 1. This simplifies ! many expressions in which we would otherwise have to have ! a separate factor of +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 / END BLOCK DATA BD2 SUBROUTINE Square (brief, fdip, iunit8, & ! inputs & mxbn, mxel, mxfel, mxnode, & ! inputs & mxstar, nfl, nodef, nodes, & ! inputs & numel, numnod, rdeth, wedge, & ! inputs & xnode, ynode, & ! modify & area, detj, & ! outputs & dxs, dys, dxsp, dysp, edgefs, & ! outputs & edgets, flen, fpflt, fpsfer, & ! outputs & farg, ncond, nodcon, sita, & ! outputs & checkn, list) ! work-space arrays ! Check, correct, and complete the geometryof the grid: IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) LOGICAL agreed, allok, brief, found, switch, vert1, vert2 LOGICAL*1 :: checkn, edgefs, edgets 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) ' / ! (1) Check that all nodes are connected to at least one ! continuum (triangular) element or fault element: 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'/ & & ' BELONG TO ANY TRIANGULAR CONTINUUM ELEMENT'/ & & ' 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 END IF ! (2) Average together the coordinates of all nodes at one "point" DO 410 i = 1, numnod checkn(i) = .FALSE. ! (means "not yet involved in averaging') 410 CONTINUE DO 490 i = 1, nfl DO 480 j1 = 1, 2 nj1 = nodef(j1, i) ! (Fault-ends are the only places that can have problems.) IF (.NOT.checkn(nj1)) THEN list(1) = nj1 checkn(nj1) = .TRUE. ! Begin list of neighbors with paired node: j2 = 5 - j1 nj2 = nodef(j2, i) list(2) = nj2 checkn(nj2) = .TRUE. ninsum = 2 ! Find shortest fault connected to either on: dx = xnode(nj1) - xnode(nj2) dy = ynode(nj1) - ynode(nj2) IF (dy > 3.14159D0) dy = dy - 6.28318D0 IF (dy < -3.14159D0) dy = dy + 6.28318D0 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 == nl1).OR.(nj2 == nl1).OR. & & (nj1 == nl2).OR.(nj2 == nl2).OR. & & (nj1 == nl3).OR.(nj2 == nl3).OR. & & (nj1 == nl4).OR.(nj2 == nl4)) THEN dx = xnode(nl1) - xnode(nl2) dy = ynode(nl1) - ynode(nl2) IF (dy > 3.14159D0) dy = dy - 6.28318D0 IF (dy < -3.14159D0) dy = dy + 6.28318D0 dy = dy * SIN(xnode(nl1)) test = SQRT(dx**2 + dy**2) short = MIN(short, test) END IF 470 CONTINUE ! COLLECT ALL CORNER NODES WITHIN 10% OF THIS toler = short / 10.0D0 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 > 3.14159D0) dy = dy - 6.28318D0 IF (dy < -3.14159D0) dy = dy + 6.28318D0 dy = dy * SIN(xnode(nj1)) r2 = dx**2 + dy**2 IF (r2 < t2) THEN ninsum = ninsum + 1 IF (ninsum > mxstar) THEN WRITE(iunit8, 421) 421 FORMAT(/' INCREASE VALUE' & & ,' OF PARAMETER MAXATP.') CALL Pause() STOP END IF list(ninsum) = k checkn(k) = .TRUE. END IF END IF 471 CONTINUE ! (Quick exit if all nodes in same place) agreed = .TRUE. DO 472 k = 2, ninsum agreed = agreed.AND. & & (xnode(list(k)) == xnode(list(1))).AND. & & (ynode(list(k)) == ynode(list(1))) 472 CONTINUE IF (agreed) GO TO 480 xsum = 0.0D0 ysum = 0.0D0 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.0D0 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 > 0.0D0) 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!') END IF END IF END IF 480 CONTINUE 490 CONTINUE ! (3) Compute derivitives of nodal ! functions at integration points; ! then check for negative areas: CALL Deriv (iunit8, mxel, mxnode, nodes, numel, & ! inputs & rdeth, xnode, ynode, & ! inputs & area, detj, & ! outputs & dxs, dys, dxsp, dysp, fpsfer, sita) ! outputs allok = .TRUE. DO 620 i = 1, numel DO 610 m = 1, 7 test = area(i) * detj(m, i) IF (test <= 0.0D0) 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. END IF 610 CONTINUE 620 CONTINUE IF (.NOT.allok) THEN CALL Pause() STOP END IF ! (4) Compute lengths of fault elements: 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 ! (5) Make a list of nodes that are on the boundary and require ! boundary conditions (nodCon); these are in counterclockwise ! order. Also make lists of element sides which contain these ! nodes: edgeTs AND edgeFs: 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 (i, j, mxel, mxfel, nfl, nodef, nodes, numel, & ! inputs & kfault, kele) ! outputs IF (kele > 0) THEN ! (ordinary interior side) edgets(j, i) = .FALSE. ELSE IF (kfault == 0) THEN ! (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. END IF IF (.NOT.checkn(n2)) THEN ncond = ncond + 1 checkn(n2) = .TRUE. END IF ELSE ! (triangular element has an exterior fault element ! adjacent to it) edgets(j, i) = .FALSE. n1 = nodes(MOD(j, 3) + 1, i) IF (nodef(2, kfault) == 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. END IF 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. END IF 808 CONTINUE END IF END IF 809 CONTINUE 810 CONTINUE IF (ncond > mxbn) THEN WRITE(iunit8, 820) ncond 820 FORMAT(/' INCREASE PARAMETER MAXBN TO',I6,' AND RECOMPILE.') CALL Pause() STOP END IF ! Stop work if no boundary nodes found (global grid): IF (ncond == 0) GO TO 899 ! 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 ! Beginning of indefinate loop which traces around the perimeter. ! Each time, it progresses by one of 3 steps: ! -2 nodes at a time along a triangle side, or ! -2 nodes at a time along a fault element side, or ! -by finding another (corner) node which shares the same location. ! Beginning of main indefinate loop: 840 node = nodcon(ndone) x = xnode(node) y = ynode(node) ! 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 == node) GO TO 846 END IF 842 CONTINUE 844 CONTINUE GO TO 850 846 n2 = nodes(MOD(j + 1, 3) + 1, i) ndone = ndone + 1 IF (ndone <= ncond) nodcon(ndone) = n2 checkn(n2) = .FALSE. nleft = nleft - 1 IF (nleft > 0) THEN GO TO 840 ELSE GO TO 870 END IF ! 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) == node) THEN n2 = nodef(2, i) GO TO 856 END IF ELSE IF (edgefs(2, i)) THEN IF (nodef(3, i) == node) THEN n2 = nodef(4, i) GO TO 856 END IF END IF 854 CONTINUE GO TO 860 856 ndone = ndone + 1 IF (ndone <= ncond) nodcon(ndone) = n2 checkn(n2) = .FALSE. nleft = nleft - 1 IF (nleft > 0) THEN GO TO 840 ELSE GO TO 870 END IF ! Else, look for another exterior corner node at same location: 860 DO 865 i = 1, numnod IF ((i /= node).AND.checkn(i)) THEN IF ( (ABS(xnode(i) - x) < 1.D-6) .AND. & & (ABS(ynode(i) - y) < 1.D-6) ) GO TO 867 END IF 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 <= ncond) nodcon(ndone) = i checkn(i) = .FALSE. nleft = nleft - 1 IF (nleft > 0) GO TO 840 ! 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(' ',2I8) 890 CONTINUE n = nodcon(1) WRITE (iunit8, 892) n 892 FORMAT(' (Note: Node ',I6,' completes the loop, but will', & & ' not be listed twice.)') END IF 899 CONTINUE ! (6) Survey fault elements and issue warnign if any element is of ! mixed type (part strike-slip, and part shallow-digging): DO 920 i = 1, nfl deld1 = fdip(1, i) - 1.570796D0 deld2 = fdip(2, i) - 1.570796D0 vert1 = ABS(deld1) <= wedge vert2 = ABS(deld2) <= wedge nvpart = 0 IF (vert1) THEN nvpart = nvpart + 1 tag1 = vertic ELSE tag1 = obliqu END IF IF (vert2) THEN nvpart = nvpart + 1 tag2 = vertic ELSE tag2 = obliqu END IF switch = ((nvpart > 0).AND.(nvpart < 2)) IF (switch) THEN dip1 = fdip(1, i) * 57.2957795D0 IF (dip1 > 90.0D0) dip1 = dip1 - 180.0D0 dip2 = fdip(2, i) * 57.2957795D0 IF (dip2 > 90.0D0) dip2 = dip2 - 180.0D0 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) <= wedge) nvpart = nvpart + 1 910 CONTINUE IF ((nvpart > 0).AND.(nvpart < 7)) THEN IF (nvpart >= 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.)') END IF END IF END IF 920 CONTINUE ! (7) Calculate fault argument (in radians, measured counterclockwise ! from +Theta) for each integration point in each fault element: ! Here use cross-product "d" of position vector "a" and "b" ! of fault element node "n1" and "n2" respectively, to calculate ! vector "d" in plate of "a" and "b" and perpendicular to "a". ! Then use "d" and "a" to decide the position vector "gold" of ! interpolation point "m" and "g" OF "m+1". Finally use the ! dot-product of "gold - pol" and "g - gold" to decide the angle ! between these two vectors, i.e., fArg. Here pol is position ! 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 (phi, theta, & ! inputs & fangle) ! output DO 900 j = 1, 2 farg(j, i) = fangle(j) fangle(j) = fangle(j) * 57.29577951D0 900 CONTINUE 1000 CONTINUE ! (8) Survey strike-slip (vertical) faults to check for conflicts in ! argument that would lock the fault: ! Loop on all fault elements (i): DO 2000 i = 1, nfl ! Loop on 2 terminal node pairs, 1-4, 2-3 (j = 1 OR 2): DO 1900 j = 1, 2 ! Dip must be within "wedge" of vertical for constraint: IF (ABS(fdip(j, i) - 1.570796D0) <= wedge) THEN nazi = j n1 = j IF(j == 1) THEN n4 = 4 ELSE n4 = 3 END IF node1 = nodef(n1, i) node4 = nodef(n4, i) ! No constraint applied where a fault ends: IF (node1 /= node4) THEN ! Endpoint pairs must be checked for duplication: ! Look for other strike-slip faults sharing this ! pair of nodes, at either end: found = .FALSE. DO 1600 l = 1, nfl IF (l /= i) THEN IF (ABS(fdip(1, l) - 1.5708D0) <= wedge) THEN IF (((node1 == nodef(1, l)).AND. & & (node4 == nodef(4, l))).OR. & & ((node1 == nodef(4, l)).AND. & & (node4 == nodef(1, l)))) THEN found = .TRUE. number = l nazl = 1 GO TO 1601 END IF END IF IF (ABS(fdip(2, l) - 1.5708D0) <= wedge) THEN IF (((node1 == nodef(2, l)).AND. & & (node4 == nodef(3, l))).OR. & & ((node1 == nodef(3, l)).AND. & & (node4 == nodef(2, l)))) THEN found = .TRUE. number = l nazl = 2 GO TO 1601 END IF END IF END IF 1600 CONTINUE ! Don't worry if this pair already checked! 1601 IF (found.AND.(number > i)) THEN ! Average arguments together (avoid cycle shifts): IF(nazi == nazl) THEN azl = farg(nazl, number) + 3.141592654D0 ELSE azl = farg(nazl, number) END IF azi = farg(nazi, i) cosz = 0.5D0 * (COS(azi) + COS(azl)) sinz = 0.5D0 * (SIN(azi) + SIN(azl)) azimut = ATan2F(sinz, cosz) farg(nazi, i) = azimut IF(nazl == nazi) THEN farg(nazl, number) = azimut - 3.141592654D0 ELSE farg(nazl, number) = azimut END IF END IF ! ^end block which looks for constraints on real nodes END IF ! ^end block which checks for distinct node numbers END IF ! ^end block which checks for dip of over 75 degrees 1900 CONTINUE ! ^end loop on 2 node pairs in fault element 2000 CONTINUE ! (9) Calculate nodal function at interpolation points ! on spherical great-circle fault: CALL FNodal (iunitt, mxfel, & ! inputs & mxnode, nfl, nodef, xnode, ynode, & ! inputs & fpflt) ! output IF (.NOT. brief) WRITE (iunit8, 9999) 9999 FORMAT (' --------------------------------------------------', & & '-----------------------------') END SUBROUTINE Square SUBROUTINE Deriv (iunitt, mxel, mxnode, & ! inputs & nodes, numel, & ! inputs & rdeth, xnode, ynode, & ! inputs & area, detj, & ! outputs & dxs, dys, dxsp, dysp, fpsfer, sita) ! outputs ! Sets up 6 vector nodal functions (fpsfer) of each spherical ! triangle finite element, at each of its 7 integration points. ! Calculates dXS and dYS, the theta-derivitive and phi-derivitive ! of each of these 6 vector nodal functions. ! Also computes area, the areas of the plane triangles. ! Also computes detJ, the local ratio of areas on the sphere ! to areas on the plane triangle. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) 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.29578D0 / 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.5D0 * rdeth * rdeth * areap IF (area(i) <= 0.0D0) THEN DO 9 k = 1, 3 xlon(k) = phi(k) * oezopi ylat(k) = 90.0D0 - 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 END IF 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 ! This part is to test if Kong's method and Bird's method ! give the same results for these derivitives: 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 DO 800 m = 1, 7 snccse = 0.0D0 sncsne = 0.0D0 cosm = 0.0D0 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) ! Bird's method: 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 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) <= 0.0D0.OR.sita(m, i) >= 3.141592654D0) THEN sitami = sita(m, i) * 57.29577951D0 WRITE(iunitt, 220) m, i, sitami 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5,' OF ELEMENT', & & I5,' IS OUT RANGE', & & E14.4) CALL Pause() STOP END IF 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 END SUBROUTINE Deriv SUBROUTINE Next (i, j, mxel, mxfel, nfl, nodef, nodes, numel, & ! inputs & kfault, kele) ! outputs ! Determine whether there are more elements adjacent to side j of ! triangular continuum element #i. ! j = 1 means the side opposite node # nodes(1,i). ! j = 2 means the side opposite node # nodes(2,i). ! j = 3 means the side opposite node # nodes(3,i). ! If a fault element is adjacent, its number is kFault; otherwise, ! kFault is set to zero. ! If another triangular continuum element is adjacent (even across ! fault element kFault!) then its number is kEle; otherwise, kEle = 0. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) LOGICAL foundf DIMENSION nodef(4, mxfel), nodes(3, mxel) ! Three node numbers along the side of interest, counterclockwise: n1 = nodes(MOD(j, 3) + 1, i) n2 = nodes(MOD(j + 1, 3) + 1, i) ! Check for adjacent fault element first: foundf = .FALSE. kfault = 0 IF (nfl > 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 == n2).AND.(m2 == n1)).OR. & & ((m3 == n2).AND.(m4 == n1))) THEN foundf = .TRUE. kfault = k GO TO 11 END IF 10 CONTINUE 11 IF (.NOT.foundf) kfault = 0 ! If there was a fault, replace 2 node numbers that we search for: IF (foundf) THEN IF (m2 == n1) THEN n1 = m3 n2 = m4 ELSE n1 = m1 n2 = m2 END IF END IF END IF ! Search for adjacent triangular continuum element: kele = 0 klow = i khigh = i ! --- Begin irregular loop, to search out nearest elements first --- 100 klow = klow - 1 IF (klow >= 1) THEN DO 110 l = 1, 3 m1 = nodes(MOD(l, 3) + 1, klow) m2 = nodes(MOD(l + 1, 3) + 1, klow) IF ((m2 == n1).AND.(m1 == n2)) THEN kele = klow RETURN END IF 110 CONTINUE END IF khigh = khigh + 1 IF (khigh <= numel) THEN DO 120 l = 1, 3 m1 = nodes(MOD(l, 3) + 1, khigh) m2 = nodes(MOD(l + 1, 3) + 1, khigh) IF ((m2 == n1).AND.(m1 == n2)) THEN kele = khigh RETURN END IF 120 CONTINUE END IF IF ((klow > 1).OR.(khigh < numel)) GO TO 100 END SUBROUTINE Next REAL*8 FUNCTION FltLen (phi1, phi2, rdeth, theta1, theta2) ! Calculates length of great-circle segment between ! point (theta1,phi1) and point (theta2,phi2), ! in physical length units (radians * planet_radius). REAL*8, INTENT(IN) :: phi1, phi2, rdeth, theta1, theta2 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 END FUNCTION FltLen SUBROUTINE ReadPm (iunit7, iunit8, names, numplt, offmax, & ! inputs & acreep, alphat, bcreep, Biot , & ! outputs & Byerly, ccreep, cfric , conduc, & ! outputs & dcreep, ecreep, everyp, ffric , gmean , & ! outputs & gradie, iconve, ipvref, & ! outputs & maxitr, okdelv, oktoqt, onekm, radio , & ! outputs & rdeth, refstr, rhoast, rhobar, rhoh2o, & ! outputs & tadiab, & ! outputs & taumax, temlim, title3, trhmax, tsurf, & ! outputs & zbasth) ! outputs ! Reads strategic and tactical input parameters from device iUnit7, ! and echoes them on device iUnit8 with annotations. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) 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) 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 user''s 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.0D0, Biot)) WRITE (iunit8, 40) Biot 40 FORMAT (' ',F10.4,' effective-pressure (Biot) coefficient,', & & ' 0.0 to 1.0') READ (iunit7, * ) Byerly Byerly = MAX(0.0D0, MIN(0.99D0, Byerly)) IF (offmax > 0.0D0) 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)') END IF CALL ReadN (iunit7, iunit8, 2, & ! inputs & acreep) ! outputs 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 (iunit7, iunit8, 2, & ! inputs & bcreep) ! outputs 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 (iunit7, iunit8, 2, & ! inputs & ccreep) ! outputs 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 (iunit7, iunit8, 2, & ! inputs & dcreep) ! outputs 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') 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) == pltref) ipvref = i 956 CONTINUE IF (ipvref == 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 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') CALL ReadN (iunit7, iunit8, 2, & ! inputs & vector) ! output taumax(1) = vector(1) taumax(2) = vector(2) ! Provide for old parameter files with only one tauMax: IF (taumax(2) <= 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) < 0.).OR.(taumax(2) < 0.)) THEN WRITE (iunit8, 107) 107 FORMAT (/' *** ERR','OR in parameter input file: ***' & & /' tauMax may not be negative.') CALL Pause() STOP END IF READ (iunit7, * ) rhoh2o WRITE (iunit8, 110) rhoh2o 110 FORMAT (' ',1P,E10.3,' density of groundwater, lakes, & oceans') CALL ReadN (iunit7, iunit8, 2, & ! inputs & rhobar) ! outputs 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 planet') CALL ReadN (iunit7, iunit8, 2, & ! inputs & alphat) ! outputs WRITE (iunit8, 160) alphat(1), alphat(2) 160 FORMAT (' ',1P,E10.2,'/',E10.2,' volumetric thermal ', & & 'expansion of rocks', & & ' (1/vol)*(d.vol/d.T). (crust/mantle)') CALL ReadN (iunit7, iunit8, 2, & ! inputs & conduc) ! outputs WRITE (iunit8, 170) conduc(1), conduc(2) 170 FORMAT (' ',1P,E10.2,'/',E10.2,' thermal conductivity, energy/', & & 'length/sec/deg. (crust/mantle)') CALL ReadN (iunit7, iunit8, 2, & ! inputs & radio) ! outputs 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 (Kelvin) scale') CALL ReadN (iunit7, iunit8, 2, & ! inputs & temlim) ! outputs WRITE (iunit8, 190) temlim(1), temlim(2) 190 FORMAT (' ', F10.0,'/',F10.0,' convecting temperature (TMax), on' & & ,' absolute (Kelvin) 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 errors allowed', & & ' due to finite stiffness'/11X, & & '(Such errors 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 errors, ', & & '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 (' --------------------------------------------------', & & '-----------------------------') END SUBROUTINE ReadPm SUBROUTINE FAngls (phi, theta, & ! inputs & fangle) ! output ! Calculate the arguments (angles counterclockwise from +Theta) ! at both ends of an arc of a great-circle. IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) DOUBLE PRECISION fpoint COMMON / sfault / fpoint DIMENSION fangle(2), fpoint(7), phi(2), theta(2) dg180 = 3.141592654D0 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.99D0 xx = s * a1 + (1.0D0 - s) * b1 yy = s * a2 + (1.0D0 - s) * b2 zz = s * a3 + (1.0D0 - 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.01D0 xx = s * a1 + (1.0D0 - s) * b1 yy = s * a2 + (1.0D0 - s) * b2 zz = s * a3 + (1.0D0 - 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 < 0.0) phai = 2.0D0 * 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) END SUBROUTINE FAngls SUBROUTINE FNodal (iunitt, mxfel, mxnode, nfl, nodef, & ! inputs & xnode, ynode, & ! inputs & fpflt) ! outputs ! Calculates vector nodal functions at interpolation point on a ! great-circle fault element: IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) 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 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 (phi, theta, & ! inputs & fpp) ! output 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 <= 0.0D0.OR.sita >= 3.141592654D0) THEN sita = sita * 57.29577951D0 WRITE(iunitt, 220) m, i, sita 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5, & & ' OF FAULT ELEMENT', & & I5,' IS OUT RANGE', & & E14.4) END IF 800 CONTINUE 900 CONTINUE END SUBROUTINE FNodal SUBROUTINE SNodal (phi, theta, & ! inputs & fpp) ! output ! Calculates vector nodal function at interpolation point on a ! great-circle side of a finite element: IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) DOUBLE PRECISION fphi, pp DIMENSION fphi(4, 7), fpp(2, 2, 2, 7), phi(2), theta(2) COMMON / fphis / fphi 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 END SUBROUTINE SNodal SUBROUTINE GetNet (iunit7, iunit8, & ! inputs & mxel, mxfel, mxnode, & ! inputs & brief, & ! outputs & continuum_LRi, & ! outputs & dqdtda, eledat, elev, & ! outputs & fault_LRi, fdip, & ! outputs & nfaken, nfl, nodef, nodes, nrealn, & ! outputs & numel, numnod, n1000, offmax, offset, & ! outputs & title1, tlnode, xnode, ynode, zmnode, & ! outputs & checke, checkf, checkn) ! work-space arrays ! Read finite element grid from UNIT "iUnit7". ! Echo the important values to a text dataset on UNIT "iUnit8". IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) CHARACTER*80 longer_line, shorter_line, title1 INTEGER continuum_LRi, fault_LRi, LRi LOGICAL allok, brief LOGICAL*1 :: checke, checkf, checkn DIMENSION checke(mxel), checkf(mxfel), checkn(mxnode), & & continuum_LRi(mxel), & & dqdtda(mxnode), eledat(3, mxel), elev(mxnode), & & fault_LRi(mxfel), & & fdip(2, mxfel), nodef(4, mxfel), & & nodes(3, mxel), offset(mxfel), tlnode(mxnode), & & xnode(mxnode), ynode(mxnode), zmnode(mxnode) DIMENSION dips(3), vector(7) 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) TRIM(title1) 3 FORMAT(/' Title of finite element grid ='/' ',A80) ! Read number of nodes. ! Input nodal locations (x,y), and (optional) elevations, heat-flow, and perhaps ! other nodal values (e.g., crust & mantle-lithosphere thicknesses). ! (One record per node). ! (Option "brief" suppresses most output.) READ (iunit7, * ) numnod, nrealn, nfaken, n1000, brief brief = .TRUE. IF (numnod /= (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 END IF IF (nrealn > n1000) THEN WRITE (iunit8, 5) nrealn, n1000 5 FORMAT (/' ERR0R: NREALN (',I6,') IS GREATER THAN' & & /' N1000 (',I6,').') CALL Pause() STOP END IF IF (numnod > 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 END IF 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',I8,' 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'/) END IF DO 90 k = 1, numnod checkn(k) = .FALSE. 90 CONTINUE DO 100 k = 1, numnod CALL ReadN (iunit7, iunit8, 7, & ! inputs & vector) ! outputs index = vector(1) + 0.5D0 IF (index > nrealn) THEN IF ((INDEX <= n1000).OR. & & (INDEX > (n1000 + nfaken))) THEN WRITE (iunit8, 91) index 91 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER: ',I6) CALL Pause() STOP END IF END IF plon = vector(2) plat = vector(3) IF (ABS(plat) > 90.0D0) THEN WRITE (iunit8, 92) index 92 FORMAT (/' ERR0R: ABS(LATITUDE) > 90 AT NODE ',I8) CALL Pause() STOP END IF IF (ABS(plat) > 89.99D0) THEN WRITE (iunit8, 93) index 93 FORMAT (/' ERR0R: NODE ',I8,' LIES ON A POLE.' & & /' THIS IS A SINGULAR POINT OF THE' & & ,' SPHERICAL COORDINATE SYSTEM.' & & /' MOVE THIS NODE, AT LEAST SLIGHTLY.') CALL Pause() STOP END IF xi = (90.0D0 - plat) * 0.017453292D0 yi = plon * 0.017453292D0 elevi = vector(4) qi = vector(5) zmi = vector(6) tli = vector(7) IF (index <= nrealn) THEN i = index ELSE i = nrealn + index - n1000 END IF checkn(i) = .TRUE. xnode(i) = xi ynode(i) = yi elev(i) = elevi dqdtda(i) = qi IF (qi < 0.0D0) THEN WRITE (iunit8, 96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL Pause() STOP END IF IF (zmi < 0.0D0) THEN WRITE (iunit8, 97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') CALL Pause() STOP END IF zmnode(i) = zmi IF (tli < 0.0D0) THEN WRITE (iunit8, 98) 98 FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS IS', & & ' NON-PHYSICAL.') CALL Pause() STOP END IF 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) END IF 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 <= nrealn) THEN INDEX = i ELSE INDEX = n1000 + i - nrealn END IF IF (.NOT.checkn(i)) WRITE(iunit8, 103)INDEX 103 FORMAT (' ',36X,I8) 104 CONTINUE CALL Pause() STOP END IF ! Read triangular elements: READ (iunit7, * ) numel IF (numel > mxel) THEN WRITE (iunit8, 108) numel 108 FORMAT(/' INCREASE PARAMETER MAXEL TO BE AT LEAST EQUAL' & & /' TO THE NUMBER OF ELEMENTS (',I8,') AND RECOMPILE.') CALL Pause() STOP END IF DO 109 k = 1, numel checke(k) = .FALSE. 109 CONTINUE IF (.NOT.brief) WRITE (iunit8, 110) numel 110 FORMAT(/' There are ',I8,' triangular continuum elements.'/ & & ' (Node numbers for each are given at corners, counter', & & 'clockwise'/ / & & ' ELEMENT C1 C2 C3') DO 200 k = 1, numel ! (Elements need not be input in order, but must all be present.) ! (ReadN method is used to allow for possible element data.) READ (iUnit7, "(A)") longer_line CALL Extract_LRi (longer_line, & ! input & LRi, shorter_line) ! output CALL ReadNFromString (shorter_line, iUnit8, 7, & ! inputs & vector) ! outputs i = NINT(vector(1)) IF ((i < 1).OR.(i > numel)) THEN WRITE (iunit8, 111) i 111 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I8) CALL Pause() STOP END IF 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) continuum_LRi(i) = LRi checke(i) = .TRUE. IF (.NOT.brief) THEN IF (LRi == 0) THEN WRITE (iUnit8, 120) i, (nodes(j, i), j = 1, 3), (eledat(j, i), j = 1, 3) 120 FORMAT (' ', I8, ':', 3I8, ES8.1, F7.1, ES8.1) ELSE WRITE (iUnit8, "(' ', I8, ':', 3I8, ES8.1, F7.1, ES8.1, ' LR', I8)") i, (nodes(j, i), j = 1, 3), (eledat(j, i), j = 1, 3), LRi END IF END IF DO 130 j = 1, 3 n = nodes(j, i) IF (n > nrealn) n = nrealn + (n - n1000) IF ((n <= 0).OR.(n > numnod)) THEN WRITE (iunit8, 125) nodes(j, i) 125 FORMAT (' NODE NUMBER ',I8,' IS ILLEGAL.') CALL Pause() STOP END IF 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,I8) 204 CONTINUE CALL Pause() STOP END IF ! Read fault elements: READ (iunit7, * ) nfl IF (nfl > mxfel) THEN WRITE (iunit8, 220)nfl 220 FORMAT (/' INCREASE PARAMETER maxFEl TO BE AT LEAST EQUAL' & & /' TO THE NUMBER OF FAULTS (',I8,') AND RECOMPILE.') CALL Pause() STOP END IF offmax = 0.0D0 DO 222 i = 1, nfl checkf(i) = .FALSE. 222 CONTINUE IF (.NOT.brief) WRITE(iunit8, 230) nfl 230 FORMAT(/ /' There are ',I8,' great-circle fault elements.') IF ((.NOT.brief).AND.(nfl > 0)) WRITE(iunit8, 231) 231 FORMAT (/' (The 4 node numbers defining each fault 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 (' ', I8, ':', 4I8, 1X, 2F6.1, 1X, F9.0) DO 300 k = 1, nFl off = 0.0D0 READ (iUnit7, "(A)") longer_line CALL Extract_LRi (longer_line, & ! input & LRi, shorter_line) ! output READ (shorter_line, * ) i, (nodef(j, k), j = 1, 4), (dips(l), l = 1, 2), off IF ((i < 1).OR.(i > nfl)) THEN WRITE (iunit8, 241) i 241 FORMAT (/' ERR0R: ILLEGAL FAULT NUMBER: ',I6) CALL Pause() STOP END IF checkf(i) = .TRUE. IF (.NOT.brief) THEN IF (LRi == 0) THEN WRITE (iunit8, 240) i, (nodef(j, i), j = 1, 4), (dips(l), l = 1, 2), off ELSE WRITE (iunit8, "(' ', I8, ':', 4I8, 1X, 2F6.1, 1X, F9.0, ' LR', I8)") i, (nodef(j, i), j = 1, 4), (dips(l), l = 1, 2), off, LRi END IF END IF DO 250 j = 1, 4 n = nodef(j, i) IF (n > nrealn) n = nrealn + (n - n1000) IF ((n <= 0).OR.(n > numnod)) THEN WRITE (iunit8, 243) nodef(j, i), i 243 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER (',I8, & & ') IN FAULT ',I8) CALL Pause() STOP END IF nodef(j, i) = n 250 CONTINUE DO 260 l = 1, 2 IF (ABS(dips(l)) > 90.0D0) 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 END IF IF (dips(l) < 0.0D0) dips(l) = 180.0D0 + dips(l) fdip(l, i) = dips(l) * 0.017453293D0 260 CONTINUE IF (off < 0.0D0) THEN WRITE (iunit8, 280) off 280 FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2, & & ' FOR FAULT ELEMENT',I8/ & & ' OFFSETS MAY NOT BE NEGATIVE.') CALL Pause() STOP END IF offset(i) = off offmax = MAX(offmax, off) fault_LRi(i) = LRi 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,I8) 304 CONTINUE CALL Pause() STOP ELSE IF (offmax > 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.') END IF END IF IF (.NOT. brief) WRITE (iunit8, 999) 999 FORMAT (' -------------------------------------------------------------------------------') END SUBROUTINE GetNet SUBROUTINE PutNet (iUnitO, & ! inputs & brief, & ! inputs & continuum_LRi, & ! inputs & dqdtda, eledat, elev, & ! inputs & fault_LRi, fdip, & ! inputs & mxel, mxfel, mxnode, n1000, & ! inputs & nfaken, nfl, nodef, nodes, & ! inputs & nrealn, numel, numnod, offset, & ! inputs & title1, tlnode, xnode, ynode, zmnode) ! inputs ! Writes finite element grid to unit "iUnitO". IMPLICIT INTEGER (I-N) IMPLICIT REAL*8 (A-H, O-Z) CHARACTER*80 title1 INTEGER continuum_LRi, fault_LRi, LRi LOGICAL anyemu, brief DIMENSION continuum_LRi(mxel), & & dqdtda(mxnode), eledat(3, mxel), elev(mxnode), & & fault_LRi(mxfel), fdip(2, mxfel), & & nodef(4, mxfel), nodes(3, mxel), np(4), & & offset(mxfel), tlnode(mxnode), & & xnode(mxnode), ynode(mxnode), zmnode(mxnode) DIMENSION dips(2) 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) WRITE (iunito, 2) numnod, nrealn, nfaken, n1000, brief 2 FORMAT(4I8,L8,' (numNod, nRealN, nFakeN, n1000, brief)') DO 100 i = 1, numnod IF (i <= nrealn) THEN iprint = i ELSE iprint = n1000 + (i - nrealn) END IF plat = 90.0D0 - xnode(i) * 57.29577951D0 plon = ynode(i) * 57.29577951D0 WRITE (iunito, 91) i, plon, plat, elev(i), dqdtda(i), zmnode(i), & & tlnode(i) 91 FORMAT (I8, 2F11.5, 4ES10.2) 100 CONTINUE WRITE (iunito, 110) numel 110 FORMAT (I10,' (numEl = number of triangular continuum elements)') anyemu = .FALSE. DO 140 i = 1, numel DO 130 k = 1, 3 IF (eledat(k, i) /= 0.0D0) 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) <= nrealn) THEN np(k) = nodes(k, i) ELSE np(k) = n1000 + (nodes(k, i) - nrealn) END IF 150 CONTINUE LRi = continuum_LRi(i) IF (anyemu) THEN IF (LRi == 0) THEN WRITE (iUnitO, 155) i, (np(k), k = 1, 3), (eledat(k, i), k = 1, 3) 155 FORMAT (I8, 3I8, ES8.1, F7.1, ES8.1) ELSE WRITE (iUnitO, "(I8, 3I8, ES8.1, F7.1, ES8.1, ' LR', I8)") i, (np(k), k = 1, 3), (eledat(k, i), k = 1, 3), LRi END IF ELSE IF (LRi == 0) THEN WRITE (iUnitO, 160) i, (np(k), k = 1, 3) 160 FORMAT (I8, 3I8) ELSE WRITE (iUnitO, "(I8, 3I8, ' LR', I8)") i, (np(k), k = 1, 3), LRi END IF END IF 200 CONTINUE WRITE (iunito, 210) nfl 210 FORMAT (I10,' (nFl = number of curvilinear fault elements)') DO 300 i = 1, nfl DO 220 k = 1, 4 IF (nodef(k, i) <= nrealn) THEN np(k) = nodef(k, i) ELSE np(k) = n1000 + (nodef(k, i) - nrealn) END IF 220 CONTINUE DO 230 k = 1, 2 dips(k) = fdip(k, i) dips(k) = dips(k) * 57.2958D0 IF (dips(k) > 90.0D0) dips(k) = dips(k) - 180.0D0 230 CONTINUE LRi = fault_LRi(i) IF (LRi == 0) THEN WRITE (iUnitO, 250) i, (np(k), k = 1, 4), (dips(k), k = 1, 2), offset(i) 250 FORMAT (I8, 4I8, 2F6.1, ES10.2) ELSE WRITE (iUnitO, "(I8, 4I8, 2F6.1, ES10.2, ' LR', I8)") i, (np(k), k = 1, 4), (dips(k), k = 1, 2), offset(i), LRi END IF 300 CONTINUE END SUBROUTINE PutNet SUBROUTINE Pause() IMPLICIT NONE WRITE(*, "(' Press [Enter]..'\)") READ(*, * ) END SUBROUTINE Pause SUBROUTINE Extract_LRi (longer_line, & ! input & LRi, shorter_line) ! output ! New routine added for Shells_v5.0+ to support multiple !"Lithospheric Rheology" (abbreviated as "LR") integer codes, ! in any line of the input .feg file which define an element !(either a triangular continuum element, or a ! linear fault element). ! CHARACTER*80, INTENT(IN) :: longer_line is the whole ! element-definition line from the .feg file. ! INTEGER, INTENT(OUT) :: LRi is the rheologic code ! (or 0, if no such code was found). ! CHARACTER*80, INTENT(OUT) :: shorter_line has the " LRi" portion removed (if any), ! so it can be interpreted by the same code as in Shells_v4.1-. IMPLICIT NONE CHARACTER*80, INTENT(IN) :: longer_line INTEGER, INTENT(OUT) :: LRi CHARACTER*80, INTENT(OUT) :: shorter_line CHARACTER*80 :: string INTEGER :: longer_length, LR_start_byte longer_length = LEN_TRIM(longer_line) LR_start_byte = INDEX(longer_line, "LR") IF (LR_start_byte > 0) THEN ! the "LR" flag was found IF (longer_length > (LR_start_byte + 1)) THEN ! some byte(s) follow the "LR" string = longer_line((LR_start_byte + 2):longer_length) READ (string, *) LRi shorter_line = longer_line(1:(LR_start_byte - 1)) ELSE ! "LR" is present, but nothing follows it; infer 0. LRi = 0 shorter_line = longer_line(1:(LR_start_byte - 1)) END IF ELSE ! no "LR" flag is present LRi = 0 shorter_line = longer_line END IF END SUBROUTINE Extract_LRi