! PROGRAM -OrbScore2-: COMPARES OUTPUT FROM -SHELLS- ! WITH DATA FROM GEODETIC NETWORKS, ! STRESS DIRECTIONS, FAULT SLIP RATES, ! SEAFLOOR SPREADING RATES, SEISMICITY, ! & UPPER-MANTLE SEISMIC ANISOTROPY, ! AND REPORTS SUMMARY SCALAR SCORES. !=========== PART OF THE "SHELLS" PACKAGE OF PROGRAMS=========== ! GIVEN A FINITE ELEMENT GRID FILE, IN THE FORMAT PRODUCED BY ! -OrbWeave- AND RENUMBERED BY -OrbNumber-, WITH NODAL DATA ! ADDED BY -OrbData-, AND NODE-VELOCITY OUTPUT FROM -SHELLS-, ! COMPUTES A VARIETY OF SCORES OF THE RESULTS. ! NOTE: Does not contain VISCOS or DIAMND, hence independent ! of changes made in May 1998, and equally compatible ! with Old_SHELLS or with improved SHELLS. ! by ! Peter Bird ! Department of Earth and Space Sciences, ! University of California, Los Angeles, California 90095-1567 ! (C) Copyright 1994, 1998, 1999, 2000, 2001, 2002, 2005, 2006 ! by Peter Bird and the Regents of ! the University of California. ! (For version data see FORMAT 1 below) ! THIS PROGRAM WAS DEVELOPED WITH SUPPORT FROM THE UNIVERSITY OF ! CALIFORNIA, THE UNITED STATES GEOLOGIC SURVEY, THE NATIONAL ! SCIENCE FOUNDATION, AND THE NATIONAL AERONAUTICS AND SPACE ! ADMINISTRATION. ! IT IS FREEWARE, AND MAY BE COPIED AND USED WITHOUT CHARGE. ! IT MAY NOT BE MODIFIED IN A WAY WHICH HIDES ITS ORIGIN ! OR REMOVES THIS MESSAGE OR THE COPYRIGHT MESSAGE. ! IT MAY NOT BE RESOLD FOR MORE THAN THE COST OF REPRODUCTION ! AND MAILING. !--------------------------------------------------------------------- ! ADDITIONAL SOFTWARE REQUIRED TO LINK ! Requires DLSARG, from IMSL (International Mathematics Subroutine ! Library). USE Numerical_Libraries, ONLY: DLSARG !CCCC USE MSIMSL, ONLY: DLSARG ! NOTE: Line above is used with COMPAQ Visual Fortran 6; ! it may not be the appropriate method for indicating linkage ! information to other compilers. !--------------------------------------------------------------------- ! DIFFERENCES between OrbScore and OrbScore2: ! (1) OrbScore is in fixed-format, due to FORTRAN 77 heritage. But, ! OrbScore2 has been (nominally, minimally) converted to free-format. ! (2) OrbScore2 uses MODULE Weighting to weight geodesy, stress, ! & anisotropy by geographic area (but not spreading ! or fault slip rates or seismicity). USE Weighting ! N.B. MODULE Weighting has the PRIVATE attribute to shield all ! of its internal variables and arrays from inadvertent ! direct access. This is critical because MODULE Weighting ! creates and stores a *different* F-E grid (a uniform subdivision ! of the icosahedron) for purposes of computing data weights. ! (3) IF (pltGEO) a new output file provides errors in predicted velocities ! of geodetic benchmarks, for plotting (e.g., with FiniteMap). ! (4) OrbScore2 includes scoring for match between upper-mantle seismic ! anisotropy (specifically, fast-direction azimuths from SKS splitting) ! and basal shear traction vectors computed from a torque report, ! in runs with ICONVE = 6. This scoring is area- and delay-weighted. ! N.B. This code, added 2006.11 for the Earth5 project, still needs ! to be generalized to other ICONVE models for basal shear traction! !--------------------------------------------------------------------- ! PARAMETER (ARRAY-SIZE) STATEMENTS ! SET THE FOLLOWING PARAMETERS AT LEAST AS LARGE AS YOUR PROBLEM: ! Numbers of pLates, boundary points (per plate) in PB2002 model of Bird [2003]: INTEGER, PARAMETER :: nPlate = 52 INTEGER, PARAMETER :: nPBnd = 1250 ! MAXNOD = MAXIMUM NUMBER OF NODES (INCLUDES BOTH "REAL"AND & "FAKE") PARAMETER (maxnod = 17000) ! MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES). PARAMETER (maxel = 27000) ! MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS); PARAMETER (maxfel = 3000) ! MAXBN = MAXIMUM NUMBER OF BOUNDARY NODES: PARAMETER (maxbn = 1200) ! MAXATP = MAXIMUM NUMBER OF NODES WHICH MAY OVERLAP AT A FAULT- ! INTERSECTION POINT. PARAMETER (maxatp = 20) ! MAXGEO = MAXIMUM NUMBER OF GEODETIC BENCHMARKS AT WHICH ! HORIZONTAL VELOCITY VECTORS ARE PROVIDED (FOR SCORING). PARAMETER (maxgeo = 6000) ! MAXSTR = MAXIMUM NUMBER OF POINTS AT WHICH MOST COMPRESSIVE ! HORIZONTAL PRINCIPAL STRESS DIRECTION IS PROVIDED ! (FOR SCORING). PARAMETER (maxstr = 10000) ! maxSKS = maximum number of points at which fast-polarization ! ("phi") azimuth and splitting time of SKS waves are ! provided (for scoring). PARAMETER (maxSKS = 10000) ! MAXDAT = MAXIMUM NUMBER OF SITES WITH DATA ON FAULT SLIP RATE: PARAMETER (maxdat = 1000) ! MAXMOR = MAXIMUM NUMBER OF POINTS AT WHICH SEAFLOOR SPREADING ! RATES ARE PROVIDED (FOR SCORING). PARAMETER (maxmor = 800) ! MAXEQS = MAXIMUM NUMBER OF EARTHQUAKES IN CATALOG USED FOR SCORING PARAMETER (maxeqs = 30000) ! MAXADJ = MAXIMUM NUMBER OF (LONGITUDE/LATITUDE/VELOCITY-DATUM/ ! PREDICTED-VELOCITY) POINTS PASSED TO -ADJUST- FOR APPLICATION ! OF A RIGID-BODY ROTATION RATE THAT MINIMIZES THE WEIGHTED ! RMS VELOCITY DIFFERENCE. PARAMETER (maxadj = 188000) ! "piO180" IS PI/180. (CONVERSION FROM DEGREES TO RADIANS): DATA piO180 / 0.017453293 / !--------------------------------------------------------------------- ! TYPE STATEMENTS ! (NOTE: THE IMPLICIT TYPING OF I-N = INTEGER, AND A-H,O-Z = REAL ! IS OBSERVED IN THIS PROGRAM. NO TYPES ARE STATED FOR SUCH NAMES.) CHARACTER*132 gpsfmt CHARACTER*80 enctit, fname, line, title1, title2, title3 CHARACTER*20 geotag, tag CHARACTER*15 frame CHARACTER*10 flag10 CHARACTER*8 c8r1, c8r2, c8rerr, c8v1, c8v2, c8verr CHARACTER*7 c7rcha, c7vcha CHARACTER*5 mortag, strtag CHARACTER*4 outcom CHARACTER*3 c3 CHARACTER*2 names, reg, strreg REAL morphi, morthe, morvel, morsig, mvmma, nlat, nlatp DOUBLE PRECISION s_dp, shrink, v ! FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: DOUBLE PRECISION points, weight ! FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: DOUBLE PRECISION fphi, fpoint, fgauss INTEGER :: subdivision INTEGER, DIMENSION(nPlate):: nBoundaryPoints LOGICAL brief, everyp, e1part, e2part, ezpart, floats, & & havenv, pltGEO, showis, sphere, vertic LOGICAL, DIMENSION(nPlate) :: slab_q ! NOTE: THE FOLLOWING ARRAYS COULD BE COMPRESSED WITH "LOGICAL*1" ! IN VS-FORTRAN: LOGICAL checke, checkf, checkn, edgets, edgefs, explod, fslips REAL nlatde REAL, DIMENSION(nPlate, nPBnd) :: pLat, pLon !--------------------------------------------------------------------- ! DIMENSION STATMENTS ! DIMENSIONS USING PARAMETER MAXNOD: DIMENSION atnode (maxnod), checkn (maxnod), dqdtda(maxnod), & & EDotnc (maxnod), EDotnm (maxnod), & & elev (maxnod), & & tlnode (maxnod), & & uvecn(3, maxnod), v (2, maxnod), & & xnode (maxnod), ynode (maxnod), zmnode(maxnod) ! DIMENSIONS USING PARAMETER MAXEL: DIMENSION area (maxel), & & cartr (3, 7, maxel), checke (maxel), detj (7, maxel), & & dxs(2, 2, 3, 7, maxel), dys(2, 2, 3, 7, maxel), & & dxsp (3, 7, maxel), dysp (3, 7, maxel), edgets(3, maxel), & & EDotec (maxel), EDotem(maxel), & & erate (3, 7, maxel), & & fpsfer(2, 2, 3, 7, maxel), & & nodes (3, maxel), & & sita (7, maxel), & & tlint(7, maxel), & & xip (7, maxel), yip (7, maxel), & & zmoho (7, maxel) ! DIMENSIONS USING PARAMETER MAXFEL: DIMENSION checkf (maxfel), edgefs (2, maxfel), & & fc (2, 2, 7, maxfel), fdip (2, maxfel), & & fimudz(7, maxfel), flen (maxfel), fpeaks(2, maxfel), & & fpflt (2, 2, 2, 7, maxfel), fslips (maxfel), & & farg (2, maxfel), ftstar(2, 7, maxfel), nodef (4, maxfel), & & offset (maxfel), ztranf (2, maxfel) ! DIMENSIONS USING PARAMETER MAXBN: DIMENSION nodcon(maxbn) ! DIMENSIONS OF (2) REFER TO CRUST/MANTLE LITHOSPHERE: DIMENSION acreep(2), alphat(2), bcreep(2), ccreep(2), conduc(2), & & dcreep(2), radio(2), rhobar(2), temlim(2) ! DIMENSIONS USING PARAMETER MAXATP: DIMENSION list (maxatp) ! DIMENSIONS USING PARAMETER nPlate: DIMENSION names(nPlate), omega(3, nPlate) ! DIMENSIONS USING PARAMETER MAXGEO: DIMENSION geotag(maxgeo), geophi(maxgeo), geothe(maxgeo), & & geovel(maxgeo), geosig(maxgeo), geoazi(maxgeo), & & geouth(maxgeo), geouph(maxgeo) ! DIMENSIONS USING PARAMETER MAXSTR: DIMENSION strTag(maxstr), strPhi(maxstr), strThe(maxstr), & & strQua(maxstr), strArg(maxstr), strReg(maxstr) ! DIMENSIONs using PARAMETER maxSTR: CHARACTER*5, DIMENSION(maxSKS) :: SKS_tag ! 5-byte ID code for datum REAL, DIMENSION(maxSKS) :: SKS_theta ! colatitude, in radians REAL, DIMENSION(maxSKS) :: SKS_phi ! colatitude, in radians REAL, DIMENSION(maxSKS) :: SKS_argument ! phi or fast direction, in radians counterclockwise from S REAL, DIMENSION(maxSKS) :: SKS_delay ! SKS splitting time, in s INTEGER, DIMENSION(maxSKS) :: SKS_iPlate ! integer identifying PB2002 plate which contains this datum REAL, DIMENSION(2, maxSKS) :: basal_shear_tractions ! (theta = S, phi = E) components, in Pa ! DIMENSIONS USING PARAMETER MAXDAT: DIMENSION delvz(2, maxdat), fname(maxdat), & & fltlat(maxdat), fltlon(maxdat), & & rlat(2, maxdat) ! DIMENSIONS USING PARAMETER MAXMOR: DIMENSION mortag(maxmor), morphi(maxmor), morthe(maxmor), & & morvel(maxmor), morsig(maxmor) ! DIMENSIONS USING PARAMETER MAXEQS: DIMENSION eqelon(maxeqs), eqnlat(maxeqs), eqmag(maxeqs) ! DIMENSIONS USING PARAMETER MAXADJ: DIMENSION DATA(2, maxadj), predic(2, maxadj), elon(maxadj), & & nlat (maxadj) ! DIMENSIONS OF FIXED SIZE: DIMENSION cartvs(3, 3), cuvec(3), equvec(3), & & ra(3), rb(3), rs(3), tempv(3), & & up(3), ut(3), vf(3), voff(3), vt(3) ! 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) ! ALLOCATABLE arrays added for OrbScore2: REAL, DIMENSION(:), ALLOCATABLE :: weights !--------------------------------------------------------------------- ! COMMON STATEMENTS ! NOTE: UN-NAMED COMMON PASSES INTEGER VARIABLES USED IN THE ! INTEGER-FUNCTION "INDEXK", TO AVOID PASSING THESE SAME ! THROUGH LONG SEQUENCES OF SUBPROGRAMS. COMMON lda, md ! NAMED COMMON BLOCKS BD1 & BD2 HOLD THE POSITIONS, ! WEIGHTS, AND NODAL FUNCTION VALUES AT THE INTEGRATION POINTS ! IN THE ELEMENTS (TRIANGULAR ELEMENTS IN BLOCK DATA BD1, AND ! FAULT ELEMENTS IN BLOCK DATA BD2). ! ENTRIES CORRESPONDING TO BD1: COMMON / s1s2s3 / points COMMON / wgtvec / weight ! ENTRIES CORRESPONDING TO BD2: COMMON / sfault / fpoint COMMON / fphis / fphi COMMON / fglist / fgauss !-------------------------------------------------------------------- ! DATA STATEMENTS ! "FLOATS" is a switch controlling whether subprogram ADJUST ! will be used to minimize the kinetic energy (RMS velocity) ! of velocity-vector lists in two contexts: ! (1) In computing RMS velocities of nodes; ! (2) In computing mismatch between geodetic and model velocities. ! For normal use, set FLOATS /.TRUE./. ! If you are SURE that your geodetic data are expressed in the ! same velocity reference frame as is used to express velocities ! in your finite element model, select FLOATS /.FALSE./ for a ! more discriminating (less forgiving) test of model quality. DATA floats / .TRUE. / ! "oezOpi" IS 180/PI (CONVERSION FROM RADIANS TO DEGREES): DATA oezOpi / 57.2957795 / ! "secPYr" IS THE NUMBER OF SECONDS IN A YEAR: DATA secPYr / 31557600. / ! "unitL" IS THE NUMBER OF MILLIMETERS PER UNIT LENGTH IN THE ! SYSTEM OF UNITS IN USE ON FORTRAN DEVICES iUnitG (FINITE ! ELEMENT GRID), AND DEVICE iUnitP (PARAMETERS), ! AND iUnitV (VELOCITY SOLUTIONS). FOR EXAMPLE, IF THESE ARE IN ! SI UNITS (METERS), THEN unitL = 1000. DATA unitL / 1000. / ! "unitT" IS THE NUMBER OF YEARS PER UNIT OF TIME IN THE ! SYSTEM OF UNITS IN USE ON FORTRAN DEVICES iUnitG (FINITE ! ELEMENT GRID), AND DEVICE iUnitP (PARAMETERS), ! AND iUnitV (VELOCITY SOLUTIONS). FOR EXAMPLE, IF THESE ARE IN ! SI UNITS (SECONDS), THEN unitT = 3.1688087E-08 DATA unitT / 3.1688087E-08 / ! "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. / ! "SUBDIP" IS THE MAXIMUM DIP (FROM HORIZONTAL, IN DEGREES) FOR A ! FAULT IN A WHOLE-EARTH MODEL (SPHERE=.TRUE.) TO BE TREATED AS ! A SUBDUCTION ZONE (IN WHICH CASE, THE FOOTWALL REQUIRES ! BOUNDARY CONDITIONS). ! IN ALL MODELS, FAULTS WITH LESS THAN THIS DIP HAVE THE DOWN-DIP ! INTEGRAL OF TRACTION LIMITED TO "TAUMAX". TAUMAX IS AN ARRAY ! OF TWO VALUES, FOR OCEANIC AND CONTINENTAL SUBDUCTION ZONES, ! RESPECTIVELY. IF SUCH LIMITS ARE NOT WANTED, THEN ! THE TAUMAX VALUES CAN BE SET TO VERY LARGE NUMBERS. DATA subdip / 19.0 / ! THE FOLLOWING ARE THE FORTRAN INPUT AND OUTPUT DEVICE NUMBERS: ! "iUnitT"= DEVICE NUMBER ASSOCIATED WITH TEXT OUTPUT, INCLUDING ! STATUS AND WARNING MESSAGES. ! (NOTE! ON SOME SYSTEMS, SYSTEM ERR0R MESSAGES ARE ALWAYS ! OUTPUT ON DEVICE 6. IF SO, THEN iUnitT SHOULD BE 6.) DATA iUnitT / 6 / ! "iUnitG"= DEVICE NUMBER ASSOCIATED WITH THE .FEG GRID INPUT FILE. DATA iUnitG / 1 / ! "iUnitP"= DEVICE NUMBER ASSOCIATED WITH THE .IN PARAMETER INPUT FILE. DATA iUnitP / 2 / ! "iUnitV"= DEVICE NUMBER ASSOCIATED WITH THE V_____.OUT VELOCITY ! SOLUTION (AT THE NODES) FILE: DATA iUnitV / 3 / ! "iUnitB"= DEVICE NUMBER ASSOCIATED WITH THE GEODETIC FILE OF ! BENCHMARK VELOCITIES. DATA iUnitB / 11 / ! "iUnitS"= DEVICE NUMBER ASSOCIATED WITH THE FILE OF MAXIMUM ! HORIZONTAL PRINCIPAL STRESS DIRECTIONS. DATA iUnitS / 12 / ! "iUnitD"= DEVICE NUMBER ASSOCIATED WITH THE GEOLOGICAL SLIP RATE FILE ! (IN .DIG FORMAT, WITH TWO TITLE LINES PER POINT) DATA iUnitD / 13 / ! "iUnitM"= DEVICE NUMBER ASSOCIATED WITH THE FILE OF MID-OCEAN ! SPREADING RATES. DATA iUnitM / 14 / ! "iUnitE"= DEVICE NUMBER ASSOCIATED WITH THE FILE OF EARTHQUAKES ! TO BE USED FOR SEISMICITY SCORING. DATA iUnitE / 15 / ! "iUnitO"= DEVICE NUMBER TO BE USED FOR OUTPUTTING THE ! .feg FILE WITH SCALAR STRAIN-RATES AT NODES ! BASED ON THE EARTHQUAKE CATALOG FROM iUnitE. DATA iUnitO / 21 / ! "iUnitZ"= DEVICE NUMBER TO BE USED FOR OUTPUTTING THE ! .feg FILE WITH SCALAR STRAIN-RATES AT NODES ! BASED ON THE MODEL (AFTER HEALING FAULTS AND SMOOTHING): DATA iUnitZ / 22 / ! "iUnitX"= Device number to be used for outputting the ! smoothed (interseismic) velocities of each node ! in v_____.out format, for use with FiniteMap. ! (This special file is only created IF(EXPLOD), ! which is true only if the "benchmarks" were ! created by Explode3 and if they match one-for- ! one with nodes of the .feg file.) DATA iUnitX / 23 / ! "iUnitY"= Device number to be used for outputting the ! error in geodetic velocity predictions, in same ! .gps format as the input dataset. ! Note that this occurs only IF (pltGEO) .AND. ! (some geodetic data are provided for scoring). DATA iUnitY / 31 / DATA pltGEO / .TRUE. / ! "iUnitK"= Device number to be used for reading file with ! upper-mantle anisotropy data, in the form of fast- ! polarization ("phi") azimuths from SKS splitting, ! together with SKS splitting times in s. DATA iUnitK / 32 / ! "iUnitQ"= Device number to be used for reading file with ! torque report (obtained by running Shells version ! of summer 2006 or later). This is used in scoring ! the anisotropy data read through iUnitK (above). DATA iUnitQ / 33 / ! "iUnitF"= Device number to be used for reading file with ! outlines of the plates (e.g., PB2002_plates.dig). ! This is used in scoring the anisotropy data read ! through iUnitK (above). DATA iUnitF / 34 / ! PB2002 plate names [Bird, 2003, G**3, Table 1]: 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' / INTEGER, PARAMETER :: ipAfri = 1 ! Index number of Africa plate in this model. ! Following rotation vectors in Cartesian (x,y,z) components, ! with units of radians per million years: ! [Bird, 2003, G**3, Table 1] DATA ((omega(i, j), i = 1, 3), j = 1, nPlate) / & & 0.002401, -0.007939, 0.013892, & ! 1 & 0.000949, -0.008643, 0.013725, & ! 2 & 0.000689, -0.006541, 0.013676, & ! 3 & 0.002042, -0.013153, 0.008856, & ! 4 & 0.008570, -0.005607, 0.017497, & ! 5 & 0.000148, -0.003070, 0.010915, & ! 6 & 0.015696, 0.002467, 0.023809, & ! 7 & 0.009349, 0.000284, 0.016252, & ! 8 & 0.000184, 0.005157, 0.001150, & ! 9 & -0.000871, -0.002268, 0.002507, & ! 10 & -0.019124, 0.030087, 0.010227, & ! 11 & 0.011506, -0.044526, 0.007197, & ! 12 & 0.001688, -0.009048, 0.012815, & ! 13 & 0.003716, -0.003791, 0.000949, & ! 14 & -0.008915, -0.026445, 0.020895, & ! 15 & -0.061175, 0.005216, -0.013755, & ! 16 & 0.070136, 0.160534, 0.094328, & ! 17 & 0.000529, -0.007235, 0.013123, & ! 18 & -0.083251, -0.002464, -0.014923, & ! 19 & 0.016256, 0.089364, 0.015035, & ! 20 & 0.008180, -0.004800, 0.016760, & ! 21 & 0.006512, 0.003176, 0.005073, & ! 22 & 0.108013, 0.299461, 0.230528, & ! 23 & 0.033318, -0.001813, 0.036441, & ! 24 & -0.013835, 0.008245, 0.015432, & ! 25 & -0.777844, 0.440872, -0.047437, & ! 26 & 0.001521, 0.007739, 0.013437, & ! 27 & 0.038223, -0.058291, 0.013679, & ! 28 & 0.001768, -0.008439, 0.009817, & ! 29 & -0.004336, 0.003769, -0.000402, & ! 30 & 0.000111, -0.006361, 0.010449, & ! 31 & 0.044913, -0.009546, 0.010601, & ! 32 & -0.055342, -0.010890, 0.006794, & ! 33 & -0.000022, -0.013417, 0.019579, & ! 34 & 0.001041, -0.008305, 0.012143, & ! 35 & -0.026223, 0.020184, 0.037208, & ! 36 & 0.000000, 0.000000, 0.000000, & ! 37 & -0.000040, -0.009291, 0.012815, & ! 38 & 0.012165, -0.012510, -0.000366, & ! 39 & -0.019183, -0.070604, 0.036798, & ! 40 & 0.000472, -0.006355, 0.009100, & ! 41 & 0.121443, -0.078836, 0.027122, & ! 42 & 0.001117, -0.007434, 0.008534, & ! 43 & -0.000833, -0.006701, 0.013323, & ! 44 & 0.001287, -0.008754, 0.014603, & ! 45 & -0.017196, 0.017186, 0.008623, & ! 46 & 0.003201, -0.010440, 0.015854, & ! 47 & 0.023380, -0.019369, -0.010465, & ! 48 & -0.009400, 0.023063, 0.008831, & ! 49 & 0.142118, 0.005616, 0.078214, & ! 50 & -0.016831, 0.018478, 0.010166, & ! 51 & -0.000836, -0.006169, 0.016274 / ! 52 !--------------------------------------------------------------------- ! STATEMENT FUNCTIONS: ! Interpolation within one fault element: fhival(s, x1, x2) = x1 + s * (x2 - x1) !--------------------------------------------------------------------- ! BEGINNING OF EXECUTABLE CODE WRITE (*, "(' OrbScore2, by Peter Bird, UCLA: version of 26 March 2007')") ! *** 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 mxel = maxel mxfel = maxfel mxbn = maxbn mxstar = maxatp mxgeo = maxgeo mxstr = maxstr mxSKS = maxSKS mxgsr = maxdat mxmor = maxmor mxadj = maxadj ! ****************************************************************** ! Mark most pLates as LACKING extensive attached slabs... DO i = 1, nPlate slab_q(i) = .FALSE. END DO ! ...except for these particular cases: slab_q( 8) = .TRUE. ! 8 = AU = Australia slab_q(14) = .TRUE. ! 14 = CL = Caroline slab_q(15) = .TRUE. ! 15 = CO = Cocos slab_q(21) = .TRUE. ! 21 = IN = India slab_q(22) = .TRUE. ! 22 = JF = Juan de Fuca slab_q(34) = .TRUE. ! 34 = NZ = Nazca slab_q(37) = .TRUE. ! 37 = PA = Pacific slab_q(39) = .TRUE. ! 39 = PS = Philippine Sea slab_q(40) = .TRUE. ! 40 = RI = Rivera slab_q(46) = .TRUE. ! 46 = SS = Solomon Sea ! N.B. This array is used when scoring upper-mantle seismic- ! anisotropy data, because in the Earth5 set of models, ! plates with slab_q(iPlate) = .TRUE. {those with ! extensive attached slabs} are assumed to have ! minimal and/or unresolvable basal shear tractions, ! and are excluded from the scoring process. ! That is, if total basal-strength traction for the ! Pacific plate is toward the Northwest, this might ! be primarily due to the attached slabs, and it does ! not preclude weak, distributed basal shear tractions ! that might point in a completely different direction! ! WRITE HEADER ON OUTPUT FILE WRITE (iUnitT, 1) 1 FORMAT ( & &' GPBOUT:'/ & &' =============================================================='/ & &' I PROGRAM -OrbScore2- I'/ & &' I A spherical-Earth, thin-plate program for scoring I'/ & &' I time-averaged (anelastic) deformation solutions I'/ & &' I computed by finite element program -SHELLS-. I'/ & &' I I'/ & &' I by I'/ & &' I Peter Bird I'/ & &' I Department of Earth and Space Sciences I'/ & &' I University of California I'/ & &' I Los Angeles, CA 90095-1567 I'/ & &' I Version of 16 November 2006 I'/ & &' ==============================================================') wedge = ABS(90. - ABS(dipmax)) * 0.017453293 slide = subdip * 0.017453293 ! INPUT FINITE ELEMENT GRID AND DATA VALUES AT NODE POINTS CALL GetNet (input, iUnitG, iUnitT, & & mxel, mxfel, mxnode, & & output, brief, dqdtda, elev, fdip, & & nfaken, nfl, nodef, nodes, nrealn, & & numel, numnod, n1000, offmax, offset, & & title1, tlnode, xnode, ynode, zmnode, & & work, checke, checkf, checkn) ! READ SCALAR PARAMETERS CALL ReadPm (input, iUnitP, iUnitT, names, nPlate, offmax, & & output, acreep, alphat, bcreep, biot , & & byerly, ccreep, cfric , conduc, & & dcreep, ecreep, everyp, ffric , gmean , & & gradie, iconve, ipvref, & & maxitr, okdelv, oktoqt, onekm, radio , & & radius, refstr, rhoast, rhobar, rhoh2o, & & tadiab, & & taumax, temlim, title3, trhmax, tsurf, & & vtimes, zbasth) ! CHECK TOPOLOGY, AND COMPUTE GEOMETRIC PROPERTIES CALL Square (input, brief, fdip, iUnitT, & & mxbn, mxel, mxfel, mxnode, & & mxstar, nfl, nodef, nodes, & & numel, numnod, radius, wedge, & & modify, xnode, ynode, & & output, area, detj, & & dxs, dys, dxsp, dysp, edgefs, & & edgets, flen, fpflt, fpsfer, & & farg, ncond, nodcon, sita, & & work, checkn, list) sphere = (ncond == 0) ! INPUT VELOCITIES AT NODE POINTS WRITE(iUnitT, 10) iUnitV 10 FORMAT(/' ATTEMPTING TO READ VELOCITIES OF NODES ON UNIT',I3/) CALL OldVel (input, iUnitT, iUnitV, mxnode, numnod, & & output, havenv, title1, title2, title3, v) IF (.NOT.havenv) THEN WRITE (iUnitT, 20) iUnitV 20 FORMAT (/' NO NODAL VELOCITIES FOUND ON UNIT ',I2) CALL PAUSE() STOP END IF !---------------------------------------------------------------- ! COMPUTE RMS VALUE OF VELOCITIES OF NODES ! (PERHAPS RMS VALUE OF SURFACE INTEGRAL WOULD BE MORE ! ELEGANT, BUT THIS IS MUCH EASIER AND ALMOST THE SAME). ! FIRST, IT IS NECESSARY TO REMOVE ANY NET ROTATION. numadj = numnod IF (mxadj < numadj) THEN WRITE (iUnitT, 220) numnod 220 FORMAT (/' IN ORDER TO REMOVE NET ROTATION FROM MODEL', & & ' VELOCITIES' & & /' (SO AS TO ACCURATELY ASSESS RMS VELOCITY)' & & /' IT IS NECESSARY TO SEND A DUMMY SET OF (0,0)' & & /' BENCHMARK VELOCITIES TO SUBPROGRAM -ADJUST-.' & & /' THIS IN TURN REQUIRES THAT PARAMETER MAXADJ' & & /' MUST BE AT LEAST EQUAL TO THE NUMBER OF' & & /' NODES, WHICH CURRENTLY IS: ',I7) CALL PAUSE() STOP END IF DO 250 i = 1, numnod nlat(i) = 90. - xnode(i) * oezopi elon(i) = ynode(i) * oezopi predic(1, i) = v(1, i) predic(2, i) = v(2, i) DATA(1, i) = 0.0 DATA(2, i) = 0.0 250 CONTINUE IF (floats) THEN ALLOCATE ( weights(numnod) ) weights = 1.0 ! whole array WRITE (iUnitT, 251) 251 FORMAT (/' CALLING -ADJUST- WITH EQUAL WEIGHT ON EACH' & & ,' NODE:') CALL Adjust (input, data, elon, iUnitT, mxadj, numnod, nlat, & & radius, weights, & & modify, predic, & & output, elonp, nlatp, rate) DEALLOCATE (weights) END IF sum = 0. DO 290 i = 1, numnod sum = sum + predic(1, i)**2 + predic(2, i)**2 290 CONTINUE rmsv = SQRT(sum / numnod) t2 = 1000. * secpyr * rmsv IF (floats) THEN rate = -rate t = ABS(rate) * 1000. * secpyr * radius WRITE (iUnitT, 295) rate, t, elonp, nlatp, rmsv, t2 295 FORMAT (/' AFTER REMOVAL OF NET ROTATION OF ',1P,E10.2, & & ' RADIANS/SECOND (',0P,F5.1,' MM/YR AT EQUATOR)' & & /' ABOUT A POLE AT ',0P,F7.2,' DEGREES EAST, ',F6.2, & & ' DEGREES NORTH,' & & /' THE RMS NODAL VELOCITY FOR THIS MODEL IS ', & & 1P,E10.2, & & ' (',0P,F7.2,' MM/YR).') ELSE WRITE (iUnitT, 296) rmsv, t2 296 FORMAT (/' THE RMS NODAL VELOCITY FOR THIS MODEL IS ', & & 1P,E10.2, & & ' (',0P,F7.2,' MM/YR).') END IF !------------------------------------------------------------------- ! ********* GEODESY SCORING SECTION ********** ! READ GEODETIC DATA FOR SCORING ! (Dislocation-in-elastic-halfspace corrections for ! temporary locking of brittle parts of faults within ! the model was added in summer 2000 by Zhen Liu and Peter Bird.) ! Note that this does NOT include corrections for ! temporary locking of subduction zones which lie ! along the boundary of a model, unless fault elements ! are used to represent the subduction shear zone.) CALL GetGEO (input, iUnitB, iUnitT, iUnitY, mxgeo, pltGEO, & & output, geotag, geophi, geothe, & & geovel, geosig, geoazi, & & gpsfmt, numgeo) ! Flag "EXPLOD" indicates that number of benchmarks ! is equal to number of nodes. This happens only when ! the "geodesy" (.gps) file was created by Explode3, and ! has "benchmark" locations for every node in the original ! .feg file, but slightly displaced away from nodes. ! In this case, processing is identical, but an extra ! output file is produced to facilitate plotting ! the smoothed (interseismic) velocities. explod = (numgeo == numnod) IF (numgeo <= 0) THEN geodes = 0.0 CLOSE (UNIT = iUnitY, DISP = "DELETE") GO TO 400 END IF WRITE (iUnitT, 301) numgeo, iUnitB 301 FORMAT (/' ',I6,' GEODETIC VELOCITY DATA WERE READ FROM', & & ' UNIT ',I2) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Find weights for geodetic benchmarks, based on associate area: IF (Is_Initialization_Needed()) THEN subdivision = 4 WRITE (*, "(/' Initializing uniform global grid at subdivision ', & & I1,' for area-weighting of data...')") subdivision CALL Initialize_Weighting (subdivision) END IF WRITE (*, "(' Computing area-weights for geodetic benchmarks...')") ALLOCATE ( weights(numgeo) ) CALL Perform_Weighting (number_of_data = numgeo, & & theta_radians = geothe, & & phi_radians = geophi, & ! INTENT(IN) & weights = weights) ! INTENT(OUT) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! FIND PREDICTED VELOCITIES AT BENCHMARKS ! (A) STEADY-STATE COMPONENTS: IF (explod) THEN WRITE (iUnitT, 302) 302 FORMAT(/' Calling FINDIT to determine element containing' & & /' each displaced node; be patient!....') END IF i = 0 303 i = i + 1 theta = geothe(i) phi = geophi(i) CALL Findit (input, mxel, mxnode, nodes, numel, xnode, ynode, & & theta, phi, & & output, iele, s1, s2, s3) IF (iele == 0) THEN ! BENCHMARK IS OUTSIDE OF GRID; DELETE IT. WRITE (iUnitT, 305) geotag(i) 305 FORMAT (' BENCHMARK ',A20,' WAS OUTSIDE THE GRID:', & & ' DELETED.') DO 310 j = i, numgeo - 1 geotag(j) = geotag(j + 1) geophi(j) = geophi(j + 1) geothe(j) = geothe(j + 1) geovel(j) = geovel(j + 1) geosig(j) = geosig(j + 1) geoazi(j) = geoazi(j + 1) 310 CONTINUE numgeo = numgeo - 1 i = i - 1 ELSE thelat = 90. - theta * oezopi thelon = phi * oezopi IF (thelon > 180.)thelon = thelon - 360. IF (thelon < -180.)thelon = thelon + 360. IF (.NOT.explod) THEN showis = .FALSE. IF (showis) THEN WRITE(iUnitT, 315)geotag(i), thelon, thelat, iele, s1, & & s2, s3 315 FORMAT(' ',A20,' (',F7.2,'E, ',F6.2,'N) = #',I5, & & ':(',F5.3,',',F5.3,',',F5.3,')') END IF END IF ! INTERPOLATE V TO (S1,S2,S3) IN ELEMENT #IELE: CALL Projec (input, iele, & & mxel, mxnode, nodes, & & s1, s2, s3, & & xnode, ynode, v, & & output, vtheta, vphi) predic(1, i) = vtheta predic(2, i) = vphi DATA(1, i) = geovel(i) * COS(geoazi(i)) DATA(2, i) = geovel(i) * SIN(geoazi(i)) nlat(i) = 90. - geothe(i) * oezopi elon(i) = geophi(i) * oezopi END IF ! END OF VARIABLE-LENGTH LOOP: IF (i < numgeo) GO TO 303 IF (explod.AND.(numgeo < numnod)) THEN ! Condition (NUMGEO==NUMNOD) that made EXPLOD true is ! no longer valid: explod = .FALSE. WRITE (iUnitT, 319) 319 FORMAT(/' This should not happen OFTEN when benchmark' & & /' locations are really node locations, slightly' & & /' displaced by Explode3.' & & /' Probably the bad benchmark(s) is/are only' & & /' slightly outside the domain of the F-E model.' & & /" Use OrbWeaver to view the Explode3'd" & & /' grid and see why some nodes/benchmarks moved' & & /' out of bounds. Then adjust the corresponding' & & /' benchmark coordinates in the false-benchmark' & & /' data file, and re-run OrbScore2.') CALL PAUSE () STOP END IF !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! FIND PREDICTED VELOCITIES AT BENCHMARKS; ! (B) Coseismic displacement rate due to brittle ! parts of faults within the model domain. ! FIND THICKNESSES OF CRUST AND MANTLE LITHOSPHERE: CALL Interp (input, zmnode, mxel, mxnode, nodes, numel, & & output, zmoho) CALL Interp (input, tlnode, mxel, mxnode, nodes, numel, & & output, tlint) ! COMPUTE TACTICAL VALUES OF LIMITS ON VISCOSITY, AND WEIGHTS FOR ! IMPOSITION OF CONSTRAINTS IN LINEAR SYSTEMS: CALL Limits (input, area, detj, iUnitT, mxel, numel, & & okdelv, radius, refstr, sphere, tlint, & & trhmax, zmoho, & & output, constr, etamax, fmumax, vismax) ! LOCATE THE BRITTLE/DUCTILE TRANSITION IN EACH FAULT ! Begin with very crude initialization: DO 322 i = 1, nfl ztranf(1, i) = MAX(zmnode(nodef(1, i)) / 2., onekm) ztranf(2, i) = MAX(tlnode(nodef(1, i)) / 2., onekm) 322 CONTINUE ! Then search for actual brittle/ductile iteratively: DO 330 imohr = 1, 3 CALL Mohr (input, acreep, alphat, bcreep, biot, byerly, & & ccreep, cfric, conduc, constr, dcreep, dqdtda, & & ecreep, elev, fdip, ffric, fmumax, & & fpflt, farg, gmean, & & mxfel, mxnode, nfl, nodef, & & offmax, offset, & & onekm, radio, rhoh2o, rhobar, & & slide, sphere, taumax, & & tlnode, tsurf, v, wedge, & & zmnode, & & modify, ztranf, & & output, fc, fimudz, fpeaks, fslips, ftstar) 330 CONTINUE ! COMPUTE BENCHMARK VELOCITIES DUE TO ELASTIC DISLOCATIONS ! FROM EARTHQUAKES IN THE BRITTLE PART OF EACH FAULT AT THE RATE ! PREDICTED BY THE FINITE ELEMENT MODEL: IF (explod) THEN WRITE (iUnitT, 340) 340 FORMAT(/' Calling COSEIS to evaluate coseismic' & & /' displacement rate for each node; be patient!....') END IF CALL Coseis (input, farg, fdip, geothe, geophi, mxnode, & & mxfel, numgeo, nfl, nodef, radius, & & v, wedge, xnode, ynode, zmnode, ztranf, & & output, geouth, geouph) ! FOR A TYPICAL YEAR WITH NO EARTHQUAKES, ! COMBINE THE PREDICTED BENCHMARK VELOCITIES ! BY SUBTRACTING THE COSEISMIC PART: DO 350 i = 1, numgeo predic(1, i) = predic(1, i) - geouth(i) predic(2, i) = predic(2, i) - geouph(i) 350 CONTINUE IF (explod) THEN WRITE (iUnitT, 351) 351 FORMAT(//' Creating file with smoothed (interseismic)' & & /' velocities in v_____.out format, to be used' & & /' as input to FiniteMap for visualizing them.') WRITE(iUnitT, 352) iUnitX 352 FORMAT(/' Attempting to create output v____.out file for' & & /' smoothed (interseismic) velocities at benchmarks' & & /' (which were slightly displaced from nodes by' & & /' Explode3) using unit ',I3/) OPEN (UNIT = iUnitX, FILE = '') line = "Explode3 of : "//trim(title1(1:66)) WRITE(iUnitX, "(A)") trim(line) line = "Interseismic: "//trim(title2(1:66)) WRITE(iUnitX, "(A)") trim(line) line = "velocities : "//trim(title3(1:66)) WRITE(iUnitX, "(A)") trim(line) WRITE(iUnitX, "(4ES18.9)")((predic(i, j), i = 1, 2), j = 1, numgeo) ! Note: NUMGEO == NUMNOD, or we wouldn't be here. CLOSE(iUnitX) WRITE (iUnitT, 353) 353 FORMAT(/' Smoothed velocities have been written.' & & /' This run of OrbScore2 will now end, as' & & /' it was not intended for actual scoring,' & & /' but just to evaluate interseismic' & & /' velocities at "benchmarks" which are' & & /' really slightly-displaced nodes.') CALL PAUSE () STOP END IF !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ADJUST PREDICTIONS BY RIGID ROTATION TO BEST-FIT GEODESY: WRITE (iUnitT, 362) 362 FORMAT (/' Calling -ADJUST- with area-weight on each benchmark' & & /' to minimize global area-integral of velocity-error-squared:') IF (floats) THEN CALL Adjust (input, data, elon, iUnitT, mxadj, numgeo, nlat, & & radius, weights, & & modify, predic, & & output, elonp, nlatp, rate) t = ABS(rate) * 1000. * secpyr * radius WRITE (iUnitT, 370) rate, t, elonp, nlatp 370 FORMAT (/' Model predictions adjusted by rotation of ', & & 1P,E10.2, & & ' radians/second (',0P,F7.2,' mm/yr at equator)' & & /' about a pole at ',0P,F7.2,' degrees East, ',F6.2, & & ' degrees North,' & & /' before comparison to geodetic data.') END IF WRITE (iUnitT, 380) 380 FORMAT (// & &' GEODETIC BENCHMARK VELOCITIES VERSUS MODEL PREDICTIONS:'// & &' BENCHMARK VELOCITY (MM/YR) AZIMUTH MODEL.V (MM/YR)' & &,' MODEL.AZ ERR0R (MM/YR) (SIGMAS) area-weight'/ & &' --------- -------- ----- ------- ------- -----', & &' -------- ----- ----- ------ ----------') 381 FORMAT (' ',A20,1P,E9.2,' (',0P,F5.1,') ',F7.1,1P,E10.2,0P, & &' (',F5.1,') ',F8.1,1P,E10.2,' (',0P,F5.1,') (',F5.1,') ',F10.5) sum0s = 0. sum1 = 0. sum1s = 0. sum2 = 0. sum2s = 0. sumn = 0. sumns = 0. DO 390 i = 1, numgeo gvmma = geovel(i) * secpyr * 1000. azim = 180. - geoazi(i) * oezopi IF (azim < 0.0) azim = azim + 360. IF (azim > 360.) azim = azim - 360. vm = SQRT(predic(1, i)**2 + predic(2, i)**2) vmmma = vm * secpyr * 1000. pazim = 180. - atan2f(predic(2, i), predic(1, i)) * oezopi IF (pazim < 0.0) pazim = pazim + 360. IF (pazim > 360.) pazim = pazim - 360. bad = SQRT((predic(1, i) - DATA(1, i))**2 + & & (predic(2, i) - DATA(2, i))**2) badmma = bad * secpyr * 1000. badsig = bad / geosig(i) WRITE (iUnitT, 381) geotag(i), geovel(i), gvmma, azim, & & vm, vmmma, pazim, bad, badmma, badsig, & & weights(i) IF (pltGEO) THEN elonde = geophi(i) / piO180 nlatde = 90.0 - (geothe(i) / piO180) nlatde = MIN(90.0, MAX(-90.0, nlatde)) vemmpa = + (predic(2, i) - DATA(2, i)) * 1000. * secpyr vnmmpa = -(predic(1, i) - DATA(1, i)) * 1000. * secpyr vesigm = 0.1 vnsigm = 0.1 correl = 0.0 frame = "[same]" tag = geotag(i) WRITE (iUnitY, gpsfmt) elonde, nlatde, vemmpa, vnmmpa, & & vesigm, vnsigm, correl, & & trim(frame), trim(tag) END IF IF (badsig > 2.) sum0s = sum0s + weights(i) sum1 = sum1 + weights(i) * bad sum1s = sum1s + weights(i) * badsig sum2 = sum2 + weights(i) * bad**2 sum2s = sum2s + weights(i) * badsig**2 sumn = MAX(sumn, bad) sumns = MAX(sumns, badsig) 390 CONTINUE sum0s = sum0s / numgeo sum1 = sum1 / numgeo sum1mm = sum1 * 1000. * secpyr sum1s = sum1s / numgeo sum2 = SQRT(sum2 / numgeo) sum2mm = sum2 * secpyr * 1000. sum2s = SQRT(sum2s / numgeo) sumnmm = sumn * 1000. * secpyr WRITE (iUnitT, 397) sum0s, sum1, sum1mm, sum1s, & & sum2, sum2mm 397 FORMAT(/' GPBJUMP:'/ & & /' SUMMARY OF GEODETIC PREDICTION ERR0RS:' & & /' FRACTION OF PREDICTIONS WRONG BY OVER 2 SIGMA: ',F5.3 & & /' MEAN VELOCITY ERR0R: ',1P,E10.2,' (',0P,F7.2,' MM/YR)' & & /' MEAN NUMBER OF SIGMAS IN ERR0R: ',F5.2, & & /' RMS VELOCITY ERR0R: ',1P,E10.2,' (',0P,F7.2,' MM/YR)' & &) IF (floats) THEN WRITE (iUnitT, 398) 398 FORMAT(' *** Note: Line above was optimized by -Adjust-' & & ,'. ***') END IF WRITE (iUnitT, 399) sum2s, sumn, sumnmm, sumns 399 FORMAT(' RMS NUMBER OF SIGMAS IN ERR0R: ',0P,F5.2 & & /' WORST VELOCITY ERR0R :',1P,E10.2,' (',0P,F7.2,' MM/YR)' & & /' LARGEST ERR0R, IN SIGMAS: ',0P,F6.2) ! SELECT ONE INDICATOR OF ERR0R (RMS ERROR, AFTER FRAME CHANGE): geodes = sum2mm IF (pltGEO) THEN CLOSE(iUnitY) WRITE(iUnitT, "(/' Wrote ERRORS.gps, for plotting', & & ' geodetic errors with FiniteMap.')") END IF DEALLOCATE (weights) ! ********* END SPACE-GEODESY SCORING SECTION ********** !------------------------------------------------------------------- ! ****** STRESS-DIRECTION SCORING SECTION ******** ! READ MOST-COMPRESSIVE HORIZONTAL PRINCIPAL STRESS AZIMUTHS. ! FOR SCORING. 400 CALL GetSTR (input, iUnitS, iUnitT, mxstr, & & output, strTag, strThe, strPhi, & & strArg, strQua, strReg, numSTR) IF (numSTR <= 0) THEN stress = 0.0 regime = 0.0 GO TO 1301 END IF WRITE (iUnitT, 401) numSTR, iUnitS 401 FORMAT (/' ',I6,' STRESS DIRECTION DATA WERE READ FROM', & & ' UNIT ',I2) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF (Is_Initialization_Needed()) THEN subdivision = 4 WRITE (*, "(/' Initializing uniform global grid at subdivision ', & & I1,' for area-weighting of data...')") subdivision CALL Initialize_Weighting (subdivision) END IF WRITE (*, "(' Computing area-weights for stress direction/regime data...')") ALLOCATE ( weights(numstr) ) weights = 1.0 ! whole array ! CALL Perform_Weighting (number_of_data = numstr, & ! & theta_radians = strthe, & ! & phi_radians = strphi, & ! INTENT(IN) ! & weights = weights) ! INTENT(OUT) ! N. B. Area-weighting turned off here when we switched to use of ! robust_interpolated_stress_for_OrbScore2, 2007.03.12 in the Earth5 project. !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !GPBfrom ! COMPUTE THE (THETA,PHI) AND THEN (X,Y,Z) OF INTEGRATION POINTS DO 450 i = 1, numel DO 420 k = 1, 3 cartvs(1, k) = COS(ynode(nodes(k, i))) * & & SIN(xnode(nodes(k, i))) cartvs(2, k) = SIN(ynode(nodes(k, i))) * & & SIN(xnode(nodes(k, i))) cartvs(3, k) = COS(xnode(nodes(k, i))) 420 CONTINUE DO 440 m = 1, 7 DO 430 j = 1, 3 tempv(j) = 0.0 DO 425 k = 1, 3 tempv(j) = tempv(j) + cartvs(j, k) * points(k, m) 425 CONTINUE 430 CONTINUE CALL Unit (modify, tempv) IF (ABS(tempv(3)) <= 0.5) THEN xip(m, i) = ACOS(tempv(3)) ELSE equpar = SQRT(tempv(1)**2 + tempv(2)**2) xip(m, i) = atan2f(equpar, tempv(3)) END IF yip(m, i) = atan2f(tempv(2), tempv(1)) cartr(1, m, i) = SIN(xip(m, i)) * COS(yip(m, i)) cartr(2, m, i) = SIN(xip(m, i)) * SIN(yip(m, i)) cartr(3, m, i) = COS(xip(m, i)) 440 CONTINUE 450 CONTINUE ! STRAIN RATES AT INTEGRATION POINTS CALL EDot (input, dxs, dys, & & fpsfer, mxel, & & mxnode, nodes, numel, radius, sita, v, & & output, erate) WRITE (iUnitT, 460) 460 FORMAT (// & &' HORIZONTAL MOST-COMPRESSIVE STRESS DIRECTIONS VERSUS ', & & 'MODEL PREDICTIONS:'/ & &/' DATUM E.LONG. N.LATT. ELEMENT POINT', & & ' QUALITY AZIMUTH MODEL.AZ ERR0R', & & ' REGIME E1 E2 E3 ERR MATCH? area-weight' & &/' ----- ------- ------- ------- ----- ------- ', & & '------- -------- -----', & & ' ------ -------- -------- -------- -------- ------ -----------') 461 FORMAT (' ',A5,1X,F7.2,1X,F7.2,1X,I7,1X,I5,1X,F7.0,1X, & & F7.0,1X,F8.0,1X,F5.0, & & 5X,A2,1P,4E9.1,1X,A4,2X,F10.5) !: ACCUMULATOR OF BAD STRESS REGIMES: nbadre = 0 ! simple integer count of bad regimes abadre = 0.0 ! area-weighted error count aregime_denominator = 0.0 !: ACCUMULATORS FOR GENERAL DATA (ALL QUALITIES): sum1 = 0.0 sum1d = 0.0 sum2 = 0.0 sum2d = 0.0 !: ACCUMULATORS FOR A,B-QUALITY DATA ONLY: n_AB_data = 0 asum1 = 0.0 asum1d = 0.0 asum2 = 0.0 asum2d = 0.0 ! NOTE: INTEGRATION POINT MUST BE WITHIN ONE-HALF ELEMENT- ! WIDTH OF DATUM, OR DATUM IS NOT USED. ! Here, the typical element side is determined by ! averaging together all element sides, ! and expressed in radians: sumsid = 0.0 DO 471 i = 1, numel DO 470 k = 1, 3 kp1 = 1 + MOD(k, 3) n1 = nodes(k, i) n2 = nodes(kp1, i) x1 = xnode(n1) x2 = xnode(n2) y1 = ynode(n1) y2 = ynode(n2) ra(1) = SIN(x1) * COS(y1) ra(2) = SIN(x1) * SIN(y1) ra(3) = COS(x1) rb(1) = SIN(x2) * COS(y2) rb(2) = SIN(x2) * SIN(y2) rb(3) = COS(x2) dot = ra(1) * rb(1) + ra(2) * rb(2) + ra(3) * rb(3) CALL Cross (input, ra, rb, & & output, cuvec) crosiz = SQRT(cuvec(1)**2 + cuvec(2)**2 + cuvec(3)**2) side = atan2f(crosiz, dot) sumsid = sumsid + side 470 CONTINUE 471 CONTINUE toler = sumsid / (6.0 * numel) numden = numstr DO 490 i = 1, numstr ! NEXT STRESS DIRECTION DATUM: thelon = oezopi * strphi(i) IF (thelon > 180.) thelon = thelon - 360. IF (thelon < -180.) thelon = thelon + 360. thelat = 90. - oezopi * strthe(i) rs(1) = SIN(strthe(i)) * COS(strphi(i)) rs(2) = SIN(strthe(i)) * SIN(strphi(i)) rs(3) = COS(strthe(i)) ! FIND CLOSEST INTEGRATION POINT IN THE GRID r2min = 9.99E29 DO 478 m = 1, 7 DO 476 j = 1, numel r2 = (rs(1) - cartr(1, m, j))**2 + & & (rs(2) - cartr(2, m, j))**2 + & & (rs(3) - cartr(3, m, j))**2 IF (r2 < r2min) THEN iele = j mip = m r2min = r2 END IF 476 CONTINUE 478 CONTINUE rmin = SQRT(r2min) IF (rmin > toler) THEN WRITE (iUnitT, 480) strtag(i), thelon, thelat, toler 480 FORMAT (' ',A5,1X,F7.2,1X,F7.2, & & ' STRESS DATUM IS MORE THAN ', & & F6.4,' RADIANS FROM NEAREST' & & ,' INTEGRATION POINT: IGNORED.') numden = numden - 1 ELSE ! GET SIGMA-1H DIRECTION FROM E-1 DIRECTION: e11h = erate(1, mip, iele) e22h = erate(2, mip, iele) e12h = erate(3, mip, iele) CALL Prince (input, e11h, e22h, e12h, & & output, e1h, e2h, u1x, u1y, u2x, u2y) ERR = -(e11h + e22h) ! (OR, = -(E1H+E2H) ) e1 = MIN(e1h, e2h, ERR) e3 = MAX(e1h, e2h, ERR) IF ((ERR > e1).AND.(ERR < e3)) THEN e2 = ERR ELSE IF ((e1h > e1).AND.(e1h < e3)) THEN e2 = e1h ELSE e2 = e2h END IF pazim = 180. - oezopi * atan2f(u1y, u1x) IF (pazim < 0.0) pazim = pazim + 180. IF (pazim < 0.0) pazim = pazim + 180. IF (pazim >= 180.) pazim = pazim - 180. IF (pazim >= 180.) pazim = pazim - 180. azim = 180. - oezopi * strArg(i) IF (azim < 0.0) azim = azim + 180. IF (azim < 0.0) azim = azim + 180. IF (azim >= 180.) azim = azim - 180. IF (azim >= 180.) azim = azim - 180. bad = ABS(pazim - azim) IF (bad > 90.) bad = ABS(bad - 180.) IF (bad > 90.) bad = ABS(bad - 180.) IF (bad > 90.) bad = ABS(bad - 180.) IF (bad > 90.) bad = ABS(bad - 180.) sum1 = sum1 + bad * strqua(i) *weights(i) sum1d = sum1d + strqua(i) * weights(i) sum2 = sum2 + strqua(i) * bad**2 * weights(i) sum2d = sum2d + strqua(i) * weights(i) ! MAINTAIN SEPARATE TOTALS FOR A,B-QUALITY ONLY: IF (strqua(i) > 2.9) THEN n_AB_data = n_AB_data + 1 asum1 = asum1 + bad * weights(i) asum1d = asum1d + weights(i) asum2 = asum2 + bad**2 * weights(i) asum2d = asum2d + weights(i) END IF ! CHECK WHETHER STRESS REGIME IS CORRECT: reg = strreg(i) outcom = 'OK? ' IF (reg == 'NF') THEN ! NORMAL FAULTING ! E1 SHOULD BE VERTICAL IF (e1 /= ERR) THEN nbadre = nbadre + 1 abadre = abadre + weights(i) outcom = 'BAD ' ELSE outcom = 'GOOD' END IF ELSE IF (reg == 'SS') THEN ! STRIKE-SLIP FAULTING ! E2 SHOULD BE VERTICAL IF (e2 /= ERR) THEN nbadre = nbadre + 1 abadre = abadre + weights(i) outcom = 'BAD ' ELSE outcom = 'GOOD' END IF ELSE IF (reg == 'TF') THEN ! THRUST FAULTING ! E3 SHOULD BE VERTICAL IF (e3 /= ERR) THEN nbadre = nbadre + 1 abadre = abadre + weights(i) outcom = 'BAD ' ELSE outcom = 'GOOD' END IF ELSE IF ((reg == 'ST').OR.(reg == 'TS')) THEN ! TRANSPRESSION ! ERR SHOULD BE POSITIVE, AND ! THE TWO HORIZONTAL PRINCIPLE STRESSES ! SHOULD HAVE OPPOSITE SIGNS IF ((ERR <= 0.).OR.((e1h * e2h) >= 0.)) THEN nbadre = nbadre + 1 abadre = abadre + weights(i) outcom = 'BAD ' ELSE outcom = 'GOOD' END IF ELSE IF ((reg == 'NS').OR.(reg == 'SN')) THEN ! TRANSTENSION ! ERR SHOULD BE NEGATIVE, AND ! THE TWO HORIZONTAL PRINCIPLE STRESSES ! SHOULD HAVE OPPOSITE SIGNS IF ((ERR >= 0.).OR.((e1h * e2h) >= 0.)) THEN nbadre = nbadre + 1 abadre = abadre + weights(i) outcom = 'BAD ' ELSE outcom = 'GOOD' END IF END IF aregime_denominator = aregime_denominator + weights(i) WRITE (iUnitT, 461) strtag(i), thelon, thelat, iele, mip, & & strqua(i), azim, pazim, bad, & & strreg(i), e1, e2, e3, ERR, outcom, weights(i) END IF 490 CONTINUE sum1 = sum1 / sum1d sum2 = SQRT(sum2 / sum2d) WRITE (iUnitT, 496) numden, sum1, sum2 496 FORMAT (/' GPBJUMP:'/ & & /' SUMMARY OF STRESS DIRECTION ERR0RS:' & & /' Number of stress directions used for scoring: ',I6,'.' & & /' Area- & quality-weighted MEAN error: ',F5.2,' degrees.' & & /' Area- & quality-weighted RMS error: ',F5.2,' degrees.') IF ((asum1d > 0.).AND.(asum2d > 0.)) THEN asum1 = asum1 / asum1d asum2 = SQRT(asum2 / asum2d) WRITE (iUnitT, 497) n_AB_data, asum1, asum2 497 FORMAT (/' OR, using only the ',I6,' A,B-quality data:' & & /' Area-weighted MEAN error: ',F5.2,' degrees.' & & /' Area-weighted RMS error: ',F5.2,' degrees.') END IF perbad = (100. * nbadre) / (1. * numden) WRITE (iUnitT, 498) nbadre, numden, perbad 498 FORMAT (' Regime:',I6,' bad stress regimes out of ',I6, & & ' data = ',F6.2,' % bad.') perbad = (100. * abadre) / aregime_denominator WRITE (iUnitT, 499) abadre, aregime_denominator, perbad 499 FORMAT (' Regime:',F10.3,' area-weighted bad stress regimes out of ',F10.3, & & ' area-weighted tests = ',F6.2,' % bad.') ! ARBITRARILY CHOOSE MEASURE OF ERR0R (MEAN ERROR, IN DEGREES): stress = sum1 DEALLOCATE (weights) ! ********* END STRESS SCORING SECTION ********** !------------------------------------------------------------------- ! ********* SLIP-RATE SCORING SECTION *********** ! READ GEOLOGICAL SLIP RATES ON FAULTS WHICH ARE REPRESENTED BY ! FAULT ELEMENTS IN THE .FEG FILE. ! SCORE STRIKE-SLIP AND RELATIVE-VERTICAL (RELATIVE-THROW) DATA ! WHERE AVAILABLE. ! THE PART OF CODE IS MODIFIED FROM THE THIN-plate SCORING PROGRAM: ! SCORE_AK.FOR ! The data format is given in: Slip_rate_format.txt. ! READ TEST DATASET OF GEOLOGIC SLIP-RATE LIMITS FOR FAULTS 1301 CALL GetSLF(input, iUnitD, iUnitT, mxgsr, & & output, fname, rlat, delvz, fltlon, fltlat, nfdata) IF (nfdata <= 0) THEN slperr = 0.0 GO TO 501 END IF ! COMPUTE SLIP RATES (IN MM/YEAR) AT SPECIFIED POINTS ! AND ASSIGN SCORES AND PROBABILITIES WRITE (iUnitT, 1310) 1310 FORMAT (/' GEOLOGIC SLIPRATE PREDICTIONS AND THEIR ERRORS' & & ,' (in mm/year):') WRITE (iUnitT, 1320) 1320 FORMAT ( & & ' Minimum Maximum Model' & &,' Dextral Minimum Maximum Model Throw ' & &/' Fault name (first 30 bytes) Dextral Dextral Dextral' & &,' Error Chance Throw Throw Throw Error Chance' & &/' ------------------------------ ------- ------- -------' & &,' ------- ------ ------- ------- ------- ------- ------') shrink = 1.0D0 ngsdat = 0 nb = 0 sumb = 0.0 sumb2 = 0.0 bigb = 0.0 DO 1390 j = 1, nfdata ! FIND FAULT ELEMENT CLOSEST TO (FLTLON,FLTLAT): ! XDATUM is Theta, S from North pole, in radians: xdatum = (90.0 - fltlat(j)) * 0.017453292 ! YDATUM is Phi, E from Greenwich meridan, in radians: ydatum = fltlon(j) * 0.017453292 IF (nfl == 0.0) THEN WRITE(iUnitT, 1331) 1331 FORMAT(' ERR0R: NO FAULT ELEMENTS IN THIS .FEG FILE' & & ' TO MATCH FAULT SLIP-RATE DATA.') CALL PAUSE() STOP END IF i = 0 s = 0.0 r2min = 9.99E+29 DO 1340 k = 1, nfl n1 = nodef(1, k) n2 = nodef(2, k) n3 = nodef(3, k) n4 = nodef(4, k) x1 = xnode(n1) x2 = xnode(n2) x3 = xnode(n3) x4 = xnode(n4) y1 = ynode(n1) y2 = ynode(n2) y3 = ynode(n3) y4 = ynode(n4) factor = SIN(0.5 * (x1 + x2)) r2flt = (x2 - x1)**2 + (factor * (y1 - y2))**2 r2n1 = (xdatum - x1)**2 + (factor * (ydatum - y1))**2 r2n2 = (xdatum - x2)**2 + (factor * (ydatum - y2))**2 IF ((r2n1 <= r2flt).OR.(r2n2 <= r2flt)) THEN ! Determine internal coordinate for this fault: IF ((r2n1 <= r2flt).AND.(r2n2 <= r2flt)) THEN st = r2n1 / (r2n1 + r2n2) smallc = SQRT(r2flt) smallb = SQRT(r2n1) smalla = SQRT(r2n2) ! A = angle at the node-1 end: cosa = (r2n1 + r2flt - r2n2) / (2.0 * smallb * smallc) a = ACOS(cosa) ! B = angle at the node-2 end: cosb = (r2n2 + r2flt - r2n1) / (2.0 * smalla * smallc) b = ACOS(cosb) ! R = offset from the fault line: r = smallc * (SIN(a) * SIN(b) / SIN(a + b)) r2 = r * r ELSE IF (r2n1 <= r2flt) THEN st = 0.0 r2 = r2n1 ELSE IF (r2n2 <= r2flt) THEN st = 1.0 r2 = r2n2 END IF ! Replace current best solution? IF (r2 <= r2min) THEN i = k s = st r2min = r2 END IF END IF 1340 CONTINUE IF (i == 0) THEN WRITE(iUnitT, 1341) j, trim(fname(j)), fltlon(j), fltlat(j) 1341 FORMAT(' ERR0R: NO FAULT ELEMENT NEAR TO DATUM ',I4 & & /' ',A & & /' AT (',F10.4,',',F10.4,')') CALL PAUSE() STOP END IF ! ASSUMING THAT I = FAULT ELEMENT IS FOUND, CONTINUE: n1 = nodef(1, i) n2 = nodef(2, i) n3 = nodef(3, i) n4 = nodef(4, i) x1 = xnode(n1) x2 = xnode(n2) x3 = xnode(n3) x4 = xnode(n4) y1 = ynode(n1) y2 = ynode(n2) y3 = ynode(n3) y4 = ynode(n4) vx1 = v(1, n1) vx2 = v(1, n2) vx3 = v(1, n3) vx4 = v(1, n4) vy1 = v(2, n1) vy2 = v(2, n2) vy3 = v(2, n3) vy4 = v(2, n4) delvx = fhival(s, vx1, vx2) - fhival(s, vx4, vx3) delvy = fhival(s, vy1, vy2) - fhival(s, vy4, vy3) s_dp = s angle = chord(farg(1, i), s_dp, farg(2, i)) unitx = COS(angle) unity = SIN(angle) unitbx = SIN(angle) unitby = -COS(angle) spread = delvx * unitbx + delvy * unitby sinist = delvx * unitx + delvy * unity dextra = -sinist ! NOTE CONVERSION OF UNITS TO MM/YEAR IN NEXT LINE rlt = dextra * unitL / unitT opens = spread closes = -spread dip = fhival(s, fdip(1, i), fdip(2, i)) vertic = (ABS(dip - 1.570796) < wedge) IF (vertic) THEN relvz = 0.0 probvz = 1.0 ELSE ! NOTE CONVERSION OF UNITS TO MM/YEAR IN NEXT LINE relvz = closes * ABS(TAN(dip)) * unitL / unitT CALL Likely (input, delvz(1, j), delvz(2, j), relvz, & & output, probvz) END IF CALL Likely (input, rlat(1, j), rlat(2, j), rlt, & & output, probrl) shrink = shrink * probrl * probvz rerr = 0.0 IF ((rlat(1, j) /= 0.0).OR.(rlat(2, j) /= 0.0)) THEN ngsdat = ngsdat + 1 IF ((rlt < rlat(1, j)).AND.(rlat(1, j) /= 0.)) THEN nb = nb + 1 rerr = rlat(1, j) - rlt sumb = sumb + rerr sumb2 = sumb2 + rerr**2 bigb = AMAX1(bigb, rerr) END IF IF ((rlt > rlat(2, j)).AND.(rlat(2, j) /= 0.)) THEN nb = nb + 1 rerr = rlt - rlat(2, j) sumb = sumb + rerr sumb2 = sumb2 + rerr**2 bigb = AMAX1(bigb, rerr) END IF END IF zerr = 0.0 IF ((delvz(1, j) /= 0.0).OR.(delvz(2, j) /= 0.0)) THEN ngsdat = ngsdat + 1 IF ((relvz < delvz(1, j)).AND.(delvz(1, j) /= 0.)) THEN nb = nb + 1 zerr = delvz(1, j) - relvz sumb = sumb + zerr sumb2 = sumb2 + zerr**2 bigb = AMAX1(bigb, zerr) END IF IF ((relvz > delvz(2, j)).AND.(delvz(2, j) /= 0.)) THEN nb = nb + 1 zerr = relvz - delvz(2, j) sumb = sumb + zerr sumb2 = sumb2 + zerr**2 bigb = AMAX1(bigb, zerr) END IF END IF ! Convert some numbers to ? or --- or "agrees" IF (rlat(1, j) == 0.0) THEN c8r1 = ' ?' ELSE WRITE(c8r1, "(F8.3)") rlat(1, j) END IF IF (rlat(2, j) == 0.0) THEN c8r2 = ' ?' ELSE WRITE(c8r2, "(F8.3)") rlat(2, j) END IF IF (delvz(1, j) == 0.0) THEN c8v1 = ' ?' ELSE WRITE(c8v1, "(F8.3)") delvz(1, j) END IF IF (delvz(2, j) == 0.0) THEN c8v2 = ' ?' ELSE WRITE(c8v2, "(F8.3)") delvz(2, j) END IF IF ((rlat(1, j) == 0.0).AND.(rlat(2, j) == 0.0)) THEN c8rerr = ' ?' ELSE IF (rerr == 0.0) THEN c8rerr = ' agrees' ELSE WRITE (c8rerr, "(F8.3)") rerr END IF IF ((delvz(1, j) == 0.0).AND.(delvz(2, j) == 0.0)) THEN c8verr = ' ?' ELSE IF (zerr == 0.0) THEN c8verr = ' agrees' ELSE WRITE (c8verr, "(F8.3)") zerr END IF IF ((rlat(1, j) == 0.0).AND.(rlat(2, j) == 0.0)) THEN c7rcha = ' ------' ELSE WRITE (c7rcha, "(F7.4)") probrl END IF IF ((delvz(1, j) == 0.0).AND.(delvz(2, j) == 0.0)) THEN c7vcha = ' ------' ELSE WRITE (c7vcha, "(F7.4)") probvz END IF IF (vertic) THEN WRITE(iUnitT, 1370) fname(j)(1:30), & & c8r1, c8r2, & & rlt, c8rerr, c7rcha 1370 FORMAT(' ',A,2A8,F8.3,A8,A7, & & ' Throw not used with strike-slip fault.') ELSE WRITE(iUnitT, 1380) fname(j)(1:30), & & c8r1, c8r2, & & rlt, c8rerr, c7rcha, & & c8v1, c8v2, & & relvz, c8verr, c7vcha 1380 FORMAT(' ',A,2A8,F8.3,A8,A7,2A8,F8.3,A8,A7) END IF 1390 CONTINUE ! SUMMARY GEOLOGIC SLIP-RATE STATISTICS FOR THIS MODEL: WRITE (iUnitT, 1392) trim(title1), trim(title2), trim(title3) 1392 FORMAT (/' FOR MODEL WITH TITLES:'/ & & ' ',A/' ',A/' ',A/) enb = FLOAT(nb) / FLOAT(ngsdat) sumb = sumb / FLOAT(ngsdat) sumb2 = SQRT(sumb2 / FLOAT(ngsdat)) WRITE (iUnitT, 1394) enb, sumb, sumb2, bigb, shrink 1394 FORMAT ( & & /' FRACTION of predictions in error =',F6.3,'.' & & /' AVERAGE error of all predictions = ',F10.3,' mm/year.' & & /' Root-Mean-Square error of all predictions =',F10.3,' mm/year.' & & /' LARGEST error of any prediction = ',F10.3,' mm/year.' & & /' Probability that model input parameters and physical equatio' & & ,'ns are consistent' & & /' with the true geology represented by' & & ,' this dataset is',1PD10.3,'.') chance = shrink**(1.D0 / ngsdat) WRITE (iUnitT, 1398) chance 1398 FORMAT( & & ' ',11('====')/ & & ' # Estimate, based on this dataset, of #'/ & & ' # probability that, in a given case, #'/ & & ' # these physical equations and parameter #'/ & & ' # values lead to a predicted slip rate #'/ & & ' # consistent with geologic data (within #'/ & & ' # its uncertainties) is ',F8.5,'. #'/ & & ' ',11('====')/) ! Use root-mean-square error (in mm/year) as summary statistic: slperr = sumb2 ! ****** END GEOLOGIC SLIP-RATE SCORING ******** !------------------------------------------------------------------- ! ****** SPREADING-RATE SCORING SECTION ******** ! READ SEAFLOOR SPREADING RATES ON MID-OCEAN RISES, FOR SCORING. 501 CALL GetMOR (input, iUnitM, iUnitT, mxmor, & & output, mortag, morphi, morthe, & & morvel, morsig, nummor) IF (nummor <= 0) THEN spread = 0.0 GO TO 600 END IF WRITE (iUnitT, 530) 530 FORMAT (// & &' SEAFLOOR SPREADING VELOCITIES VERSUS MODEL PREDICTIONS:'// & &' RIDGE E.LON. N.LAT. VELOCITY (MM/YR) MODEL.V (MM/YR)', & &' ERR0R (MM/YR) (SIGMAS)'/ & &' ----- ------ ------ -------- ----- ------- -----', & &' ----- ----- ------') 531 FORMAT (' ',A5,F7.2,F7.2,1P,E9.2,' (',0P,F5.1,')',1P,E10.2,0P, & &' (',F5.1,')',1P,E9.2,' (',0P,F5.1,') (',F4.1,')') sum0s = 0. sum1 = 0. sum1s = 0. sum2 = 0. sum2s = 0. sumn = 0. sumns = 0. numden = nummor DO 590 i = 1, nummor ! FOR EACH MID-OCEAN RISE SPREADING RATE, thelon = oezopi * morphi(i) thelat = 90.0 - oezopi * morthe(i) ! FIND CLOSEST NORMAL FAULT rs(1) = SIN(morthe(i)) * COS(morphi(i)) rs(2) = SIN(morthe(i)) * SIN(morphi(i)) rs(3) = COS(morthe(i)) rmin = 9.99E29 DO 550 j = 1, nfl dip = 0.5 * (fdip(1, j) + fdip(2, j)) ! DIP MUST BE GREATER THAN 37.5 DEGREES: IF (ABS(dip - 1.570796) < 0.916294) THEN ! DIP MUST BE LESS THAN 64.0 DEGREES: IF (ABS(dip - 1.570796) > 0.448585) THEN ! POSITIONS OF END-NODES OF FAULT ra(1) = SIN(xnode(nodef(1, j))) * & & COS(ynode(nodef(1, j))) ra(2) = SIN(xnode(nodef(1, j))) * & & SIN(ynode(nodef(1, j))) ra(3) = COS(xnode(nodef(1, j))) rb(1) = SIN(xnode(nodef(2, j))) * & & COS(ynode(nodef(2, j))) rb(2) = SIN(xnode(nodef(2, j))) * & & SIN(ynode(nodef(2, j))) rb(3) = COS(xnode(nodef(2, j))) ! VECTOR ALONG FAULT FROM NODE 1 TO 2 vf(1) = rb(1) - ra(1) vf(2) = rb(2) - ra(2) vf(3) = rb(3) - ra(3) flong = SQRT(vf(1)**2 + vf(2)**2 + vf(3)**2) ! VECTOR FROM NODE 1 TO DATUM voff(1) = rs(1) - ra(1) voff(2) = rs(2) - ra(2) voff(3) = rs(3) - ra(3) ! DETERMINE INTERNAL COORDINATE (0..1) OF DATUM frac = (voff(1) * vf(1) + voff(2) * vf(2) + voff(3) * vf(3)) / & & (flong**2) IF (frac < 0.0) THEN frac = 0. r = SQRT(voff(1)**2 + voff(2)**2 + voff(3)**2) ELSE IF (frac > 1.) THEN frac = 1. r = SQRT((rs(1) - rb(1))**2 + & & (rs(2) - rb(2))**2 + & & (rs(3) - rb(3))**2) ELSE CALL Cross (input, vf, voff, output, vt) vtl = SQRT(vt(1)**2 + vt(2)**2 + vt(3)**2) r = vtl / flong END IF IF (r < rmin) THEN rmin = r iele = j fsave = frac END IF END IF END IF 550 CONTINUE IF (rmin > toler) THEN WRITE (iUnitT, 551) i, mortag(i), thelon, thelat, toler 551 FORMAT (' #',I5,' ON RISE ',A5,' AT ',F7.2, & & 'E, ',F6.2,'N IS', & & ' MORE THAN ',F5.3,' RADIANS', & & ' FROM NORMAL FAULT: IGNORED.') numden = numden - 1 ELSE ! COMPUTE SPREADING RATE (IELE,FSAVE) dvat = v(1, nodef(4, iele)) - v(1, nodef(1, iele)) dvap = v(2, nodef(4, iele)) - v(2, nodef(1, iele)) dvbt = v(1, nodef(3, iele)) - v(1, nodef(2, iele)) dvbp = v(2, nodef(3, iele)) - v(2, nodef(2, iele)) dvt = dvat + fsave * (dvbt - dvat) dvp = dvap + fsave * (dvbp - dvap) ! (FIND THETA AND PHI UNIT VECTORS AT THIS POINT) up(1) = -rs(2) up(2) = rs(1) up(3) = 0.0 CALL Unit (modify, up) CALL Cross (input, up, rs, output, ut) voff(1) = dvt * ut(1) + dvp * up(1) voff(2) = dvt * ut(2) + dvp * up(2) voff(3) = dvt * ut(3) + dvp * up(3) ra(1) = SIN(xnode(nodef(1, iele))) * & & COS(ynode(nodef(1, iele))) ra(2) = SIN(xnode(nodef(1, iele))) * & & SIN(ynode(nodef(1, iele))) ra(3) = COS(xnode(nodef(1, iele))) rb(1) = SIN(xnode(nodef(2, iele))) * & & COS(ynode(nodef(2, iele))) rb(2) = SIN(xnode(nodef(2, iele))) * & & SIN(ynode(nodef(2, iele))) rb(3) = COS(xnode(nodef(2, iele))) ! VECTOR ALONG FAULT FROM NODE 1 TO 2 vf(1) = rb(1) - ra(1) vf(2) = rb(2) - ra(2) vf(3) = rb(3) - ra(3) flong = SQRT(vf(1)**2 + vf(2)**2 + vf(3)**2) CALL Cross (input, vf, voff, output, vt) dot = vt(1) * rs(1) + vt(2) * rs(2) + vt(3) * rs(3) spread = dot / flong bad = ABS(spread - morvel(i)) badsig = bad / morsig(i) mvmma = morvel(i) * secpyr * 1000. sprmma = spread * secpyr * 1000. badmma = bad * secpyr * 1000. WRITE (iUnitT, 531) mortag(i), thelon, thelat, & & morvel(i), mvmma, & & spread, sprmma, bad, badmma, badsig IF (badsig > 2.) sum0s = sum0s + 1 sum1 = sum1 + bad sum1s = sum1s + badsig sum2 = sum2 + bad**2 sum2s = sum2s + badsig**2 sumn = MAX(sumn, bad) sumns = MAX(sumns, badsig) END IF 590 CONTINUE sum0s = sum0s / numden sum1 = sum1 / numden sum1mm = sum1 * 1000. * secpyr sum1s = sum1s / numden sum2 = SQRT(sum2 / numden) sum2mm = sum2 * secpyr * 1000. sum2s = SQRT(sum2s / numden) sumnmm = 1000. * secpyr * sumn WRITE (iUnitT, 599) numden, sum0s, sum1, sum1mm, sum1s, & & sum2, sum2mm, sum2s, sumn, sumnmm, sumns 599 FORMAT(/' SUMMARY OF SPREADING PREDICTION ERR0RS:' & & /' NUMBER OF SPREADING RATES USED FOR SCORING: ',I6 & & /' FRACTION OF PREDICTIONS WRONG BY OVER 2 SIGMA: ',F5.3 & & /' MEAN VELOCITY ERR0R: ',1P,E10.2,' (',0P,F7.2,' MM/YR)' & & /' MEAN NUMBER OF SIGMAS IN ERR0R: ',F7.2, & & /' RMS VELOCITY ERR0R : ',1P,E10.2,' (',0P,F7.2,' MM/YR)' & & /' RMS NUMBER OF SIGMAS IN ERR0R: ',0P,F5.2, & & /' WORST VELOCITY ERR0R :',1P,E10.2,' (',0P,F7.2,' MM/YR)' & & /' WORST ERR0R, IN SIGMAS: ',0P,F7.2) ! ARBITRARILY SELECT AN ERR0R MEASURE (MEAN ERROR, IN MM/YEAR) spread = sum1mm ! *********** END SPREADING-RATE SCORING ************ !---------------------------------------------------------------- ! *********** BEGIN SEISMICITY SCORING ************ 600 numeqs = 0 time1 = 9999. * 3.15576E7 time2 = 0. ! ATTEMPT TO READ CATALOG FILE (WHICH MAY NOT BE PRESENT) ! IN FORMAT PRODUCED BY MY SEISMICITY.F90 PROGRAM WRITE(iUnitT, 601) iUnitE 601 FORMAT(/' ATTEMPTING TO READ SEISMIC CATALOG FILE ON UNIT',I3/) 602 READ(iUnitE, 603, iostat = ios) flag10, & & jyear, month, kday, & & ihour, minute, isecon, itenth, & & eqelon(numeqs + 1), eqnlat(numeqs + 1), & & idepth, eqmag(numeqs + 1) ! CHANGED BY Z.L '.' IS CHANGED AS 1X IN FORMAT 602 603 FORMAT(A10, & & I4,1X,I2,1X,I2, 1X, & & I2,1X,I2,1X,I2,1X,I1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) IF (ios /= 0) THEN ! PROBLEM WITH LAST READ (E.G., NO FILE, OR END-OF-FILE) GO TO 610 ELSE numeqs = numeqs + 1 IF (numeqs > maxeqs) THEN WRITE(iUnitT, 604) 604 FORMAT(/' INCREASE PARAMETER MAXEQS TO EQUAL', & & ' THE NUMBER OF EARTHQUAKES IN INPUT FILE.') CALL PAUSE() STOP END IF tsec = jyear * 3.15576E7 + (month - 1) * 2.6298E6 + (kday - 1) * 86400. + & & (ihour - 1) * 3600. + (minute - 1) * 60. + isecon + itenth / 10. ! TIME1 AND TIME2 ARE EARLIEST AND LAST EARTHQUAKE TIMES, ! IN SECONDS. time1 = MIN(time1, tsec) time2 = MAX(time2, tsec) GO TO 602 END IF 610 seismi = 0.0 ! (INITIALIZING ERROR MEASURE IN CASE OF NO CATALOG) IF (numeqs <= 0) GO TO 700 WRITE(iUnitT, 611) numeqs 611 FORMAT(//' -----------------------------------------------' & & /' BEGINNING SEIMICITY SCORING SECTION' & & //' READ CATALOG FILE, FOUND ',I6,' EARTHQUAKES.') ! SMOOTH DELTA-FUNCTION SEISMICITY WITH A GAUSSIAN FILTER, ! AND SCALE TO CONVERT TO A SCALAR STRAIN-RATE; ! EVALUATE THIS STRAIN-RATE AT EACH NODE (FOR PLOTTING). ! ======================= sigma = 60. * 1000. ! ======================= ! SIGMA IS THE SMOOTHING DISTANCE, IN METERS. ! NOTE THAT IT SHOULD ALLOW FOR EPICENTER MISLOCATION AND ! ALSO FOR FINITE SOURCE DIMENSIONS. ! HOWEVER, IT IS NOT USED TO ACHIEVE A REGIONALLY-SMOOTH ! SCALAR STRAIN-RATE FIELD; THIS WILL BE DONE BELOW. sigma2 = sigma**2 vp = 8100. rho = 3300. ! (SI VALUES APPROPRIATE FOR UPPER MANTLE/OCEANIC LITHOSPHERE) vs = vp / SQRT(3.) simu = rho * vs**2 ! SIMU IS THE ELASTIC SHEAR MODULUS, IN PA. thick = 20000. ! THICK IS THE THICKNESS OF THE SEISMOGENIC LAYER, IN M. deltat = time2 - time1 ! DELTAT IS THE LENGTH OF THE CATALOG, IN SECONDS. factor = 1. / (3.14159 * simu * sigma2 * thick * deltat) ! FACTOR IS ABOUT 1.E-35 IN SI; POSSIBLE UNDERFLOW TROUBLE! DO 615 i = 1, numnod EDotnc(i) = 0. uvecn(1, i) = SIN(xnode(i)) * COS(ynode(i)) uvecn(2, i) = SIN(xnode(i)) * SIN(ynode(i)) uvecn(3, i) = COS(xnode(i)) 615 CONTINUE DO 630 j = 1, numeqs IF ((flag10(1:3) == 'ISC').OR. & & (flag10(1:3) == 'IGN')) THEN ! ISC OR IGN CATALOG, SO MAGNITUDES ARE BODY-WAVE ! MAGNITUDES OR MB: eqmb = eqmag(j) eqms = (eqmb - 2.5) / 0.63 ! CONVERSION FROM MB TO SURFACE-WAVE MAGNITUDE (EQMS) ! FROM RICHTER'S ELEMENTARY SEISMOLOGY BOOK. eqm0 = 10.**(1.5 * eqms + 9.14) ! CONVERSION FROM MS TO MOMENT (EQM0) FROM ! EKSTROM & DZIEWONSKI (1998, NATURE, 332, 319-323), ! WHICH UPDATES FORMULA OF ! KANAMORI AND ANDERSON (BSSA, VOL. 65, P. 1073), ! IN WHICH THE CONSTANT WAS +8.5. ! THE FORMULA GIVES EQM0 IN NEWTON-METERS (SI). ELSE IF (flag10(1:10) == 'HarvardCMT') THEN ! MAGNITUDES ARE MOMENT-MAGNITUDES eqmw = eqmag(j) eqm0 = 10.**(1.5 * eqmw + 9.1) ! PER DEFINITION OF MOMENT-MAGNITUDE IN ! HANKS AND KANAMORI (1979, J.G.R., 2348-2350). ELSE IF (flag10(1:2) == 'NZ') THEN ! CHANGED BY Z.L. NEW ZEALAND LOCAL MAGNITUDE ML IS CALCULATED BASED ON ! FOMULA GIVEN BY HAINES IN BSSA,71,275-294,1981. SYSTEMATIC DIFFERENCE ! BETWEEN ML AND BODY MAGNITUDE mb EXISTS AND IS DEPTH-DEPENDENT. SEE AL ! DAVID HARTE AND DAVID VERE-JONES, N.Z.J.GEOL.GEOPHYS.,42,237-253,1999 ! VERY ROUGHLY, mb-ML ~0.2 FOR EVENTS SHALLOWER THAN 40KM BASED ON HISTO ! GIVEN IN HARTE AND VERE-JONES' PAPER. SO WE SIMPLY ADD 0.2 FROM ML TO ! IT INTO mb. THEN CONVERT mb INTO MS AND MW USING ABOVE FORMULA eqmb = eqmag(j) + 0.2 eqms = (eqmb - 2.5) / 0.63 eqm0 = 10.**(1.5 * eqms + 9.14) ELSE WRITE (iUnitT, 617) flag10 617 FORMAT(/' ILLEGAL CATALOG FLAG: ',A) CALL PAUSE() STOP ' ' END IF equvec(1) = COS(eqnlat(j) * 0.01745329) * & & COS(eqelon(j) * 0.01745329) equvec(2) = COS(eqnlat(j) * 0.01745329) * & & SIN(eqelon(j) * 0.01745329) equvec(3) = SIN(eqnlat(j) * 0.01745329) DO 620 i = 1, numnod dot = equvec(1) * uvecn(1, i) + & & equvec(2) * uvecn(2, i) + & & equvec(3) * uvecn(3, i) CALL Cross (input, equvec, uvecn(1, i), & & output, cuvec) crosiz = SQRT(cuvec(1)**2 + cuvec(2)**2 + cuvec(3)**2) arc = ATAN2(crosiz, dot) r2 = (arc * radius)**2 arg = -r2 / sigma2 IF (arg > -80.) THEN ! (AVOID UNDERFLOWS FOR DISTANT NODES) EDotnc(i) = EDotnc(i) + factor * eqm0 * EXP(arg) END IF 620 CONTINUE 630 CONTINUE ! AVERAGE CATALOG STRAIN-RATE TO ELEMENT CENTERS ! AND IMPOSE REASONABLE MINIMUM OF 0.1%/5 GA = 6.E-21 ! (TO PREVENT CORRELATION TEST OF LOGARITHMS FROM ! BEING DOMINATED BY ARTIFICIAL, UNREASONABLE LOW-END VALUES!) DO 640 i = 1, numel EDotec(i) = (EDotnc(nodes(1, i)) + & & EDotnc(nodes(2, i)) + & & EDotnc(nodes(3, i))) / 3. EDotec(i) = MAX(EDotec(i), 6.E-21) 640 CONTINUE ! END SEISMIC-CATALOG SCALAR STRAIN-RATE FIELD (ROUGH!) ! ------------------------------------------------------------ ! BEGIN FINITE-ELEMENT-MODEL SCALAR STRAIN-RATE FIELD ! MODIFY VELOCITY SOLUTION BY AVERAGING TOGETHER THE VELOCITIES ! OF ALL NODES THAT SHARE A COMMON LOCATION, THUS FORCING THE ! VELOCITY FIELD TO BE CONTINUOUS (AND ELIMINATING INFINITE ! STRAIN-RATES IN FAULT ELEMENTS). DO 642 i = 1, numnod checkn(i) = .FALSE. ! (MEANS "NOT YET INVOLVED IN AVERAGING') 642 CONTINUE DO 670 i = 1, nfl DO 660 j1 = 1, 2 nj1 = nodef(j1, i) ! (FAULT ENDS ARE THE ONLY PLACES THAT CAN HAVE DELTA.V'S) 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 ONE dx = xnode(nj1) - xnode(nj2) dy = ynode(nj1) - ynode(nj2) IF (dy > 3.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 dy = dy * SIN(xnode(nj1)) short = SQRT(dx**2 + dy**2) DO 646 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.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 dy = dy * SIN(xnode(nl1)) test = SQRT(dx**2 + dy**2) short = MIN(short, test) END IF 646 CONTINUE ! COLLECT ALL CORNER NODES WITHIN 10% OF THIS toler = short / 10. t2 = toler**2 DO 650 k = 1, numnod IF (.NOT.checkn(k)) THEN dx = xnode(nj1) - xnode(k) dy = ynode(nj1) - ynode(k) IF (dy > 3.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 dy = dy * SIN(xnode(nj1)) r2 = dx**2 + dy**2 IF (r2 < t2) THEN ninsum = ninsum + 1 IF (ninsum > mxstar) THEN WRITE(iUnitT, 648) 648 FORMAT(/' INCREASE VALUE' & & ,' OF PARAMETER MAXATP.') CALL PAUSE() STOP END IF list(ninsum) = k checkn(k) = .TRUE. END IF END IF 650 CONTINUE ! AVERAGING ACTUALLY BEGINS HERE, BASED ON LIST: xsum = 0. ysum = 0. DO 652 k = 1, ninsum xsum = xsum + v(1, list(k)) ysum = ysum + v(2, list(k)) 652 CONTINUE xmean = xsum / ninsum ymean = ysum / ninsum DO 656 k = 1, ninsum v(1, list(k)) = xmean v(2, list(k)) = ymean 656 CONTINUE END IF 660 CONTINUE 670 CONTINUE ! COMPUTE MODEL STRAIN-RATES AFTER REMOVING DELTA.V'S: CALL EDot (input, dxs, dys, & & fpsfer, mxel, & & mxnode, nodes, numel, radius, sita, v, & & output, erate) ! CONVERT ELEMENT-CENTER VALUES TO SCALARS: ! ((E3*-E1*)/2) AT INTEGRATION POINT 1 ! (NOTE: * MEANS STRAIN-RATE PARTITIONED(?) IF E2 /= 0) DO 672 i = 1, numel exx = erate(1, 1, i) eyy = erate(2, 1, i) exy = erate(3, 1, i) CALL Prince (input, exx, eyy, exy, & & output, e1, e2, u1x, u1y, u2x, u2y) ez = -(exx + eyy) IF ((e2 == 0.).AND.(e1 == 0)) THEN EDotem(i) = 1.E-20 ELSE IF ((e2 * ez) > 0.) THEN ! E1 HAS THE UNIQUE SIGN AND IS PARTITIONED e1part = .TRUE. e2part = .FALSE. ezpart = .FALSE. ELSE IF ((e1 * ez) > 0.) THEN ! E2 HAS THE UNIQUE SIGN AND IS PARTITIONED e1part = .FALSE. e2part = .TRUE. ezpart = .FALSE. ELSE ! EZZ HAS THE UNIQUE SIGN AND IS PARTITIONED e1part = .FALSE. e2part = .FALSE. ezpart = .TRUE. END IF ! STRIKE-SLIP RATE (E2ME1) IF (e1part) THEN e2me1 = 2. * ABS(e2) ELSE IF (e2part) THEN e2me1 = 2. * ABS(e1) ELSE e2me1 = ABS(e2 - e1) END IF ! THRUST-FAULTING RATE (EZME1) IF (e1part) THEN ezme1 = 2. * ABS(ez) ELSE IF (ezpart) THEN ezme1 = 2. * ABS(e1) ELSE ezme1 = ABS(ez - e1) END IF ! NORMAL FAULTING RATE (E2MEZ) IF (e2part) THEN e2mez = 2. * ABS(ez) ELSE IF (ezpart) THEN e2mez = 2. * ABS(e2) ELSE e2mez = ABS(e2 - ez) END IF elarge = MAX(e2me1, ezme1, e2mez) EDotem(i) = MAX(elarge / 2., 1.E-20) END IF 672 CONTINUE ! REPEATEDLY SMOOTH SCALAR STRAIN-RATE FIELDS ! BY AVERAGING ELEMENT VALUES ONTO NODES ! (TREATING FAULTS AS HEALED), ! AND THEN INTERPOLATING FROM NODES TO ELEMENT CENTERS. ! THE SAME SMOOTHING OPERATOR IS APPLIED TO BOTH ! CATALOG AND FINITE-ELEMENT-MODEL SCALAR STRAIN-RATE FIELDS! ! ==================================== nbland = 16 ! ==================================== nbland = MAX(nbland, 1) WRITE(iUnitT, 673)nbland 673 FORMAT(/' USING ',I2,' SMOOTHING CYCLES ON BOTH STRAIN-RATES') DO 684 nb = 1, nbland ! (A) EXTRAPOLATE FROM ELEMENTS TO NODES (ACROSS FAULTS) DO 674 i = 1, numnod EDotnc(i) = 0. EDotnm(i) = 0. atnode(i) = 0. 674 CONTINUE DO 676 i = 1, numel DO 675 j = 1, 3 node = nodes(j, i) EDotnc(node) = EDotnc(node) + EDotec(i) EDotnm(node) = EDotnm(node) + EDotem(i) atnode(node) = atnode(node) + 1. 675 CONTINUE 676 CONTINUE DO 678 i = 1, nfl DO 677 j = 1, 2 node1 = nodef(j, i) node2 = nodef(5 - j, i) t1 = EDotnc(node1) t2 = EDotnc(node2) EDotnc(node1) = EDotnc(node1) + t2 EDotnc(node2) = EDotnc(node2) + t1 t1 = EDotnm(node1) t2 = EDotnm(node2) EDotnm(node1) = EDotnm(node1) + t2 EDotnm(node2) = EDotnm(node2) + t1 t1 = atnode(node1) t2 = atnode(node2) atnode(node1) = atnode(node1) + t2 atnode(node2) = atnode(node2) + t1 677 CONTINUE 678 CONTINUE DO 679 i = 1, numnod EDotnc(i) = EDotnc(i) / MAX(atnode(i), 1.) EDotnm(i) = EDotnm(i) / MAX(atnode(i), 1.) 679 CONTINUE ! (B) INTERPOLATE FROM NODES TO ELEMENT CENTERS DO 680 i = 1, numel EDotec(i) = (EDotnc(nodes(1, i)) + & & EDotnc(nodes(2, i)) + & & EDotnc(nodes(3, i))) / 3. EDotem(i) = (EDotnm(nodes(1, i)) + & & EDotnm(nodes(2, i)) + & & EDotnm(nodes(3, i))) / 3. 680 CONTINUE 684 CONTINUE ! ----- OUTPUT SMOOTHED SCALAR STRAIN-RATE FIELDS ! TO TWO .FEG FILES: ---------------------- WRITE(c3, "(I3)") nbland WRITE(iUnitT, 686) sigma, iUnitO 686 FORMAT(/' EQ MAGNITUDES (mb) HAVE BEEN CONVERTED TO MOMENTS,' & & /' THE MOMENTS DISTRIBUTED WITH A GAUSSIAN FILTER OF' & & /' CHARACTERISTIC DISTANCE ',1P,E10.3,' m;' & & /' THEN THE SCALAR STRAIN-RATE WAS DIFFUSIVELY SMOOTHED;' & & /' USE THE .feg FILE WRITTEN TO UNIT ',I3 & & ,' TO PLOT THE RESULTS.') enctit = 'STRAIN-RATES FROM SEISMIC CATALOG, ' & & //'diffusively smoothed (n ='//c3//')' CALL PutNet (input, iUnitO, & & .TRUE., EDotnc, EDotnc, fdip, & & mxel, mxfel, mxnode, n1000, & & nfaken, nfl, nodef, nodes, & & nrealn, numel, numnod, offset, & & enctit, EDotnc, xnode, ynode, EDotnc) WRITE(iUnitT, 688) iUnitZ 688 FORMAT(/' FINITE-ELEMENT MODEL VELOCITIES WERE MADE CONTINUOUS' & & /' BY HEALING ALL ACTIVE FAULTS;' & & /' THEN THE SCALAR STRAIN-RATE WAS DIFFUSIVELY SMOOTHED;' & & /' USE THE .feg FILE WRITTEN TO UNIT ',I3 & & ,' TO PLOT THE RESULTS.') enctit = 'STRAIN-RATES FROM HEALED-FAULT FE MODEL, ' & & //'diffusively smoothed (n ='//c3//')' CALL PutNet (input, iUnitZ, & & .TRUE., EDotnm, EDotnm, fdip, & & mxel, mxfel, mxnode, n1000, & & nfaken, nfl, nodef, nodes, & & nrealn, numel, numnod, offset, & & enctit, EDotnm, xnode, ynode, EDotnm) ! ----------------------------------------------- ! COMPUTE CORRELATION COEFFICIENT BETWEEN: ! *LOG10(SMOOTHED SCALAR STRAIN-RATE FROM SEISMIC CATALOG) ! *LOG10(SMOOTHED SCALAR STRAIN-RATE FROM FAULT-HEALED MODEL) ! CONVERT TO COMMON LOGARITHMS: DO 690 i = 1, numel EDotec(i) = ALOG10(MAX(EDotec(i), 1.E-30)) EDotem(i) = ALOG10(MAX(EDotem(i), 1.E-30)) 690 CONTINUE ! FIND AVERAGES OF LOG10'S OF SCALAR STRAIN-RATES; CATALOG & MODEL avec = 0. avem = 0. sum = 0. DO 693 i = 1, numel avec = avec + area(i) * EDotec(i) avem = avem + area(i) * EDotem(i) sum = sum + area(i) 693 CONTINUE avec = avec / sum avem = avem / sum varc = 0. varm = 0. crossp = 0. DO 694 i = 1, numel varc = varc + area(i) * (EDotec(i) - avec)**2 varm = varm + area(i) * (EDotem(i) - avem)**2 crossp = crossp + area(i) * (EDotec(i) - avec) * (EDotem(i) - avem) 694 CONTINUE correl = crossp / SQRT(varc * varm) WRITE(iUnitT, 695) correl 695 FORMAT(/' CORRELATION COEFFICIENT BETWEEN:' & & /' * LOG10(SMOOTHED SCALAR STRAIN-RATE ' & & ,'FROM SEISMIC CATALOG)' & & /' * LOG10(SMOOTHED SCALAR STRAIN-RATE ' & & ,'FROM FAULT-HEALED MODEL)' & & /' IS: ',F7.4) ! ARBITRARILY SELECT AN ERR0R MEASURE (CORRELATION COEFFICIENT) seismi = correl ! *********** END SEISMICITY SCORING ************ !---------------------------------------------------------------- !GPB_begin_new_code ! ****** UPPER-MANTLE SEISMIC ANISOTROPY SCORING SECTION ******** ! N.B. The basic assumption is that upper-mantle seismic anisotropy ! (expressed as SKS fast-polarization {"phi"} azimuths and ! splitting times in s} is due to simple shear in the ! asthenosphere, and thus is comparable to directions of ! basal shear tractions inferred from a torque-report file ! (but only for slab-less plates, with slab_q(iPlate) = .FALSE.). ! Bear in mind that anisotropy could also be due to several other causes: ! (1) Most-extensional direction of Neogene finite strain in lithosphere. ! (2) Most-extensional direction of ancient finite strain in lithosphere. ! (3) Azimuth of aligned magma-filled dikes (or water-filled cracks). ! and in these cases the data would NOT be comparable as assumed! ! READ fast-polarization azimuths ("phi") and splitting times ! of SKS waves (e.g., compilation of Matt Fouch). 700 CALL GetSKS (iUnitK, iUnitT, mxSKS, & ! input & numSKS, SKS_tag, SKS_theta, SKS_phi, & & SKS_argument, SKS_delay) ! output IF (numSKS <= 0) THEN anisotropy = 0.0 ! summary error measure GO TO 900 END IF WRITE (iUnitT, 701) numSKS, iUnitK 701 FORMAT (/' ',I6,' Fast SKS azimuths and delay times were read from unit ',I2) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF (Is_Initialization_Needed()) THEN subdivision = 4 WRITE (*, "(/' Initializing uniform global grid at subdivision ', & & I1,' for area-weighting of data...')") subdivision CALL Initialize_Weighting (subdivision) END IF WRITE (*, "(' Computing area-weights for SKS azimuth/delay data...')") ALLOCATE ( weights(numSKS) ) CALL Perform_Weighting (number_of_data = numSKS, & & theta_radians = SKS_theta, & & phi_radians = SKS_phi, & ! INTENT(IN) & weights = weights) ! INTENT(OUT) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Read plate outlines, for assigning sites to the correct plate: WRITE (*, 2) iUnitF 2 FORMAT (/' Attempting to read outlines of plates from unit' & & ,I3/) CALL GetPBx (iUnitF, iUnitT, names, nPBnd, nPlate, & ! input & nBoundaryPoints, pLat, pLon) ! output !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (*, "(' Assigning SKS azimuth/delay data to plates...')") DO 720 i = 1, numSKS CALL ThetaPhi_2_pLate (iUnitT, & & names, nPBnd, nBoundaryPoints, & & nPlate, omega, pLat, pLon, & & SKS_theta(i), SKS_phi(i), & ! INTENT(IN) & iPlate) ! INTENT(OUT) SKS_iPlate(i) = iPlate 720 CONTINUE !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (*, "(' Finding basal shear azimuths to compare with SKS azimuths...')") CALL Tractor(iunitQ, iUnitT, nPlate, numSKS, & & slab_q, SKS_phi, SKS_theta, SKS_iPlate, & ! INTENT(IN) & basal_shear_tractions) ! INTENT(OUT) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Prepare to compute azimuths of fastest-horizontal-extension principal strain-rate ! axes, as alternate predictor of fast-polarization azimuths (phi) of SKS data: !GPBto ! COMPUTE THE (THETA,PHI) AND THEN (X,Y,Z) OF INTEGRATION POINTS DO 750 i = 1, numel DO 725 k = 1, 3 cartvs(1, k) = COS(ynode(nodes(k, i))) * & & SIN(xnode(nodes(k, i))) cartvs(2, k) = SIN(ynode(nodes(k, i))) * & & SIN(xnode(nodes(k, i))) cartvs(3, k) = COS(xnode(nodes(k, i))) 725 CONTINUE DO 740 m = 1, 7 DO 730 j = 1, 3 tempv(j) = 0.0 DO 728 k = 1, 3 tempv(j) = tempv(j) + cartvs(j, k) * points(k, m) 728 CONTINUE 730 CONTINUE CALL Unit (modify, tempv) IF (ABS(tempv(3)) <= 0.5) THEN xip(m, i) = ACOS(tempv(3)) ELSE equpar = SQRT(tempv(1)**2 + tempv(2)**2) xip(m, i) = atan2f(equpar, tempv(3)) END IF yip(m, i) = atan2f(tempv(2), tempv(1)) cartr(1, m, i) = SIN(xip(m, i)) * COS(yip(m, i)) cartr(2, m, i) = SIN(xip(m, i)) * SIN(yip(m, i)) cartr(3, m, i) = COS(xip(m, i)) 740 CONTINUE 750 CONTINUE ! STRAIN RATES AT INTEGRATION POINTS CALL EDot (input, dxs, dys, & & fpsfer, mxel, & & mxnode, nodes, numel, radius, sita, v, & & output, erate) ! NOTE: INTEGRATION POINT MUST BE WITHIN ONE-HALF ELEMENT- ! WIDTH OF DATUM, OR DATUM IS NOT USED. ! Here, the typical element side is determined by ! averaging together all element sides, ! and expressed in radians: sumsid = 0.0 DO 756 i = 1, numel DO 753 k = 1, 3 kp1 = 1 + MOD(k, 3) n1 = nodes(k, i) n2 = nodes(kp1, i) x1 = xnode(n1) x2 = xnode(n2) y1 = ynode(n1) y2 = ynode(n2) ra(1) = SIN(x1) * COS(y1) ra(2) = SIN(x1) * SIN(y1) ra(3) = COS(x1) rb(1) = SIN(x2) * COS(y2) rb(2) = SIN(x2) * SIN(y2) rb(3) = COS(x2) dot = ra(1) * rb(1) + ra(2) * rb(2) + ra(3) * rb(3) CALL Cross (input, ra, rb, & & output, cuvec) crosiz = SQRT(cuvec(1)**2 + cuvec(2)**2 + cuvec(3)**2) side = atan2f(crosiz, dot) sumsid = sumsid + side 753 CONTINUE 756 CONTINUE toler = sumsid / (6.0 * numel) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (iUnitT, 760) 760 FORMAT (//' Fast SKS-polarization azimuths versus model predictions (in slab-less plates):'/ & &/' Datum E.Long. N.Latt. Plate Delay,s Azimuth Model.Az1,2 Error Area-weight Slabless plate?' & &/' ----- ------- ------- ----- ------- ------- ----------- ----- ----------- ---------------') 761 FORMAT ( ' ',A5,1X,F7.2,1X,F7.2,4X,A2,1X,F7.2,1X,F7.0,1X,F5.0,',',F5.0,1X,F5.0,1X,F11.4,1X,' .TRUE. ') 762 FORMAT ( ' ',A5,1X,F7.2,1X,F7.2,4X,A2,1X,F7.2,1X,F7.0,1X,' ??, ??',1X,' --',9X,'---',1X,' .FALSE. ') !: Initialize accumulators: nSKS_bad = 0 ! SKS data misfit by more than 20 degrees (L0 numerator) nSKS_used = 0 ! Count of data in slabless plates (L0 denominator) sum1 = 0.0 ! Numerator for L1 (mean size of error) sum1d = 0.0 ! Denominator for L1 sum2 = 0.0 ! Numerator for L2 (RMS size of error) sum2d = 0.0 ! Denominator for L2 DO 790 i = 1, numSKS ! Next SKS datum: thelon = oezopi * SKS_phi(i) IF (thelon > 180.) thelon = thelon - 360. IF (thelon < -180.) thelon = thelon + 360. thelat = 90. - oezopi * SKS_theta(i) !Cartesian unit vector: rs(1) = SIN(SKS_theta(i)) * COS(SKS_phi(i)) rs(2) = SIN(SKS_theta(i)) * SIN(SKS_phi(i)) rs(3) = COS(SKS_theta(i)) !Datum: Azimuth of fast polarization of SKS, in degrees: azimuth = 180. - oezopi * SKS_argument(i) IF (azimuth < 0.0) azimuth = azimuth + 180. IF (azimuth < 0.0) azimuth = azimuth + 180. IF (azimuth >= 180.) azimuth = azimuth - 180. IF (azimuth >= 180.) azimuth = azimuth - 180. IF (slab_q(SKS_iPlate(i))) THEN !No scoring; just echo the datum, without providing any model predictions... WRITE (iUnitT, 762) SKS_tag(i), thelon, thelat, names(SKS_iPlate(i)), SKS_delay(i), azimuth ELSE ! slabless-plate; proceed with scoring! !Predictor1: Azimuth of simple shear in asthenosphere, in degrees: bAzim = 180. - oezopi * ATan2f(basal_shear_tractions(2, i), basal_shear_tractions(1, i)) ! (y = phi = E, x = theta = S) IF (bAzim < 0.0) bAzim = bAzim + 180. IF (bAzim < 0.0) bAzim = bAzim + 180. IF (bAzim >= 180.) bAzim = bAzim - 180. IF (bAzim >= 180.) bAzim = bAzim - 180. !Predictor2: Azimuth of e3 (fastest-extension) horizontal principal strain-rate in lithosphere, in degrees: !(a) Find closest integration point in the grid... r2min = 9.99E29 DO 770 m = 1, 7 DO 765 j = 1, numel r2 = (rs(1) - cartr(1, m, j))**2 + & & (rs(2) - cartr(2, m, j))**2 + & & (rs(3) - cartr(3, m, j))**2 IF (r2 < r2min) THEN iele = j mip = m r2min = r2 END IF 765 CONTINUE 770 CONTINUE rmin = SQRT(r2min) IF (rmin > toler) THEN WRITE (iUnitT, 780) strtag(i), thelon, thelat, toler 780 FORMAT (' ',A5,1X,F7.2,1X,F7.2, & & ' SKS DATUM IS MORE THAN ', & & F6.4,' RADIANS FROM NEAREST' & & ,' INTEGRATION POINT: IGNORED.') e3Azim = bAzim ! so prediction errors will be the same, and bAzim error won't be underbid. ELSE !(b) Get e_1H (fastest-horizontal-shortening) direction: e11h = erate(1, mip, iele) e22h = erate(2, mip, iele) e12h = erate(3, mip, iele) CALL Prince (input, e11h, e22h, e12h, & & output, e1h, e2h, u1x, u1y, u2x, u2y) ERR = -(e11h + e22h) !(or, = -(E1H+E2H) ) e1 = MIN(e1h, e2h, ERR) e3 = MAX(e1h, e2h, ERR) IF ((ERR > e1).AND.(ERR < e3)) THEN e2 = ERR ELSE IF ((e1h > e1).AND.(e1h < e3)) THEN e2 = e1h ELSE e2 = e2h END IF pazim = 180. - oezopi * atan2f(u1y, u1x) IF (pazim < 0.0) pazim = pazim + 180. IF (pazim < 0.0) pazim = pazim + 180. IF (pazim >= 180.) pazim = pazim - 180. IF (pazim >= 180.) pazim = pazim - 180. e3Azim = pazim + 90. END IF IF (e3Azim < 0.0) e3Azim = e3Azim + 180. IF (e3Azim < 0.0) e3Azim = e3Azim + 180. IF (e3Azim >= 180.) e3Azim = e3Azim - 180. IF (e3Azim >= 180.) e3Azim = e3Azim - 180. !Prediction error of Predictor1 (simple shear in asthenosphere): simple_bad = ABS(bAzim - azimuth) IF (simple_bad > 90.) simple_bad = ABS(simple_bad - 180.) IF (simple_bad > 90.) simple_bad = ABS(simple_bad - 180.) IF (simple_bad > 90.) simple_bad = ABS(simple_bad - 180.) IF (simple_bad > 90.) simple_bad = ABS(simple_bad - 180.) !Prediction error of Predictor2 (pure-shear in lithosphere): pure_bad = ABS(e3Azim - azimuth) IF (pure_bad > 90.) pure_bad = ABS(pure_bad - 180.) IF (pure_bad > 90.) pure_bad = ABS(pure_bad - 180.) IF (pure_bad > 90.) pure_bad = ABS(pure_bad - 180.) IF (pure_bad > 90.) pure_bad = ABS(pure_bad - 180.) !Select smaller of two prediction errors: bad = MIN(simple_bad, pure_bad) IF (bad > 20.0) nSKS_bad = nSKS_bad + 1 ! SKS data misfit by more than 20 degrees (L0 numerator) nSKS_used = nSKS_used + 1 ! Count of data in slabless plates (L0 denominator) sum1 = sum1 + bad * SKS_delay(i) * weights(i) ! Numerator for L1 (mean size of error) sum1d = sum1d + SKS_delay(i) * weights(i) ! Denominator for L1 sum2 = sum2 + bad**2 * SKS_delay(i) * weights(i) ! Numerator for L2 (RMS size of error) sum2d = sum2d + SKS_delay(i) * weights(i) ! Denominator for L2 WRITE (iUnitT, 761) SKS_tag(i), thelon, thelat, names(SKS_iPlate(i)), SKS_delay(i), azimuth, bAzim, e3Azim, bad, weights(i) END IF ! slab_q(SKS_iPlate(i)), or NOT <== proceed with scoring 790 CONTINUE perbad = (100. * nSKS_bad) / (1. * nSKS_used) sum1 = sum1 / sum1d sum2 = SQRT(sum2 / sum2d) WRITE (iUnitT, 796) nSKS_used, perbad, sum1, sum2 796 FORMAT (/' GPBJUMP:'/ & & /' SUMMARY OF SEISMIC ANISOTROPY ERR0RS:' & & /' Number of SKS data in slabless plates, used for scoring: ',I6,'.' & & /' Percentage of directions mis-predicted by more than 20 degrees: ',F6.2 & & /' Area- & delay-weighted MEAN error: ',F5.2,' degrees.' & & /' Area- & delay-weighted RMS error: ',F5.2,' degrees.') ! ARBITRARILY CHOOSE MEASURE OF ERR0R (MEAN ERROR, IN DEGREES): anisotropy = sum1 DEALLOCATE (weights) ! ********* END ANISOTROPY SCORING SECTION ********** !GPB_end_new_code !---------------------------------------------------------------- ! ****** SUMMARY OF SCORES ******** 900 WRITE (iUnitT, 901) IF (geodes /= 0.0) WRITE (iUnitT, 910) geodes IF (spread /= 0.0) WRITE (iUnitT, 920) spread IF (stress /= 0.0) WRITE (iUnitT, 930) stress IF (slperr /= 0.0) WRITE (iUnitT, 940) slperr IF (seismi /= 0.0) WRITE (iUnitT, 950) seismi IF (anisotropy /= 0.0) WRITE (iUnitT, 970) anisotropy 901 FORMAT (' ------------------------------------------------' & & //' REVIEW OF SUMMARY MEASURES OF ERROR OR QUALITY:') 910 FORMAT (' GEODETIC VELOCITY (ERROR): ',F7.2) 920 FORMAT (' SEAFLOOR SPREADING RATE (ERROR): ',F7.2) 930 FORMAT (' STRESS DIRECTION (ERROR): ',F7.2) 940 FORMAT (' FAULT SLIP RATE (ERROR): ',F7.2) 950 FORMAT (' (SMOOTHED) SEISMICITY CORRELATION (QUALITY): ',F7.4) 970 FORMAT (' SEISMIC ANISOTROPY (ERROR): ',F7.2) WRITE (iUnitT, 999) trim(title1), trim(title2), trim(title3) 999 FORMAT (/' FOR MODEL:' & & /' ',A & & /' ',A & & /' ',A/) CALL PAUSE() STOP END !GPBEND SUBROUTINE Adjust (input, data, elon, iUnitT, mxadj, n, nlat, & & radius, weights, & & modify, predic, & & output, elonp, nlatp, rate) ! COMPARES PREDICTED (HORIZONTAL) VELOCITIES "predic" ! (WHOSE FIRST COMPONENT IS THE SOUTHWARD OR THETA VELOCITY, ! AND WHOSE SECOND COMPONENT IS THE EASTWARD OR PHI VELOCITY) ! WITH GEODETIC HORIZONTAL VELOCITIES "data" (IN THE SAME FORMAT) ! which have assigned (dimensionless) "weights" of mean value 1.0, ! AT A SET OF POINTS i=1, ... ,"n" <= "mxadj" ! (DEFINED BY "elon(i)" AND "nlat(i)", IN DEGREES EAST & NORTH) ! AND FINDS THE ROTATION CORRECTION WHICH MINIMIZES THE RMS ! OF THE MAGNITUDE OF THE VELOCITY DIFFERENCES. ! IT THEN ADDS THIS CORRECTION TO "predic", AND REPORTS ! THE POLE IN USER-FRIENDLY COORDINATES "elonp" AND "nlatp" ! (IN DEGREES EAST AND NORTH) ! AND GIVES THE ROTATION RATE "rate" IN RADIANS/SECOND. ! SEE NASA NOTES FOR 23 SEPTEMBER 1994. REAL elon, nlat, elonp, nlatp DOUBLE PRECISION coef, equat, right, solut DIMENSION coef(3, 3), DATA(2, mxadj), elon(mxadj), & & nlat(mxadj), predic(2, mxadj), right(3), solut(3) REAL, DIMENSION(mxadj), INTENT(IN) :: weights DATA piO180 / 0.0174533 / DATA oezopi / 57.2958 / IF (n > mxadj) THEN WRITE (iUnitT, 1) mxadj, n 1 FORMAT (//' ERR0R: PARAMETER MXADJ (NOW ',I6,')' & & /' MUST BE INCREASED TO AT LEAST ',I6) CALL PAUSE() STOP END IF coef(1, 1) = 0.0D0 coef(1, 2) = 0.0D0 coef(1, 3) = 0.0D0 coef(2, 2) = 0.0D0 coef(2, 3) = 0.0D0 coef(3, 3) = 0.0D0 right(1) = 0.0D0 right(2) = 0.0D0 right(3) = 0.0D0 DO 100 i = 1, n theta = (90. - nlat(i)) * piO180 phi = elon(i) * piO180 coslat = SIN(theta) coslon = COS(phi) sinlat = COS(theta) sinlon = SIN(phi) rx = radius * coslat * coslon ry = radius * coslat * sinlon rz = radius * sinlat rzu = sinlat a = -rz * rzu * sinlon - ry * coslat b = -rz * coslon c = + rz * rzu * coslon + rx * coslat d = -rz * sinlon e = -ry * rzu * coslon + rx * rzu * sinlon f = + ry * sinlon + rx * coslon coef(1, 1) = coef(1, 1) + weights(i) * ( a * a + b * b ) coef(1, 2) = coef(1, 2) + weights(i) * ( a * c + b * d ) coef(1, 3) = coef(1, 3) + weights(i) * ( a * e + b * f ) coef(2, 2) = coef(2, 2) + weights(i) * ( c * c + d * d ) coef(2, 3) = coef(2, 3) + weights(i) * ( c * e + d * f ) coef(3, 3) = coef(3, 3) + weights(i) * ( e * e + f * f ) right(1) = right(1) + weights(i) * & & (a * (DATA(1, i) - predic(1, i)) + & & b * (DATA(2, i) - predic(2, i))) right(2) = right(2) + weights(i) * & & (c * (DATA(1, i) - predic(1, i)) + & & d * (DATA(2, i) - predic(2, i))) right(3) = right(3) + weights(i) * & & (e * (DATA(1, i) - predic(1, i)) + & & f * (DATA(2, i) - predic(2, i))) 100 CONTINUE coef(2, 1) = coef(1, 2) coef(3, 1) = coef(1, 3) coef(3, 2) = coef(2, 3) !CCCC DO 110 I=1,3 !CCCC WRITE (6,101) (COEF(I,J),J=1,3),RIGHT(I) !C101 FORMAT (/' ',1P,3E11.3,' X ',E11.3) !C110 CONTINUE ! USING A ROUTINE FROM IMSL = ! International Mathematics Sub routine Library: CALL Dlsarg (3, coef, 3, right, 1, solut) rate = SQRT(solut(1)**2 + solut(2)**2 + solut(3)**2) equat = SQRT(solut(1)**2 + solut(2)**2) elonp = oezopi * ATAN2(solut(2), solut(1)) nlatp = oezopi * ATAN2(solut(3), equat) DO 200 i = 1, n theta = (90. - nlat(i)) * piO180 phi = elon(i) * piO180 coslat = SIN(theta) coslon = COS(phi) sinlat = COS(theta) sinlon = SIN(phi) rx = radius * coslat * coslon ry = radius * coslat * sinlon rz = radius * sinlat rzu = sinlat vx = solut(2) * rz - solut(3) * ry vy = solut(3) * rx - solut(1) * rz vz = solut(1) * ry - solut(2) * rx tx = rzu * coslon ty = rzu * sinlon tz = -coslat vtheta = vx * tx + vy * ty + vz * tz px = -sinlon py = coslon pz = 0. vphi = vx * px + vy * py predic(1, i) = predic(1, i) + vtheta predic(2, i) = predic(2, i) + vphi 200 CONTINUE RETURN END DOUBLE PRECISION FUNCTION Atanpv (y, x) ! RETURNS PRINCIPAL VALUE OF INVERSE TANGENT OF (Y,X) IMPLICIT DOUBLE PRECISION (a - h, o - z) IF (x >= 0.) THEN xx = x yy = y ELSE xx = -x yy = -y END IF IF(yy == 0.0.OR.xx == 0.0) THEN IF(yy == 0.0) Atanpv = 0.0 IF(xx == 0.0) Atanpv = 1.570796327D0 ELSE Atanpv = ATAN2(yy, xx) END IF RETURN END 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 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 bd2 ! DEFINE "FPHI" (NODAL FUNCTIONS) AND "FGAUSS" (GAUSSIAN INTEGRATION ! WEIGHTS) AT THE 7 INTEGRATION POINTS IN EACH FAULT ELEMENT, ! DEFINED BY INTERNAL COORDINATE "FPOINT(M=1,...,7)" ! WHICH CONTAINS THE RELATIVE POSITION ! (FRACTIONAL LENGTH) OF 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 SEVEN 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 OR +1 OR -1 FOR THE TWO SIDES OF THE FAULT. DATA fphi / & & 0.9745539D0, 0.0254461D0, -0.0254461D0, -0.9745539D0, & & 0.8707656D0, 0.1292344D0, -0.1292344D0, -0.8707656D0, & & 0.7029226D0, 0.2970774D0, -0.2970774D0, -0.7029226D0, & & 0.5000000D0, 0.5000000D0, -0.5000000D0, -0.5000000D0, & & 0.2970774D0, 0.7029226D0, -0.7029226D0, -0.2970774D0, & & 0.1292344D0, 0.8707656D0, -0.8707656D0, -0.1292344D0, & & 0.0254461D0, 0.9745539D0, -0.9745539D0, -0.0254461D0 / END REAL FUNCTION Chord (angle1, s, angle2) ! RETURNS AN ANGLE OBTAINED BY INTERPOLATION BETWEEN ANGLE1 ! AND ANGLE2. THE INTERPOLATION METHOD IS NOT SENSITIVE TO ! POSSIBLE CYCLE SHIFTS (OF 2*N*PI) BETWEEN ANGLE1 AND ANGLE2. ! UNIT VECTORS ARE CONSTRUCTED FOR ANGLE1 AND ANGLE2, AND A ! LINEAR CHORD IS DRAWN BETWEEN THEIR TIPS. ! DOUBLE PRECISION S IS THE INTERNAL COORDINATE ALONG THE CHORD; ! IT IS DIMENSIONLESS, WITH VALUE 0.0 AT ANGLE1 AND 1.0 AT ! ANGLE2. (THE USER MAY INPUT S VALUES OUTSIDE THIS RANGE ! TO GET RESULTS OUTSIDE THE (SMALLER) ANGLE BETWEEN ANGLE1 AND ! ANGLE2, IF DESIRED.) THE ANGLE RETURNED IS THAT FROM THE ! ORIGIN TO THIS CHORD POINT. ! THIS ALGORITHM SHOULD WORK EQUALLY WELL FOR ANGLES MEASURED ! EITHER CLOCKWISE OR COUNTERCLOCKWISE FROM ANY REFERENCE, AS ! LONG AS THE USAGE IS CONSISTENT. ! BOTH THE INPUT ANGLES AND THE RESULT "CHORD" ARE IN RADIANS. DOUBLE PRECISION s REAL angle1, angle2, uvec1, uvec2 DIMENSION uvec1(2), uvec2(2), vecs(2) uvec1(1) = COS(angle1) uvec1(2) = SIN(angle1) uvec2(1) = COS(angle2) uvec2(2) = SIN(angle2) vecs(1) = (1.0D0 - s) * uvec1(1) + s * uvec2(1) vecs(2) = (1.0D0 - s) * uvec1(2) + s * uvec2(2) Chord = Atan2f(vecs(2), vecs(1)) RETURN END SUBROUTINE Aura (input, u4, theta4, x4, l4, dbot4, dtop4, & & output, v4) ! COMPUTES DISPLACEMENTS IN AN INFINITE UNIFORM ELASTIC HALFSPACE ! OF POISSON'S RATIO 0.25 CAUSED BY UNIFORM SLIP OVER THE SURFACE ! OF A BURIED RECTANGULAR DISLOCATION WITH HORIZONTAL AND PLUNGING ! SIDES. TO CALCULATE DISPLACEMENTS CAUSED BY MORE GENERAL SLIP, ! USE THIS ROUTINE MANY TIMES IN THE KERNEL OF AN INTEGRAL, ! ASSIGNING THE LOCAL AVERAGE SLIP OVER EACH VERY SMALL RECTANGLE. ! EQUATIONS BASED ON MANSINHA AND SMYLIE (1971), ! BULL.SEISM.SOC.AMER., 61, 1433-1440. IMPLICIT DOUBLE PRECISION (a - z) REAL output INTEGER i, input DIMENSION u(3), u4(3), v(3), v4(3), x(3), x4(3) DO 1 i = 1, 3 u(i) = u4(i) x(i) = x4(i) 1 CONTINUE theta = theta4 l = l4 dbot = dbot4 dtop = dtop4 ! U(3) IS THE DISLOCATION VECTOR (I.E. THE TOTAL SLIP VECTOR) ! OF THE HANGING WALL RELATIVE TO THE FOOTWALL, ! IN A CARTESIAN COORDINATE SYSTEM WHERE X1 IS PARALLEL TO THE ! FAULT STRIKE, X2 IS ALSO HORIZONTAL AND POINTING TO HANGING ! WALL, AND X3 IS DOWN (TAKE CARE TO MAKE THESE RIGHT-HANDED). ! THETA IS THE FAULT DIP ANGLE (IN RADIANS) FROM +X2 AXIS. ! X(3) GIVES THE OBSERVATION POINT IN THESE COORDINATES. ! L IS THE HALF-LENGTH OF THE RECTANGLE, ALONG STRIKE. ! DBOT AND DTOP ARE OBLIQUE DISTANCES IN THE FAULT PLANE (AND ! DOWN THE DIP) FROM THE SURFACE TO THE BOTTOM AND TOP OF THE ! RECTANGULAR DISLOCATION, RESPECTIVELY. ! NOTE: FAULT DOES NOT BREAK THE SURFACE, SO DTOP > 0., ! ALTHOUGH IT MAY BE SMALL. ! V(3) IS THE ANSWER: THE 3-COMPONENT DISPLACEMENT AT THE ! OBSERVATION POINT IN SAME UNITS AS U AND IN X1,X2,X3 COORDINATES. sinth = SIN(theta) costh = COS(theta) tanth = sinth / costh ! NOTE: FAULT MAY NOT BE VERTICAL IN THIS ROUTINE. secth = 1. / costh r2 = x(2) * sinth - x(3) * costh r3 = x(2) * costh + x(3) * sinth q2 = x(2) * sinth + x(3) * costh q3 = -x(2) * costh + x(3) * sinth DO 10 i = 1, 3 10 v(i) = 0.D0 DO 100 i = 1, 4 factor = 1.00D0 IF (i == 2.OR.i == 3) factor = -1.00D0 xi1 = l IF(i >= 3) xi1 = -l xi = dbot IF(i == 2.OR.i == 4) xi = dtop xi2 = xi * costh xi3 = xi * sinth rsquar = (x(1) - xi1)**2 + (x(2) - xi2)**2 + (x(3) - xi3)**2 qsquar = (x(1) - xi1)**2 + (x(2) - xi2)**2 + (x(3) + xi3)**2 k = SQRT(qsquar - (q3 + xi)**2) h = SQRT(qsquar - (x(1) - xi1)**2) ucap1 = u(1) ucap = SQRT(u(2)**2 + u(3)**2) IF(u(2) < 0.)ucap = -ucap q = SQRT(qsquar) r = SQRT(rsquar) q = MAX(q, 1.00001D0 * ABS(x(1) - xi1)) r = MAX(r, 1.00001D0 * ABS(x(1) - xi1)) term1 = (x(1) - xi1) * (2. * r2 / (r * (r + r3 - xi)) & & - (4. * q2 - 2. * x(3) * costh) / (q * (q + q3 + xi)) & & - 3. * tanth / (q + x(3) + xi3) & & + 4. * q2 * x(3) * sinth / q**3 & & - 4. * q2 * q3 * x(3) * sinth * (2. * q + q3 + xi) / & & (q**3 * (q + q3 + xi)**2)) & & - 6. * tanth**2 * Atanpv(((k - q2 * costh) * (q - k) + & & (q3 + xi) * k * sinth), ((x(1) - xi1) * (q3 + xi) * costh)) & & + 3. * Atanpv((x(1) - xi1) * (r3 - xi), (r2 * r)) & & - 3. * Atanpv((x(1) - xi1) * (q3 + xi), (q2 * q)) term2 = sinth * (3. * tanth * secth * log(q + x(3) + xi3) - & & log(r + r3 - xi) - (1. + 3. * tanth**2) * log & & (q + q3 + xi)) + 2. * r2**2 * sinth / (r * (r + r3 - xi)) & & + 2. * r2 * costh / r - 2. * sinth * & & (2. * x(3) * (q2 * costh - q3 * sinth) + q2 * (q2 + x(2) * sinth)) & & / (q * (q + q3 + xi)) - 3. * tanth * (x(2) - xi2) / & & (q + x(3) + xi3) + 2. * (q2 * costh - q3 * sinth - & & x(3) * sinth**2) / q + 4. * q2 * x(3) * sinth * & & ((x(2) - xi2) + q3 * costh) / q**3 - 4. * q2**2 * q3 * x(3) & & * sinth**2 * (2. * q + q3 + xi) / (q**3 * (q + q3 + xi)**2) ! TERM3=COSTH*(LOG(R+R3-XI)+(1.+3.*TANTH**2)*LOG( ! 1 Q+Q3+XI)-3.*TANTH*SECTH*LOG(Q+X(3)+XI3)) ! 2 +2.*R2*SINTH/R+2.*SINTH*(Q2+X(2)*SINTH)/Q ! 3 -2.*R2**2*COSTH/(R*(R+R3-XI))+ ! 4 (4.*Q2*X(3)*SINTH**2-2.*(Q2+X(2)*SINTH)*(X(3)+Q3*SINTH))/ ! 5 (Q*(Q+Q3+XI))+ ! 6 4.*Q2*X(3)*SINTH*(X(3)+XI3-Q3*SINTH)/Q**3 ! 7 -4.*Q2**2*Q3*X(3)*COSTH*SINTH*(2.*Q+Q3+XI)/ ! 8 (Q**3*(Q+Q3+XI)**2) term4 = (x(2) - xi2) * sinth * (2. / r + 4. / q - 4. * xi3 * x(3) / q**3 & & - 3. / (q + x(3) + xi3)) - costh * (3. * log(q + x(3) + xi3) & & + 2. * (x(3) - xi3) / r + 4. * (x(3) - xi3) / q + 4. * xi3 * x(3) * (x(3) + xi3) & & / q**3) + 3. * (log(q + x(3) + xi3) - sinth * log(q + & & q3 + xi)) / costh + 6. * x(3) * (costh / q - q2 * sinth / & & (q * (q + q3 + xi))) term5 = sinth * (-log(r + x(1) - xi1) + log(q + x(1) - xi1) & & + 4. * xi3 * x(3) / (q * (q + x(1) - xi1)) + 3. * (x(1) - xi1) / & & (q + x(3) + xi3) + (x(2) - xi2)**2 * (2. / (r * (r + x(1) - xi1)) & & + 4. / (q * (q + x(1) - xi1)) - 4. * xi3 * x(3) * (2. * q + x(1) - xi1) / & & (q**3 * (q + x(1) - xi1)**2))) & & - costh * ((x(2) - xi2) * (2. * (x(3) - xi3) / (r * (r + x(1) - xi1)) & & + 4. * (x(3) - xi3) / (q * (q + x(1) - xi1)) + 4. * xi3 * x(3) * & & (x(3) + xi3) * (2. * q + x(1) - xi1) / (q**3 * (q + x(1) - xi1)**2)) & & + 6. * Atanpv((x(1) - xi1) * (x(2) - xi2), ((h + x(3) + xi3) * (q + h))) & & - 3. * Atanpv((x(1) - xi1) * (r3 - xi), (r2 * r)) + & & 6. * Atanpv((x(1) - xi1) * (q3 + xi), (q2 * q))) + & & 6. * (secth * Atanpv(((k - q2 * costh) * (q - k) + (q3 + xi) * k * sinth) & & , ((x(1) - xi1) * (q3 + xi) * costh)) & & + x(3) * (((sinth**2 - costh**2) * (q3 + xi) + 2. * q2 * & & costh * sinth) / (q * (q + x(1) - xi1)) + (x(1) - xi1) * sinth**2 & & / (q * (q + q3 + xi)))) ! TERM6=SINTH*((X(2)-XI2)*(2.*(X(3)-XI3)/(R*(R+X(1)-XI1)) ! 1 +4.*(X(3)-XI3)/(Q*(Q+X(1)-XI1))-4.*XI3*X(3)* ! 2 (X(3)+XI3)*(2.*Q+X(1)-XI1)/(Q**3*(Q+X(1)-XI1)**2)) ! 3 -6.*ATANPV((X(1)-XI1)*(X(2)-XI2),((H+X(3)+XI3)*(Q+H))) ! 4 +3.*ATANPV((X(1)-XI1)*(R3-XI),(R2*R))- ! 5 6.*ATANPV((X(1)-XI1)*(Q3+XI),(Q2*Q))) ! 6 +COSTH*(LOG(R+X(1)-XI1)-LOG(Q+X(1)-XI1)- ! 7 2.*(X(3)-XI3)**2/(R*(R+X(1)-XI1)) ! 8 -4.*((X(3)+XI3)**2-XI3*X(3))/(Q*(Q+X(1)-XI1)) ! 9 -4.*XI3*X(3)*(X(3)+XI3)**2*(2.*Q+X(1)-XI1)/ ! A (Q**3*(Q+X(1)-XI1)**2)) ! B +6.*X(3)*(COSTH*SINTH*(2.*(Q3+XI)/(Q*(Q+X(1)-XI1)) ! C +(X(1)-XI1)/(Q*(Q+Q3+XI))) ! D -Q2*(SINTH**2-COSTH**2)/(Q*(Q+X(1)-XI1))) deld = dbot - dtop v(1) = v(1) + factor * (term1 * ucap1 + term4 * ucap) / 37.69911184D0 v(2) = v(2) + factor * (term2 * ucap1 + term5 * ucap) / 37.69911184D0 ! V(3)=V(3)+FACTOR*(TERM3*UCAP1+TERM6*UCAP)/37.69911184D0 100 CONTINUE v4(1) = v(1) v4(2) = v(2) ! V4(3)=V(3) RETURN END SUBROUTINE Change (input, argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & radius, & & slip, & & wedge, & & ztop, zbot, & & output, duthet, duphi) ! This routine is a "driver" or "wrapper" for the ! Mansinha & Smylie routine HALO (vertical dislocation patch) ! or AURA (dipping dislocation patch) ! which permits them to be called from spherical-planet ! programs because it converts coordinates from ! (theta, phi) = (colatitude, longitude) in radians ! to the (Z1,Z2,Z3) Cartesian units of a locally flat planet. ! Also, output (horizontal) vector (DUTHET,DUPHI) is expressed ! in the spherical-planet (theta,phi) system, so components are ! (South, East). Units are the same as input vector SLIP, ! which may be either a relative displacement or ! a relative velocity of the two sides of the fault, ! FTHETA and FPHI should give (theta, phi) coordinates of the ! midpoint of the trace of the plane containing the ! rectangular dislocation patch. ! BTHETA and BPHI should give (theta, phi) coordinates of the ! benchmark or test point at which the displacement (-rate) ! is desired. ! LF is the half-length (from center to end) of the dislocation ! patch, measured along a horizontal strike line. ! Units are the same as RADIUS, ZTOP, ZBOT (below). ! ARGUME is the argument of the trace of the dislocation patch, ! measured in radians counterclockwise from +theta (from South). ! DIPF is fault dip in radians measured clockwise ! (initially, down) from horizontal on the right side of the ! fault (when looking in direction ARGUME). ! SLIP must be already rotated into fault coordinates: ! SLIP(1) is the component parallel to the fault, ! and it is positive for sinistral sense of slip. ! SLIP(2) is the horizontal component in the direction ! perpendicular to the trace of the fault, ! and it is positive for divergence across the trace. ! SLIP(3) is the relative vertical component, and it ! is positive when the right side of the fault ! (when looking along direction ARGUME) is down. ! WEDGE is a tolerance (in radians) for dip; if DIPF is ! within WEDGE of Pi/2, then fault is considered vertical ! and routine HALO is called; otherwise, AURA is called. ! ZTOP and ZBOT are (positive) depths to top and bottom of ! the dislocation patch, respectively. Units are the same ! as RADIUS, which gives the size of the planet. ! NOTE: If distance from "F" point to "B" point exceeds ! 50*ZBOT, result is approximated as zero and ! HALO/AURA are never called. ! Following must be DP to agree with HALO, AURA: DOUBLE PRECISION dbot, dtop, lf, theta, u, ucap, x, z ! Following are DP for precision in determination of ! variables FAR and ALPHA. (At short distances, ! unit vectors must be DP or location precision is ! lost and benchmark can end up on wrong side of fault!) DOUBLE PRECISION buvec, crossp, dot1, dot2, dotpro, dvec, & & fuvec, tvec, uphi, utheta INTEGER input LOGICAL vert REAL output DIMENSION buvec(3), dvec(3), fuvec(3), slip(3), tvec(3), & & u(3), ucap(3), uphi(3), utheta(3), x(3), z(3) ! Compute position relative to fault trace; ! this is where we convert from a spherical-Earth ! to a flat-Earth model, implicitly using a gnomonic ! projection (which preserves distance and azimuth ! from the center of the trace of the dislocation plane). fuvec(1) = SIN(dble(ftheta)) * COS(dble(fphi)) fuvec(2) = SIN(dble(ftheta)) * SIN(dble(fphi)) fuvec(3) = COS(dble(ftheta)) buvec(1) = SIN(dble(btheta)) * COS(dble(bphi)) buvec(2) = SIN(dble(btheta)) * SIN(dble(bphi)) buvec(3) = COS(dble(btheta)) dotpro = fuvec(1) * buvec(1) + fuvec(2) * buvec(2) + fuvec(3) * buvec(3) tvec(1) = fuvec(2) * buvec(3) - fuvec(3) * buvec(2) tvec(2) = fuvec(3) * buvec(1) - fuvec(1) * buvec(3) tvec(3) = fuvec(1) * buvec(2) - fuvec(2) * buvec(1) crossp = SQRT(tvec(1)**2 + tvec(2)**2 + tvec(3)**2) far = radius * ATAN2(crossp, dotpro) IF (far > (50.D0 * zbot)) THEN duthet = 0.0 duphi = 0.0 ELSE dvec(1) = buvec(1) - fuvec(1) dvec(2) = buvec(2) - fuvec(2) dvec(3) = buvec(3) - fuvec(3) utheta(1) = COS(dble(ftheta)) * COS(dble(fphi)) utheta(2) = COS(dble(ftheta)) * SIN(dble(fphi)) utheta(3) = -SIN(dble(ftheta)) dot1 = dvec(1) * utheta(1) + dvec(2) * utheta(2) + dvec(3) * utheta(3) uphi(1) = -SIN(dble(fphi)) uphi(2) = COS(dble(fphi)) uphi(3) = 0.0D0 dot2 = dvec(1) * uphi(1) + dvec(2) * uphi(2) + dvec(3) * uphi(3) alpha = ATAN2(dot2, dot1) z(1) = far * COS(alpha - argume) z(2) = -far * SIN(alpha - argume) ! Note: Z(3) = 0 because we believe it is more important ! to convey that the benchmark is on the free surface ! than it is to convey the precise angular relationship ! to the dislocation. However, some day one might test ! a positive Z(3), in proportion to square of FAR. z(3) = 0.0 vert = ABS(dipf - 1.5708) <= wedge IF (vert) THEN tazimf = argume theta = dipf ! Z values computed above are still good ucap(1) = slip(1) ucap(2) = 0.0 ucap(3) = 0.0 dtop = ztop dbot = zbot ! - - - - - - - - - - - - - - - - - - - - - - - CALL Halo (input, ucap, z, lf, dbot, dtop, & & output, u) ! - - - - - - - - - - - - - - - - - - - - - - - ELSE ! non-vertical, dipping fault: IF (dipf <= 1.5708) THEN tazimf = argume theta = dipf DO 10 i = 1, 3 x(i) = z(i) ucap(i) = slip(i) 10 CONTINUE ELSE ! Look at this fault from the other end/side, ! so that dip THETA will be less than Pi/2: tazimf = argume + 3.141593 theta = 3.141593 - dipf x(1) = -z(1) x(2) = -z(2) x(3) = z(3) ucap(1) = slip(1) ucap(2) = slip(2) ucap(3) = -slip(3) ! Note: Reversal of Z1, Z2 axes (but not Z3) ! during name-change to X1,X2,X3 axes (for AURA) ! combined with change of "moving" block ! leaves "sinistral" and "opening" components ! unchanged, while "relative vertical" component ! reverses. END IF dtop = ztop / SIN(theta) dbot = zbot / SIN(theta) ! Kludge to insure that DTOP is not zero ! because AURA gives erroneous results in that case: dtop = MAX(dtop, 0.001D0 * dbot) ! - - - - - - - - - - - - - - - - - - - - - - - CALL Aura (input, ucap, theta, x, lf, dbot, dtop, & & output, u) ! - - - - - - - - - - - - - - - - - - - - - - - END IF sinaz = SIN(tazimf) cosaz = COS(tazimf) duthet = + cosaz * u(1) + sinaz * u(2) duphi = + sinaz * u(1) - cosaz * u(2) END IF RETURN END SUBROUTINE Coseis (input, farg, fdip, geothe, geophi, & & mxnode, mxfel, numgeo, nfl, nodef, radius, & & v, wedge, xnode, ynode, zmnode, ztranf, & & output, utheta, uphi) ! Computes the horizontal-plane velocities ! (UTHETA,UPHI) for each of NUMGEO benchmark sites, due to ! earthquake dislocations on the brittle parts of ! finite-element faults within the model domain of SHELLS. ! (The usual reason for doing this is so the coseismic part ! of the benchmark velocities can be subtracted from the ! velocities predicted by SHELLS before comparing to geodetic ! data in an aseismic year.) ! GEOTHE and GEOPHI should be the theta and phi coordinates of ! the benchmarks, in radians. Theta is colatitude, measured ! Southward from the North pole. Phi is latitude, measured ! Eastward from the prime meridian. The same convention ! should be used for node positions in XNODE (= theta) and ! YNODE (= phi); the use of X for theta and Y for phi is ! a memnonic convention adopted when flat-planet pLatES ! was rewritten as spherical-planet SHELLS. ! Similarly, the 1st (X, or theta) component of V should be ! the Southward velocity, while the 2nd (Y, or phi) component ! of V should be the Eastward velocity at the nodes. ! Most of the actual work is done by the CALL CHANGE ! statement; this routine is only reponsible for cutting ! each fault into several segments (for greater accuracy) ! and summing the corrections due to each segment of ! each fault. ! The following parameter NSEGME determines how many segments ! each fault element is divided into; higher values may ! give more accuracy for benchmarks close to faults with ! spatially-varying slip rates. PARAMETER (nsegme = 3) DOUBLE PRECISION lf, v(2, mxnode) INTEGER nodef(4, nfl) REAL Chord, fhival, Fltlen, output, & & radius, s, smid, vx1, vx2, vx3, vx4, & & vy1, vy2, vy3, vy4, wedge, x1, x2, x3, x4, & & xbegin, xend, y1, y2, y3, y4, ybegin, yend REAL farg(2, mxfel), fdip(2, nfl), & & geothe(numgeo), geophi(numgeo), moho, & & utheta(numgeo), uphi(numgeo), & & xnode(mxnode), ynode(mxnode), zmnode(mxnode), ztranf(2, nfl) REAL slip(3) !--------------------------------------------------------------------- ! STATEMENT FUNCTIONS: ! Interpolation within one fault element: fhival(s, x1, x2) = x1 + s * (x2 - x1) !--------------------------------------------------------------------- ! Kludge: Convert parameter to variable: nseg = nsegme ! Initialize to zero before beginning sums: DO 10 j = 1, numgeo utheta(j) = 0.0 uphi(j) = 0.0 10 CONTINUE ! Sum contributions from fault elements: DO 100 i = 1, nfl n1 = nodef(1, i) n2 = nodef(2, i) n3 = nodef(3, i) n4 = nodef(4, i) x1 = xnode(n1) x2 = xnode(n2) x3 = xnode(n3) x4 = xnode(n4) y1 = ynode(n1) y2 = ynode(n2) y3 = ynode(n3) y4 = ynode(n4) vx1 = v(1, n1) vx2 = v(1, n2) vx3 = v(1, n3) vx4 = v(1, n4) vy1 = v(2, n1) vy2 = v(2, n2) vy3 = v(2, n3) vy4 = v(2, n4) xbegin = x1 ybegin = y1 ! Sum over NSEG segments within each fault: DO 90 m = 1, nseg smid = (REAL(m) - 0.5) / REAL(nseg) CALL OnArc(input, smid, x1, y1, x2, y2, & & output, xmid, ymid) s = REAL(m) / REAL(nseg) IF(m == nseg) THEN xend = x2 yend = y2 ELSE CALL OnArc(input, s, x1, y1, x2, y2, & & output, xend, yend) END IF al = Fltlen(ybegin, yend, radius, xbegin, xend) lf = 0.5D0 * al ! FARG is the direction from node 1 to node 2, ! in radians counterclockwise from the ! +theta (South) direction, and so is ARGUME: argume = Chord(farg(1, i), smid * 1.0D0, farg(2, i)) dip = fhival(smid, fdip(1, i), fdip(2, i)) ! We consider that the node-1,2 side of the fault ! moves while the node-3,4 side is fixed: delvx = fhival(smid, vx1, vx2) - fhival(smid, vx4, vx3) delvy = fhival(smid, vy1, vy2) - fhival(smid, vy4, vy3) ! UNIT is a horizontal unit vector pointing ! along the fault from the node-1 end to the node-2 end; ! this will be the same as the Z1 direction of HALO if ! the fault is vertical, and equal or opposite to the ! X1 direction of AURA if the fault is dipping: unitx = COS(argume) unity = SIN(argume) ! CROSS is a horizontal unit vector pointing ! across the fault, 90 degrees clockwise from UNIT ! (when viewed from outside the planet); ! this will be the same as the Z2 direction of HALO ! if the fault is exactly vertical or ! equal or opposite to the X2 direction of AURA ! if the fault is dipping. crossx = + unity crossy = -unitx slip(1) = delvx * unitx + delvy * unity ! Note: positive for sinistral slip slip(2) = delvx * crossx + delvy * crossy ! Note: positive for opening across trace ! Vertical component is positive when down: IF (ABS(1.570796 - dip) > wedge) THEN slip(3) = slip(2) * TAN(dip) ELSE slip(3) = 0.0 END IF ! Consider both crustal and mantle patches: moho = fhival(smid, zmnode(n1), zmnode(n2)) DO 80 ilayer = 1, 2 height = ztranf(ilayer, i) IF (height > 0.0) THEN IF (ilayer == 1) THEN ztop = 0.0 ELSE ztop = moho END IF zbot = ztop + height ! Compute effects of patch at all benchmarks: DO 70 j = 1, numgeo ! ------------------------------------------- CALL Change(input, argume, & & geothe(j), geophi(j), & & dip, lf, & & xmid, ymid, & & radius, & & slip, & & wedge, & & ztop, zbot, & & output, duthet, duphi) ! ------------------------------------------- utheta(j) = utheta(j) + duthet uphi(j) = uphi(j) + duphi 70 CONTINUE END IF 80 CONTINUE ! Prepare to loop to next segment: xbegin = xend ybegin = yend 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE Cross (input, a, b, & & output, c) ! COMPUTES VECTOR CROSS PRODUCT C = A x B DIMENSION a(3), b(3), c(3) c(1) = a(2) * b(3) - a(3) * b(2) c(2) = a(3) * b(1) - a(1) * b(3) c(3) = a(1) * b(2) - a(2) * b(1) RETURN END SUBROUTINE Deriv3 (input, iUnitT, mxel, mxnode, & & nodes, numel, & & radius, xnode, ynode, & & output, area, detj, & & dxs, dys, dxsp, dysp, fpsfer, sita) ! 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 TRIANGLES. 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) COMMON / s1s2s3 / points 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) = radius * radius * (0.5 * areap) 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 ARE SAME FOR ! CALCULATING DEERIVATIVE 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.0 sncsne = 0.0 cosm = 0.0 DO 200 j = 1, 3 snccse = snccse + points(j, m) * SIN(theta(j)) * COS(phi(j)) sncsne = sncsne + points(j, m) * SIN(theta(j)) * SIN(phi(j)) cosm = cosm + points(j, m) * COS(theta(j)) 200 CONTINUE xyzp = SQRT(snccse * snccse + sncsne * sncsne + cosm * cosm) snccse = snccse / xyzp sncsne = sncsne / xyzp cosm = cosm / xyzp sitaj = ACOS(cosm) ty = sncsne tx = snccse phaij = Atan2f(ty, tx) csccse = COS(sitaj) * COS(phaij) cscsne = COS(sitaj) * SIN(phaij) ! 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.0.OR.sita(m, i) >= 3.141592654) THEN sitami = sita(m, i) * 57.29577951 WRITE(iUnitT, 220) m, i, sitami 220 FORMAT(' COLATITUDE 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 RETURN END SUBROUTINE EDot (input, dxs, dys, & & fpsfer, mxel, & & mxnode, nodes, numel, radius, sita, v, & & output, erate) ! COMPUTE STRAIN-RATE COMPONENTS EDotXX, EDotYY, AND ! EDotXY (TENSOR FORM; EQUAL TO ! (1/2) * ((DVX/DY)+(DVY/DX)) ! AT THE INTEGRATION POINTS OF TRIANGULAR CONTINUUM ELEMENTS. DOUBLE PRECISION points, v DIMENSION dxs(2, 2, 3, 7, mxel), dys(2, 2, 3, 7, mxel), & & erate(3, 7, mxel), & & fpsfer(2, 2, 3, 7, mxel), & & nodes(3, mxel), points(3, 7), & & sita(7, mxel), v(2, mxnode) COMMON / s1s2s3 / points DO 1000 m = 1, 7 DO 900 i = 1, numel exx = 0. eyy = 0. exy = 0. DO 800 j = 1, 3 node = nodes(j, i) vx = v(1, node) vy = v(2, node) dy11 = dys(1, 1, j, m, i) / SIN(sita(m, i)) dy21 = dys(2, 1, j, m, i) / SIN(sita(m, i)) dy12 = dys(1, 2, j, m, i) / SIN(sita(m, i)) dy22 = dys(2, 2, j, m, i) / SIN(sita(m, i)) fp11 = fpsfer(1, 1, j, m, i) / TAN(sita(m, i)) fp21 = fpsfer(2, 1, j, m, i) / TAN(sita(m, i)) fp12 = fpsfer(1, 2, j, m, i) / TAN(sita(m, i)) fp22 = fpsfer(2, 2, j, m, i) / TAN(sita(m, i)) exx = exx + vx * dxs(1, 1, j, m, i) + vy * dxs(2, 1, j, m, i) eyy = eyy + vx * dy12 + vy * dy22 + vx * fp11 + vy * fp21 exy = exy + vx * dy11 + vy * dy21 & & + vx * dxs(1, 2, j, m, i) + vy * dxs(2, 2, j, m, i) & & - vx * fp12 - vy * fp22 800 CONTINUE erate(1, m, i) = exx / radius erate(2, m, i) = eyy / radius erate(3, m, i) = 0.5 * exy / radius 900 CONTINUE 1000 CONTINUE RETURN END SUBROUTINE Fangls (input, phi, theta, & & output, fangle) ! CALCULATE THE ARGUMENTS (ANGLES COUNTERCLOCKWISE FROM +THETA) ! AT BOTH ENDS OF AN ARC OF A GREAT CIRCLE. DOUBLE PRECISION fpoint COMMON / sfault / fpoint DIMENSION fangle(2), fpoint(7), phi(2), theta(2) dg180 = 3.141592654 a1 = SIN(theta(1)) * COS(phi(1)) a2 = SIN(theta(1)) * SIN(phi(1)) a3 = COS(theta(1)) b1 = SIN(theta(2)) * COS(phi(2)) b2 = SIN(theta(2)) * SIN(phi(2)) b3 = COS(theta(2)) s = 0.99 xx = s * a1 + (1.0 - s) * b1 yy = s * a2 + (1.0 - s) * b2 zz = s * a3 + (1.0 - s) * b3 xval = SQRT(xx * xx + yy * yy + zz * zz) xx = xx / xval yy = yy / xval zz = zz / xval dx = xx - a1 dy = yy - a2 dz = zz - a3 sita = theta(1) phai = phi(1) s1 = COS(sita) * COS(phai) s2 = COS(sita) * SIN(phai) s3 = -SIN(sita) p1 = -SIN(phai) p2 = COS(phai) dxx = dx * s1 + dy * s2 + dz * s3 dyy = dx * p1 + dy * p2 fangle(1) = Atan2f(dyy, dxx) s = 0.01 xx = s * a1 + (1.0 - s) * b1 yy = s * a2 + (1.0 - s) * b2 zz = s * a3 + (1.0 - s) * b3 xval = SQRT(xx * xx + yy * yy + zz * zz) xx = xx / xval yy = yy / xval zz = zz / xval dx = b1 - xx dy = b2 - yy dz = b3 - zz sita = ACOS(zz) phai = Atan2f(yy, xx) IF(phai < 0.0) phai = 2.0 * dg180 + phai s1 = COS(sita) * COS(phai) s2 = COS(sita) * SIN(phai) s3 = -SIN(sita) p1 = -SIN(phai) p2 = COS(phai) dxx = dx * s1 + dy * s2 + dz * s3 dyy = dx * p1 + dy * p2 fangle(2) = Atan2f(dyy, dxx) RETURN END REAL FUNCTION Fltlen (phi1, phi2, radius, theta1, theta2) ! CALCULATES LENGTH OF GREAT CIRCLE SEGMENT BETWEEN ! POINT (THETA1,PHI1) AND POINT (THETA2,PHI2), ! IN PHYSICAL LENGTH UNITS (RADIANS*PLANET RADIUS). DOUBLE PRECISION ab, sqrads, thmean LOGICAL plane thmean = 0.5D0 * theta1 + 0.5D0 * theta2 sqrads = ((phi2 - phi1) * SIN(thmean))**2 + (theta2 - theta1)**2 ! Use plane-geometry method for short distances (under 1 degree): plane = (sqrads < 3.D-4) IF (plane) THEN ab = SQRT(sqrads) ELSE ab = SIN(theta1) * SIN(theta2) * COS(phi1) * COS(phi2) + & & SIN(theta1) * SIN(theta2) * SIN(phi1) * SIN(phi2) + & & COS(theta1) * COS(theta2) ab = ACOS(ab) END IF Fltlen = ab * radius RETURN END SUBROUTINE Findit (input, mxel, mxnode, nodes, numel, xnode, ynode, & & theta, phi, & & output, iele, s1, s2, s3) ! DETERMINE WHICH ELEMENT (IELE) CONTAINS THE SURFACE POINT ! (THETA,PHI), AND ALSO ITS INTERNAL TRIANGULAR COORDINATES ! (S1,S2,S3) IN THE PLANE TRIANGLE: S1+S2+S3=1. ! IF THE POINT IS NOT IN ANY ELEMENT, IELE=0 IS RETURNED. DIMENSION nodes(3, mxel), xnode(mxnode), ynode(mxnode) DIMENSION perp(3), pole(3), ra(3), rb(3), rn(3, 3), rp(3), rs(3) DATA eps / 1.E-4 / ! CARTESIAN (X,Y,Z) RADIUS VECTOR TO POINT ON SPHERE rs(1) = SIN(theta) * COS(phi) rs(2) = SIN(theta) * SIN(phi) rs(3) = COS(theta) DO 100 i = 1, numel DO 10 k = 1, 3 ! FIND CARTESIAN RADIUS VECTOR FOR EACH NODE rn(1, k) = SIN(xnode(nodes(k, i))) * COS(ynode(nodes(k, i))) rn(2, k) = SIN(xnode(nodes(k, i))) * SIN(ynode(nodes(k, i))) rn(3, k) = COS(xnode(nodes(k, i))) 10 CONTINUE ! FIND NORMAL VECTOR OF PLANE ELEMENT ra(1) = rn(1, 2) - rn(1, 1) ra(2) = rn(2, 2) - rn(2, 1) ra(3) = rn(3, 2) - rn(3, 1) rb(1) = rn(1, 3) - rn(1, 2) rb(2) = rn(2, 3) - rn(2, 2) rb(3) = rn(3, 3) - rn(3, 2) CALL Cross (input, ra, rb, output, perp) CALL Unit (modify, perp) ! SHORTEN/LENGTHEN RADIUS VECTOR SO IT POINTS TO PLANE ELEMENT dot = rs(1) * perp(1) + rs(2) * perp(2) + rs(3) * perp(3) IF (dot <= 0.0) GO TO 100 dotp = rn(1, 1) * perp(1) + rn(2, 1) * perp(2) + rn(3, 1) * perp(3) growth = dotp / dot rp(1) = rs(1) * growth rp(2) = rs(2) * growth rp(3) = rs(3) * growth ! COMPUTE AND TEST S1, THEN S2, THEN S3... CALL Cross (input, rn(1, 2), rn(1, 3), output, pole) s1num = rp(1) * pole(1) + rp(2) * pole(2) + rp(3) * pole(3) s1den = rn(1, 1) * pole(1) + rn(2, 1) * pole(2) + rn(3, 1) * pole(3) IF (s1den /= 0.0) THEN s1 = s1num / s1den ELSE GO TO 100 END IF IF ((s1 >= -eps).AND.(s1 <= (1.00 + eps))) THEN CALL Cross (input, rn(1, 3), rn(1, 1), output, pole) s2num = rp(1) * pole(1) + rp(2) * pole(2) + rp(3) * pole(3) s2den = rn(1, 2) * pole(1) + rn(2, 2) * pole(2) + rn(3, 2) * pole(3) IF (s2den /= 0.0) THEN s2 = s2num / s2den ELSE GO TO 100 END IF IF ((s2 >= -eps).AND.(s2 <= (1.00 + eps))) THEN s3 = 1.00 - s1 - s2 IF ((s3 >= -eps).AND.(s3 <= (1.00 + eps))) THEN iele = i GO TO 101 END IF END IF END IF 100 CONTINUE iele = 0 101 RETURN END SUBROUTINE Fnodal (input, iUnitT, mxfel, mxnode, nfl, nodef, & & xnode, ynode, & & output, fpflt) ! CALCULATES VECTOR NODAL FUNCTIONS AT INTERPOLATION POINT ON A ! GREAT CIRCLE FAULT ELEMENT 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 ! SNCCOP: SIN(SITA)*COS(PHAI) AT INTERPOLATION POINT ! SNCSNP SIN(SITA)*SIN(PHAI) AT INTERPOLATION POINT DO 900 i = 1, nfl n1 = nodef(1, i) n2 = nodef(2, i) theta(1) = xnode(n1) theta(2) = xnode(n2) phi(1) = ynode(n1) phi(2) = ynode(n2) CALL Snodal (input, phi, theta, & & output, fpp) DO 800 m = 1, 7 DO 500 j = 1, 2 DO 400 k = 1, 2 DO 300 l = 1, 2 fpflt(l, k, j, m, i) = fpp(l, k, j, m) 300 CONTINUE 400 CONTINUE 500 CONTINUE snccop = fphi(1, m) * SIN(theta(1)) * COS(phi(1)) + & & fphi(2, m) * SIN(theta(2)) * COS(phi(2)) sncsnp = fphi(1, m) * SIN(theta(1)) * SIN(phi(1)) + & & fphi(2, m) * SIN(theta(2)) * SIN(phi(2)) cosm = fphi(1, m) * COS(theta(1)) + fphi(2, m) * COS(theta(2)) pp = SQRT(snccop * snccop + sncsnp * sncsnp + cosm * cosm) cosm = cosm / pp sita = ACOS(cosm) IF(sita <= 0.0.OR.sita >= 3.141592654) THEN sita = sita * 57.29577951 WRITE(iUnitT, 220) m, i, sita 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5, & & ' OF FAULT ELEMENT', & & I5,' IS OUT RANGE', & & E14.4) END IF 800 CONTINUE 900 CONTINUE RETURN END SUBROUTINE GetGEO (input, iUnitB, iUnitT, iUnitY, mxgeo, pltGEO, & & output, geotag, geophi, geothe, & & geovel, geosig, geoazi, & & gpsfmt, numgeo) ! READS RELATIVE HORIZONTAL VELOCITY VECTORS OF BENCHMARKS, ! DETERMINED BY GEODESY. ! The input file format expected is: ! Peter Bird's .gps format of 2002.08.01: ! 3 lines of text/format/headers, ! then one line per benchmark, giving ! coordinates, E & N velocity components in mm/a, ! E & N standard deviations of velocity components in mm/a, ! correlation between E and N components of velocity, ! reference frame, and benchmark identification. ! THIS SUBPROGRAM WILL RETURN: ! GEOPHI = EAST LONGITUDE (RADIANS) ! GEOTHE = CO-LATITUDE (RADIANS) ! GEOVEL = VELOCITY MAGNITUDE, M/S ! GEOAZI = AZIMUTH OF VELOCITY, IN RADIANS COUNTERCLOCKWISE ! FROM SOUTH. ! GEOSIG = STANDARD DEVIATION OF VELOCITY (ASSUMES CIRCULAR ! UNCERTAINTY ELLIPSE, AND NO CORRELATION OF ERR0RS) ! GEOTAG = BENCHMARK NAME CHARACTER*15 frame CHARACTER*20 geotag, tag CHARACTER*132 gpsfmt, header LOGICAL pltGEO REAL nlatdg DIMENSION geotag(mxgeo), geophi(mxgeo), geothe(mxgeo), & & geovel(mxgeo), geosig(mxgeo), geoazi(mxgeo) ! IF FIRST READ FAILS (NO FILE CONNECTED), RETURNS NUMGEO=0 ! "piO180" IS PI/180. (CONVERSION FROM DEGREES TO RADIANS): DATA piO180 / 0.017453293 / ! "SECPYR" IS THE NUMBER OF SECONDS IN A YEAR: DATA secpyr / 31557600. / WRITE(iUnitT, 1) iUnitB 1 FORMAT(/' ATTEMPTING TO READ GEODETIC BENCHMARK VELOCITIES' & & ,' (.gps file) FROM UNIT',I3/) numgeo = 0 IF (pltGEO) OPEN(UNIT = iUnitY, FILE = "ERRORS.gps") ! read, but do not use, the line of text in line 1 READ (iUnitB, * , iostat = ios) IF (ios /= 0) RETURN IF (pltGEO) WRITE(iUnitY, 5) 5 FORMAT("ERROR in geodetic velocity predictions") ! read and store the FORMAT in line 2 READ (iUnitB, "(A)", iostat = ios) gpsfmt IF (pltGEO) WRITE(iUnitY, "(A)") trim(gpsfmt) ! read, but do not use, the column headers in line 3 READ (iUnitB, "(A)", iostat = ios) header IF (ios /= 0) RETURN IF (pltGEO) WRITE(iUnitY, "(A)") trim(header) 10 READ (iUnitB, gpsfmt, END = 100) elondg, nlatdg, vemmpa, vnmmpa, & & vesigm, vnsigm, correl, frame, tag numgeo = numgeo + 1 geophi(numgeo) = elondg * piO180 geothe(numgeo) = (90. - nlatdg) * piO180 geovel(numgeo) = SQRT(vemmpa**2 + vnmmpa**2) / (1000. * secpyr) geoazi(numgeo) = atan2f(vemmpa, -vnmmpa) geosig(numgeo) = SQRT(vesigm**2 + vnsigm**2) / (1000. * secpyr) IF (geosig(numgeo) <= 0.0) THEN WRITE (iUnitT, 19) numgeo 19 FORMAT (/' ERROR: Standard deviation of velocity' & & /' is <= 0.0 for benchmark ',I5 & & /' All values must be positive.') CALL PAUSE () STOP END IF geotag(numgeo) = trim(tag) !CCCC GPBhere IF (numgeo < mxgeo) GO TO 10 WRITE (iUnitT, 90) mxgeo 90 FORMAT (/' SUBPROGRAM GetGEO READ ONLY',I6,' BENCHMARKS.' & & /' INCREASE PARAMETER MAXGEO AND RECOMPILE.') CALL PAUSE() STOP 100 RETURN END SUBROUTINE GetMOR (input, iUnitM, iUnitT, mxmor, & & output, mortag, morphi, morthe, & & morvel, morsig, nummor) ! READS TOTAL SEAFLOOR SPREADING RATES AT MID-OCEAN RISES. ! CURRENT FORMAT MATCHES FILE "SPREADIN.NUVEL1". ! MORTAG = BENCHMARK NAME ! MORPHI = EAST LONGITUDE (RADIANS) ! MORTHE = CO-LATITUDE (RADIANS) ! MORVEL = VELOCITY MAGNITUDE, M/S ! MORSIG = STANDARD DEVIATION OF VELOCITY ! IF FIRST READ FAILS (NO FILE CONNECTED), RETURNS NUMMOR=0 CHARACTER*5 mortag, tag REAL morphi, morthe, morvel, morsig, nlat DIMENSION mortag(mxmor), morphi(mxmor), morthe(mxmor), & & morvel(mxmor), morsig(mxmor) ! "piO180" IS PI/180. (CONVERSION FROM DEGREES TO RADIANS): DATA piO180 / 0.017453293 / ! "SECPYR" IS THE NUMBER OF SECONDS IN A YEAR: DATA secpyr / 31557600. / WRITE(iUnitT, 1) iUnitM 1 FORMAT(/' ATTEMPTING TO READ SEAFLOOR SPREADING RATES FROM' & & ,' UNIT',I3/) nummor = 0 READ (iUnitM, * , iostat = ios) IF (ios /= 0) RETURN READ (iUnitM, * ) READ (iUnitM, * ) 10 READ (iUnitM, 11, END = 100) tag, nlat, elon, v, sigma 11 FORMAT (A5,F7.2,F8.2,F6.1,F5.1) nummor = nummor + 1 mortag(nummor) = tag morphi(nummor) = elon * piO180 morthe(nummor) = (90. - nlat) * piO180 morvel(nummor) = v / (1000. * secpyr) morsig(nummor) = sigma / (secpyr * 1000.) IF (nummor < mxmor) GO TO 10 WRITE (iUnitT, 90) mxmor 90 FORMAT (/' SUBPROGRAM GetMOR READ ONLY',I6,' BENCHMARKS.' & & /' INCREASE PARAMETER MAXMOR AND RECOMPILE.') CALL PAUSE() STOP 100 RETURN END SUBROUTINE GetNet (input, iunit7, iunit8, & & mxel, mxfel, mxnode, & & output, brief, dqdtda, elev, fdip, & & nfaken, nfl, nodef, nodes, nrealn, & & numel, numnod, n1000, offmax, offset, & & title1, tlnode, xnode, ynode, zmnode, & & work, checke, checkf, checkn) ! READ FINITE ELEMENT GRID FROM UNIT "IUNIT7". ! ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8". CHARACTER*80 title1 LOGICAL allok, brief ! NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL checke, checkf, checkn DIMENSION checke(mxel), checkf(mxfel), checkn(mxnode), & & dqdtda(mxnode), elev(mxnode), & & fdip(2, mxfel), nodef(4, mxfel), & & nodes(3, mxel), offset(mxfel), tlnode(mxnode), & & xnode(mxnode), ynode(mxnode), zmnode(mxnode) DIMENSION dips(3), vecton(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) title1 3 FORMAT(/' TITLE OF FINITE ELEMENT GRID ='/' ',A80) ! READ NUMBER OF NODES. ! INPUT NODAL LOCATIONS (X,Y), ELEVATIONS, HEAT-FLOW, AND ISOSTATIC ! GRAVITY ANOMALIES (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',I5,' NODES IN THE GRID') WRITE (iunit8, 50) 50 FORMAT (/ & & 77X,' MANTLE'/ & & 77X,' CRUSTAL LITHOSPHERE'/ & & ' NODE E-LONGITUDE N-LATITUDE', & & ' THETA PHI ELEVATION', & & ' HEAT-FLOW THICKNESS THICKNESS'/) END IF DO 90 k = 1, numnod checkn(k) = .FALSE. 90 CONTINUE DO 100 k = 1, numnod CALL ReadN (input, iunit7, iunit8, 7, & & output, vecton) INDEX = vecton(1) + 0.5 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 = vecton(2) pLat = vecton(3) xi = (90.0 - pLat) * 0.017453292 yi = pLon * 0.017453292 elevi = vecton(4) qi = vecton(5) zmi = vecton(6) tli = vecton(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.) THEN WRITE (iunit8, 96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL PAUSE() STOP END IF IF (zmi < 0.) THEN WRITE (iunit8, 97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') CALL PAUSE() STOP END IF zmnode(i) = zmi IF (tli < 0.) 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,I6) 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 (',I6,') 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 ',I6,' TRIANGULAR CONTINUUM ELEMENTS.'/ & & ' (NODE NUMBERS FOR EACH ARE GIVEN CORNERS, COUNTER', & & 'CLOCKWISE'/ / & & ' ELEMENT C1 C2 C3') DO 200 k = 1, numel ! (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) READ (iunit7, * ) i, (nodes(j, i), j = 1, 3) IF ((i < 1).OR.(i > numel)) THEN WRITE (iunit8, 111) i 111 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I6) CALL PAUSE() STOP END IF checke(i) = .TRUE. IF (.NOT.brief) WRITE (iunit8, 120) i, (nodes(j, i), j = 1, 3) 120 FORMAT (' ',I6,':',3I10) 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 ',I6,' 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,I6) 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 (',I6,') AND RECOMPILE.') CALL PAUSE() STOP END IF offmax = 0. DO 222 i = 1, nfl checkf(i) = .FALSE. 222 CONTINUE IF (.NOT.brief) WRITE(iunit8, 230) nfl 230 FORMAT(/ /' THERE ARE ',I6,' GREAT CIRCLE FAULT ELEMENTS.') IF ((.NOT.brief).AND.(nfl > 0)) WRITE(iunit8, 231) 231 FORMAT (/' (THE 4 NODE NUMBERS DEFINING EACH ELEMENT MUST BE', & & ' IN A COUNTERCLOCKWISE ORDER:'/ & & ' N1, AND N2 ARE IN LEFT-TO-RIGHT SEQUENCE ON THE', & & ' NEAR SIDE,'/ & & ' THEN N3 IS OPPOSITE N2, N4 IS OPPOSITE N1 '/, & & ' (FAULT DIPS ARE GIVEN AT N1, N2, ', & & ' IN DEGREES FROM HORIZONTAL;'/ & & ' POSITIVE DIPS ARE TOWARD N1, AND N2, RESPECTIVELY, '/ & & ' WHILE NEGATIVE DIPS ARE TOWARD N4, AND N3.)'/ & & ' OFFSET IS THE TOTAL PAST SLIP OF THE FAULT.'/ / & & ' ELEMENT N1 N2 N3 N4 DIP1 DIP2', & & ' OFFSET'/) 240 FORMAT (' ',I6,':',4I5,1X,2F6.1,1X,F9.0) DO 300 k = 1, nfl off = 0. READ(iunit7, * ) i, (nodef(j, k), j = 1, 4), (dips(l), l = 1, 2), off IF ((i < 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) WRITE (iunit8, 240) i, (nodef(j, i), j = 1, 4), & & (dips(l), l = 1, 2), off DO 250 j = 1, 4 n = nodef(j, i) IF (n > 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 (',I6, & & ') IN FAULT ',I6) CALL PAUSE() STOP END IF nodef(j, i) = n 250 CONTINUE DO 260 l = 1, 2 IF (ABS(dips(l)) > 90.) THEN WRITE(iunit8, 252) dips(l) 252 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN', & & ' RANGE OF -90. TO +90. DEGREES.'/ & & ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE', & & ' HORIZONAL;'/ & & ' A + PREFIX (OR NONE) INDICATES A DIP', & & ' TOWARD THE N1-N2 SIDE;'/ & & ' A - PREFIX INDICATES A DIP TOWARD', & & ' THE N4-N3 SIDE.)') CALL PAUSE() STOP END IF IF (dips(l) < 0.) dips(l) = 180. + dips(l) fdip(l, i) = dips(l) * 0.017453293 260 CONTINUE IF (off < 0.) THEN WRITE (iunit8, 280) off 280 FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2, & & ' FOR FAULT ELEMENT',I6/ & & ' OFFSETS MAY NOT BE NEGATIVE.') CALL PAUSE() STOP END IF offset(i) = off offmax = MAX(offmax, off) 300 CONTINUE allok = .TRUE. DO 301 i = 1, nfl allok = allok.AND.checkf(i) 301 CONTINUE IF (.NOT.allok) THEN WRITE(iunit8, 302) 302 FORMAT(' THE FOLLOWING FAULTS WERE NEVER READ:') DO 304 i = 1, nfl IF (.NOT.checkf(i)) WRITE(iunit8, 303)i 303 FORMAT(' ',36X,I6) 304 CONTINUE CALL PAUSE() STOP ELSE IF (offmax > 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 (' --------------------------------------------------', & & '-----------------------------') RETURN END SUBROUTINE GetPBx (iUnitM, iUnitT, names, nPBnd, nPlate, & ! INTENT(IN) & nBoundaryPoints, pLat, pLon) ! INTENT(OUT) ! Sets up arrays defining the pLates in the PB2002 model of: ! Bird [2003; G**3]. ! (The rotation vectors of the pLates are contained in DATA ! statements in the main PROGRAM.) ! The digitized boundaries of the pLates (continuous closed curves, ! always circling counterclockwise, and redundantly descriging ! each plate boundary *twice*, from each side) ! are read here, from an input file such as 'PB2002_plates.dig', ! on Fortran input device iUnitM. ! The convention for identifying the pLates is a 2-character symbol. ! See array NAMES in the main PROGRAM. !------------------------------------------------------- CHARACTER*2 names, symbol CHARACTER*3 stars DIMENSION names (nPlate) , nBoundaryPoints(nPlate) DIMENSION pLat (nPlate, nPBnd), pLon (nPlate, nPBnd) !------------------------------------------------------ nread = 0 100 READ (iUnitM, 101, END = 201, iostat = ios) symbol IF((nread == 0).AND.(ios /= 0)) THEN WRITE(*, "(' ERR','OR:File not found, or file empty.')") CALL PAUSE() STOP END IF 101 FORMAT (A2) DO 120 l = 1, nPlate IF(symbol == names(l)) THEN ip = l GO TO 140 END IF 120 CONTINUE WRITE (iUnitT, 121) iUnitM 121 FORMAT (/' ERR0R: BAD plate NAME ON INPUT DEVICE ',I3) CALL PAUSE() STOP 140 nread = nread + 1 IF (nread > nPlate) THEN WRITE (iUnitT, 141) WRITE (*, 141) 141 FORMAT(/' Increase nPlate and recompile.') CALL PAUSE() STOP END IF i = 0 142 READ (iUnitM, 145, END = 201) stars 145 FORMAT (A3) IF (stars == '***') THEN nBoundaryPoints(ip) = i GO TO 100 END IF BACKSPACE iUnitM i = i + 1 IF (i > nPBnd) THEN WRITE (iUnitT, 146) WRITE (*, 146) 146 FORMAT(/' Increase nPBnd and recompile.') CALL PAUSE() STOP END IF READ (iUnitM, * ) pLon(ip, i), pLat(ip, i) pLon(ip, i) = pLon(ip, i) * 0.017453293 pLat(ip, i) = pLat(ip, i) * 0.017453293 GO TO 142 201 IF(nread < nPlate) THEN WRITE(iUnitT, "(' ERR','OR: Expecting ',I3,' pLates but' & & /' read outlines of only ',I3)")nPlate, nread WRITE(*, "(' ERR','OR: Expecting ',I3,' pLates but' & & /' read outlines of only ',I3)")nPlate, nread CALL PAUSE() STOP END IF RETURN END SUBROUTINE GetSLF(input, iUnitD, iUnitT, mxgsr, & & output, fname, rlat, delvz, fltlon, fltlat, & & nfdata) CHARACTER*8 c8r1, c8r2, c8v1, c8v2 CHARACTER*80 fname LOGICAL dummy DIMENSION delvz(2, mxgsr), fltlat(mxgsr), fltlon(mxgsr), & & fname(mxgsr), rlat(2, mxgsr) ! DETECT WHETHER THE FILE EXISTS nfdata = 0 WRITE (iUnitT, 101) iUnitD 101 FORMAT (/' ATTEMPTING TO READ FAULT SLIP RATES FROM UNIT',I3/) READ (iUnitD, * , iostat = ios) IF (ios == 0) THEN REWIND(iUnitD) ELSE RETURN END IF DO 106 i = 1, mxgsr fname(nfdata + 1) = ' ' READ(iUnitD, "(A)", END = 111) fname(nfdata + 1) ! NOTE: The line with four slip-rate numbers (in mm/year) ! is marked with initial letter "F". ! This is to prevent program -Projector- ! from interpreting such lines as (x, y) data. ! When we read such lines here, we assign the letter "F" ! to LOGICAL DUMMY (and never use it). READ(iUnitD, * , END = 111) dummy, & & rlat(1, nfdata + 1), rlat(2, nfdata + 1), & & delvz(1, nfdata + 1), delvz(2, nfdata + 1) READ(iUnitD, * , END = 111) fltlon(nfdata + 1), fltlat(nfdata + 1) READ(iUnitD, * , END = 111) nfdata = nfdata + 1 106 CONTINUE WRITE (iUnitT, 107) mxgsr 107 FORMAT (/' ERR0R: PARAMETER MAXDAT = ',I6,' IS NOT LARGE'/ & & ' ENOUGH FOR THE NUMBER OF GEOLOGIC SLIP-RATE DATA.') CALL PAUSE() STOP 111 IF (nfdata == 0) THEN WRITE (iUnitT, 112) iUnitD 112 FORMAT (/' NO FAULT SLIP RATE DATA COULD BE READ FROM UNIT ' & & ,I2) ELSE WRITE (iUnitT, 113) nfdata, iUnitD 113 FORMAT (/' ',I4,' DATA ON FAULT SLIP RATES', & & ' WERE READ FROM UNIT ',I2,':'/) WRITE (iUnitT, 114) 114 FORMAT ( & & ' ', & & ' -in millimeters/year- -in degrees-'/ & & ' ', & & ' Min. Max. Min. Max. East North'/ & & ' Fault ', & & 'Dextral Dextral Throw Throw Longitude Latitude'/ & & ' --------------------------- ', & & '------- ------- ------- ------- --------- --------') DO 150 i = 1, nfdata IF (rlat(1, i) == 0.0) THEN c8r1 = ' ?' ELSE WRITE (c8r1, "(F8.2)") rlat(1, i) END IF IF (rlat(2, i) == 0.0) THEN c8r2 = ' ?' ELSE WRITE (c8r2, "(F8.2)") rlat(2, i) END IF IF (delvz(1, i) == 0.0) THEN c8v1 = ' ?' ELSE WRITE (c8v1, "(F8.2)") delvz(1, i) END IF IF (delvz(2, i) == 0.0) THEN c8v2 = ' ?' ELSE WRITE (c8v2, "(F8.2)") delvz(2, i) END IF WRITE(iUnitT, 145)fname(i)(1:27), & & c8r1, c8r2, & & c8v1, c8v2, & & fltlon(i), fltlat(i) 145 FORMAT(' ',A,4A8,F10.3,F9.3) 150 CONTINUE WRITE (iUnitT, 170) 170 FORMAT (' NOTE: Throw-rates are always positive in', & & ' cases of convergence (thrusting),'/ & & ' and negative in cases of divergence', & & ' (normal faulting).') END IF RETURN END SUBROUTINE GetSKS (iUnitK, iUnitT, mxSKS, & ! input & numSKS, SKS_tag, SKS_theta, SKS_phi, & & SKS_argument, SKS_delay) ! output ! READs data on upper-mantle seismic anisotropy, from SKS splitting. ! Current FORMAT matches file "Fouch_upper-mantle_SKS_splitting-2004.dat": ! If first READ fails (no file connected), returns numSKS = 0 IMPLICIT NONE INTEGER, INTENT(IN) :: iUnitK ! Fortran device number to read from INTEGER, INTENT(IN) :: iUnitT ! Fortran devide number for error messages INTEGER, INTENT(IN) :: mxSKS ! maximum number of SKS data that can be stored INTEGER, INTENT(OUT) :: numSKS ! number of SKS data actually read CHARACTER*5, DIMENSION(mxSKS), INTENT(OUT) :: SKS_tag ! 5-byte ID code for datum REAL, DIMENSION(mxSKS), INTENT(OUT) :: SKS_theta ! colatitude, in radians REAL, DIMENSION(mxSKS), INTENT(OUT) :: SKS_phi ! colatitude, in radians REAL, DIMENSION(mxSKS), INTENT(OUT) :: SKS_argument ! phi or fast direction, in radians counterclockwise from S REAL, DIMENSION(mxSKS), INTENT(OUT) :: SKS_delay ! SKS splitting time, in s CHARACTER*5 tag INTEGER :: azimuth, ios LOGICAL :: problem REAL :: dt, latitude, longitude, piO180 ! "piO180" is pi/180. (conversion from degrees to radians): DATA piO180 / 0.017453293 / WRITE(iUnitT, 1) iUnitK 1 FORMAT(/' Attempting to read fast SKS azimuths' & & ,' and delays from unit',I3/) numSKS = 0 ! Do some test READs to see if a file is connected: problem = .FALSE. READ (iUnitK, "(A)", IOSTAT = ios) tag ! title line? problem = problem.OR.(ios /= 0) READ (iUnitK, "(A)", IOSTAT = ios) tag ! header line? problem = problem.OR.(ios /= 0) IF (problem) RETURN ! Now proceed to read actual datum lines: ! MAIN LOOP 10 READ (iUnitK, 32, END = 100) tag, latitude, longitude, azimuth, dt 32 FORMAT (A5,1X,2F10.4,I8,F8.2) numSKS = numSKS + 1 SKS_tag(numSKS) = tag SKS_phi(numSKS) = longitude * piO180 SKS_theta(numSKS) = (90.0 - latitude) * piO180 SKS_argument(numSKS) = (180.0 - azimuth) * piO180 SKS_delay(numSKS) = dt IF (numSKS < mxSKS) GO TO 10 WRITE (iUnitT, 90) mxSKS 90 FORMAT (/' Subprogram GetSTR READ only',I6,' data points.' & & /' Increase PARAMETER maxSKS and recompile.') CALL Pause() STOP 100 RETURN END SUBROUTINE GetSTR (input, iUnitS, iUnitT, mxSTR, & & output, strTag, strThe, strPhi, & & strArg, strQua, strReg, numSTR) ! READS MOST-COMPRESSIVE HORIZONTAL PRINCIPAL STRESS AZIMUTHS. ! CURRENT FORMAT MATCHES FILE "WSM1092.LRECL116": ! strTag = SHORT 5-BYTE IDENTIFICATION ! strPhi = EAST LONGITUDE (RADIANS) ! strThe = CO-LATITUDE (RADIANS) ! strArg = AZIMUTH OF SIGMA.1, IN RADIANS COUNTERCLOCKWISE ! FROM SOUTH. ! strQua = RELATIVE QUALITY, CONVERTED FROM LETTER TO REAL WEIGHT. ! strReg = 2-LETTER ABBREVIATION FOR STRESS REGIME. ! IF FIRST READ FAILS (NO FILE CONNECTED), RETURNS NUMSTR=0 CHARACTER*5 strTag, tag CHARACTER*2 strReg, regime CHARACTER*1 qualit REAL nlat INTEGER azimut DIMENSION strTag(mxstr), strPhi(mxstr), strThe(mxstr), & & strQua(mxstr), strArg(mxstr), strReg(mxstr) ! "piO180" IS PI/180. (CONVERSION FROM DEGREES TO RADIANS): DATA piO180 / 0.017453293 / WRITE(iUnitT, 1) iUnitS 1 FORMAT(/' ATTEMPTING TO READ STRESS DIRECTIONS' & & ,' FROM UNIT',I3/) numstr = 0 ! DO ONE TEST READ TO SEE IF A FILE IS CONNECTED: READ (iUnitS, 11, iostat = ios) tag, nlat, elon, azimut, qualit, regime IF (ios == 0) THEN ! BACKSPACE 11 ! WRONG UNIT# FOR iUnitS? CHANGED BY Z.L BACKSPACE 12 ELSE RETURN END IF ! MAIN LOOP 10 READ (iUnitS, 11, END = 100) tag, nlat, elon, azimut, qualit, regime 11 FORMAT (A5,2F9.3,7X,I3,17X,A1,5X,A2) IF ((azimut.LT.-180).OR.(azimut.GT.360)) THEN WRITE (*,"(' ERROR: Bad azimuth in stress data: ',I6)") azimut WRITE (*, 12) tag, nlat, elon, azimut, qualit, regime 12 FORMAT (' ',A5,2F9.3,7X,I3,17X,A1,5X,A2) STOP END IF numSTR = numSTR + 1 strTag(numstr) = tag strPhi(numstr) = elon * piO180 strThe(numstr) = (90. - nlat) * piO180 strArg(numstr) = (180. - azimut) * piO180 IF (qualit == 'A') THEN strQua(numstr) = 4. ELSE IF (qualit == 'B') THEN strQua(numstr) = 3. ELSE IF (qualit == 'C') THEN strQua(numstr) = 2. ELSE IF (qualit == 'D') THEN strQua(numstr) = 1. ELSE strQua(numstr) = 0. END IF strReg(numstr) = regime IF (numSTR < mxSTR) GO TO 10 WRITE (iUnitT, 90) mxSTR 90 FORMAT (/' SUBPROGRAM GetSTR READ ONLY',I6,' AZIMUTHS.' & & /' INCREASE PARAMETER MAXSTR AND RECOMPILE.') CALL PAUSE() STOP 100 RETURN END SUBROUTINE Halo (input, ucap, z, l, dbot, dtop, & & output, u) ! COMPUTES DISPLACEMENTS IN A UNIFORM ELASTIC HALF-SPACE ! OF POISSON'S RATION 0.25 DUE TO UNIFORM SLIP ON A ! RECTANGULAR DISLOCATION IN THE VERTICAL PLANE. ! EQUATIONS BASED ON MANSINHA AND SMYLIE (1967), ! J. GEOPHYS. RES., 72, 4731-4743. IMPLICIT DOUBLE PRECISION (a - z) REAL output INTEGER i, input DIMENSION ucap(3), z(3), u(3) ! THE COORDINATE SYSTEM IS CARTESIAN AND RIGHT-HANDED, ! WITH Z3 DOWN, Z1 PARALLEL TO THE FAULT TRACE, AND Z2 ! PERPENDICULAR. THE DISLOCATION SURFACE IS AT Z2=0, ! DTOP <= Z3 <= DBOT (NOTE THAT DTOP > 0.), AND ! -L <= Z1 <= L. ! UCAP CONTAINS THE THREE COMPONENTS OF THE SLIP OF THE ! POSITIVE-Z2 SIDE RELATIVE TO THE OTHER SIDE; UCAP(2) ! SHOULD THEREFORE BE ZERO (BUT IT IS NOT USED HERE). ! U CONTAINS THE THREE COMPONENTS OF THE DISPLACEMENT ! AT THE OBSERVATION POINT Z. DO 10 i = 1, 3 10 u(i) = 0. DO 100 i = 1, 4 factor = 1. IF(i == 2.OR.i == 3) factor = -1. zeta1 = l IF(i == 3.OR.i == 4) zeta1 = -l zeta3 = dbot IF(i == 2.OR.i == 4) zeta3 = dtop r = SQRT((z(1) - zeta1)**2 + z(2)**2 + (z(3) - zeta3)**2) q = SQRT((z(1) - zeta1)**2 + z(2)**2 + (z(3) + zeta3)**2) term1 = z(2) * (z(1) - zeta1) * (2. / (r * (r + z(3) - zeta3)) & & - (5. * q + 8. * zeta3) / (2. * q * (q + z(3) + zeta3)**2) & & + 4. * z(3) * zeta3 * (2. * q + z(3) + zeta3) / & & (q**3 * (q + z(3) + zeta3)**2)) & & + 3. * Atanpv((z(1) - zeta1) * (z(3) - zeta3), (r * z(2))) & & - 3. * Atanpv((z(1) - zeta1) * (z(3) + zeta3), (q * z(2))) term2 = -log(r + z(3) - zeta3) & & + 0.5 * log(q + z(3) + zeta3) & & - (5. * z(3) - 3. * zeta3) / (2. * (q + z(3) + zeta3)) & & - 4. * z(3) * zeta3 / (q * (q + z(3) + zeta3)) & & + z(2)**2 * (2. / (r * (r + z(3) - zeta3)) & & - (5. * q + 8. * zeta3) / (2. * q * (q + z(3) + zeta3)**2) & & + 4. * z(3) * zeta3 * (2. * q + z(3) + zeta3) / & & (q**3 * (q + z(3) + zeta3)**2)) ! TERM3=Z(2)*(2./R-2./Q+4.*Z(3)*ZETA3/Q**3 ! 1 +3./(Q+Z(3)+ZETA3)+2.*(Z(3)+3.*ZETA3)/ ! 2 (Q*(Q+Z(3)+ZETA3))) ! TERM4=Z(2)*(2./R+4./Q-4.*Z(3)*ZETA3/Q**3 ! 1 -6.*Z(3)/(Q*(Q+Z(3)+ZETA3))) ! TERM5= -LOG(R+Z(1)-ZETA1)+LOG(Q+Z(1)-ZETA1) ! 1 +(6.*Z(3)**2+10.*Z(3)*ZETA3)/(Q*(Q+Z(1)-ZETA1)) ! 2 +(6.*Z(3)*(Z(1)-ZETA1))/(Q*(Q+Z(3)+ZETA3)) ! 3 +Z(2)**2*(2./(R*(R+Z(1)-ZETA1)) ! 4 +4./(Q*(Q+Z(1)-ZETA1)) ! 5 -4.*Z(3)*ZETA3*(2.*Q+Z(1)-ZETA1)/(Q**3*(Q+Z(1)-ZETA1)**2)) ! TERM6=Z(2)*(2.*(Z(3)-ZETA3)/(R*(R+Z(1)-ZETA1)) ! 1 -2.*(Z(3)+2.*ZETA3)/(Q*(Q+Z(1)-ZETA1)) ! 2 -4.*Z(3)*ZETA3*(Z(3)+ZETA3)*(2.*Q+Z(1)-ZETA1)/ ! 3 (Q**3*(Q+Z(1)-ZETA1)**2)) ! 4 +3.*ATANPV((Z(1)-ZETA1)*(Z(3)-ZETA3),(R*Z(2))) ! 5 -3.*ATANPV((Z(1)-ZETA1)*(Z(3)+ZETA3),(Q*Z(2))) u(1) = u(1) + factor * term1 * ucap(1) / 37.7 u(2) = u(2) + factor * term2 * ucap(1) / 37.7 ! U(1)=U(1)+FACTOR*TERM4*UCAP(3)/37.7 ! U(2)=U(2)+FACTOR*TERM5*UCAP(3)/37.7 ! U(3)=U(3)+FACTOR*(TERM3*UCAP(1)+TERM6*UCAP(3))/37.7 100 CONTINUE RETURN END SUBROUTINE Interp (input, fatnod, mxel, mxnode, nodes, numel, & & output, fatip) ! INTERPOLATE A SCALAR FUNCTION KNOWN AT THE NODES ("FATNOD") ! TO VALUES AT THE 7 INTEGRATION POINTS IN EACH TRIANGULAR ! CONTINUUM ELEMENT. DOUBLE PRECISION points COMMON / s1s2s3 / points DIMENSION points(3, 7) DIMENSION fatnod(mxnode), fatip(7, mxel), nodes(3, mxel) DO 100 m = 1, 7 DO 90 i = 1, numel fatip(m, i) = points(1, m) * fatnod(nodes(1, i)) + & & points(2, m) * fatnod(nodes(2, i)) + & & points(3, m) * fatnod(nodes(3, i)) 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE Likely (input, vmin, vmax, v, & & output, prob) ! COMPUTE THE PROBABILITY THAT VELOCITY V IS CONSISTENT WITH LIMIT ! VMIN AND VMAX (BOTH OF WHICH CONTAIN RANDOM ERR0R). LOGICAL boxcar, post, point DATA badata / 0.20 / , sdmodv / 0.10 / , sdage / 0.15 / , e / 2.71828 / , & & tpim1h / 0.398942 / prob = 1. IF ((vmin == 0.) .AND. (vmax == 0.)) RETURN v1 = vmin IF (vmin == 0.) v1 = -1.E38 v2 = vmax IF (vmax == 0.) v2 = + 1.E38 vcent = (v1 + v2) / 2. spread = ABS((v2 - v1) / (v1 + v2)) boxcar = (spread >= 1.3 * sdage).OR.(v1 == -1.E38).OR.(v2 == 1.E38) point = spread == 0. post = (.NOT.boxcar).AND.(.NOT.point) vu = v * (1. + 2. * sdmodv) vl = v * (1. - 2. * sdmodv) IF(v >= 0.)GO TO 4 temp = vu vu = vl vl = temp 4 v1u = v1 * (1. + 2. * sdage) v1l = v1 * (1. - 2. * sdage) IF(v1 >= 0.)GO TO 6 temp = v1u v1u = v1l v1l = temp 6 v2u = v2 * (1. + 2. * sdage) v2l = v2 * (1. - 2. * sdage) IF(v2 >= 0.)GO TO 8 temp = v2u v2u = v2l v2l = temp 8 IF(vu > v1l)GO TO 10 prob = badata RETURN 10 IF(vl < v2u) GO TO 20 prob = badata RETURN 20 IF((vl < v1u).OR.(vu > v2l)) GO TO 30 prob = 1. RETURN 30 prob = 0. b = -3.1 db = 0.10 d1 = ABS(v1) * sdage * 1.41421 d2 = ABS(v2) * sdage * 1.41421 dw = db * sdmodv * ABS(v) w = v - 31. * dw DO 40 i = 1, 61 b = b + db pv = tpim1h * e**(-0.5 * b**2) w = w + dw IF(point) pgf = e**(-MIN(99., 0.5 * ((w - v1) / (sdage * v1))**2)) IF(post)pgf = e**(-MIN(99., (0.832 * (w - vcent) / (spread * vcent))**2)) IF(boxcar)pgf = MIN(0.5 * (1. + erf((w - v1) / d1)), & & 1., 0.5 * (1. + erf((v2 - w) / d2))) pgf = AMAX1(pgf, 1.E-50) pg = badata + (1. - badata) * pgf 40 prob = prob + pv * pg * db RETURN END SUBROUTINE Limits (input, area, detj, iUnitT, mxel, numel, & & okdelv, radius, refstr, sphere, & & tlint, trhmax, zmoho, & & output, constr, etamax, fmumax, vismax) ! COMPUTE AREA, MEAN THICKNESS, AND OTHER DIMENSIONAL PARAMETERS ! OF THE plate, THEN DETERMINE VALUES OF STIFFNESS LIMITS NEEDED ! TO KEEP VELOCITY ERR0RS DOWN TO ORDER "OKDELV" AT SHEAR STRESS ! LEVEL "REFSTR". LOGICAL sphere DOUBLE PRECISION weight COMMON / wgtvec / weight DIMENSION weight(7) DIMENSION area(mxel), detj(7, mxel), tlint(7, mxel), zmoho(7, mxel) totala = 0. totalv = 0. DO 20 m = 1, 7 DO 10 i = 1, numel da = area(i) * detj(m, i) * weight(m) totala = totala + da totalv = totalv + da * (zmoho(m, i) + tlint(m, i)) 10 CONTINUE 20 CONTINUE thick = totalv / totala IF (sphere) THEN side = radius nfault = 1 ELSE side = SQRT(totala) nfault = 4 END IF constr = nfault * refstr * thick / okdelv etamax = refstr * thick / (side * okdelv) etamax = MIN(etamax, trhmax / okdelv) fmumax = nfault * refstr / okdelv vismax = 0.25 * refstr * side / okdelv WRITE (iUnitT, 50) totala, totalv, thick, side, constr, etamax, & & fmumax, vismax 50 FORMAT (/ /' SUBPROGRAM -LIMITS- PERFORMS DIMENSIONAL ANALYSIS'/ & & ' AND ESTIMATES NECESSARY STIFFNESS LIMITS TO BALANCE'/1P, & & ' THE CONFLICTING OBJECTIVES OF ACCURACY AND PRECISION:'/ / & & ' AREA OF MODEL = ',E10.3,' LENGTH**2'/ & & ' VOLUME OF MODEL = ',E10.3,' LENGTH**3'/ & & ' TYPICAL THICKNESS = ',E10.3,' LENGTH'/ & & ' TYPICAL WIDTH = ',E10.3,' LENGTH'/ & & ' CONSTR (CONSTRAINT WEIGHT) = ',E10.3,' FORCE-SEC/LENGTH**2'/ & & ' ETAMAX (MAX. BASAL COUPLING) = ',E10.3,' FORCE-SEC/LENGTH**3'/ & & ' FMUMAX (MAX. FAULT STIFFNESS) = ',E10.3,' FORCE-SEC/LENGTH**3'/ & & ' VISMAX (MAX. BLOCK VISCOSITY) = ',E10.3,' FORCE-SEC/LENGTH**2') RETURN END SUBROUTINE Mohr (input, acreep, alphat, bcreep, biot, byerly, & & ccreep, cfric, conduc, constr, dcreep, dqdtda, & & ecreep, elev, fdip, ffric, fmumax, & & fpflt, farg, gmean, & & mxfel, mxnode, nfl, nodef, & & offmax, offset, & & onekm, radio, rhoh2o, rhobar, & & slide, sphere, taumax, & & tlnode, tsurf, v, wedge, & & zmnode, & & modify, ztranf, & & output, fc, fimudz, fpeaks, fslips, ftstar) ! THIS SUBPROGRAM CONTAINS THE NONLINEAR RHEOLOGY OF THE FAULTS. ! FOR EACH OF 7 INTEGRATION POINTS ALONG THE LENGTH OF EACH FAULT ! ELEMENT, IT: ! (1) COMPUTES THE SLIP-RATE VECTOR ON THE FAULT SURFACE, ! (2) DETERMINES THE SHEAR STRESS ON THE FAULT SURFACE BY MOHR/ ! COULOMB/NAVIER THEORY (THIS STRESS IS PROPORTIONAL TO DEPTH, ! SO THE CALCULATION IS ACTUALLY DONE AT UNIT DEPTH AND THEN ! SCALED), ! (3) PROCEEDS DOWN THE DIP OF THE FAULT, CHECKING TEMPERATURE, ! STRAIN-RATE, AND PRESSURE TO SEE IF FRICTIONAL OR CREEP ! SHEAR STRESS IS LOWER, ! (4) REPORTS THE VERTICAL INTEGRAL OF "MU" (THE RATIO OF SHEAR ! STRESS TO SLIP RATE) DOWN THE FAULT AS "FIMUDZ". ! (NOTE THAT THE INTEGRAL IS VERTICAL, NOT ON A SLANT, EVEN THOUGH ! CONDITIONS ARE EVALUATED ALONG A SLANT PATH.) ! (5) FOR DIPPING, OBLIQUE-SLIP FAULT ONLY, ALSO REPORTS RECOMMENDED ! TACTICAL VALUES FOR THE MATRIX "FC" AND THE VECTOR "FTSTAR" ! WHICH JOINTLY DESCRIBE A LINEARIZED RHEOLOGY STIFFER THAN ! THE ACTUAL NONLINEAR RHEOLOGY. ! (6) "ZTRANF" HOLDS THE LATEST ESTIMATES OF THE DEPTHS ! TO THE BRITTLE/DUCTILE TRANSITION, AT THE FAULT MIDPOINT; ! FIRST VALUE REFERS TO CRUST, AND SECOND TO MANTLE LITHOSPHERE. ! EACH VALUE IS MEASURED FROM THE TOP OF THE LAYER. ! (7) LOGICAL VARIABLE "FSLIPS" INDICATES WHETHER THE FAULT IS ! SLIPPING AT ITS MIDPOINT. OTHERWISE, IT IS IN THE ARTIFICIAL ! LINEARIZED REGIME, WITH STIFFNESS "FMUMAX". ! (8) "FPEAKS" GIVES THE PEAK SHEAR STRESS AT THE MIDPOINT OF EACH ! FAULT, EVALUATED AT THE BRITTLE/DUCTILE TRANSITION. ! (9) FAULTS WITH DIP LESS THAN "SLIDE" (RADIANS) ARE LIMITED ! TO A MAXIMUM DOWN-DIP INTEGRAL SHEAR TRACTION OF "TAUMAX". ! NOTE THAT PORE PRESSURES ARE CONSIDERED IN THE CALCULATION OF ! FRICTIONAL STRENGTH: ! *NORMAL PORE PRESSURES REDUCE THE EFFECTIVE NORMAL STRESS ON THE ! FAULT SURFACE BY THE AMOUNT ! -BIOT*GMEAN*RHOH20*Z ! *IF (OFFMAX > 0.), THEN THE REMAINING EFFECTIVE FRICTIONAL STRENGTH ! OF THE FAULT IS MULTIPLIED BY THE REDUCING FACTOR ! *(1.-BYERLY*OFFSET(I)/OFFMAX). ! THIS IS ALSO A PORE PRESSURE EFFECT, BECAUSE BYERLY'S MODEL IS ! THAT GOUGE LAYERS HAVE THICKNESS IN PROPORTION TO OFFSET, AND ! THAT THEY SUPPORT NON-DARCY STATIC PORE PRESSURE GRADIENTS WHICH ! REDUCE THE EFFECTIVE FRICTION OF THE FAULT. ! FOLLOWING PARAMETER GIVES NUMBER OF STEPS IN VERTICAL INTEGRAL ! OF CREEP SHEAR STRESS ON DUCTILE PARTS OF FAULTS: PARAMETER (nstep = 20) ! HIGHER VALUES OBVIOUSLY COST MORE. ON THE OTHER HAND, SMALL VALUES ! DO NOT MERELY APPROXIMATE THE CREEP LAW; THEY ALSO INTRODUCE ! SOME RANDOM ERR0R WHICH CAN PUT A FLOOR ON CONVERGENCE. ! NOTE: IN VS-FORTRAN, FOLLOWING TYPE COULD BE LOGICAL*1: LOGICAL fslips LOGICAL locked, puress, sloped, sphere DOUBLE PRECISION dpt1, v DOUBLE PRECISION fphi COMMON / fphis / fphi REAL mantle, normal ! DIMENSIONS PER COMMON BLOCK: DIMENSION fphi(4, 7) ! DIMENSIONS OF INTERNAL CONVENIENCE ARRAYS: DIMENSION dlepdz(2), dsfdz(2), rho(2), sheart(2), tmean(2), ztrans(2) ! DIMENSIONS OF EXTERNAL ARGUMENTS ARRAYS: DIMENSION acreep(2), alphat(2), bcreep(2), ccreep(2), conduc(2), & & dcreep(2), dqdtda(mxnode), elev(mxnode), & & fc(2, 2, 7, mxfel), fdip(2, mxfel), & & fimudz(7, mxfel), fpeaks(2, mxfel), & & fpflt(2, 2, 2, 7, mxfel), fslips(mxfel), & & farg(2, mxfel), ftstar(2, 7, mxfel), nodef(4, mxfel), & & offset(mxfel), radio(2), rhobar(2), & & taumax(2), tlnode(mxnode), & & v(2, mxnode), zmnode(mxnode), ztranf(2, mxfel) ! FOLLOWING TWO NUMBERS ARE "VERY SMALL" AND "VERY LARGE", BUT NOT ! SO EXTREME AS TO CAUSE UNDERFLOW OR OVERFLOW. THEY MAY NEED TO ! BE ADJUSTED, DEPENDING ON THE COMPUTER AND COMPILER BEING USED. DATA tiny / 2.E-38 / DATA huge / 1.E+38 / cgamma = (1. + SIN(ATAN(cfric))) / (1. - SIN(ATAN(cfric))) DO 100 i = 1, nfl IF (offmax <= 0.) THEN fric = ffric ELSE fric = ffric * (1. - byerly * offset(i) / offmax) END IF n1 = nodef(1, i) n2 = nodef(2, i) n3 = nodef(3, i) n4 = nodef(4, i) ! IS THIS A PURELY STRIKE-SLIP FAULT ELEMENT? puress = (ABS(fdip(1, i) - 1.570796) <= wedge).AND. & & (ABS(fdip(2, i) - 1.570796) <= wedge) ! IF SO, COMPUTE ESTIMATE OF RELATIVE NORMAL STRESS ! (RELATIVE TO VERTICAL STRESS) BY USING AMOUNT OF DIVERGENCE ! BETWEEN AVERAGE OF NODE N1 AND N2 AND AVERAGE OF NODE N3 ! AND NODE N4 (IN SPITE OF CONSTRAINT EQUATION): IF (puress) THEN angle = Chord(farg(1, i), 0.5D0, farg(2, i)) unitbx = SIN(angle) unitby = -COS(angle) delvx = v(1, n1) * fpflt(1, 1, 1, 4, i) + v(2, n1) * fpflt(2, 1, 1, 4, i) & & + v(1, n2) * fpflt(1, 1, 2, 4, i) + v(2, n2) * fpflt(2, 1, 2, 4, i) & & - v(1, n3) * fpflt(1, 1, 2, 4, i) - v(2, n3) * fpflt(2, 1, 2, 4, i) & & - v(1, n4) * fpflt(1, 1, 1, 4, i) - v(2, n4) * fpflt(2, 1, 1, 4, i) delvy = v(1, n1) * fpflt(1, 2, 1, 4, i) + v(2, n1) * fpflt(2, 2, 1, 4, i) & & + v(1, n2) * fpflt(1, 2, 2, 4, i) + v(2, n2) * fpflt(2, 2, 2, 4, i) & & - v(1, n3) * fpflt(1, 2, 2, 4, i) - v(2, n3) * fpflt(2, 2, 2, 4, i) & & - v(1, n4) * fpflt(1, 2, 1, 4, i) - v(2, n4) * fpflt(2, 2, 1, 4, i) spread = delvx * unitbx + delvy * unitby deltau = constr * spread tlan = 0.5 * (tlnode(n1) + tlnode(n2)) zman = 0.5 * (zmnode(n1) + zmnode(n2)) IF ((tlan <= 0.).OR.(ztranf(2, i) <= 0.)) THEN ! CRUST ALONE RESISTS CONVERGENCE: dpmax = -2. * deltau / ztranf(1, i) ddpndz = dpmax / ztranf(1, i) ELSE ! MANTLE LITHOSPHERE HELPS TO RESIST CONVERGENCE: ddpndz = -deltau / & & (0.5 * ztranf(1, i)**2 + ztranf(2, i) * zman + & & 0.5 * ztranf(2, i)**2) END IF ! DDPNDZ IS THE GRADIENT OF EXCESS NORMAL PRESSURE (IN ! EXCESS OF VERTICAL PRESSURE) WITH DEPTH ON THIS FAULT; ! CHECK THAT IT LIES WITHIN FRICTIONAL LIMITS OF BLOCKS: q = 0.25 * (dqdtda(n1) + dqdtda(n2) + & & dqdtda(n3) + dqdtda(n4)) ttrans = tsurf + ztranf(1, i) * q / conduc(1) - & & ztranf(1, i)**2 * radio(1) / (2. * conduc(1)) tmeanc = (tsurf + ttrans) / 2. rhoc = rhobar(1) * (1. - alphat(1) * tmeanc) dlepdc = gmean * (rhoc - rhoh2o * biot) thrust = dlepdc * cgamma normal = dlepdc / cgamma ddpndz = MAX(ddpndz, normal - dlepdc) ddpndz = MIN(ddpndz, thrust - dlepdc) ELSE ! DIFFERENT LOGIC WILL BE USED; THIS PARAMETER IS NOT ! REALLY NEEDED. ZERO IT JUST TO BE CAREFUL. ddpndz = 0. END IF DO 90 m = 1, 7 ! ELEVATION: elevat = elev(n1) * fphi(1, m) + elev(n2) * fphi(2, m) ! HEAT-FLOW: q = dqdtda(n1) * fphi(1, m) + dqdtda(n2) * fphi(2, m) ! CRUSTAL THICKNESS: crust = zmnode(n1) * fphi(1, m) + zmnode(n2) * fphi(2, m) ! MANTLE LITHOSPHERE THICKNESS: mantle = tlnode(n1) * fphi(1, m) + tlnode(n2) * fphi(2, m) mantle = MAX(mantle, 0.) ! MOHO TEMPERATURE: tmoho = tsurf + crust * q / conduc(1) - & & crust**2 * radio(1) / (2. * conduc(1)) ! TEMPERATURE AT BASE OF plate: tasth = tmoho + mantle * (q - crust * radio(1)) / conduc(2) - & & mantle**2 * radio(2) / (2. * conduc(2)) ! MEAN TEMPERATURES: tmean(1) = (tsurf + tmoho) / 2. tmean(2) = (tmoho + tasth) / 2. ! MEAN DENSITIES: rho(1) = rhobar(1) * (1. - alphat(1) * tmean(1)) rho(2) = rhobar(2) * (1. - alphat(2) * tmean(2)) ! DERIVITIVES OF LITHOSTATIC EFFECTIVE PRESSURE WRT DEPTH dlepdz(1) = gmean * (rho(1) - rhoh2o * biot) epmoho = dlepdz(1) * crust dlepdz(2) = gmean * (rho(2) - rhoh2o * biot) ! ANGLE IS THE FAULT STRIKE, IN RADIANS CCLKWS FROM +X. angle = Chord(farg(1, i), fphi(2, m), farg(2, i)) ! UNITA IS A UNIT VECTOR ALONG THE FAULT, FROM N1 TO N2. unitax = COS(angle) unitay = SIN(angle) ! UNITB IS A PERPENDICULAR UNIT VECTOR, POINTING OUT ! TOWARD THE N4-N3 SIDE. unitbx = -unitay unitby = + unitax ! RELATIVE VELOCITIES ARE FOR N1-2 SIDE RELATIVE TO ! THE N4-3 SIDE: delvx = v(1, n1) * fpflt(1, 1, 1, m, i) + v(2, n1) * fpflt(2, 1, 1, m, i) & & + v(1, n2) * fpflt(1, 1, 2, m, i) + v(2, n2) * fpflt(2, 1, 2, m, i) & & - v(1, n3) * fpflt(1, 1, 2, m, i) - v(2, n3) * fpflt(2, 1, 2, m, i) & & - v(1, n4) * fpflt(1, 1, 1, m, i) - v(2, n4) * fpflt(2, 1, 1, m, i) delvy = v(1, n1) * fpflt(1, 2, 1, m, i) + v(2, n1) * fpflt(2, 2, 1, m, i) & & + v(1, n2) * fpflt(1, 2, 2, m, i) + v(2, n2) * fpflt(2, 2, 2, m, i) & & - v(1, n3) * fpflt(1, 2, 2, m, i) - v(2, n3) * fpflt(2, 2, 2, m, i) & & - v(1, n4) * fpflt(1, 2, 1, m, i) - v(2, n4) * fpflt(2, 2, 1, m, i) ! SINISTRAL STRIKE-SLIP RATE COMPONENT: sinist = delvx * unitax + delvy * unitay ! CONVERGENCE RATE COMPONENT (IN HORIZONTAL PLANE): CLOSE = delvx * unitbx + delvy * unitby ! DIP OF THE FAULT (FROM HORIZONTAL ON THE N1-N2 SIDE): dip = fdip(1, i) * fphi(1, m) + fdip(2, i) * fphi(2, m) sloped = ABS(dip - 1.570796) > wedge IF (.NOT.sloped) THEN ! CASE OF A NEAR-VERTICAL FAULT: dsfdz(1) = (dlepdz(1) + ddpndz) * fric sfmoho = dsfdz(1) * crust dsfdz(2) = (dlepdz(2) + ddpndz) * fric slip = ABS(sinist) locked = .FALSE. ELSE ! CASE OF A SHALLOW-DIPPING FAULT: ! VUPDIP IS THE UP-DIP VELOCITY COMPONENT, IN THE ! FAULT PLANE, OF THE BLOCK ON THE N1-N3 SIDE. vupdip = CLOSE / COS(dip) ! RAKE ANGLE IS MEASURED COUNTERCLOCKWISE IN ! FAULT PLANE FROM HORIZONTAL & PARALLEL TO ANGLE. rake = Atan2f(vupdip, sinist) ! DERIVITIVE OF EFFECTIVE NORMAL PRESSURE ! WITH RESPECT TO SHEAR TRACTION ON FAULT: depdst = TAN(dip) * SIN(rake) ! (NOTICE THAT WHEN SENSE OF DIP REVERSES, SIGN ! CHANGE CAUSED BY TAN(DIP) IS CANCELLED BY SIGN ! CHANGE CAUSED BY SIN(RAKE).) ! ACCORDING TO THEORY, THE EQUATION TO SOLVE IS: ! D(SHEAR_TRACTION)/DZ = ! "FRIC"*("DLEPDZ"+"DEPDST"*D(SHEAR_TRACTION)/DZ) ! THIS MAY HAVE A PHYSICAL SOLUTION (ONE WITH ! POSITIVE SHEAR_TRACTION). IF NOT, THE ! FAULT IS LOCKED. locked = (fric * depdst) >= 1.00 IF (locked) THEN dsfdz(1) = huge dsfdz(2) = huge ELSE dsfdz(1) = fric * dlepdz(1) / (1.00 - fric * depdst) sfmoho = dsfdz(1) * crust dsfdz(2) = fric * dlepdz(2) / (1.00 - fric * depdst) END IF slip = SQRT(sinist**2 + vupdip**2) END IF slip = MAX(slip, 1.E8 * tiny) ! LOCATE PLASTIC/CREEP TRANSITION(S) ! BY ITERATED HALVING OF DOMAIN: IF (mantle > 0.) THEN limit = 2 ELSE limit = 1 ztrans(2) = 0. sheart(2) = 0. END IF DO 60 layer = 1, limit topz = 0. IF (layer == 1) THEN basez = crust sf0 = 0. t0 = tsurf q0 = q z0 = 0. ELSE basez = mantle sf0 = sfmoho t0 = tmoho q0 = q - crust * radio(1) z0 = crust END IF DO 50 kiter = 1, 15 z = 0.5 * (topz + basez) zabs = z + z0 shearf = z * dsfdz(layer) + sf0 shearp = MIN(shearf, dcreep(layer)) t = t0 + q0 * z / conduc(layer) - (radio(layer) / & & (2. * conduc(layer))) * z**2 IF (zabs <= (15. * onekm)) THEN t90pc = 0.5 * zabs ELSE IF (zabs < (45. * onekm)) THEN t90pc = (405. / 8.) * onekm + & & (-7.) * zabs + & & (13. / 40.) * onekm * (zabs / onekm)**2 + & & (-1. / 300.) * onekm * (zabs / onekm)**3 ELSE t90pc = 2. * zabs END IF ! SEE TURCOTTE ET AL (1980) JGR, 85, B11, 6224-6230 strain = slip / t90pc shearc = acreep(layer) * (strain**ecreep) * & & EXP((bcreep(layer) + ccreep(layer) * z) / t) IF (shearc < shearp) THEN basez = z ELSE topz = z END IF 50 CONTINUE ztrans(layer) = 0.5 * (topz + basez) sheart(layer) = ztrans(layer) * dsfdz(layer) + sf0 60 CONTINUE ! PLASTIC PART OF VERTICAL INTEGRAL(S) OF TRACTION: ! (A) CRUST: IF (sheart(1) <= dcreep(1)) THEN vitdz = 0.5 * sheart(1) * ztrans(1) ELSE zp = ztrans(1) * dcreep(1) / sheart(1) vitdz = dcreep(1) * (ztrans(1) - 0.5 * zp) END IF ! (B) MANTLE LITHOSPHERE: IF ((mantle > 0.).AND.(sheart(2) > sfmoho)) THEN IF (sheart(2) <= dcreep(2)) THEN vitdz = vitdz + 0.5 * (sfmoho + sheart(2)) * ztrans(2) ELSE zp = ztrans(2) * (dcreep(2) - sfmoho) / & & (sheart(2) - sfmoho) zp = MAX(zp, 0.) vitdz = vitdz + 0.5 * (sfmoho + sheart(2)) * zp + & & dcreep(2) * (ztrans(2) - zp) END IF END IF ! ADD CREEP PART(S) OF INTEGRAL, USING PARABOLIC RULE sum = 0. DO 80 layer = 1, limit IF (layer == 1) THEN thick = crust t0 = tsurf q0 = q zabs = 0. ELSE thick = mantle t0 = tmoho q0 = q - crust * radio(1) zabs = crust END IF dz = (thick - ztrans(layer)) / nstep oldsc = sheart(layer) oldsc = MIN(oldsc, dcreep(layer)) z0 = ztrans(layer) DO 70 j = 1, nstep zhalf = z0 + 0.5 * dz zfull = z0 + dz azhalf = zhalf + zabs azfull = zfull + zabs thalf = t0 + q0 * zhalf / conduc(layer) - & & (radio(layer) / & & (2. * conduc(layer))) * zhalf**2 tfull = t0 + q0 * zfull / conduc(layer) - & & (radio(layer) / & & (2. * conduc(layer))) * zfull**2 IF (azhalf <= (15. * onekm)) THEN whalf = 0.5 * azhalf ELSE IF (azhalf < (45. * onekm)) THEN whalf = (405. / 8.) * onekm + & & (-7.) * azhalf + & & (13. / 40.) * onekm * (azhalf / onekm)**2 + & & (-1. / 300.) * onekm * (azhalf / onekm)**3 ELSE whalf = 2. * azhalf END IF IF (azfull <= (15. * onekm)) THEN wfull = 0.5 * azfull ELSE IF (azfull < (45. * onekm)) THEN wfull = (405. / 8.) * onekm + & & (-7.) * azfull + & & (13. / 40.) * onekm * (azfull / onekm)**2 + & & (-1. / 300.) * onekm * (azfull / onekm)**3 ELSE wfull = 2. * azhalf END IF ! SEE TURCOTTE ET AL (1980) JGR, 85, B11, 6224-6230 ehalf = slip / whalf efull = slip / wfull schalf = acreep(layer) * (ehalf**ecreep) * & & EXP((bcreep(layer) + ccreep(layer) * zhalf) & & / thalf) schalf = MIN(schalf, dcreep(layer)) scfull = acreep(layer) * (efull**ecreep) * & & EXP((bcreep(layer) + ccreep(layer) * zfull) & & / tfull) scfull = MIN(scfull, dcreep(layer)) sum = sum + dz * (0.1666667 * oldsc + & & 0.6666667 * schalf + & & 0.1666666 * scfull) z0 = zfull oldsc = scfull 70 CONTINUE 80 CONTINUE vitdz = vitdz + sum ! LIMIT SHEAR TRACTION ON SUBDUCTION ZONES ONLY: dippy = MIN(dip, 3.141592654 - dip) IF (dippy <= slide) THEN IF (elevat < 0.0) THEN ! APPLY OCEANIC SUBDUCTION ZONE LIMIT: vitdz = MIN(vitdz, taumax(1) * SIN(dip)) ELSE ! APPLY CONTINENTAL SUBDUCTION ZONE LIMIT: vitdz = MIN(vitdz, taumax(2) * SIN(dip)) END IF END IF dpt1 = (1.D0 * vitdz) / slip vimudz = MIN(dpt1, 1.D38) fimudz(m, i) = MIN(vimudz, fmumax * (crust + mantle)) ! DIPPING, OBLIQUE-SLIP INTEGRATION ! POINTS ARE ALSO CHARACTERIZED ! BY "FC" AND "FTSTAR": IF (sloped) THEN ts = sinist * fimudz(m, i) tu = vupdip * fimudz(m, i) IF (locked) THEN fc(1, 1, m, i) = fimudz(m, i) fc(1, 2, m, i) = 0. fc(2, 1, m, i) = 0. fc(2, 2, m, i) = fimudz(m, i) ELSE sinr = SIN(rake) cosr = COS(rake) tand = TAN(dip) ! *** IMPORTANT NOTE: *** ! THE FOLLOWING 7 STATEMENTS ARE -NOT- THE ! RESULT OF THEORY, BUT A TACTICAL CHOICE ! WHICH ATTEMPTS TO COMPROMISE BETWEEN ! STABILITY OF THE LINEAR SYSTEM, STABILITY ! OF THE ITERATION, AND EFFICIENCY. ! THEY MAY BE CHANGED IF THE PROGRAM DOES ! NOT CONVERGE SATISFACTORILY! tune = 2. fc(1, 1, m, i) = fimudz(m, i) * & & (1. - tune * sinr * cosr**2 * tand) fc(1, 2, m, i) = fimudz(m, i) * & & (tune * cosr**3 * tand) fc(2, 1, m, i) = fimudz(m, i) * & & (-tune * sinr**2 * cosr * tand) fc(2, 2, m, i) = fimudz(m, i) * & & (1. + tune * sinr * cosr**2 * tand) ! (OFTEN, FC(1,2) IS THE BIGGEST TERM. ! IN SOME CASES, DIAGONALS BECOME NEGATIVE. ! FOR STABILITY, BE SURE THAT THE FC ! MATRIX REMAINS POSITIVE DEFINITE: fc(1, 1, m, i) = MAX(fc(1, 1, m, i), ABS(fc(1, 2, m, i))) fc(2, 2, m, i) = MAX(fc(2, 2, m, i), ABS(fc(1, 2, m, i))) END IF ftstar(1, m, i) = ts - fc(1, 1, m, i) * sinist - & & fc(1, 2, m, i) * vupdip ftstar(2, m, i) = tu - fc(2, 1, m, i) * sinist - & & fc(2, 2, m, i) * vupdip END IF ! PROVIDE INTERESTING DIAGNOSTIC DATA AT MIDPOINTS ONLY: IF (m == 4) THEN fslips(i) = (.NOT.locked).AND. & & (fimudz(m, i) < (0.99 * fmumax * (crust + mantle))) ztranf(1, i) = ztrans(1) fpeaks(1, i) = MIN(sheart(1), dcreep(1)) ztranf(2, i) = ztrans(2) fpeaks(2, i) = MIN(sheart(2), dcreep(2)) END IF 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE Next (input, i, j, mxel, mxfel, nfl, nodef, nodes, numel, & & output, kfault, kele) ! 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. 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 RETURN END SUBROUTINE OldVel (input, iUnitT, iUnitV, mxnode, numnod, & & output, havenv, title1, title2, title3, v) ! READ OLD VELOCITY SOLUTION FROM UNIT iUnitV, OR ELSE SET LOGICAL ! VARIABLE 'HAVENV" TO FALSE. ! WRITES 3 TITLE LINES TO iUnitT. LOGICAL havenv CHARACTER*80 title1, title2, title3 DOUBLE PRECISION v DIMENSION v(2, mxnode) READ (iUnitV, '(A80)', END = 100, ERR = 100) title1 READ (iUnitV, '(A80)', END = 100, ERR = 100) title2 READ (iUnitV, '(A80)', END = 100, ERR = 100) title3 READ (iUnitV, * , END = 100, ERR = 100) ((v(j, i), j = 1, 2), i = 1, numnod) havenv = .TRUE. WRITE (iUnitT, 50) iUnitV, title1, title2, title3 50 FORMAT (/ /' VELOCITY SOLUTION WAS', & & ' READ FROM UNIT',I3,'; TITLES WERE:'/3(/' ',A80)) GO TO 900 ! ------------------(THIS SECTION EXECUTED ONLY IF READ FAILS)--------- 100 WRITE (iUnitT, 110) iUnitV 110 FORMAT (/ /' NO FURTHER VELOCITY SOLUTIONS FOUND ON UNIT', & & I3) havenv = .FALSE. ! --------------------------------------------------------------------- 900 WRITE (iUnitT, 999) 999 FORMAT (' --------------------------------------------------', & & '-----------------------------') RETURN END SUBROUTINE OnArc (input, s, theta1, phi1, theta2, phi2, & & output, theta, phi) ! Computes coordinates (THETA,PHI) = ! (colatitude, longitude) in radians ! for a point which lies on the great circle arc from ! (THETA1,PHI1) to (THETA2,PHI2). ! Input parameter S expresses the fractional distance, ! so input S = 0.0 returns (THETA1,PHI1) and ! input S = 1.0 returns (THETA2,PHI2). REAL tvec(3), uvec(3), uvec1(3), uvec2(3) ! These are all unit vectors (in the unit sphere) ! in a Cartesian coordinate system with ! X = (0E, 0N), Y = (90E, 0N), Z = North pole. uvec1(1) = SIN(theta1) * COS(phi1) uvec1(2) = SIN(theta1) * SIN(phi1) uvec1(3) = COS(theta1) uvec2(1) = SIN(theta2) * COS(phi2) uvec2(2) = SIN(theta2) * SIN(phi2) uvec2(3) = COS(theta2) comple = 1.0 - s tvec(1) = comple * uvec1(1) + s * uvec2(1) tvec(2) = comple * uvec1(2) + s * uvec2(2) tvec(3) = comple * uvec1(3) + s * uvec2(3) size = SQRT(tvec(1)**2 + tvec(2)**2 + tvec(3)**2) uvec(1) = tvec(1) / size uvec(2) = tvec(2) / size uvec(3) = tvec(3) / size equat = SQRT(uvec(1)**2 + uvec(2)**2) theta = ATAN2(equat, uvec(3)) phi = Atan2f(uvec(2), uvec(1)) RETURN END SUBROUTINE Pause() ! This routine is called after any error message, ! and also after successful program completion, ! in order to let the user have time to read the message ! when the job is run interactively. ! To run the job as a batch job (e.g., from a Windows .bat file), ! comment out the READ so that this routine does not stop ! the flow of the program. IMPLICIT none WRITE(*,"(' Press [Enter]...')") ! ------------------------------------- !CCCC Following line could be commented out so that jobs won't pause !CCCC on batch-processing computer systems. !CCCC Remove the CCCCC for interactive computer systems! !CCCC READ(*,*) ! ------------------------------------- RETURN END SUBROUTINE Prince (input, e11, e22, e12, & & output, e1, e2, u1x, u1y, u2x, u2y) ! FIND PRINCIPAL VALUES (E1,E2) OF THE SYMMETRIC 2X2 TENSOR E11 E12 ! E12 E22 ! AND ALSO THE ASSOCIATED EIGENVECTORS #1=(U1X,U1Y),#2=(U2X,U2Y). ! THE CONVENTION IS THAT E1 <= E2. r = SQRT(((e11 - e22) / 2.)**2 + e12**2) c = (e11 + e22) / 2. e1 = c - r e2 = c + r theta = Atan2f(e11 - e1, -e12) u1x = COS(theta) u1y = SIN(theta) u2x = u1y u2y = -u1x RETURN END SUBROUTINE Projec (input, iele, & & mxel, mxnode, nodes, & & s1, s2, s3, & & xnode, ynode, v, & & output, vtheta, vphi) ! COMPUTES HORIZONTAL VELOCITY (VTHETA,VPHI) AT POINT (S1,S2,S3) ! IN TRIANGULAR ELEMENT IELE. DOUBLE PRECISION v DIMENSION nodes(3, mxel), xnode(mxnode), ynode(mxnode) DIMENSION rn(3, 3), rp(3), rs(3), & & utheta(3, 3), uphi(3, 3), ut(3), up(3), & & v(2, mxnode), vcart(3, 3), vp(3), vs(3) DO 10 k = 1, 3 ! FIND CARTESIAN RADIUS VECTOR FOR EACH NODE rn(1, k) = SIN(xnode(nodes(k, iele))) * COS(ynode(nodes(k, iele))) rn(2, k) = SIN(xnode(nodes(k, iele))) * SIN(ynode(nodes(k, iele))) rn(3, k) = COS(xnode(nodes(k, iele))) ! FIND THETA-POINTING CARTESIAN UNIT VECTOR AT EACH NODE utheta(1, k) = COS(xnode(nodes(k, iele))) * & & COS(ynode(nodes(k, iele))) utheta(2, k) = COS(xnode(nodes(k, iele))) * & & SIN(ynode(nodes(k, iele))) utheta(3, k) = -SIN(xnode(nodes(k, iele))) ! FIND PHI-POINTING CARTESIAN UNIT VECTOR AT EACH NODE uphi(1, k) = -SIN(ynode(nodes(k, iele))) uphi(2, k) = COS(ynode(nodes(k, iele))) uphi(3, k) = 0.0 ! CREATE CARTESIAN VELOCITY AT EACH NODE vcart(1, k) = v(1, nodes(k, iele)) * utheta(1, k) + & & v(2, nodes(k, iele)) * uphi(1, k) vcart(2, k) = v(1, nodes(k, iele)) * utheta(2, k) + & & v(2, nodes(k, iele)) * uphi(2, k) vcart(3, k) = v(1, nodes(k, iele)) * utheta(3, k) + & & v(2, nodes(k, iele)) * uphi(3, k) 10 CONTINUE ! INTERPOLATE VELOCITY IN PLANE vp(1) = s1 * vcart(1, 1) + s2 * vcart(1, 2) + s3 * vcart(1, 3) vp(2) = s1 * vcart(2, 1) + s2 * vcart(2, 2) + s3 * vcart(2, 3) vp(3) = s1 * vcart(3, 1) + s2 * vcart(3, 2) + s3 * vcart(3, 3) ! INTERPOLATE POSITION IN PLANE rp(1) = s1 * rn(1, 1) + s2 * rn(1, 2) + s3 * rn(1, 3) rp(2) = s1 * rn(2, 1) + s2 * rn(2, 2) + s3 * rn(2, 3) rp(3) = s1 * rn(3, 1) + s2 * rn(3, 2) + s3 * rn(3, 3) size = SQRT(rp(1)**2 + rp(2)**2 + rp(3)**2) growth = 1. / size ! PROJECT POSITION TO SPHERE rs(1) = rp(1) * growth rs(2) = rp(2) * growth rs(3) = rp(3) * growth ! PROJECT VELOCITY TO VELOCITY ON SPHERE vs(1) = vp(1) * growth vs(2) = vp(2) * growth vs(3) = vp(3) * growth ! FIND THETA AND PHI UNIT VECTORS AT THIS POINT up(1) = -rs(2) up(2) = rs(1) up(3) = 0.0 CALL Unit (modify, up) CALL Cross (input, up, rs, output, ut) ! FIND HORIZONTAL COMPONENTS OF VELOCITY vtheta = vs(1) * ut(1) + vs(2) * ut(2) + vs(3) * ut(3) vphi = vs(1) * up(1) + vs(2) * up(2) + vs(3) * up(3) RETURN END SUBROUTINE PutNet (input, iUnitO, & & brief, dqdtda, elev, fdip, & & mxel, mxfel, mxnode, n1000, & & nfaken, nfl, nodef, nodes, & & nrealn, numel, numnod, offset, & & title1, tlnode, xnode, ynode, zmnode) ! WRITES FINITE ELEMENT GRID TO UNIT "iUnitO". CHARACTER*80 title1 LOGICAL brief DIMENSION dqdtda(mxnode), elev(mxnode), & & fdip(2, mxfel), np(4), & & nodef(4, mxfel), nodes(3, mxel), & & offset(mxfel), tlnode(mxnode), & & xnode(mxnode), ynode(mxnode), zmnode(mxnode) DIMENSION dips(2) 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. - xnode(i) * 57.29577951 pLon = ynode(i) * 57.29577951 WRITE (iUnitO, 91) i, pLon, pLat, elev(i), dqdtda(i), zmnode(i), & & tlnode(i) 91 FORMAT (I5,0P,2F11.5,1P,4E10.2) 100 CONTINUE WRITE (iUnitO, 110) numel 110 FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') 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 WRITE (iUnitO, 160) i, (np(k), k = 1, 3) 160 FORMAT (I5,3I5) 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.2958 IF (dips(k) > 90.01) dips(k) = dips(k) - 180. 230 CONTINUE WRITE (iUnitO, 250) i, (np(k), k = 1, 4), & & (dips(k), k = 1, 2), offset(i) 250 FORMAT (I5,4I6,2F5.0,1P,E10.2) 300 CONTINUE RETURN END SUBROUTINE ReadN (input, iUnitP, iUnitT, n, & & output, vecton) ! 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. CHARACTER*1 c CHARACTER*80 line LOGICAL anyin, dotted, expon, signed DIMENSION vecton(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, * ) (vecton(i), i = 1, number) IF (number < n) THEN DO 99 i = number + 1, n vecton(i) = 0. 99 CONTINUE END IF END IF RETURN END SUBROUTINE ReadPm (input, iunit7, iunit8, names, nPlate, offmax, & & output, acreep, alphat, bcreep, biot , & & byerly, ccreep, cfric , conduc, & & dcreep, ecreep, everyp, ffric , gmean , & & gradie, iconve, ipvref, & & maxitr, okdelv, oktoqt, onekm, radio , & & radius, refstr, rhoast, rhobar, rhoh2o, & & tadiab, & & taumax, temlim, title3, trhmax, tsurf, & & vtimes, zbasth) ! READS STRATEGIC AND TACTICAL INPUT PARAMETERS FROM DEVICE IUNIT7, ! FOLLOWED BY PLOT-CONTROL PARAMETERS, ! AND ECHOES THEM ON DEVICE IUNIT8 WITH ANNOTATIONS. CHARACTER*2 names, pltref CHARACTER*80 title3 LOGICAL everyp DIMENSION acreep(2), alphat(2), bcreep(2), ccreep(2), conduc(2), & & dcreep(2), names(nPlate), radio(2), & & rhobar(2), temlim(2), tempv(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 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 > 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)') END IF CALL ReadN (input, iunit7, iunit8, 2, & & output, acreep) WRITE (iunit8, 50) acreep(1), acreep(2) 50 FORMAT (' ',1P, E10.2,'/',E10.2,' A FOR CREEP = ', & & 'PRE-EXPONENTIAL SHEAR', & & ' STRESS CONSTANT FOR CREEP. (CRUST/MANTLE)') CALL ReadN (input, iunit7, iunit8, 2, & & output, bcreep) WRITE (iunit8, 60) bcreep(1), bcreep(2) 60 FORMAT (' ', F10.0,'/',F10.0,' B FOR CREEP =(ACTIVATION ', & & 'ENERGY)/R/N', & & ' IN K. (CRUST/MANTLE)') CALL ReadN (input, iunit7, iunit8, 2, & & output, ccreep) WRITE (iunit8, 70) ccreep(1), ccreep(2) 70 FORMAT (' ',1P, E10.2,'/',E10.2,' C FOR CREEP = DERIVATIVE OF B', & & ' WITH RESPECT TO DEPTH. (CRUST/MANTLE)') CALL ReadN (input, iunit7, iunit8, 2, & & output, dcreep) WRITE (iunit8, 80) dcreep(1), dcreep(2) 80 FORMAT (' ',1P, E10.2,'/',E10.2,' D FOR CREEP = MAXIMUM SHEAR ', & & 'STRESS UNDER ANY CONDITIONS. (CRUST/MANTLE)') READ (iunit7, * ) ecreep WRITE (iunit8, 90) ecreep 90 FORMAT (' ', F10.6,' E FOR CREEP = STRAIN-RATE EXPONENT FOR', & & ' CREEP (1/N). (SAME FOR CRUST AND MANTLE!)') READ (iunit7, * ) tadiab, gradie WRITE (iunit8, 92) tadiab, gradie 92 FORMAT (' ',F10.0,1P,E10.2,' INTERCEPT AND GRADIENT OF UPPER' & & ,' MANTLE ADIABAT (K, K/M)') READ (iunit7, * ) zbasth WRITE (iunit8, 94) zbasth 94 FORMAT (' ',1P,E10.2,' DEPTH OF BASE OF ASTHENOSPHERE') 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, nPlate IF (names(i) == pltref) ipvref = i 956 CONTINUE IF (ipvref == 0) THEN WRITE (iunit8, 958) (names(i), i = 1, nPlate) 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:' & & /' ',10X,' 0 = NONE.' & & /' ',10X,' 1 = HAGER AND OCONNELL (1979) MODEL II' & & /' ',10X,' 2 = BAUMGARDNER (1988) FIGURE 7A-F, *10' & & /' ',10X,' 3 = NUVEL-1 (DE METS ET AL., 1990)' & & /' ',10X,' 4 = NUVEL-1A DRAGS CONTINENTS, NO OCEAN DRAG' & & /' ',10X,' 5 = DRAG ON BASE OF SUBDUCTION FOREARC ONLY' & & ) IF (iconve > 0) THEN BACKSPACE iunit7 CALL ReadN (input, iunit7, iunit8, 2, & & output, tempv) IF (tempv(2) > 0) THEN vtimes = tempv(2) WRITE (iunit8, 98) vtimes 98 FORMAT (' ',F10.4,' SPEED FACTOR FOR CONVECTION', & & ' MODEL IDENTIFIED ABOVE') ELSE WRITE (iunit8, 99) 99 FORMAT (' UNINTERPRETABLE VALUE FOR VTIMES; SET ' & & ,'TO 1.0') vtimes = 1.0 END IF ELSE vtimes = 1.0 END IF READ (iunit7, * ) trhmax WRITE (iunit8, 101) trhmax 101 FORMAT (' ',1P,E10.2,' LIMIT ON HORIZONTAL TRACTIONS', & & ' APPLIED TO BASE OF plate') READ (iunit7, * ) taumax WRITE (iunit8, 102) taumax 102 FORMAT (' ',1P,E10.2,' MAX. DOWN-DIP INTEGRAL OF SHEAR TRACTION', & & ' IN SUBDUCTION ZONES (FAULTS WITH DIP < SUBDIP)') READ (iunit7, * ) rhoh2o WRITE (iunit8, 110) rhoh2o 110 FORMAT (' ',1P,E10.3,' DENSITY OF GROUNDWATER, LAKES, & OCEANS') CALL ReadN (input, iunit7, iunit8, 2, & & output, rhobar) WRITE (iunit8, 120) rhobar(1), rhobar(2) 120 FORMAT (' ',1P,E10.3,'/',E10.3,' MEAN DENSITY,', & & ' CORRECTED TO 0 DEGREES KELVIN. (CRUST/MANTLE)') READ (iunit7, * ) rhoast WRITE (iunit8, 130) rhoast 130 FORMAT (' ',1P,E10.3,' DENSITY OF ASTHENOSPHERE') READ (iunit7, * ) gmean WRITE (iunit8, 140) gmean 140 FORMAT (' ',1P,E10.3,' MEAN GRAVITATIONAL ACCELERATION', & & ' (LENGTH/SEC**2)') READ (iunit7, * ) onekm WRITE (iunit8, 150) onekm 150 FORMAT (' ',1P,E10.3,' NUMBER OF LENGTH UNITS NEEDED TO', & & ' MAKE 1 KILOMETER'/11X, & & ' (E.G., 1000. IN SI, 1.E5 IN CGS)') READ (iunit7, * ) radius WRITE (iunit8, 155) radius 155 FORMAT (' ',1P,E10.3,' RADIUS OF THE PLANET') CALL ReadN (input, iunit7, iunit8, 2, & & output, alphat) WRITE (iunit8, 160) alphat(1), alphat(2) 160 FORMAT (' ',1P,E10.2,'/',E10.2,' THERMAL ', & & 'EXPANSION', & & ' (1/VOL)*(D.VOL/D.T). (CRUST/MANTLE)') CALL ReadN (input, iunit7, iunit8, 2, & & output, conduc) WRITE (iunit8, 170) conduc(1), conduc(2) 170 FORMAT (' ',1P,E10.2,'/',E10.2,' CONDUCTIVITY, ENERGY/', & & 'LENGTH/SEC/DEG. (CRUST/MANTLE)') CALL ReadN (input, iunit7, iunit8, 2, & & output, radio) WRITE (iunit8, 180) radio(1), radio(2) 180 FORMAT (' ',1P,E10.2,'/',E10.2,' RADIOACTIVE HEATING', & & ' ENERGY/VOLUME/SEC. (CRUST/MANTLE)') READ (iunit7, * ) tsurf WRITE (iunit8, 185) tsurf 185 FORMAT (' ', F10.0,' SURFACE TEMPERATURE, ON', & & ' ABSOLUTE SCALE') CALL ReadN (input, iunit7, iunit8, 2, & & output, temlim) WRITE (iunit8, 190) temlim(1), temlim(2) 190 FORMAT (' ', F10.0,'/',F10.0,' CONVECTING TEMPERATURE (TMAX),' & & ,' KELVIN. (CRUST/MANTLE)') WRITE (iunit8, 199) 199 FORMAT (/' ========== TACTICAL PARAMETERS (HOW TO REACH ', & & 'THE SOLUTION) =========='/) READ (iunit7, * ) maxitr WRITE (iunit8, 200) maxitr 200 FORMAT (' ',I10,' MAXIMUM ITERATIONS WITHIN VELOCITY SOLUTION') READ (iunit7, * ) oktoqt WRITE (iunit8, 210) oktoqt 210 FORMAT (' ',F10.6,' ACCEPTABLE FRACTIONAL CHANGE IN VELOCITY ', & & '(STOPS ITERATION EARLY)') READ (iunit7, * ) refstr WRITE (iunit8, 220) refstr 220 FORMAT (' ',1P,E10.2,' EXPECTED MEAN VALUE OF SHEAR STRESS IN', & & ' CRUST'/' ',10X, & & ' (USED TO INITIALIZE AND SET STIFFNESS LIMITS)') READ (iunit7, * ) okdelv WRITE (iunit8, 230) okdelv 230 FORMAT (' ',1P,E10.2,' MAGNITUDE OF VELOCITY ERR0RS ALLOWED', & & ' DUE TO FINITE STIFFNESS'/11X, & & '(SUCH ERR0RS MAY APPEAR IN SUCH FORMS AS:'/11X, & & ' 1. FICTICIOUS BASAL SLIP OF CRUST OVER MANTLE'/11X, & & ' 2. ERRONEOUS CONVERGENCE/DIVERGENCE AT VERTICAL FAULTS'/ & & 11X, & & ' 3. VELOCITY EFFECT OF FICTICIOUS VISCOUS COMPLIANCES'/11X, & & ' HOWEVER, VALUES WHICH ARE TOO SMALL WILL CAUSE ILL-CONDITIONED' & & /11X, & & ' LINEAR SYSTEMS AND STRESS ERR0RS, ', & & 'AND MAY PREVENT CONVERGENCE!)') READ (iunit7, * ) everyp WRITE (iunit8, 240) everyp 240 FORMAT (' ',L10,' SHOULD NODAL VELOCITIES BE OUTPUT EVERY STE', & & 'P? (FOR CONVERGENCE STUDY)') WRITE (iunit8, 999) 999 FORMAT (' --------------------------------------------------', & & '-----------------------------') RETURN END REAL FUNCTION Atan2f (y, x) ! CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS IF ((y /= 0.).OR.(x /= 0.)) THEN Atan2f = ATAN2(y, x) ELSE Atan2f = 0. END IF RETURN END SUBROUTINE Snodal (input, phi, theta, & & output, fpp) ! CALCULATES VECTOR NODAL FUNCTION AT INTERPOLATION POINT ON A ! GREAT CIRCLE OF A SIDE OF A FINITE ELEMENT DOUBLE PRECISION fphi, pp DIMENSION fphi(4, 7), fpp(2, 2, 2, 7), phi(2), theta(2) COMMON / fphis / fphi ! SNCCOP: SIN(SITA)*COS(PHAI) AT INTERPOLATION POINT ! SNCSNP SIN(SITA)*SIN(PHAI) AT INTERPOLATION POINT x1 = SIN(theta(1)) * COS(phi(1)) y1 = SIN(theta(1)) * SIN(phi(1)) z1 = COS(theta(1)) x2 = SIN(theta(2)) * COS(phi(2)) y2 = SIN(theta(2)) * SIN(phi(2)) z2 = COS(theta(2)) xn = x1 + x2 yn = y1 + y2 zn = z1 + z2 xyzn = SQRT(xn * xn + yn * yn + zn * zn) xn = xn / xyzn yn = yn / xyzn zn = zn / xyzn dd = x1 * xn + y1 * yn + z1 * zn DO 800 m = 1, 7 xx = fphi(1, m) * x1 + fphi(2, m) * x2 yy = fphi(1, m) * y1 + fphi(2, m) * y2 zz = fphi(1, m) * z1 + fphi(2, m) * z2 pp = SQRT(xx * xx + yy * yy + zz * zz) xx = xx / pp yy = yy / pp zz = zz / pp sitaj = ACOS(zz) phaij = Atan2f(yy, xx) rn = xx * xn + yy * yn + zz * zn ppm = rn / dd cscs = COS(sitaj) * COS(phaij) cssn = COS(sitaj) * SIN(phaij) snsn = SIN(sitaj) * SIN(phaij) snc = SIN(sitaj) snp = SIN(phaij) csp = COS(phaij) DO 500 j = 1, 2 fp = fphi(j, m) * ppm fpp(1, 1, j, m) = ( COS(theta(j)) * COS(phi(j)) * cscs & & + COS(theta(j)) * SIN(phi(j)) * cssn & & + SIN(theta(j)) * snc) * fp fpp(2, 1, j, m) = (-SIN(phi(j)) * cscs + COS(phi(j)) * cssn) * fp fpp(1, 2, j, m) = (-COS(theta(j)) * COS(phi(j)) * snp & & + COS(theta(j)) * SIN(phi(j)) * csp) * fp fpp(2, 2, j, m) = ( SIN(phi(j)) * snp & & + COS(phi(j)) * csp) * fp 500 CONTINUE 800 CONTINUE RETURN END SUBROUTINE Square (input, brief, fdip, iunit8, & & mxbn, mxel, mxfel, mxnode, & & mxstar, nfl, nodef, nodes, & & numel, numnod, radius, wedge, & & modify, xnode, ynode, & & output, area, detj, & & dxs, dys, dxsp, dysp, edgefs, & & edgets, flen, fpflt, fpsfer, & & farg, ncond, nodcon, sita, & & work, checkn, list) ! CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID LOGICAL agreed, allok, brief, found, switch, vert1, vert2 ! NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL 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 ONE dx = xnode(nj1) - xnode(nj2) dy = ynode(nj1) - ynode(nj2) IF (dy > 3.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 dy = dy * SIN(xnode(nj1)) short = SQRT(dx**2 + dy**2) DO 470 k = 1, nfl nl1 = nodef(1, k) nl2 = nodef(2, k) nl3 = nodef(3, k) nl4 = nodef(4, k) IF ((nj1 == 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.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 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. 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.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 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. 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 > 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!') 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 Deriv3 (input, iUnitT, mxel, mxnode, nodes, numel, & & radius, xnode, ynode, & & output, area, detj, & & dxs, dys, dxsp, dysp, fpsfer, sita) allok = .TRUE. DO 620 i = 1, numel DO 610 m = 1, 7 test = area(i) * detj(m, i) IF (test <= 0.) THEN WRITE(iunit8, 605) m, i 605 FORMAT(/' EXCESSIVELY DISTORTED ELEMENT LEADS TO ' & & ,'NEGATIVE AREA AT POINT ',I1,' IN ELEMENT ', & & I5) WRITE(iunit8, 606) area(i), detj(m, i) 606 FORMAT('AREA = ',1P,E12.4,' DETJ: ',0P,F12.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, radius, 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 A 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 (input, i, j, mxel, mxfel, nfl, nodef, nodes, numel, & & output, kfault, kele) 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.E-6) .AND. & & (ABS(ynode(i) - y) < 1.E-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(' ',2I6) 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 WARNING IF ANY ELEMENT IS OF ! MIXED TYPE (PART STRIKE-SLIP, AND PART SHALLOW-DIPPING: DO 920 i = 1, nfl deld1 = fdip(1, i) - 1.570796 deld2 = fdip(2, i) - 1.570796 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.2957795 IF (dip1 > 90.) dip1 = dip1 - 180. dip2 = fdip(2, i) * 57.2957795 IF (dip2 > 90.) dip2 = dip2 - 180. WRITE (iunit8, 905) i, dip1, tag1, dip2, tag2 905 FORMAT(/ /' CAUTION:'/ & & ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ & & ' ',F7.2,' DEGREES ',A21/ & & ' ',F7.2,' DEGREES ',A21/ & & ' WHICH MAKES IT MIXED-MODE.'/ & & ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ & & ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ & & ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', & & ' FAULT ELEMENT TO THE NEXT.)') ELSE nvpart = 0 DO 910 m = 1, 7 deld = deld1 * fphi(1, m) + deld2 * fphi(2, m) IF (ABS(deld) <= 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 PLANE 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 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) alfad = flen(i) / radius CALL Fangls(input, phi, theta, & & output, fangle) DO 900 j = 1, 2 farg(j, i) = fangle(j) fangle(j) = fangle(j) * 57.29577951 900 CONTINUE 1000 CONTINUE ! (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.570796) <= 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.5708) <= 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.5708) <= 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.141592654 ELSE azl = farg(nazl, number) END IF azi = farg(nazi, i) cosz = 0.5 * (COS(azi) + COS(azl)) sinz = 0.5 * (SIN(azi) + SIN(azl)) azimut = Atan2f(sinz, cosz) farg(nazi, i) = azimut IF(nazl == nazi) THEN farg(nazl, number) = azimut - 3.141592654 ELSE farg(nazl, number) = azimut END IF IF (ABS(azi - azl) > 0.02) THEN dazi = azi * 57.2957795 dazl = azl * 57.2957795 daz = azimut * 57.2957795 np1 = node1 np4 = node4 delon = 57.2957795 * ynode(node1) dnlat = 90. - 57.2957795 * xnode(node1) ! WRITE STATEMENT HERE WAS REMOVED; AT SCORING TIME IS TOO LATE TO FIX! 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 (input, iUnitT, mxfel, & & mxnode, nfl, nodef, xnode, ynode, & & output, fpflt) IF (.NOT. brief) WRITE (iunit8, 9999) 9999 FORMAT (' --------------------------------------------------', & & '-----------------------------') RETURN END SUBROUTINE ThetaPhi_2_pLate (iUnitT, & & names, nPBnd, nBoundaryPoints, & & nPlate, omega, pLat, pLon, & & theta, phi, & ! INTENT(IN) & iPlate) ! INTENT(OUT) ! Assigns an integer plate ID# to one surface point ! with coordinates (theta = colatitude, in radians, ! phi = longitude, in radians). ! NOTE: This subprogram should NOT be called to assign ! F-E nodes to pLates if those nodes might belong to ! fault elements along plate edges! Results would be purely ! random, and would not consider topological connectivity ! to an adjacent plate. Instead, use subprogram ASSIGN ! from Shells (summer 2006). INTEGER, intent(in) :: iUnitT CHARACTER*2, DIMENSION(nPlate) :: names INTEGER, intent(in) :: nPBnd INTEGER, DIMENSION(nPlate), intent(in) :: nBoundaryPoints INTEGER, intent(in) :: nPlate REAL, DIMENSION(3, nPlate), intent(in) :: omega REAL, DIMENSION(nPlate, nPBnd), intent(in):: pLat, pLon REAL, intent(in) :: theta, phi INTEGER, intent(out) :: iPlate LOGICAL :: gotout REAL, DIMENSION(3) :: alongv, crossv ! PB2002 model of Bird [2003; G**3]; ! Already has plate NAMES and OMEGA vectors in ! main program (DATA statements); ! must also have digitised plate ! outlines in arrays pLat and pLon, ! presumably read from file "PB2002_plates.dig". ! That is, this routine will not read any file. xo = COS(phi) * SIN(theta) yo = SIN(phi) * SIN(theta) zo = COS(theta) oxyz = xo * xo + yo * yo + zo * zo oxyz = SQRT(oxyz) xo = xo / oxyz yo = yo / oxyz zo = zo / oxyz npoint = 0 angle = 0.0 iPlate = 0 DO 500 ip = 1, nPlate tangl = 0.0 nend = nBoundaryPoints(ip) DO 300 j = 1, nend j2 = j + 1 IF(j == nend) THEN j2 = 1 END IF a1 = COS(pLon(ip, j)) * COS(pLat(ip, j)) a2 = SIN(pLon(ip, j)) * COS(pLat(ip, j)) a3 = SIN(pLat(ip, j)) b1 = COS(pLon(ip, j2)) * COS(pLat(ip, j2)) b2 = SIN(pLon(ip, j2)) * COS(pLat(ip, j2)) b3 = SIN(pLat(ip, j2)) ao = xo * a1 + yo * a2 + zo * a3 bo = xo * b1 + yo * b2 + zo * b3 a1 = a1 / ao a2 = a2 / ao a3 = a3 / ao b1 = b1 / bo b2 = b2 / bo b3 = b3 / bo a1 = a1 - xo a2 = a2 - yo a3 = a3 - zo b1 = b1 - xo b2 = b2 - yo b3 = b3 - zo aa = SQRT(a1 * a1 + a2 * a2 + a3 * a3) bb = SQRT(b1 * b1 + b2 * b2 + b3 * b3) ab1 = a2 * b3 - a3 * b2 ab2 = a3 * b1 - a1 * b3 ab3 = a1 * b2 - a2 * b1 stheta = (ab1 * xo + ab2 * yo + ab3 * zo) / (aa * bb) ! prevent stupid abends due to imprecision: stheta = MAX(-1., MIN(1., stheta)) tangl = tangl + ASIN(stheta) 300 CONTINUE dangle = tangl - 3.1416 IF(dangle >= 0.0001) THEN npoint = npoint + 1 iPlate = ip END IF 500 CONTINUE IF(npoint >= 3) THEN xpoint = 90.0 - theta * 57.29577951 ypoint = phi * 57.29577951 WRITE(iUnitT, 505) xpoint, ypoint 505 FORMAT(' POINT ',2F10.3,' WAS FOUND IN MORE THAN TWO' & & ,'pLatES; SOMETHING IS WRONG !!!!') CALL PAUSE() STOP END IF IF (iPlate == 0) THEN xpoint = 90.0 - theta * 57.29577951 ypoint = phi * 57.29577951 WRITE(iUnitT, 600) xpoint, ypoint CALL PAUSE() STOP 600 FORMAT('THE POINT ', 2F10.2,' DOES NOT BELONG TO' & & ,' ANY plate !!!!') CALL PAUSE() STOP END IF END SUBROUTINE Tractor(iUnitQ, iUnitT, nPlates, nPoints, & & slab_q, phi_list, theta_list, whichP, & ! INTENT(IN) & basal_shear_tractions) ! INTENT(OUT) ! Requests file name of an existing torque report ! (including traction pole vectors for each plate, ! created by a previous experiment with SHELLS, ! usually one that had TRHMAX = 0. and extra internal ! velocity boundary conditions for each slabless plate). ! Reads this file, extracts the traction pole vectors, ! and uses them to precompute basal_shear_tractions shear tractions ! at each point specified in (theta_list, phi_list) = ! (colatitude, longitude) in radians. These points must ! already be identified as to plate affinity, in whichP. ! For further clarification of "traction pole vectors" ! see subprogram -TWIST- in SHELLS. ! Derived from subprogram TRACT of Shells (summer 2006), ! but modified to free-format, and variable names changed ! to drop previous implication that the list of points ! to be processed is a list of node locations. IMPLICIT NONE INTEGER, INTENT(IN) :: iUnitQ, iUnitT, nPlates, nPoints LOGICAL, DIMENSION(nPlates), INTENT(IN) :: slab_q REAL, DIMENSION(nPoints), INTENT(IN) :: theta_list, phi_list INTEGER, DIMENSION(nPoints), INTENT(IN) :: whichP REAL, DIMENSION(2, nPoints), INTENT(OUT):: basal_shear_tractions CHARACTER*132 line, trqfil INTEGER i, ios, iPlate, j LOGICAL, DIMENSION(:), allocatable :: tpread REAL equat, lat, length, lon, t, tequat REAL, DIMENSION(3) :: tvec, uphi, utheta, uvec REAL, DIMENSION(:, :), allocatable :: tpvecs ALLOCATE (tpvecs(3, nPlates)) ALLOCATE (tpread(nPlates)) ! Zero whole array; advisable because some plates may not ! appear in report. DO 30 j = 1, nPlates DO 20 i = 1, 3 tpvecs(i, j) = 0. 20 CONTINUE tpread(j) = .FALSE. 30 CONTINUE WRITE(iUnitT, 35) iUnitQ 35 FORMAT(/' ATTEMPTING TO READ TORQUE REPORT ON UNIT',I3/) ! Waste first 6 lines (titles & 2 blanks & header) of torque file. DO 40 i = 1, 6 READ (iUnitQ, "(A)") line 40 CONTINUE ! Loop on plates in report (up to nPlates for whole-Earth model): DO 100 j = 1, nPlates READ(iUnitQ, * , iostat = ios) IF (ios == -1) GO TO 101 READ(iUnitQ, "(8X,I6)", iostat = ios) iPlate IF (ios == -1) GO TO 101 ! Waste 23 more lines of each plate report DO 50 i = 1, 23 READ(iUnitQ, * , iostat = ios) IF (ios == -1) GO TO 101 50 CONTINUE READ(iUnitQ, "(56X,ES10.3,2F10.2)") t, lon, lat ! T is magnitude, in Pa, at location 90 deg. from (LON, LAT). tpvecs(1, iPlate) = t * COS(lat / 57.296) * COS(lon / 57.296) tpvecs(2, iPlate) = t * COS(lat / 57.296) * SIN(lon / 57.296) tpvecs(3, iPlate) = t * SIN(lat / 57.296) tpread(iPlate) = .TRUE. ! Waste 14 lines to get past the "=======" at the bottom of ! each torque report: DO 60 i = 1, 14 READ(iUnitQ, * , iostat = ios) IF (ios == -1) GO TO 101 60 CONTINUE 100 CONTINUE 101 CLOSE(iUnitQ) DO 200 i = 1, nPoints iPlate = whichP(i) IF (slab_q(iPlate)) THEN ! no need for inferred basal_shear_tractions-strength traction: basal_shear_tractions(1, i) = 0.0D0 basal_shear_tractions(2, i) = 0.0D0 ELSE IF (tpread(iPlate)) THEN ! Uvec is unit vector to node location: uvec(1) = SIN(theta_list(i)) * COS(phi_list(i)) uvec(2) = SIN(theta_list(i)) * SIN(phi_list(i)) uvec(3) = COS(theta_list(i)) ! Tvec is cross-product with traction pole vector: tvec(1) = tpvecs(2, iPlate) * uvec(3) - & & tpvecs(3, iPlate) * uvec(2) tvec(2) = tpvecs(3, iPlate) * uvec(1) - & & tpvecs(1, iPlate) * uvec(3) tvec(3) = tpvecs(1, iPlate) * uvec(2) - & & tpvecs(2, iPlate) * uvec(1) t = SQRT(tvec(1)**2 + tvec(2)**2 + tvec(3)**2) ! Unit vectors at this site (NOT a pole): uphi(1) = -uvec(2) uphi(2) = uvec(1) equat = SIN(theta_list(i)) uphi(1) = uphi(1) / equat uphi(2) = uphi(2) / equat uphi(3) = 0.0 tequat = uvec(3) utheta(3) = -equat utheta(1) = tequat * uvec(1) / equat utheta(2) = tequat * uvec(2) / equat length = SQRT(utheta(1)**2 + utheta(2)**2 + & & utheta(3)**2) utheta(1) = utheta(1) / length utheta(2) = utheta(2) / length utheta(3) = utheta(3) / length ! Horizontal components of shear traction: basal_shear_tractions(1, i) = tvec(1) * utheta(1) + tvec(2) * utheta(2) + & & tvec(3) * utheta(3) basal_shear_tractions(2, i) = tvec(1) * uphi(1) + tvec(2) * uphi(2) + & & tvec(3) * uphi(3) ELSE basal_shear_tractions(1, i) = 0.0D0 basal_shear_tractions(2, i) = 0.0D0 END IF END IF 200 CONTINUE DEALLOCATE (tpread) DEALLOCATE (tpvecs) RETURN END SUBROUTINE Unit (modify, anyvec) ! CONVERTS ANY 3-COMPONENT VECTOR TO A UNIT VECTOR DIMENSION anyvec(3) r2 = anyvec(1) * anyvec(1) + anyvec(2) * anyvec(2) + & & anyvec(3) * anyvec(3) IF (r2 > 0.) THEN size = 1. / SQRT(r2) anyvec(1) = anyvec(1) * size anyvec(2) = anyvec(2) * size anyvec(3) = anyvec(3) * size ELSE anyvec(1) = 1. anyvec(2) = 0. anyvec(3) = 0. END IF RETURN END