C-------------------------------------------------------------------- C PROGRAM OrbData5 C C An update (2005.06, 2006.08) of the existing program OrbData. C Note that the "5" in "OrbData5" alludes to both "Earth5" C and to "2005". C Both versions are intended to fill in nodal data C (elevation, heat flow, crustal thickness, mantle C lithosphere thickness) in finite-element grids (*.feg) C intended for use with program Shells. C C The difference is that: C -OrbData computed crustal thickness and mantle lithosphere C thickness from the elevation and heat flow, C under assumptions of local isostasy and C steady-state geotherms. C -OrbData5 takes crustal thickness from the CRUST2 model, C and obtains mantle lithosphere thickness either C from seismic travel-time anomalies (in continents) or C from an ad-hoc analytic function of age (in oceans). C -A global grid of sea-floor ages (or warning flags), C such as age_1p5.grd derived from Mueller et al. [1997a], C was optional in OrbData but is required by OrbData5. C C The necessary 2 new degrees of freedom at each node are: C *Chemical density anomaly of the entire lithosphere; C *Extra curvature of the geotherm (throughout the lithosphere) C due to transient cooling (common) or heating (rare). C These extra nodal data are appended to each nodal line C in the .feg file, where they will be read and used C (if present) by the latest versions of OrbNumber, C OrbWin, Shells, and FiniteMap. C Chemical density anomaly is subject to some fairly C arbitrary limits in size, so that seafloor trenches C will NOT be in isostatic balance, and will still imply C anomalous low pressure (due to slab pull) at the base C of the model lithosphere column. C In other places, computed chemical density anomalies may C be real, OR may be artifacts of errors in heat flow, C crustal thickness (from CRUST2) or mantle lithosphere C thickness (from seismic travel-time anomalies); C therefore, be VERY CAREFUL when presenting and interpreting C them; they will often be artifacts! C Extra curvature of the geotherm is subject to the upper bound C that the geotherm should be monotonically increasing with C depth, so that it is not allowed to overshoot into the C temperature range of the asthenosphere. This may result C in mantle lithosphere thicknesses which are locally less C than would be dictated by the guiding principles of OrbData5. C C Some featues are common to both OrbData5 and OrbData: C -Only zero values of elevation and heat flow will be replaced; C non-zero values are assumed to result from hand-editing C and will be preserved (unless heat flow falls below the C lower limit). C -All values of crustal thickness and mantle lithosphere C thickness are replaced in every run. Existing values C are simply ignored. C C Written in Fortran 77 (fixed-format) syntax, C but the need for dynamic memory allocation (ALLOCATE) C for the (up to) 5 huge arrays of elevation, C heat-flow, seafloor age, crustal thickenss, and C vertical S wave travel-time anomaly C (between the Moho and 400 km depth) C has now forced this code into Fortran 90. C Therefore, it now has ALLOCATE in main, C REAL, DIMENSION(:,:) in both main and ASSIGN, C an interface block for ASSIGN, and matching C Fortran 90 parameter declarations in ASSIGN. C C By Peter Bird, Dept. of Earth & Space Sciences, C University of California, Los Angeles, CA 90095-1567. C (For version date, see 1st WRITE below.) C C Comments below are inherited from OrbData (TEXT IN ALL-CAPS) C and occasionally edited for OrbData5 (lower-case text): C------------------------------------------------------------------ C C READS A FINITE-ELEMENT GRID PRODUCED BY ORBWEAVE or OrbWin C (in Shells mode) C AND A PARAMETER FILE (FORMATTED FOR SHELLS) WITH C THERMAL CONDUCTIVITIES, DENSITIES, ET CETERA; C THEN, FILLS IN OR ADJUSTS NODAL DATA: C -ELEVATION (+ ABOVE SEA LEVEL, - BELOW); C -HEAT-FLOW (ASSUMING STEADY-STATE GEOTHERM); C -CRUSTAL THICKNESS; C -MANTLE LITHOSPHERE THICKNESS (NOT INCLUDING CRUST). C C ALL UNITS ARE SI: m, W/m**2 (not mW/m**2), m, m. C C IN THE OUTPUT DATA, THE BASE OF THE MANTLE LITHOSPHERE C WILL BE AN ISOTHERMAL SURFACE. THE ISOTHERM VALUE C IS CHOSEN TO LIE ON THE ASTHENOSPHERE ADIABAT C (EVALUATED AT AN ARBITRARY DEPTH OF 100 km). C C IF A NODAL ELEVATION IS NON-ZERO, IT WILL ALWAYS BE LEFT C UNCHANGED, TO PRESERVE EFFECTS OF HAND-EDITING. C IF ANY ELEVATION(S) IS/ARE ZERO, NEW VALUE(S) WILL BE C INTERPOLATED FROM A GRIDDED DATASET. C C IF A NODAL HEAT-FLOW IS NON-ZERO, IT WILL ALWAYS BE LEFT C UNCHANGED, TO PRESERVE EFFECTS OF HAND-EDITING. C IF ANY HEAT-FLOW(S) IS/ARE ZERO, NEW VALUE(S) WILL BE C INTERPOLATED FROM A GRIDDED DATASET. C C HEAT-FLOW and lithosphere thickness C IN OCEANS IS BASED ON AGE OF SEAFLOOR, C WHERE KNOWN, FROM A GRIDDED-AGE DATASET. C C HEAT-FLOW (FROM ANY SOURCE) IS SUBJECT TO A MINIMUM VALUE. C C--------------------------------------------------------------------- INTERFACE SUBROUTINE ASSIGN (INPUT,AARRAY,AX1,ADX,AX2,NAX,ADY,AY2,NAY, + ALPHAT,CLIMIT,CONDUC, + CARRAY,CX1,CDX,CX2,NCX,CDY,CY2,NCY, + Delta_rho_limit, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT, + ONEKM, + PLON,PLAT, + QLIM0,DQL_DE,QLIM1, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + SARRAY,SX1,SDX,SX2,NSX,SDY,SY2,NSY, + TASTHK,TEMLIM,TSURF, + MODIFY,ELEVAT,HEATFL, + OUTPUT,THICKC,THICKM,chemical_Delta_rho, + cooling_curvature) INTEGER, INTENT(IN) :: INPUT,IUNITL,IUNITT,MODIFY,NAX,NAY, & NCX,NCY,NEX,NEY,NQX,NQY,NSX,NSY REAL, INTENT(IN) :: AX1,ADX,ADY,AX2,AY2,CLIMIT, & CX1,CDX,CX2,CDY,CY2, & Delta_rho_limit, & EX1,EDX,EX2,EDY,EY2, & GMEAN,HCMAX,HLMAX, & ONEKM,OUTPUT,PLON,PLAT, & QLIM0,DQL_DE,QLIM1,QX1,QDX,QX2,QDY,QY2, & RHOAST,RHOH2O, & SX1,SDX,SX2,SDY,SY2, & TASTHK,TSURF REAL, DIMENSION(2), INTENT(IN) :: ALPHAT,CONDUC,RADIO,RHOBAR, & TEMLIM REAL, DIMENSION(:,:), INTENT(IN) :: AARRAY,CARRAY,EARRAY, & QARRAY,SARRAY REAL, INTENT(INOUT) :: ELEVAT,HEATFL REAL, INTENT(OUT) :: THICKC,THICKM,chemical_Delta_rho, & cooling_curvature END SUBROUTINE ASSIGN END INTERFACE C--------------------------------------------------------------------- C PARAMETER (ARRAY-SIZE) STATEMENTS C C SET THE FOLLOWING PARAMETERS AT LEAST AS LARGE AS YOUR PROBLEM: C C MAXNOD = MAXIMUM NUMBER OF NODES C [N.B. The use of "fake" nodes as in program FAULTS has not been C supported in the Shells family; here all nodes are "real".] PARAMETER (MAXNOD=20000) C C MAXBN = MAXIMUM NUMBER OF BOUNDARY NODES C [N.B. The use of "fake" nodes as in program FAULTS has not been C supported in the Shells family; here all nodes are "real".] PARAMETER (MAXBN=6000) C C MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES). PARAMETER (MAXEL=30000) C C MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS); PARAMETER (MAXFEL=6000) C C MAXATP = MAXIMUM NUMBER OF NODES WHICH MAY OVERLAP AT A FAULT- C INTERSECTION POINT. PARAMETER (MAXATP=10) C C PARAMETER GIVING THE (EXACT) NUMBER OF PB2002 PLATES: PARAMETER (NUMPLT=52) C C--------------------------------------------------------------------- C TYPE STATEMENTS C C (NOTE: THE IMPLICIT TYPING OF I-N = INTEGER, AND A-H,O-Z = REAL C IS OBSERVED IN THIS PROGRAM. NO TYPES ARE STATED FOR SUCH NAMES.) C CHARACTER*1 C1 CHARACTER*2 NAMES CHARACTER*80 TITLE1,TITLE3 C LOGICAL BRIEF,DOTHIS,EVERYP,NEEDE,NEEDQ C REAL LRANGE,NLAT C C NOTE: THE FOLLOWING ARRAYS COULD BE COMPRESSED WITH "LOGICAL*1" C IN VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN,EDGETS,EDGEFS C C THE FOLLOWING TO AGREE WITH BLOCK DATA BD1: DOUBLE PRECISION POINTS,WEIGHT C THE FOLLOWING TO AGREE WITH BLOCK DATA BD2: DOUBLE PRECISION FPHI,FPOINT,FGAUSS C C--------------------------------------------------------------------- C DIMENSION STATMENTS C C DIMENSIONs using dynamic memory allocation: REAL, DIMENSION(:,:), ALLOCATABLE :: EARRAY,QARRAY,AARRAY, + CARRAY,SARRAY C C DIMENSIONS USING PARAMETER MAXNOD: DIMENSION CHECKN(MAXNOD), chemical_Delta_rho_list(MAXNOD), + cooling_curvature_list(MAXNOD), DQDTDA(MAXNOD), + ELEV (MAXNOD), + TLNODE(MAXNOD), + XNODE (MAXNOD), YNODE (MAXNOD), ZMNODE(MAXNOD) C C DIMENSIONS USING PARAMETER MAXBN: DIMENSION NODCON (MAXBN) C C DIMENSIONS USING PARAMETER MAXEL: DIMENSION AREA (MAXEL), CHECKE (MAXEL), + DETJ (7,MAXEL), + DXS(2,2,3,7,MAXEL), DYS (2,2,3,7,MAXEL), + DXSP (3,7,MAXEL), DYSP (3,7,MAXEL), + EDGETS (3,MAXEL), ELEDAT (3,MAXEL), + FPSFER(2,2,3,7,MAXEL), + NODES (3,MAXEL), + SITA (7,MAXEL) C C DIMENSIONS USING PARAMETER MAXFEL: DIMENSION CHECKF (MAXFEL), EDGEFS (2,MAXFEL), + FDIP (2,MAXFEL), + FLEN (MAXFEL), + FPFLT(2,2,2,7,MAXFEL), + FARG (2,MAXFEL), NODEF (4,MAXFEL), + OFFSET (MAXFEL) C C DIMENSIONS USING PARAMETER MAXATP: DIMENSION LIST (MAXATP) C C DIMENSIONS USING PARAMETER NUMPLT: DIMENSION NAMES (NUMPLT) C C FIXED DIMENSIONS: DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),RADIO(2), + RHOBAR(2),TAUMAX(2),TEMLIM(2) C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: DIMENSION POINTS(3,7),WEIGHT(7) C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: DIMENSION FPHI(4,7),FPOINT(7),FGAUSS(7) C C--------------------------------------------------------------------- C COMMON STATEMENTS C C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD1: COMMON /S1S2S3/ POINTS COMMON /WGTVEC/ WEIGHT C FOLLOWING STATEMENT TO AGREE WITH BLOCK DATA BD2: COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS C------------------------------------------------------------------- C DATA STATEMENTS C C "QLIM0" is the lower limit on heat-flow for points C with an elevation of zero. DATA QLIM0 /0.000/ C [N.B. Former limit, in program OrbData, was:} C DATA QLIM0 /0.045/ C C "DQL_DE" is the derivitive d(QLIMIT)/d(elevation), C which adjusts the minimum heat-flow for elevation. DATA DQL_DE /0.00E-06/ C C [N.B. Former limit, in program OrbData, was:} C DATA DQL_DE /1.43E-06/ C C "QLIM1 is the upper limit on heat-flow for all points. DATA QLIM1 /0.300/ C C "CLIMIT" IS A LOWER LIMIT ON CRUSTAL THICKNESS: DATA CLIMIT /6570./ C Changed from 5000. to 6570. on 2005.05.24, to agree with CRUST2.grd. C C "HCMAX" IS AN UPPER LIMIT ON CRUSTAL THICKNESS: DATA HCMAX /75000./ C Already agreed with CRUST2.grd; no need to change! C C "HLMAX" IS AN UPPER LIMIT ON TOTAL LITHOSPHERE THICKNESS: DATA HLMAX /400000./ C This limit prevents "lithosphere" from extending into the C transition zone, where new mineral phases appear, C and the input physical parameters would not be valid. C C "Delta_rho_limit" is the maximum permitted size of chemical C density anomalies throughout the lithosphere: DATA Delta_rho_limit /50./ C C NOTE THAT ALL OF THE ABOVE ARE IN SI UNITS (W/m**2, m, m, m, C kg/m**3, etc.) C C "IUNITT" = FORTRAN DEVICE NUMBER FOR TEXT MESSAGES: DATA IUNITT /6/ C C "IUNITL" = FORTRAN DEVICE NUMBER FOR -ASSIGN- LOG FILE: DATA IUNITL /13/ C C plate names (in alphabetical order): DATA names / 'AF','AM','AN', + 'AP','AR','AS', + 'AT','AU','BH', + 'BR','BS','BU', + 'CA','CL','CO', + 'CR','EA','EU', + 'FT','GP','IN', + 'JF','JZ','KE', + 'MA','MN','MO', + 'MS','NA','NB', + 'ND','NH','NI', + 'NZ','OK','ON', + 'PA','PM','PS', + 'RI','SA','SB', + 'SC','SL','SO', + 'SS','SU','SW', + 'TI','TO','WL', + 'YA' / C C "DIPMAX" IS THE MAXIMUM DIP (FROM HORIZONTAL, IN DEGREES) FOR A C FAULT ELEMENT TO BE TREATED AS A DIP-SLIP FAULT, WITH TWO DEGREES C OF FREEDOM PER NODE-PAIR. AT STEEPER DIPS, THE DEGREE OF FREEDOM C CORRESPONDING TO OPENING OR CONVERGENCE OF THE OPPOSITE SIDES IS C ELIMINATED BY A CONSTRAINT EQUATION, AND THE FAULT IS TREATED AS C A VERTICAL STRIKE-SLIP FAULT. THIS ARBITRARY LIMIT IS NECESSARY C BECAUSE THE EQUATIONS FOR DIP-SLIP FAULTS BECOME SINGULAR AS THE C DIP APPROACHES 90 DEGREES. IN PRACTICE, IT IS BEST TO SPECIFY DIPS C AS EITHER (1) VERTICAL, OR (2) CLEARLY LESS THAN "DIPMAX", WITHIN C EACH FAULT ELEMENT. IF THE DIP VARIES WITHIN AN ELEMENT IN SUCH A C WAY THAT IT PASSES THROUGH THIS LIMIT WITHIN THE ELEMENT, THEN C THE REPRESENTATION OF THAT FAULT ELEMENT IN THE EQUATIONS MAY C BE INACCURATE. DATA DIPMAX /75./ C--------------------------------------------------------------------- C C STATEMENT FUNCTIONS: EASTOF(X)=MOD((X + 720.),360.) C C--------------------------------------------------------------------- C C BEGINNING OF EXECUTABLE CODE C C *** KLUDGE ALERT *** C CONVERSION OF PARAMETERS (CONSTANTS) TO VARIABLES SHOULD LOGICALLY C HAVE NO EFFECT, BUT IN FACT HELPS TO SUPPRESS SOME SPURIOUS C MESSAGES FROM THE IBM VS-FORTRAN COMPILER. MXNODE=MAXNOD MXBN =MAXBN MXEL =MAXEL MXFEL =MAXFEL MXSTAR=MAXATP C WEDGE=ABS(90.-ABS(DIPMAX))*0.017453293 C--------------------------------------------------------------------- C C INTRODUCTION / HEADER SECTION: C WRITE (IUNITT,1) 1 FORMAT ( + ' -------------------------------------------------------------' + /' Program OrbData5 (version of 24 August 2006)' +//' by Peter Bird, Dept. of Earth & Space Sciences,' + /' University of California, Los Angeles, CA 90095-1567.' +//' Reads a finite element grid produced by OrbWeave or OrbWin' + /' (run in Shells mode)' + /' and a parameter file (formatted for Shells) with' + /' thermal conductivities, densities, et cetera;' + /' Then, fills in or computes nodal data:' + /' -elevation (+ above sea level, - below);' + /' -heat-flow;' + /' -crustal thickness;' + /' -mantle lithosphere thickness (NOT including crust);' + /' -chemical density anomaly (of whole lithosphere);' + /' -geotherm curvature due to cooling (whole lithosphere).') CALL Pause() C WRITE (IUNITT, 2) 2 FORMAT( +//' All units are SI: m, W/m**2 (NOT mW/m**2), m, m, kg/m**3,' + /' and C/m**2.' + /' Note that the change in the quadratic coefficient of the' + /' geotherm caused by cooling_curvature is:' + /' Delta_GEOTH3 = Delta_GEOTH7 = -0.5 * cooling_curvature.' +//' In the output data, the base of the mantle lithosphere' + /' will be an isothermal surface. The isotherm value' + /' is chosen to lie on the asthenosphere adiabat' + /' (evaluated at an arbitrary depth of 100 km).') CALL Pause() C WRITE (IUNITT,3) QLIM0, DQL_DE, QLIM1 3 FORMAT ( + /' If an elevation is non-zero, this elevation will be left' + /' unchanged, to preserve effects of hand-editing.' + /' If any elevation is zero, a new values will be' + /' interpolated from a gridded dataset, such as ETOPO5.grd.' + /' If a nodal heat-flow is non-zero, this heat-flow will be left' + /' unchanged, to preserve effects of hand-editing.' + /' If any heat-flow is zero, a new values will be' + /' interpolated from a gridded dataset.' + /' Ocean-floor heat-flow will be based on age of seafloor,' + /' where known (< 200 Ma), from a gridded-age dataset,' + /' such as age_1p5.grd derived from Mueller et al. [1997a].' +//' Heat-flow (from any source) is subject to a minimum value of' + /' (',F5.3,' + ',ES10.3,' * elevation)' + /' and a maximum value of ',F5.3,'.' + /' ------------------------------------------------------------') CALL Pause() C C INITIATE LOG FILE FOR -ASSIGN- OUTPUT: C WRITE (IUNITT,91) IUNITL 91 FORMAT (/ /' Attempting to create detailed log file on unit', + I3/) WRITE (IUNITL, "('Log file of a run of OrbData5:')") C C Echo the limits compiled-in-place: C WRITE (IUNITT,5) QLIM0,DQL_DE, QLIM1, CLIMIT, HCMAX, HLMAX, + Delta_rho_limit WRITE (IUNITL,5) QLIM0,DQL_DE, QLIM1, CLIMIT, HCMAX, HLMAX, + Delta_rho_limit 5 FORMAT(/' The following limits apply in this run:' +/' Lower limit on heat-flow = ',F5.3,'+',ES10.3,' * elevation' +/' Upper limit on heat-flow = ',F5.3 +/' Lower limit on crustal thickness = ',F7.0 +/' Upper limit on crustal thickness = ',F7.0 +/' Upper limit on total lithosphere thickness = ',F7.0 +/' Upper limit of chemical density variation = ',F7.0) C WRITE (IUNITL,92) 92 FORMAT (/'Table of results, possibly interrupted by' + /'warning messages from subprogram -ASSIGN-:') C C READ PARAMETERS ON UNIT 1 C CALL READPM (INPUT, 1, 6, NAMES, NUMPLT, OFFMAX, + OUTPUT,ACREEP, ALPHAT, BCREEP, BIOT , + BYERLY, CCREEP, CFRIC , CONDUC, + DCREEP, ECREEP, EVERYP, FFRIC , GMEAN , + GRADIE, ICONVE, IPVREF, + MAXITR, OKDELV, OKTOQT, ONEKM, RADIO , + RDETH, REFSTR, RHOAST, RHOBAR, RHOH2O, + TADIAB, + TAUMAX, TEMLIM, TITLE3, TRHMAX, TSURF, + ZBASTH) C C LITHOSPHERE/ASTHENOSPHERE TEMPERATURE, IN KELVIN, IS DETERMINED C FROM THE ASTHENOSPHERE ADIABAT IN THE PARAMETER FILE, C EVALUATED (ARBITRARILY) AT 100 KM: C TASTHK = TADIAB + GRADIE * 100.E3 C C READ FINITE ELEMENT GRID ON UNIT 2 C CALL GETNET (INPUT, 2,IUNITT, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + NFAKEN,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,OFFMAX,OFFSET, + TITLE1,XNODE,YNODE, + WORK,CHECKE,CHECKF,CHECKN) C [N.B. This version ignores any crustal thicknesses, C mantle lithosphere thicknesses, chemical density anomalies, C or cooling curvatures that might be in the .feg file.] WRITE (IUNITT,40) C 40 FORMAT (/' Successfully read FE grid, now verifying topology...') CALL SQUARE (INPUT,BRIEF,FDIP,IUNITT, + MXBN,MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES, + NUMEL,NUMNOD,RDETH,WEDGE, + MODIFY,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,EDGEFS, + EDGETS,FLEN,FPFLT,FPSFER, + FARG,NCOND,NODCON,SITA, + WORK,CHECKN,LIST) WRITE (IUNITT,50) 50 FORMAT (/' Topology of finite element grid verified.'/) C C READ IN ELEVATION ARRAY ON UNIT 3, IF NEEDED C NEEDE=.FALSE. DO 60 I=1,NUMNOD IF (ELEV(I).EQ.0.0) NEEDE=.TRUE. 60 CONTINUE IF (NEEDE) THEN WRITE(*,61) 61 FORMAT(/ /' Attempting to read gridded elevations:'/) READ (3,*) EX1,EDX,EX2 READ (3,*) EY1,EDY,EY2 NEX=(EX2-EX1)/EDX + 1.5 NEY=(EY2-EY1)/EDY + 1.5 ALLOCATE ( EARRAY(NEY, NEX) ) READ (3,*) ((EARRAY(IROW,JCOL),JCOL=1,NEX),IROW=1,NEY) ELSE WRITE (*,69) 69 FORMAT (/' All nodes have non-zero elevation.' & /' No elevation grid is needed.') ENDIF C C READ IN HEAT-FLOW ARRAY ON UNIT 4, IF NEEDED C NEEDQ=.FALSE. DO 70 I=1,NUMNOD IF (DQDTDA(I).EQ.0.0) NEEDQ=.TRUE. 70 CONTINUE C IF (NEEDQ) THEN WRITE(*,71) 71 FORMAT(/ /' Attempting to read gridded heat-flow:'/) READ (4,*) QX1,QDX,QX2 READ (4,*) QY1,QDY,QY2 NQX=(QX2-QX1)/QDX + 1.5 NQY=(QY2-QY1)/QDY + 1.5 ALLOCATE ( QARRAY(NQY, NQX) ) READ (4,*) ((QARRAY(IROW,JCOL),JCOL=1,NQX),IROW=1,NQY) ELSE WRITE (*,79) 79 FORMAT (/' All nodes have non-zero heat-flow.' & /' No heat-flow grid is needed.') C C IMPOSE HEAT-FLOW LIMITS AT NODES, TO PREVENT UNREASONABLY C THICK AND STIFF LITHOSPHERE ANYWHERE C DO 80 I=1,NUMNOD QLIMIT=QLIM0+DQL_DE*ELEV(I) IF (DQDTDA(I).NE.0.0) DQDTDA(I)=MAX(DQDTDA(I),QLIMIT) DQDTDA(I)=MIN(DQDTDA(I),QLIM1) 80 CONTINUE END IF C C READ DATASET OF GRIDDED AGES, ON UNIT 7: C (NOTE: AGE > 200 Ma MEANS "UNKNOWN" OR "CONTINENTAL") WRITE(*,73) 73 FORMAT(/ /' Attempting to read gridded ages of seafloor:'/) READ (7,*) AX1,ADX,AX2 READ (7,*) AY1,ADY,AY2 NAX=(AX2-AX1)/ADX + 1.5 NAY=(AY2-AY1)/ADY + 1.5 ALLOCATE ( AARRAY(NAY, NAX) ) READ (7,*) ((AARRAY(IROW,JCOL),JCOL=1,NAX),IROW=1,NAY) C C Read in array of crustal thickness, such as CRUST2.grd, on unit 11: C WRITE(*,81) 81 FORMAT(/ /' Attempting to read gridded crustal thickness:'/) READ (11,*) CX1,CDX,CX2 READ (11,*) CY1,CDY,CY2 NCX=(CX2-CX1)/CDX + 1.5 NCY=(CY2-CY1)/CDY + 1.5 ALLOCATE ( CARRAY(NCY, NCX) ) READ (11,*) ((CARRAY(IROW,JCOL),JCOL=1,NCX),IROW=1,NCY) C C Must convert CRUST2.grd from km to m units: DO IROW = 1, NCY DO JCOL = 1, NCX CARRAY(IROW,JCOL) = CARRAY(IROW,JCOL) * 1000. END DO END DO C C Read in array of vertical S wave travel time anomalies C (from the Moho to 400 km depth), C such as delta_ts.grd, on unit 12: C WRITE(*,82) 82 FORMAT(/ /' Attempting to read gridded S-wave travel-time' & /' anomalies (in the upper mantle):'/) READ (12,*) SX1,SDX,SX2 READ (12,*) SY1,SDY,SY2 NSX=(SX2-SX1)/SDX + 1.5 NSY=(SY2-SY1)/SDY + 1.5 ALLOCATE ( SARRAY(NSY, NSX) ) READ (12,*) ((SARRAY(IROW,JCOL),JCOL=1,NSX),IROW=1,NSY) C C PROCESS ALL NODES EQUALLY: C WRITE (IUNITT,600) 600 FORMAT (/' Computing layer thicknesses at all nodes...' + /' 0 nodes completed...') WRITE (IUNITL,601) 601 FORMAT (' NODE LONGITUDE LATITUDE ELEV DQDTDA', & ' ZMNODE TLNODE', & ' chemical_Delta_rho cooling_curvature') DO 680 I=1,NUMNOD PLON=YNODE(I)*57.2957795 PLAT=90.-XNODE(I)*57.2957795 ELEVAT=ELEV(I) HEATFL=DQDTDA(I) CALL ASSIGN (INPUT, + AARRAY,AX1,ADX,AX2,NAX,ADY,AY2,NAY, + ALPHAT,CLIMIT,CONDUC, + CARRAY,CX1,CDX,CX2,NCX,CDY,CY2,NCY, + Delta_rho_limit, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT, + ONEKM, + PLON,PLAT, + QLIM0,DQL_DE,QLIM1, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + SARRAY,SX1,SDX,SX2,NSX,SDY,SY2,NSY, + TASTHK,TEMLIM,TSURF, + MODIFY,ELEVAT,HEATFL, + OUTPUT,THICKC,THICKM,chemical_Delta_rho, + cooling_curvature) ELEV(I)=ELEVAT DQDTDA(I)=HEATFL ZMNODE(I)=MAX(THICKC,CLIMIT) TLNODE(I)=MAX(THICKM,0.) chemical_Delta_rho_list(I) = chemical_Delta_rho cooling_curvature_list(I) = cooling_curvature WRITE (IUNITL,678) I,PLON,PLAT, + ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I), + chemical_Delta_rho_list(I), + cooling_curvature_list(I) 678 FORMAT (' ',I5,2F10.4,4ES10.2,ES19.2,ES18.2) WRITE (IUNITT,679) I 679 FORMAT ('+',I6,' nodes completed...') 680 CONTINUE WRITE (IUNITT,700) 700 FORMAT (/ /' CHECK THE LOG FILE CAREFULLY FOR PROBLEMS!') C C OUTPUT THE MODIFIED .FEG FILE: C CALL PUTNET (INPUT, 14, + BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE, + chemical_Delta_rho_list, + cooling_curvature_list) STOP END C C C REAL FUNCTION ATAN2F (Y,X) C C CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS C IF ((Y.NE.0.).OR.(X.NE.0.)) THEN ATAN2F=ATAN2(Y,X) ELSE ATAN2F=0. ENDIF RETURN END C C C SUBROUTINE READN (INPUT,IUNITP,IUNITT,N, + OUTPUT,VECTOR) C C A UTILITY ROUTINE DESIGNED TO PERMIT -FAULTS- INPUT FILES C TO ALSO BE USED BY -PLATES-, WHICH EXPECTS MORE NUMBERS C IN SOME RECORDS. C THIS ROUTINE ATTEMPTS TO READ 'N' FLOATING-POINT VALUES C (USING * FORMAT) FROM THE NEXT RECORD ON DEVICE 'IUNITP'. C IF ANYTHING GOES WRONG, THE MISSING VALUES ARE SET TO ZERO. C CHARACTER*1 C CHARACTER*132 LINE LOGICAL ANYIN,DOTTED,EXPON,SIGNED DIMENSION VECTOR(N) C READ (IUNITP,1) LINE 1 FORMAT (A) C NUMBER=0 ANYIN=.FALSE. EXPON=.FALSE. SIGNED=.FALSE. DOTTED=.FALSE. DO 10 I=1,132 C=LINE(I:I) IF ((C.EQ.' ').OR.(C.EQ.',').OR.(C.EQ.'/')) THEN SIGNED=.FALSE. EXPON=.FALSE. DOTTED=.FALSE. IF (ANYIN) THEN NUMBER=NUMBER+1 ANYIN=.FALSE. ENDIF ELSE IF ((C.EQ.'+').OR.(C.EQ.'-')) THEN IF (SIGNED) THEN GO TO 50 ELSE SIGNED=.TRUE. ENDIF ELSE IF ((C.EQ.'E').OR.(C.EQ.'D').OR. + (C.EQ.'e').OR.(C.EQ.'d')) THEN IF (EXPON) THEN GO TO 50 ELSE EXPON=.TRUE. SIGNED=.FALSE. DOTTED=.TRUE. ENDIF ELSE IF (C.EQ.'.') THEN IF (DOTTED) THEN GO TO 50 ELSE DOTTED=.TRUE. ENDIF ELSE IF ((C.EQ.'0').OR.(C.EQ.'1').OR.(C.EQ.'2').OR. + (C.EQ.'3').OR.(C.EQ.'4').OR.(C.EQ.'5').OR. + (C.EQ.'6').OR.(C.EQ.'7').OR.(C.EQ.'8').OR. + (C.EQ.'9')) THEN SIGNED=.TRUE. ANYIN=.TRUE. ELSE GO TO 50 ENDIF 10 CONTINUE IF (ANYIN) NUMBER=NUMBER+1 C 50 IF (NUMBER.EQ.0) THEN WRITE (IUNITT,91) N,LINE 91 FORMAT (/' ERR0R: A LINE OF PARAMETER INPUT WHICH', + ' WAS SUPPOSED TO CONTAIN 1-',I2,' NUMBERS'/ + ' COULD NOT BE INTERPRETED. LINE FOLLOWS:'/ + ' ',A80) CALL PAUSE() STOP ELSE NUMBER=MIN(NUMBER,N) BACKSPACE IUNITP READ (IUNITP,*) (VECTOR(I),I=1,NUMBER) IF (NUMBER.LT.N) THEN DO 99 I=NUMBER+1,N VECTOR(I)=0. 99 CONTINUE ENDIF ENDIF RETURN END C C C SUBROUTINE SQUEEZ (INPUT,ALPHAT,density_anomaly_kgpm3,ELEVAT, + GEOTH1,GEOTH2,GEOTH3,GEOTH4, + GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN, + IUNITT,ONEKM,RHOAST,RHOBAR,RHOH2O, + TEMLIM,ZM,ZSTOP, + OUTPUT,TAUZZ,SIGZZB) C C CALCULATES "TAUZZ", THE VERTICAL INTEGRAL THROUGH THE PLATE C OF THE VERTICAL STANDARDIZED STRESS ANOMALY, WHICH IS C RELATIVE TO A COLUMN OF MANTLE AT ASTHENOSPHERE TEMPERATURE C WITH A 5 KM CRUST AND A 2.7 KM OCEAN ON TOP, LIKE A MID-OCEAN C RISE. THE INTEGRAL IS FROM EITHER THE LAND SURFACE OR THE C SEA SURFACE, DOWN TO A DEPTH OF "ZSTOP" BELOW THE TOP OF C THE CRUST. C IF "ZSTOP" EXCEEDS MOHO DEPTH "ZM", THEN PROPERTIES OF THE MANTLE C WILL BE USED IN THE LOWER PART OF THE INTEGRAL. C ALSO RETURNS "SIGZZB", THE STANDARDIZED VERTICAL STRESS ANOMALY C AT DEPTH "ZSTOP" BELOW THE SOLID ROCK SURFACE. C NOTE: THIS VERSION IS DIFFERENT FROM THE VERSION FOUND IN THE LARAMY C PROGRAM PACKAGE. FIRST, IT ACTS ON ONLY A SINGLE POINT. C SECOND, IT INFERS SUB-PLATE NORMAL-STRESS ANOMALIES FROM C THE GIVEN TOPOGRAPHY, INSTEAD OF FROM MODEL STRUCTURE. C PARAMETER (NDREF=300) LOGICAL CALLED C INTERNAL ARRAYS: DIMENSION DREF(NDREF),PREF(0:NDREF) C ARGUMENT ARRAYS: DIMENSION ALPHAT(2),RHOBAR(2),TEMLIM(2) SAVE CALLED,DREF,PREF DATA CALLED /.FALSE./ C C STATEMENT FUNCTIONS: TEMPC(H)=MIN(TEMLIM(1),GEOTH1+GEOTH2*H+GEOTH3*H**2 + +GEOTH4*H**3) TEMPM(H)=MIN(TEMLIM(2),GEOTH5+GEOTH6*H+GEOTH7*H**2 + +GEOTH8*H**3) C C CREATE REFERENCE TEMPERATURE & DENSITY PROFILES TO DEPTH OF NDREF KM C IF (.NOT.CALLED) THEN RHOTOP=RHOBAR(1)*(1.-ALPHAT(1)*GEOTH1) DREF(1)=RHOH2O DREF(2)=RHOH2O DREF(3)=0.7*RHOH2O+0.3*RHOTOP DREF(4)=RHOTOP DREF(5)=RHOTOP DREF(6)=RHOTOP DREF(7)=RHOTOP DREF(8)=0.7*RHOTOP+0.3*RHOAST DO 50 J=9,NDREF DREF(J)=RHOAST 50 CONTINUE PREF(0)=0. DO 100 I=1,NDREF PREF(I)=PREF(I-1)+DREF(I)*GMEAN*ONEKM 100 CONTINUE ENDIF C C ROUTINE PROCESSING (ON EVERY CALL): C IF (ELEVAT.GT.0.) THEN C LAND ZTOP= -ELEVAT ZBASE=ZSTOP-ELEVAT DENSE1=RHOBAR(1)*(1.-GEOTH1*ALPHAT(1))+ + density_anomaly_kgpm3 H=0. LAYER1=1 ELSE C OCEAN ZTOP=0. ZBASE=ZSTOP+(-ELEVAT) DENSE1=RHOH2O H=ELEVAT LAYER1=0 ENDIF LASTDR=ZBASE/ONEKM IF (ZBASE.GT.ONEKM*LASTDR) LASTDR=LASTDR+1 IF (LASTDR.GT.NDREF) THEN WRITE(IUNITT,110) LASTDR 110 FORMAT(' IN SUBPROGRAM SQUEEZ, PARAMETER NDREF '/ + ' MUST BE INCREASED TO AT LEAST ',I10) CALL PAUSE() STOP ENDIF NSTEP=(ZBASE-ZTOP)/ONEKM OLDSZZ=0. OLDPR=0. SIGZZ=0. TAUZZ=0. Z=ZTOP DO 200 I=1,NSTEP Z=Z+ONEKM H=H+ONEKM IF (H.GT.0.) THEN IF (H.LE.ZM) THEN T=TEMPC(H) DENSE2=RHOBAR(1)*(1.-T*ALPHAT(1))+ + density_anomaly_kgpm3 LAYER2=1 ELSE T=TEMPM(H-ZM) DENSE2=RHOBAR(2)*(1.-T*ALPHAT(2))+ + density_anomaly_kgpm3 LAYER2=2 ENDIF ELSE DENSE2=RHOH2O LAYER2=0 ENDIF IF ((LAYER1.EQ.0).AND.(LAYER2.EQ.1)) THEN FRAC2=H/ONEKM FRAC1=1.-FRAC2 ELSE IF ((LAYER1.EQ.1).AND.(LAYER2.EQ.2)) THEN FRAC2=(H-ZM)/ONEKM FRAC1=1.-FRAC2 ELSE FRAC1=0.5 FRAC2=0.5 ENDIF DENSE=FRAC1*DENSE1+FRAC2*DENSE2 IF (Z.GT.0.) THEN N1=Z/ONEKM N2=N1+1 FRAC=Z/ONEKM-N1 PR=PREF(N1)+FRAC*(PREF(N2)-PREF(N1)) ELSE PR=0. ENDIF SIGZZ=SIGZZ-DENSE*GMEAN*ONEKM+(PR-OLDPR) TAUZZ=TAUZZ+0.5*(SIGZZ+OLDSZZ)*ONEKM DENSE1=DENSE2 OLDSZZ=SIGZZ OLDPR=PR LAYER1=LAYER2 200 CONTINUE RESID=ZBASE-Z H=ZSTOP Z=ZBASE IF (ZSTOP.LE.ZM) THEN T=TEMPC(H) DENSE2=RHOBAR(1)*(1.-T*ALPHAT(1))+ + density_anomaly_kgpm3 ELSE T=TEMPM(H-ZM) DENSE2=RHOBAR(2)*(1.-T*ALPHAT(2))+ + density_anomaly_kgpm3 ENDIF DENSE=0.5*(DENSE1+DENSE2) IF (Z.GT.0.) THEN N1=Z/ONEKM N2=N1+1 FRAC=Z/ONEKM-N1 PR=PREF(N1)+FRAC*(PREF(N2)-PREF(N1)) ELSE PR=0. ENDIF SIGZZB=SIGZZ-DENSE*GMEAN*RESID+(PR-OLDPR) TAUZZ=TAUZZ+0.5*(SIGZZB+OLDSZZ)*RESID CALLED=.TRUE. RETURN END C C C SUBROUTINE ASSIGN (INPUT,AARRAY,AX1,ADX,AX2,NAX,ADY,AY2,NAY, + ALPHAT,CLIMIT,CONDUC, + CARRAY,CX1,CDX,CX2,NCX,CDY,CY2,NCY, + Delta_rho_limit, + EARRAY,EX1,EDX,EX2,NEX,EDY,EY2,NEY, + GMEAN,HCMAX,HLMAX, + IUNITL,IUNITT, + ONEKM, + PLON,PLAT, + QLIM0,DQL_DE,QLIM1, + QARRAY,QX1,QDX,QX2,NQX,QDY,QY2,NQY, + RADIO,RHOAST,RHOBAR,RHOH2O, + SARRAY,SX1,SDX,SX2,NSX,SDY,SY2,NSY, + TASTHK,TEMLIM,TSURF, + MODIFY,ELEVAT,HEATFL, + OUTPUT,THICKC,THICKM,chemical_Delta_rho, + cooling_curvature) C C DETERMINES ELEVATION (IF 0.0), HEAT FLOW (IF 0.0), AND C VALUES OF CRUSTAL AND MANTLE LITHOSPHERE THICKNESS, AT ONE POINT: C (PLON = EAST LONGITUDE, PLAT = NORTH LATITUDE), BOTH IN DEGREES. C (UNLIKE SUBPROGRAM -ASSIGN- IN ALGRIDDR, THIS VERSION NEVER C DOES REGIONAL AVERAGING. NEITHER DOES IT INTERPOLATE BY C KRIGING; SIMPLE INTERPOLATION WITHIN RECTANGULAR GRIDS ARE C USED INSTEAD.) C C This version of ASSIGN was completely re-written for OrbData5. C Instead of solving for crustal thickness from isostasy, it C reads crustal thickness from a grid of values. C Instead of computing lithosphere thickness from heat flow, C it computes it from seafloor age or from vertical S-wave C travel-time anomaly under continents. Then extra degrees of C freedom "chemical_Delta_rho" and "cooling_curvature" are C used to create a self-consistent structure. C C Note that the result may *FAIL* to be isostatic with mid-ocean C ridges if that would require chemical_Delta_rho to go beyond C the prescribed amplitude. Also, if the requested lithosphere C thickness is more than ~twice that expected based on heat-flow C (and a steady-state assumption), then the lithosphere thickness C returned will be limited so as to prevent geotherm overshoot C and negative temperature gradients. C C Revision by Peter Bird, UCLA, June 2005 C INTEGER, INTENT(IN) :: INPUT,IUNITL,IUNITT,MODIFY,NAX,NAY, & NCX,NCY,NEX,NEY,NQX,NQY,NSX,NSY REAL, INTENT(IN) :: AX1,ADX,ADY,AX2,AY2, & CLIMIT, & CX1,CDX,CX2,CDY,CY2, & Delta_rho_limit, & EX1,EDX,EX2,EDY,EY2, & GMEAN,HCMAX,HLMAX, & ONEKM,OUTPUT,PLON,PLAT, & QLIM0,DQL_DE,QLIM1,QX1,QDX,QX2,QDY,QY2, & RHOAST,RHOH2O, & SX1,SDX,SX2,SDY,SY2, & TASTHK,TSURF REAL, DIMENSION(2), INTENT(IN) :: ALPHAT,CONDUC,RADIO,RHOBAR, & TEMLIM REAL, DIMENSION(:,:), INTENT(IN) :: AARRAY,CARRAY,EARRAY, & QARRAY,SARRAY REAL, INTENT(INOUT) :: ELEVAT,HEATFL REAL, INTENT(OUT) :: THICKC,THICKM,chemical_Delta_rho, & cooling_curvature C LOGICAL BADP,BADT,DIVERG,OUTSID,NEEDE,NEEDQ, + WARNC1, WARNC2, WARNM1, WARNL2, WAYOUT REAL m_star C--------------------------------------------------------------------- OUTSID=.FALSE. WAYOUT=.FALSE. WARNC1=.FALSE. WARNC2=.FALSE. WARNM1=.FALSE. WARNL2=.FALSE. C C DETERMINE THE ELEVATION, IF NEEDED C NEEDE=(ELEVAT.EQ.0.0) IF (NEEDE) THEN IR1=((EY2-PLAT)/EDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NEY-1) IR2=IR1+1 FR=((EY2-EDY*(IR1-1))-PLAT)/EDY C IF ((PLON.LT.EX1).OR.(PLON.GT.EX2)) THEN IF (((PLON+360.).GE.EX1).AND. + ((PLON+360.).LE.EX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.EX1).AND. + ((PLON-360.).LE.EX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-EX1)/EDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NEX-1) IC2=IC1+1 FC=(USELON-(EX1+EDX*(IC1-1)))/EDX C OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=EARRAY(IR1,IC1)+FC*(EARRAY(IR1,IC2)-EARRAY(IR1,IC1)) BOT=EARRAY(IR2,IC1)+FC*(EARRAY(IR2,IC2)-EARRAY(IR2,IC1)) ELEVAT=TOP+FR*(BOT-TOP) ENDIF C C Determine age of sea-floor, which is always needed: C IR1=((AY2-PLAT)/ADY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NAY-1) IR2=IR1+1 FR=((AY2-ADY*(IR1-1))-PLAT)/ADY C IF ((PLON.LT.AX1).OR.(PLON.GT.AX2)) THEN IF (((PLON+360.).GE.AX1).AND. + ((PLON+360.).LE.AX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.AX1).AND. + ((PLON-360.).LE.AX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-AX1)/ADX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NAX-1) IC2=IC1+1 FC=(USELON-(AX1+ADX*(IC1-1)))/ADX C OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=AARRAY(IR1,IC1)+FC*(AARRAY(IR1,IC2)-AARRAY(IR1,IC1)) BOT=AARRAY(IR2,IC1)+FC*(AARRAY(IR2,IC2)-AARRAY(IR2,IC1)) AGEMA=TOP+FR*(BOT-TOP) C C DETERMINE HEAT FLOW, IF NEEDED C NEEDQ=(HEATFL.EQ.0.0) IF (NEEDQ) THEN IR1=((QY2-PLAT)/QDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NQY-1) IR2=IR1+1 FR=((QY2-QDY*(IR1-1))-PLAT)/QDY C IF ((PLON.LT.QX1).OR.(PLON.GT.QX2)) THEN IF (((PLON+360.).GE.QX1).AND. + ((PLON+360.).LE.QX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.QX1).AND. + ((PLON-360.).LE.QX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-QX1)/QDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NQX-1) IC2=IC1+1 FC=(USELON-(QX1+QDX*(IC1-1)))/QDX C OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=QARRAY(IR1,IC1)+FC*(QARRAY(IR1,IC2)-QARRAY(IR1,IC1)) BOT=QARRAY(IR2,IC1)+FC*(QARRAY(IR2,IC2)-QARRAY(IR2,IC1)) HEATFL=TOP+FR*(BOT-TOP) C C CHECK FOR SEAFLOOR AGE OVERRIDING HEAT-FLOW GRID VALUE: C IF (AGEMA.LT.200.0) THEN C SEAFLOOR AGE IS VALID (NOT "UNKNOWN" OR "CONTINENTAL") IF (AGEMA.LE.0.0) THEN HEATFL=QLIM1 ELSE C Carol A. Stein & Seth Stein [1992] C A model for the global variation in oceanic C depth and heat flow with lithospheric age, C Nature, v. 359, 10 September, p. 123-129. C According to their preferred GDH1 model: IF (AGEMA.LE.55.0) THEN HEATFL=0.510/SQRT(AGEMA) ELSE HEATFL=0.048+0.096*EXP(-0.0278*AGEMA) END IF END IF END IF C END IF C C APPLY LIMITS ON Q: C QLIMIT=QLIM0+DQL_DE*ELEVAT HEATFL=MAX(HEATFL,QLIMIT) HEATFL=MIN(HEATFL,QLIM1) C C Obtain crustal thickness from grid CARRAY: C IR1=((CY2-PLAT)/CDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NCY-1) IR2=IR1+1 FR=((CY2-CDY*(IR1-1))-PLAT)/CDY IF ((PLON.LT.CX1).OR.(PLON.GT.CX2)) THEN IF (((PLON+360.).GE.CX1).AND. + ((PLON+360.).LE.CX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.CX1).AND. + ((PLON-360.).LE.CX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-CX1)/CDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NCX-1) IC2=IC1+1 FC=(USELON-(CX1+CDX*(IC1-1)))/CDX OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=CARRAY(IR1,IC1)+FC*(CARRAY(IR1,IC2)-CARRAY(IR1,IC1)) BOT=CARRAY(IR2,IC1)+FC*(CARRAY(IR2,IC2)-CARRAY(IR2,IC1)) THICKC=TOP+FR*(BOT-TOP) THICKC=MAX(THICKC,CLIMIT) THICKC=MIN(THICKC,HCMAX) C C GIVE WARNING IF ANY NECESSARY DATA COULD NOT BE READ C IF (OUTSID) THEN WRITE (IUNITL,1) PLON,PLAT 1 FORMAT(' WARNING: (',F8.3,'E, ',F7.3,'N) IS ', + 'SLIGHTLY OUTSIDE DATA GRID(S).') END IF IF (WAYOUT) THEN WRITE (IUNITL,2) PLON,PLAT 2 FORMAT(' ERROR: (',F8.3,'E, ',F7.3,'N) IS ', + 'FAR OUTSIDE DATA GRID(S).') CCCCC WRITE (IUNITT,2) PLON,PLAT CCCCC CALL Pause() CCCCC STOP END IF C C Method of determining lithosphere thickness depends on whether C node has known (and valid) sea-floor age: C IF (AGEMA.LT.200.0) THEN C C This node is in oceanic lithosphere. C C (1)Determine age-dependent model Q: C C Carol A. Stein & Seth Stein [1992] C A model for the global variation in oceanic C depth and heat flow with lithospheric age, C Nature, v. 359, 10 September, p. 123-129. C According to their preferred GDH1 model: IF (AGEMA.LE.55.0) THEN Q_GDH1 = 0.510 / SQRT(AGEMA) ELSE Q_GDH1 = 0.048 + 0.096 * EXP(-0.0278 * AGEMA) END IF C C (2)Find total lithosphere thickness expected according to C a steady-state conduction model (like old OrbData; see C spreadsheet S20RTS_delta_ts_vs_age.xls for calibration). C Q_GDH1 = MAX(Q_GDH1, QLIMIT) Q_GDH1 = MIN(Q_GDH1, QLIM1) Q_radioactivity = 0.007 Delta_T = TASTHK - TSURF h_Earth3 = (Delta_T * CONDUC(2))/(Q_GDH1 - Q_radioactivity) C C (3)Take geometric mean of this "Earth3" or "old OrbData" C thickness with the constant plate thickness of Stein & C Stein [1992]: C h_plate = 95000. h_Earth5 = SQRT(h_Earth3 * h_plate) C C (4)Mantle lithosphere excludes crust: C THICKM = h_Earth5 - THICKC C C (5)Apply lower and upper limits: C THICKM = MAX(THICKM, 0.) THICKM = MIN(THICKM, (HLMAX - THICKC)) C ELSE C C This node is in continental (or unknown) lithosphere. C C (1)Read vertical S-wave travel-time anomaly (in the upper C mantle, from Moho to 400 km depth) from array SARRAY: C IR1=((SY2-PLAT)/SDY)+1.00001 IR1=MAX(IR1,1) IR1=MIN(IR1,NSY-1) IR2=IR1+1 FR=((SY2-SDY*(IR1-1))-PLAT)/SDY IF ((PLON.LT.SX1).OR.(PLON.GT.SX2)) THEN IF (((PLON+360.).GE.SX1).AND. + ((PLON+360.).LE.SX2)) THEN USELON=PLON+360.0 ELSE IF (((PLON-360.).GE.SX1).AND. + ((PLON-360.).LE.SX2)) THEN USELON=PLON-360.0 ELSE USELON=PLON END IF ELSE USELON=PLON END IF IC1=((USELON-SX1)/SDX)+1.00001 IC1=MAX(IC1,1) IC1=MIN(IC1,NSX-1) IC2=IC1+1 FC=(USELON-(SX1+SDX*(IC1-1)))/SDX OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR. + (FC.LT.-0.01).OR.(FC.GT.1.01) WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR. + (FC.LT.-1.01).OR.(FC.GT.2.01) FR=MIN(1.,MAX(0.,FR)) FC=MIN(1.,MAX(0.,FC)) TOP=SARRAY(IR1,IC1)+FC*(SARRAY(IR1,IC2)-SARRAY(IR1,IC1)) BOT=SARRAY(IR2,IC1)+FC*(SARRAY(IR2,IC2)-SARRAY(IR2,IC1)) Delta_ts = TOP + FR * (BOT - TOP) C C (2)Scale to mantle lithosphere thickness, based on C calibration performed in oceanic lithosphere of known age C (see diary.doc entry for 2005.06.02): C slope = (1.494) / (85000.) C (in units of seconds per meter of lithosphere) ridge_Delta_ts = 1.636 C (seconds; see diary.doc entry for 2005.06.02) C THICKM = (ridge_Delta_ts - Delta_ts) / slope C C (3)Apply lower and upper limits: C THICKM = MAX(THICKM, 0.) THICKM = MIN(THICKM, (HLMAX - THICKC)) C END IF C (End of branch on oceanic vs. continental/unknown lithosphere) C C Compute steady-state GEOTHerm, as in old OrbData: C GEOTH1=TSURF GEOTH2=HEATFL/CONDUC(1) GEOTH3= -RADIO(1)/(2.*CONDUC(1)) GEOTH4=0. TMOHO=GEOTH1+GEOTH2*THICKC+GEOTH3*THICKC**2 GEOTH5=TMOHO QRED=HEATFL-THICKC*RADIO(1) GEOTH6=QRED/CONDUC(2) GEOTH7= -RADIO(2)/(2.*CONDUC(2)) GEOTH8=0. C C Check for Moho hotter than asthenosphere, and C reduce heat flow if necessary. C N.B. This typically occurs where spreading ridges C with rule-based high Q intersect continental C margins, where CRUST2 shows thick crust! C TEST=GEOTH1+GEOTH2*THICKC+GEOTH3*THICKC**2 IF (TEST.GT.TASTHK) THEN TERR0R=TEST-TASTHK DELTAQ=-TERR0R*CONDUC(1)/THICKC HEATFL=HEATFL+DELTAQ QLIMIT=QLIM0+DQL_DE*ELEVAT HEATFL=MAX(HEATFL,QLIMIT) HEATFL=MIN(HEATFL,QLIM1) GEOTH2=HEATFL/CONDUC(1) TMOHO=TASTHK GEOTH5=TMOHO QRED=HEATFL-THICKC*RADIO(1) GEOTH6=QRED/CONDUC(2) END IF C C Compute trial value of cooling_curvature C (referring to GEOTHerm without this effect): C TEST=GEOTH5+GEOTH6*THICKM+GEOTH7*THICKM**2 TERR0R=TEST-TASTHK total_lithosphere = THICKC + THICKM delta_quadratic = -TERR0R / total_lithosphere**2 cooling_curvature = -2. * delta_quadratic C C (Note that we are postponing the addition of C delta_quadratic to GEOTH3 and GEOTH7 until this C value passes some tests!) C C Test whether this value would cause any temperature C maximum within the lithosphere (forbidden!)?: C IF (cooling_curvature > 0.0) THEN t_geoth1 = TSURF t_geoth2 = HEATFL / CONDUC(1) t_geoth3= delta_quadratic - RADIO(1) / (2. * CONDUC(1)) t_TMoho = t_geoth1 + t_geoth2 * THICKC + & t_geoth3 * THICKC**2 t_geoth5 = t_TMoho t_QRed = HEATFL - THICKC * RADIO(1) - & cooling_curvature * THICKC * CONDUC(2) t_geoth6 = QRed / CONDUC(2) t_geoth7= delta_quadratic - RADIO(2) / (2. * CONDUC(2)) C C Check for temperature maximum within mantle lithosphere: c0_of_mantle_gradient = t_geoth6 c1_of_mantle_gradient = 2. * t_geoth7 z_of_maximum = c0_of_mantle_gradient / & (-c1_of_mantle_gradient) IF ((z_of_maximum < THICKM) .AND. & (z_of_maximum > 0.0)) THEN C C Must take corrective action; reduce THICKM C and also change cooling_curvature, C delta_quadratic, etc.. For algebra, see Peter C Bird's notes of 2005.06.15. C z_star = (GEOTH7 * THICKC**2 + GEOTH5 - & GEOTH6 * THICKC - TASTHK) / & (GEOTH7 * THICKC - 0.5 * GEOTH6) m_star = z_star - THICKC THICKM = MAX(MIN(m_star, THICKM), 0.) total_lithosphere = THICKC + THICKM TERR0R = (GEOTH5 + GEOTH6 * THICKM + & GEOTH7 * THICKM**2) - TASTHK delta_quadratic = -TERR0R / total_lithosphere**2 cooling_curvature = -2. * delta_quadratic C END IF END IF C C Build geotherm again with final cooling_curvature: C GEOTH1 = TSURF GEOTH2 = HEATFL / CONDUC(1) GEOTH3= delta_quadratic - RADIO(1) / (2. * CONDUC(1)) GEOTH4 = 0. TMOHO = GEOTH1 + GEOTH2 * THICKC + & GEOTH3 * THICKC**2 GEOTH5 = TMOHO QRED = HEATFL - THICKC * RADIO(1) - & cooling_curvature * THICKC * CONDUC(2) GEOTH6 = QRED / CONDUC(2) GEOTH7= delta_quadratic - RADIO(2) / (2. * CONDUC(2)) GEOTH8 = 0. C chemical_Delta_rho = 0.0 C (Try this case first; adjust density_anomaly below.) C CALL SQUEEZ (INPUT,ALPHAT,chemical_Delta_rho,ELEVAT, + GEOTH1,GEOTH2,GEOTH3,GEOTH4, + GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN,IUNITT, + ONEKM,RHOAST,RHOBAR,RHOH2O, + TEMLIM,THICKC,THICKC+THICKM, + OUTPUT,TAUZZ,SIGZZB) C C Trial value of chemical_density_anomaly: C chemical_Delta_rho = SIGZZB / (GMEAN * total_lithosphere) C C Apply limits: C chemical_Delta_rho = MIN(chemical_Delta_rho, Delta_rho_limit) chemical_Delta_rho = MAX(chemical_Delta_rho, -Delta_rho_limit) C C Repeat isostasy test: C CALL SQUEEZ (INPUT,ALPHAT,chemical_Delta_rho,ELEVAT, + GEOTH1,GEOTH2,GEOTH3,GEOTH4, + GEOTH5,GEOTH6,GEOTH7,GEOTH8, + GMEAN,IUNITT, + ONEKM,RHOAST,RHOBAR,RHOH2O, + TEMLIM,THICKC,THICKC+THICKM, + OUTPUT,TAUZZ,SIGZZB) C C TEST FOR SUCCESSFUL CALCULATION: C C ISOSTASY IS BAD IF BASAL PRESSURE ANOMALY IS LARGE: BADP=(ABS(SIGZZB).GT.4.E7) C C GEOTHERM MIGHT NOT CONNECT TO ASTHENOSPHERE ADIABAT: TEST=GEOTH5+GEOTH6*THICKM+GEOTH7*THICKM**2 TERR0R=TEST-TASTHK BADT=(ABS(TERR0R).GT.20.) C IF (BADP) THEN WRITE (IUNITL,102) PLON,PLAT 102 FORMAT(/' ERROR: (',F8.3,'E, ',F7.3,'N) IS NOT ', + 'ISOSTATIC WITH RISE:') END IF C IF (BADT) THEN WRITE (IUNITL,103) PLON,PLAT 103 FORMAT(/' ERROR: (',F8.3,'E, ',F7.3,'N) DOES NOT ' + ,'CONNECT TO ASTHENOSPHERE ADIABAT') END IF C RETURN END C C C BLOCK DATA BD1 C C DEFINE "WEIGHT" (GAUSSIAN INTEGRATION WEIGHTS) OF THE C SEVEN INTEGRATION POINTS IN EACH ELEMENT, DEFINED BY INTERNAL C COORDINATES "POINTS(3,7)", WHERE POINTS(1-3,M)=S1-S3 OF C INTEGRATION POINT NUMBER M. C BECAUSE ALL OF THESE ARRAYS ARE FUNCTIONS OF INTERNAL C COORDINATES, THEY ARE NOT AFFECTED BY SCALING OR DEFORMATION OF C THE ELEMENTS. C DOUBLE PRECISION POINTS,WEIGHT COMMON /S1S2S3/ POINTS COMMON /WGTVEC/ WEIGHT DIMENSION POINTS(3,7),WEIGHT(7) C C C "POINTS" CONTAINS THE INTERNAL COORDINATES (S1,S2,S3) OF THE 7 C GAUSSIAN INTEGRATION POINTS (FOR AREA INTEGRALS) OF THE C TRIANGULAR ELEMENTS. "POINTS" IS ALSO THE SET OF NODAL FUNCTIONS C FOR UNPROJECTED SCALAR QUANTITIES WITHIN AN ELEMENT. DATA POINTS / + 0.3333333333333333D0, 0.3333333333333333D0, 0.3333333333333333D0, + 0.0597158733333333D0, 0.4701420633333333D0, 0.4701420633333333D0, + 0.4701420633333333D0, 0.0597158733333333D0, 0.4701420633333333D0, + 0.4701420633333333D0, 0.4701420633333333D0, 0.0597158733333333D0, + 0.7974269866666667D0, 0.1012865066666667D0, 0.1012865066666667D0, + 0.1012865066666667D0, 0.7974269866666667D0, 0.1012865066666667D0, + 0.1012865066666667D0, 0.1012865066666667D0, 0.7974269866666667D0/ C C "WEIGHT" IS THE GAUSSIAN WEIGHT (FOR AREA INTEGRALS) OF THE 7 C INTEGRATION POINTS IN EACH TRIANGULAR ELEMENT. DATA WEIGHT / 0.2250000000000000D0, + 0.1323941500000000D0, 0.1323941500000000D0, 0.1323941500000000D0, + 0.1259391833333333D0, 0.1259391833333333D0, 0.1259391833333333D0/ C END C C C BLOCK DATA BD2 C C DEFINE "FPHI" (NODAL FUNCTIONS) AND "FGAUSS" (GAUSSIAN INTEGRATION C WEIGHTS) AT THE 7 INTEGRATION POINTS IN EACH FAULT ELEMENT, C DEFINED BY INTERNAL COORDINATE "FPOINT(M=1,...,7)" C WHICH CONTAINS THE RELATIVE POSITION C (FRACTIONAL LENGTH) OF THE INTEGRATION POINTS. C BECAUSE ALL OF THESE ARRAYS ARE FUNCTIONS OF INTERNAL C COORDINATES, THEY ARE NOT AFFECTED BY SCALING OR DEFORMATION OF C THE ELEMENTS. C DOUBLE PRECISION FPHI,FPOINT,FGAUSS COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS DIMENSION FPHI(4,7),FPOINT(7),FGAUSS(7) C C "FPOINT" CONTAINS THE SEVEN INTEGRATION POINT LOCATIONS FOR THE FAULT C ELEMENTS. EACH VALUE GIVES A POSITION AS A FRACTION OF TOTAL LENGTH C MEASURED FROM NODE1 TO NODE2 (OF ARRAY "NODEF"). DATA FPOINT/ 1 0.0254461D0, 2 0.1292344D0, 3 0.2970774D0, 4 0.5000000D0, 5 0.7029226D0, 6 0.8707656D0, 7 0.9745539D0 / C C "FGAUSS" CONTAINS THE SEVEN CORRESPONDING WEIGHT FACTORS FOR USE IN C LINE INTEGRALS. DATA FGAUSS/ 1 0.0647425D0, 2 0.1398527D0, 3 0.1909150D0, 4 0.2089796D0, 5 0.1909150D0, 6 0.1398527D0, 7 0.0647425D0/ C C "FPHI" CONTAINS THE VALUES OF THE 4 NODAL FUNCTIONS (ONE PER NODE) C AT EACH OF THESE 7 INTEGRATION POINTS IN THE FAULT ELEMENT. C A SPECIAL CONVENTION IS THAT THE NODAL FUNCTION OF NODE 3 C IS THE NEGATIVE OF THAT FOR NODE 2, WHILE THE NODAL FUNCTION C FOR NODE 4 IS THE NEGATIVE OF THAT FOR NODE 1. THIS SIMPLIFIES C MANY EXPRESSIONS IN WHICH WE WOULD OTHERWISE HAVE TO HAVE C A SEPARATE FACTOR OR +1 OR -1 FOR THE TWO SIDES OF THE FAULT. DATA FPHI/ + 0.9745539D0, 0.0254461D0, -0.0254461D0, -0.9745539D0, + 0.8707656D0, 0.1292344D0, -0.1292344D0, -0.8707656D0, + 0.7029226D0, 0.2970774D0, -0.2970774D0, -0.7029226D0, + 0.5000000D0, 0.5000000D0, -0.5000000D0, -0.5000000D0, + 0.2970774D0, 0.7029226D0, -0.7029226D0, -0.2970774D0, + 0.1292344D0, 0.8707656D0, -0.8707656D0, -0.1292344D0, + 0.0254461D0, 0.9745539D0, -0.9745539D0, -0.0254461D0/ C END C C C SUBROUTINE SQUARE (INPUT,BRIEF,FDIP,IUNIT8, + MXBN,MXEL,MXFEL,MXNODE, + MXSTAR,NFL,NODEF,NODES, + NUMEL,NUMNOD,RDETH,WEDGE, + MODIFY,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,EDGEFS, + EDGETS,FLEN,FPFLT,FPSFER, + FARG,NCOND,NODCON,SITA, + WORK,CHECKN,LIST) C C CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID C LOGICAL AGREED,ALLOK,BRIEF,FOUND,SWITCH,VERT1,VERT2 C C NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL CHECKN,EDGEFS,EDGETS C CHARACTER*21 OBLIQU,TAG1,TAG2,VERTIC DOUBLE PRECISION FPHI,FPOINT,FGAUSS COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS DIMENSION FANGLE(2),FPHI(4,7),FPOINT(7),FGAUSS(7),PHI(2),THETA(2) DIMENSION AREA(MXEL),CHECKN(MXNODE), + DETJ(7,MXEL), + DXS(2,2,3,7,MXEL),DYS(2,2,3,7,MXEL), + DXSP(3,7,MXEL),DYSP(3,7,MXEL), + EDGEFS(2,MXFEL),EDGETS(3,MXEL),FDIP(2,MXFEL), + FLEN(MXFEL), + FPFLT(2,2,2,7,MXFEL), + FPSFER(2,2,3,7,MXEL),FARG(2,MXFEL), + LIST(MXSTAR),NODCON(MXBN), + NODEF(4,MXFEL),NODES(3,MXEL), + SITA(7,MXEL),XNODE(MXNODE),YNODE(MXNODE) DATA OBLIQU /'(DIP SLIP IS ALLOWED)'/ DATA VERTIC /'(STRIKE-SLIP ONLY) '/ C C (1) CHECK THAT ALL NODES ARE CONNECTED TO AT LEAST ONE C CONTINUUM (TRIANGULAR) ELEMENT OR FAULT ELEMENT; C DO 110 I=1,NUMNOD CHECKN(I)=.FALSE. 110 CONTINUE DO 130 I=1,NUMEL DO 120 J=1,3 CHECKN(NODES(J,I))=.TRUE. 120 CONTINUE 130 CONTINUE DO 136 I=1,NFL DO 134 J=1,4 CHECKN(NODEF(J,I))=.TRUE. 134 CONTINUE 136 CONTINUE ALLOK=.TRUE. DO 140 I=1,NUMNOD ALLOK=ALLOK.AND.CHECKN(I) 140 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,150) 150 FORMAT(' BAD GRID TOPOLOGY: FOLLOWING REAL NODES DO NOT'/ 1 ' BELONG TO ANY TRIANGULAR CONTINUUM ELEMENT'/ 2 ' OR FAULT ELEMENT:') DO 160 I=1,NUMNOD IF (.NOT.CHECKN(I)) WRITE (IUNIT8,155) I 155 FORMAT (' ',43X,I6) 160 CONTINUE CALL PAUSE() STOP ENDIF C C (2) AVERAGE TOGETHER THE COORDINATES OF ALL NODES AT ONE "POINT" C DO 410 I=1,NUMNOD CHECKN(I)=.FALSE. C (MEANS "NOT YET INVOLVED IN AVERAGING') 410 CONTINUE DO 490 I=1,NFL DO 480 J1=1,2 NJ1=NODEF(J1,I) C (FAULT ENDS ARE THE ONLY PLACES THAT CAN HAVE PROBLEMS) IF (.NOT.CHECKN(NJ1)) THEN LIST(1)=NJ1 CHECKN(NJ1)=.TRUE. C BEGIN LIST OF NEIGHBORS WITH PAIRED NODE J2=5-J1 NJ2=NODEF(J2,I) LIST(2)=NJ2 CHECKN(NJ2)=.TRUE. NINSUM=2 C FIND SHORTEST FAULT CONNECTED TO EITHER ONE DX=XNODE(NJ1)-XNODE(NJ2) DY=YNODE(NJ1)-YNODE(NJ2) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NJ1)) SHORT=SQRT(DX**2+DY**2) DO 470 K=1,NFL NL1=NODEF(1,K) NL2=NODEF(2,K) NL3=NODEF(3,K) NL4=NODEF(4,K) IF ((NJ1.EQ.NL1).OR.(NJ2.EQ.NL1).OR. + (NJ1.EQ.NL2).OR.(NJ2.EQ.NL2).OR. + (NJ1.EQ.NL3).OR.(NJ2.EQ.NL3).OR. + (NJ1.EQ.NL4).OR.(NJ2.EQ.NL4)) THEN DX=XNODE(NL1)-XNODE(NL2) DY=YNODE(NL1)-YNODE(NL2) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NL1)) TEST=SQRT(DX**2+DY**2) SHORT=MIN(SHORT,TEST) ENDIF 470 CONTINUE C COLLECT ALL CORNER NODES WITHIN 10% OF THIS TOLER=SHORT/10. T2=TOLER**2 DO 471 K=1,NUMNOD IF (.NOT.CHECKN(K)) THEN DX=XNODE(NJ1)-XNODE(K) DY=YNODE(NJ1)-YNODE(K) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NJ1)) R2=DX**2+DY**2 IF (R2.LT.T2) THEN NINSUM=NINSUM+1 IF (NINSUM.GT.MXSTAR) THEN WRITE(IUNIT8,421) 421 FORMAT(/' INCREASE VALUE' + ,' OF PARAMETER MAXATP.') CALL PAUSE() STOP ENDIF LIST(NINSUM)=K CHECKN(K)=.TRUE. ENDIF ENDIF 471 CONTINUE C (QUICK EXIT IF ALL NODES IN SAME PLACE) AGREED=.TRUE. DO 472 K=2,NINSUM AGREED=AGREED.AND. + (XNODE(LIST(K)).EQ.XNODE(LIST(1))).AND. + (YNODE(LIST(K)).EQ.YNODE(LIST(1))) 472 CONTINUE IF (AGREED) GO TO 480 XSUM=0. YSUM=0. DO 473 K=1,NINSUM XSUM=XSUM+XNODE(LIST(K)) YSUM=YSUM+YNODE(LIST(K)) 473 CONTINUE XMEAN=XSUM/NINSUM YMEAN=YSUM/NINSUM RMAX=0. DO 474 K=1,NINSUM R=SQRT((XNODE(LIST(K))-XMEAN)**2+ + (YNODE(LIST(K))-YMEAN)**2) RMAX=MAX(RMAX,R) 474 CONTINUE DO 475 K=1,NINSUM XNODE(LIST(K))=XMEAN YNODE(LIST(K))=YMEAN 475 CONTINUE IF (.NOT.BRIEF) THEN IF (RMAX.GT.0.) THEN WRITE(IUNIT8,476) NINSUM, + (LIST(N),N=1,NINSUM) 476 FORMAT(/ + ' AVERAGING TOGETHER THE POSITIONS OF', + ' THESE ',I6,' NODES:',(/' ',12I6)) WRITE (IUNIT8,477) RMAX 477 FORMAT (' MAXIMUM CORRECTION TO ', + 'ANY POSITION IS',1P,E10.2/ + ' YOU ARE RESPONSIBLE FOR ', + ' DECIDING WHETHER THIS IS A', + ' SERIOUS ERR0R!') ENDIF ENDIF ENDIF 480 CONTINUE 490 CONTINUE C C (3) COMPUTE DERIVITIVES OF NODAL C FUNCTIONS AT INTEGRATION POINTS; C THEN CHECK FOR NEGATIVE AREAS C CALL DERIV (INPUT,IUNIT8,MXEL,MXNODE,NODES,NUMEL, + RDETH,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,FPSFER,SITA) ALLOK=.TRUE. DO 620 I=1,NUMEL DO 610 M=1,7 TEST=AREA(I)*DETJ(M,I) IF (TEST.LE.0.) THEN WRITE(IUNIT8,605) M,I 605 FORMAT(/' EXCESSIVELY DISTORTED ELEMENT LEADS TO ' + ,'NEGATIVE AREA AT POINT ',I1,' IN ELEMENT ', + I5) WRITE(IUNIT8,606) AREA(I),DETJ(M,I) 606 FORMAT('AREA = ',1P,E14.4,' DETJ: ',0P,F14.6) ALLOK=.FALSE. ENDIF 610 CONTINUE 620 CONTINUE IF (.NOT.ALLOK) THEN CALL PAUSE() STOP ENDIF C C (4) COMPUTE LENGTHS OF FAULT ELEMENTS. C DO 750 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) THETA1=XNODE(N1) THETA2=XNODE(N2) PHI1 =YNODE(N1) PHI2 =YNODE(N2) FLEN(I)=FLTLEN (PHI1,PHI2,RDETH,THETA1,THETA2) 750 CONTINUE C C (5) MAKE A LIST OF NODES THAT ARE ON THE BOUNDARY AND REQUIRE C BOUNDARY CONDITIONS (NODCON); THESE ARE IN COUNTERCLOCKWISE C ORDER. ALSO MAKE A LISTS OF ELEMENT SIDES WHICH CONTAIN THESE C NODES: EDGETS AND EDGEFS. C NCOND=0 DO 801 I=1,NUMNOD CHECKN(I)=.FALSE. 801 CONTINUE DO 802 I=1,NFL EDGEFS(1,I)=.FALSE. EDGEFS(2,I)=.FALSE. 802 CONTINUE DO 810 I=1,NUMEL DO 809 J=1,3 CALL NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) IF (KELE.GT.0) THEN C (ORDINARY INTERIOR SIDE) EDGETS(J,I)=.FALSE. ELSE IF (KFAULT.EQ.0) THEN C (EXTERIOR SIDE) EDGETS(J,I)=.TRUE. N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J+1,3)+1,I) IF (.NOT.CHECKN(N1)) THEN NCOND=NCOND+1 CHECKN(N1)=.TRUE. ENDIF IF (.NOT.CHECKN(N2)) THEN NCOND=NCOND+1 CHECKN(N2)=.TRUE. ENDIF ELSE C (TRIANGULAR ELEMENT HAS AN EXTERIOR FAULT ELEMENT C ADJACENT TO IT) EDGETS(J,I)=.FALSE. N1=NODES(MOD(J, 3)+1,I) IF (NODEF(2,KFAULT).EQ.N1) THEN EDGEFS(2,KFAULT)=.TRUE. DO 806 K=3,4 N=NODEF(K,KFAULT) IF (.NOT.CHECKN(N)) THEN NCOND=NCOND+1 CHECKN(N)=.TRUE. ENDIF 806 CONTINUE ELSE EDGEFS(1,KFAULT)=.TRUE. DO 808 K=1,2 N=NODEF(K,KFAULT) IF (.NOT.CHECKN(N)) THEN NCOND=NCOND+1 CHECKN(N)=.TRUE. ENDIF 808 CONTINUE ENDIF ENDIF 809 CONTINUE 810 CONTINUE IF (NCOND.GT.MXBN) THEN WRITE(IUNIT8,820) NCOND 820 FORMAT(/' INCREASE PARAMETER MAXBN TO',I6,' AND RECOMPILE.') CALL PAUSE() STOP ENDIF C C STOP WORK IF NO BOUNDARY NODES FOUND (GLOBAL GRID) C IF (NCOND.EQ.0) GO TO 899 C C BEGIN CIRCUIT WITH LOWEST-NUMBERED BOUNDARY NODE DO 830 I=1,NUMNOD IF (CHECKN(I)) GO TO 831 830 CONTINUE 831 NODCON(1)=I NDONE=1 NLEFT=NCOND C BEGINNING OF INDEFINATE LOOP WHICH TRACES AROUND THE PERIMETER. C EACH TIME, IT PROGRESSES BY ONE OF 3 STEPS: C -2 NODES AT A TIME ALONG A TRIANGLE SIDE, OR C -2 NODES AT A TIME ALONG A FAULT ELEMENT SIDE, OR C -BY FINDING ANOTHER (CORNER) NODE WHICH SHARES THE SAME LOCATION. C BEGINNING OF MAIN INDEFINATE LOOP: 840 NODE=NODCON(NDONE) X=XNODE(NODE) Y=YNODE(NODE) C LOOK FOR AN ADJACENT TRIANGULAR ELEMENT USING THIS NODE. DO 844 I=1,NUMEL DO 842 J=1,3 IF (EDGETS(J,I)) THEN N1=NODES(MOD(J,3)+1,I) IF (N1.EQ.NODE) GO TO 846 ENDIF 842 CONTINUE 844 CONTINUE GO TO 850 846 N2=NODES(MOD(J+1,3)+1,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR AN ADJACENT FAULT ELEMENT USING THIS NODE. 850 DO 854 I=1,NFL IF (EDGEFS(1,I)) THEN IF (NODEF(1,I).EQ.NODE) THEN N2=NODEF(2,I) GO TO 856 ENDIF ELSE IF (EDGEFS(2,I)) THEN IF (NODEF(3,I).EQ.NODE) THEN N2=NODEF(4,I) GO TO 856 ENDIF ENDIF 854 CONTINUE GO TO 860 856 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR ANOTHER EXTERIOR CORNER NODE AT SAME LOCATION 860 DO 865 I=1,NUMNOD IF ((I.NE.NODE).AND.CHECKN(I)) THEN IF ( (ABS(XNODE(I)-X).LT.1.E-6) .AND. + (ABS(YNODE(I)-Y).LT.1.E-6) ) GO TO 867 ENDIF 865 CONTINUE WRITE(IUNIT8,866) NODE 866 FORMAT(' BAD GRID TOPOLOGY: WHILE TRACING PERIMETER,'/ + ' COULD NOT FIND ANY WAY TO CONTINUE FROM NODE ',I6/ + ' EITHER THROUGH SHARED BOUNDARY ELEMENTS, OR'/ + ' THROUGH OTHER BOUNDARY NODES SHARING THE SAME ', + 'POSITION.') CALL PAUSE() STOP 867 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=I CHECKN(I)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) GO TO 840 C END OF INDEFINATE LOOP WHICH TRACES AROUND PERIMETER. 870 IF (.NOT.BRIEF) THEN WRITE(IUNIT8,880) 880 FORMAT(/ /' HERE FOLLOWS A LIST, IN CONSECUTIVE ORDER,'/ + ' OF THE NODES WHICH DEFINE THE PERIMETER'/ + ' OF THE MODEL; THESE NODES REQUIRE BOUNDARY', + ' CONDITIONS:'/' BC# NODE') DO 890 I=1,NCOND N=NODCON(I) WRITE(IUNIT8,882) I, N 882 FORMAT(' ',2I6) 890 CONTINUE N=NODCON(1) WRITE (IUNIT8,892) N 892 FORMAT(' (NOTE: NODE ',I6,' COMPLETES THE LOOP, BUT WILL', + ' NOT BE LISTED TWICE.)') ENDIF 899 CONTINUE C C (6) SURVEY FAULT ELEMENTS AND ISSUE WARNING IF ANY ELEMENT IS OF C MIXED TYPE (PART STRIKE-SLIP, AND PART SHALLOW-DIPPING: C DO 920 I=1,NFL DELD1=FDIP(1,I)-1.570796 DELD2=FDIP(2,I)-1.570796 VERT1=ABS(DELD1).LE.WEDGE VERT2=ABS(DELD2).LE.WEDGE NVPART=0 IF (VERT1) THEN NVPART=NVPART+1 TAG1=VERTIC ELSE TAG1=OBLIQU ENDIF IF (VERT2) THEN NVPART=NVPART+1 TAG2=VERTIC ELSE TAG2=OBLIQU ENDIF SWITCH=((NVPART.GT.0).AND.(NVPART.LT.2)) IF (SWITCH) THEN DIP1=FDIP(1,I)*57.2957795 IF (DIP1.GT.90.) DIP1=DIP1-180. DIP2=FDIP(2,I)*57.2957795 IF (DIP2.GT.90.) DIP2=DIP2-180. WRITE (IUNIT8,905) I,DIP1,TAG1,DIP2,TAG2 905 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',F7.2,' DEGREES ',A21/ + ' ',F7.2,' DEGREES ',A21/ + ' WHICH MAKES IT MIXED-MODE.'/ + ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ + ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ + ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', + ' FAULT ELEMENT TO THE NEXT.)') ELSE NVPART=0 DO 910 M=1,7 DELD=DELD1*FPHI(1,M)+DELD2*FPHI(2,M) IF (ABS(DELD).LE.WEDGE) NVPART=NVPART+1 910 CONTINUE IF ((NVPART.GT.0).AND.(NVPART.LT.7)) THEN IF (NVPART.GE.4) THEN WRITE (IUNIT8,912) I,DIP1,DIP2 912 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',F7.2,' DEGREES, AND'/ + ' ',F7.2,' DEGREES'/ + ' WHICH APPEAR TO MAKE IT STRIKE-SLIP.'/ + ' HOWEVER, THESE VALUES ARE SUCH THAT DIP-SLIP'/ + ' IS PERMITTED AT ONE OR MORE INTEGRATION POINTS.'/ + ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ + ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ + ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', + ' FAULT ELEMENT TO THE NEXT.)') ELSE WRITE (IUNIT8,914) I,DIP1,DIP2 914 FORMAT(/ /' CAUTION:'/ + ' FAULT ELEMENT ',I6,' HAS DIPS OF '/ + ' ',F7.2,' DEGREES, AND'/ + ' ',F7.2,' DEGREES'/ + ' WHICH APPEAR TO MAKE IT FREE-SLIPPING.'/ + ' HOWEVER, THESE VALUES ARE SUCH THAT DIP-SLIP'/ + ' IS PROHIBITED AT ONE OR MORE INTEGRATION POINTS.'/ + ' SUCH ELEMENTS ARE INACCURATE AND NOT RECOMMENDED.'/ + ' PREFERABLY EACH ELEMENT SHOULD BE OF A SINGLE TYPE.'/ + ' (REMEMBER, DIP NEED NOT BE CONTINUOUS FROM ONE', + ' FAULT ELEMENT TO THE NEXT.)') ENDIF ENDIF ENDIF 920 CONTINUE C C (7) CALCULATE FAULT ARGUMENT (IN RADIANS, MEASURED COUNTERCLOCKWISE C FROM +THETA) FOR EACH INTEGRATION POINT IN EACH FAULT ELEMENT. C C HERE USE CROSS PRODUCT "D" OF POSITION VECTOR "A" AND "B" C OF FAULT ELEMENT NODE "N1" AND "N2" RESPECTIVELY, TO CALCULATE C VECTOR "D" IN PLANE OF "A" AND "B" AND PERPENDICULAR TO "A". C THEN USE "D" AND "A" TO DECIDE THE POSITION VECTOR "GOLD" OF C INTERPOLATION POINT "M" AND "G" OF "M+1". FINALLY USE THE C DOT PRODUCT OF OF "GOLD-POL" AND "G-GOLD" TO DECIDE THE ANGLE C BETWEEN THESE TWO VECTORS, I.E., FARG. HERE POL IS POSITION C VECTOR OF NORTH POLE. DO 1000 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) THETA(1)=XNODE(N1) THETA(2)=XNODE(N2) PHI(1) =YNODE(N1) PHI(2) =YNODE(N2) CALL FANGLS (INPUT,PHI,THETA, + OUTPUT,FANGLE) DO 900 J=1,2 FARG(J,I)=FANGLE(J) FANGLE(J)=FANGLE(J)*57.29577951 900 CONTINUE 1000 CONTINUE C C (8) SURVEY STRIKE-SLIP (VERTICAL) FAULTS TO CHECK FOR CONFLICTS IN C ARGUMENT THAT WOULD LOCK THE FAULT. C C LOOP ON ALL FAULT ELEMENTS (I): DO 2000 I=1,NFL C LOOP ON 2 TERMINAL NODE PAIRS, 1-4, 2-3 (J = 1 OR 2): DO 1900 J=1,2 C DIP MUST BE WITHIN "WEDGE" OF VERTICAL FOR CONSTRAINT: IF (ABS(FDIP(J,I)-1.570796).LE.WEDGE) THEN NAZI=J N1=J IF(J.EQ.1) THEN N4=4 ELSE N4=3 ENDIF NODE1=NODEF(N1,I) NODE4=NODEF(N4,I) C NO CONSTRAINT APPLIED WHERE A FAULT ENDS: IF (NODE1.NE.NODE4) THEN C ENDPOINT PAIRS MUST BE CHECKED FOR DUPLICATION: C LOOK FOR OTHER STRIKE-SLIP FAULTS SHARING THIS C PAIR OF NODES, AT EITHER END: FOUND=.FALSE. DO 1600 L=1,NFL IF (L.NE.I) THEN IF (ABS(FDIP(1,L)-1.5708).LE.WEDGE) THEN IF (((NODE1.EQ.NODEF(1,L)).AND. + (NODE4.EQ.NODEF(4,L))).OR. + ((NODE1.EQ.NODEF(4,L)).AND. + (NODE4.EQ.NODEF(1,L)))) THEN FOUND=.TRUE. NUMBER=L NAZL=1 GO TO 1601 ENDIF ENDIF IF (ABS(FDIP(2,L)-1.5708).LE.WEDGE) THEN IF (((NODE1.EQ.NODEF(2,L)).AND. + (NODE4.EQ.NODEF(3,L))).OR. + ((NODE1.EQ.NODEF(3,L)).AND. + (NODE4.EQ.NODEF(2,L)))) THEN FOUND=.TRUE. NUMBER=L NAZL=2 GO TO 1601 ENDIF ENDIF ENDIF 1600 CONTINUE C DON'T WORRY IF THIS PAIR ALREADY CHECKED! 1601 IF (FOUND.AND.(NUMBER.GT.I)) THEN C AVERAGE ARGUMENTS TOGETHER (AVOID CYCLE SHIFTS): IF(NAZI.EQ.NAZL) THEN AZL=FARG(NAZL,NUMBER)+3.141592654 ELSE AZL=FARG(NAZL,NUMBER) ENDIF AZI=FARG(NAZI,I) COSZ=0.5*(COS(AZI)+COS(AZL)) SINZ=0.5*(SIN(AZI)+SIN(AZL)) AZIMUT=ATAN2F(SINZ,COSZ) FARG(NAZI,I)=AZIMUT IF(NAZL.EQ.NAZI) THEN FARG(NAZL,NUMBER)=AZIMUT-3.141592654 ELSE FARG(NAZL,NUMBER)=AZIMUT ENDIF ENDIF C ^END BLOCK WHICH LOOKS FOR CONSTRAINTS ON REAL NODES ENDIF C ^END BLOCK WHICH CHECKS FOR DISTINCT NODE NUMBERS ENDIF C ^END BLOCK WHICH CHECKS FOR DIP OF OVER 75 DEGREES 1900 CONTINUE C ^END LOOP ON 2 NODE PAIRS IN FAULT ELEMENT 2000 CONTINUE C C (9) CALCULATE NODAL FUNCTION AT INTERPOLATION POINTS C ON SPHERICAL GREAT CIRCLE FAULT. C CALL FNODAL (INPUT,IUNITT,MXFEL, + MXNODE,NFL,NODEF,XNODE,YNODE, + OUTPUT,FPFLT) C IF (.NOT. BRIEF) WRITE (IUNIT8,9999) 9999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE DERIV (INPUT,IUNITT,MXEL,MXNODE, + NODES,NUMEL, + RDETH,XNODE,YNODE, + OUTPUT,AREA,DETJ, + DXS,DYS,DXSP,DYSP,FPSFER,SITA) C C SETS UP 6 VECTOR NODAL FUNCTIONS (FPSFER) OF EACH SPHERICAL C TRIANGLE FINITE ELEMENT, AT EACH OF ITS 7 INTEGRATION POINTS. C CALCULATES DXS AND DYS, THE THETA-DERIVITIVE AND PHI-DERIVITIVE C OF EACH OF THESE 6 VECTOR NODAL FUNCTIONS. C ALSO COMPUTES AREA, THE AREAS OF THE PLANE TRIANGLES. C ALSO COMPUTES DETJ, THE LOCAL RATIO OF AREAS ON THE SPHERE C TO AREAS ON THE PLANE TRIANGLES. C DOUBLE PRECISION POINTS DOUBLE PRECISION FFF,SKKC,SKKE,SNCSNE,SNCCSE,CSCCSE,CSCSNE DOUBLE PRECISION XA,XB,XC,YA,YB,YC,ZA,ZB,ZC,XYZP DIMENSION XNODE(MXNODE),YNODE(MXNODE),NODES(3,MXEL),AREA(MXEL) DIMENSION DETJ(7,MXEL) DIMENSION DXS(2,2,3,7,MXEL),DYS(2,2,3,7,MXEL) DIMENSION DXSP(3,7,MXEL),DYSP(3,7,MXEL),POINTS(3,7) DIMENSION PHI(3),THETA(3),SKKC(3),SKKE(3),FFF(3), + SITA(7,MXEL),FPSFER(2,2,3,7,MXEL) DIMENSION XLON(3),YLAT(3) COMMON /S1S2S3/ POINTS DATA OEZOPI /57.29578/ C DO 900 I=1,NUMEL DO 100 J=1,3 THETA(J)=XNODE(NODES(J,I)) PHI(J) =YNODE(NODES(J,I)) 100 CONTINUE X21=SIN(THETA(2))*COS(PHI(2))-SIN(THETA(1))*COS(PHI(1)) X31=SIN(THETA(3))*COS(PHI(3))-SIN(THETA(1))*COS(PHI(1)) Y21=SIN(THETA(2))*SIN(PHI(2))-SIN(THETA(1))*SIN(PHI(1)) Y31=SIN(THETA(3))*SIN(PHI(3))-SIN(THETA(1))*SIN(PHI(1)) Z21=COS(THETA(2))-COS(THETA(1)) Z31=COS(THETA(3))-COS(THETA(1)) A=Y21*Z31-Y31*Z21 B=Z21*X31-Z31*X21 C=X21*Y31-X31*Y21 AREAP=SQRT(A*A+B*B+C*C) AREA(I)=0.5*RDETH*RDETH*AREAP IF (AREA(I).LE.0.0) THEN DO 9 K=1,3 XLON(K)=PHI(K)*OEZOPI YLAT(K)=90.-THETA(K)*OEZOPI 9 CONTINUE WRITE (IUNITT,10) AREA(I), I, + (NODES(K,I), XLON(K), YLAT(K), K=1,3) 10 FORMAT (/' ERR0R: NON-POSITIVE AREA OF ',1P,E10.2, + ' FOR ELEMENT ',I5/ + ' NODES LONGITUDE LATITUDE', + 3(/' ',I7,0P,2F10.3)) CALL PAUSE() STOP ENDIF PNX=A/AREAP PNY=B/AREAP PNZ=C/AREAP DD1=SIN(THETA(1))*COS(PHI(1))*PNX DD2=SIN(THETA(1))*SIN(PHI(1))*PNY DD3=COS(THETA(1))*PNZ DD=DD1+DD2+DD3 C C THIS PART IS TO TEST IF KONG'S METHOD AND BIRD'S METHOD ARE SAME FOR C CALCULATING DEERIVATIVE C XA=SIN(THETA(1))*COS(PHI(1)) XB=SIN(THETA(2))*COS(PHI(2)) XC=SIN(THETA(3))*COS(PHI(3)) YA=SIN(THETA(1))*SIN(PHI(1)) YB=SIN(THETA(2))*SIN(PHI(2)) YC=SIN(THETA(3))*SIN(PHI(3)) ZA=COS(THETA(1)) ZB=COS(THETA(2)) ZC=COS(THETA(3)) CKA=(YB*ZC-ZB*YC)*XA+(ZB*XC-XB*ZC)*YA+(XB*YC-YB*XC)*ZA C C C DO 800 M=1,7 SNCCSE=0.0 SNCSNE=0.0 COSM=0.0 DO 200 J=1,3 SNCCSE=SNCCSE+POINTS(J,M)*SIN(THETA(J))*COS(PHI(J)) SNCSNE=SNCSNE+POINTS(J,M)*SIN(THETA(J))*SIN(PHI(J)) COSM=COSM+POINTS(J,M)*COS(THETA(J)) 200 CONTINUE XYZP=SQRT(SNCCSE*SNCCSE+SNCSNE*SNCSNE+COSM*COSM) SNCCSE=SNCCSE/XYZP SNCSNE=SNCSNE/XYZP COSM=COSM/XYZP SITAJ=ACOS(COSM) TY=SNCSNE TX=SNCCSE PHAIJ=ATAN2F(TY,TX) CSCCSE=COS(SITAJ)*COS(PHAIJ) CSCSNE=COS(SITAJ)*SIN(PHAIJ) C C BIRD'S METHOD C FFF(1)=((YB*ZC-ZB*YC)*SNCCSE+(ZB*XC-XB*ZC)*SNCSNE + +(XB*YC-YB*XC)*COSM)/CKA FFF(2)=((YC*ZA-ZC*YA)*SNCCSE+(ZC*XA-XC*ZA)*SNCSNE + +(XC*YA-YC*XA)*COSM)/CKA FFF(3)=((YA*ZB-ZA*YB)*SNCCSE+(ZA*XB-XA*ZB)*SNCSNE + +(XA*YB-YA*XB)*COSM)/CKA SKKC(1)=((YB*ZC-ZB*YC)*CSCCSE + +(ZB*XC-XB*ZC)*CSCSNE + -(XB*YC-YB*XC)*SIN(SITAJ))/CKA SKKC(2)=((YC*ZA-ZC*YA)*CSCCSE + +(ZC*XA-XC*ZA)*CSCSNE + -(XC*YA-YC*XA)*SIN(SITAJ))/CKA SKKC(3)=((YA*ZB-ZA*YB)*CSCCSE + +(ZA*XB-XA*ZB)*CSCSNE + -(XA*YB-YA*XB)*SIN(SITAJ))/CKA SKKE(1)=(-(YB*ZC-ZB*YC)*SNCSNE + +(ZB*XC-XB*ZC)*SNCCSE)/CKA SKKE(2)=(-(YC*ZA-ZC*YA)*SNCSNE + +(ZC*XA-XC*ZA)*SNCCSE)/CKA SKKE(3)=(-(YA*ZB-ZA*YB)*SNCSNE + +(ZA*XB-XA*ZB)*SNCCSE)/CKA C C C SITA(M,I)=SITAJ RR1=SIN(SITAJ)*COS(PHAIJ) RR2=SIN(SITAJ)*SIN(PHAIJ) RR3=COS(SITAJ) RN=RR1*PNX+RR2*PNY+RR3*PNZ PP=DD/RN DPDC=(COS(SITAJ)*COS(PHAIJ)*PNX+COS(SITAJ)*SIN(PHAIJ)*PNY + -SIN(SITAJ)*PNZ) DPDE=(-SIN(SITAJ)*SIN(PHAIJ)*PNX+ + SIN(SITAJ)*COS(PHAIJ)*PNY) DDPN=PP/RN DPDC=-DDPN*DPDC DPDE=-DDPN*DPDE IF(SITA(M,I).LE.0.0.OR.SITA(M,I).GE.3.141592654) THEN SITAMI=SITA(M,I)*57.29577951 WRITE(IUNITT,220) M,I,SITAMI 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5,' OF ELEMENT', + I5,' IS OUT RANGE', + E14.4) CALL PAUSE() STOP ENDIF DO 500 J=1,3 DXSP(J,M,I)= DPDC*FFF(J)+PP*SKKC(J) DYSP(J,M,I)= DPDE*FFF(J)+PP*SKKE(J) CSCS=COS(THETA(J))*COS(PHI(J)) CSSN=COS(THETA(J))*SIN(PHI(J)) SNC=SIN(THETA(J)) SNE=SIN(PHI(J)) CSE=COS(PHI(J)) FPSFER(1,1,J,M,I)= CSCS*CSCCSE+CSSN*CSCSNE + +SNC*SIN(SITAJ) FPSFER(2,1,J,M,I)=-SNE*CSCCSE+CSE*CSCSNE FPSFER(1,2,J,M,I)=-CSCS*SIN(PHAIJ)+CSSN*COS(PHAIJ) FPSFER(2,2,J,M,I)=SNE*SIN(PHAIJ)+CSE*COS(PHAIJ) DXS(1,1,J,M,I)=(-CSCS*SNCCSE-CSSN*SNCSNE + +SNC*COS(SITAJ))*FFF(J) + +FPSFER(1,1,J,M,I)*SKKC(J) DXS(2,1,J,M,I)=(SNE*SNCCSE-CSE*SNCSNE)*FFF(J) + +FPSFER(2,1,J,M,I)*SKKC(J) DYS(1,1,J,M,I)=(-CSCS*CSCSNE+CSSN*CSCCSE)*FFF(J) + +FPSFER(1,1,J,M,I)*SKKE(J) DYS(2,1,J,M,I)=(SNE*CSCSNE+CSE*CSCCSE)*FFF(J) + +FPSFER(2,1,J,M,I)*SKKE(J) DXS(1,2,J,M,I)=FPSFER(1,2,J,M,I)*SKKC(J) DXS(2,2,J,M,I)=FPSFER(2,2,J,M,I)*SKKC(J) DYS(1,2,J,M,I)=(-CSCS*COS(PHAIJ)-CSSN*SIN(PHAIJ)) + *FFF(J) + +FPSFER(1,2,J,M,I)*SKKE(J) DYS(2,2,J,M,I)=(SNE*COS(PHAIJ)-CSE*SIN(PHAIJ)) + *FFF(J) + +FPSFER(2,2,J,M,I)*SKKE(J) FPSFER(1,1,J,M,I)=FPSFER(1,1,J,M,I)*FFF(J) FPSFER(2,1,J,M,I)=FPSFER(2,1,J,M,I)*FFF(J) FPSFER(1,2,J,M,I)=FPSFER(1,2,J,M,I)*FFF(J) FPSFER(2,2,J,M,I)=FPSFER(2,2,J,M,I)*FFF(J) 500 CONTINUE PFQ=FFF(1)+FFF(2)+FFF(3) DETJ(M,I)=RN**3/(DD*DD) 800 CONTINUE 900 CONTINUE RETURN END C C C SUBROUTINE NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) C C DETERMINE WHETHER THERE ARE MORE ELEMENTS ADJACENT TO SIDE J OF C TRIANGULAR CONTINUUM ELEMENT #I. C J = 1 MEANS THE SIDE OPPOSITE NODE # NODES(1,I). C J = 2 MEANS THE SIDE OPPOSITE NODE # NODES(2,I). C J = 3 MEANS THE SIDE OPPOSITE NODE # NODES(3,I). C IF A FAULT ELEMENT IS ADJACENT, ITS NUMBER IS KFAULT; OTHERWISE, C KFAULT IS SET TO ZERO. C IF ANOTHER TRIANGULAR CONTINUUM ELEMENT IS ADJACENT (EVEN ACROSS C FAULT ELEMENT KFAULT!) THEN ITS NUMBER IS KELE; OTHERWISE, KELE = 0. C LOGICAL FOUNDF DIMENSION NODEF(4,MXFEL),NODES(3,MXEL) C C THREE NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE: N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J+1,3)+1,I) C CHECK FOR ADJACENT FAULT ELEMENT FIRST: FOUNDF=.FALSE. KFAULT=0 IF (NFL.GT.0) THEN DO 10 K=1,NFL M1=NODEF(1,K) M2=NODEF(2,K) M3=NODEF(3,K) M4=NODEF(4,K) IF (((M1.EQ.N2).AND.(M2.EQ.N1)).OR. + ((M3.EQ.N2).AND.(M4.EQ.N1))) THEN FOUNDF=.TRUE. KFAULT=K GO TO 11 ENDIF 10 CONTINUE 11 IF (.NOT.FOUNDF) KFAULT=0 C IF THERE WAS A FAULT, REPLACE 2 NODE NUMBERS THAT WE SEARCH FOR: IF (FOUNDF) THEN IF (M2.EQ.N1) THEN N1=M3 N2=M4 ELSE N1=M1 N2=M2 ENDIF ENDIF ENDIF C SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT: KELE=0 KLOW=I KHIGH=I C --- BEGIN IRREGULAR LOOP, TO SEARCH OUT NEAREST ELEMENTS FIRST --- 100 KLOW=KLOW-1 IF (KLOW.GE.1) THEN DO 110 L=1,3 M1=NODES(MOD(L, 3)+1,KLOW) M2=NODES(MOD(L+1,3)+1,KLOW) IF ((M2.EQ.N1).AND.(M1.EQ.N2)) THEN KELE=KLOW RETURN ENDIF 110 CONTINUE ENDIF KHIGH=KHIGH+1 IF (KHIGH.LE.NUMEL) THEN DO 120 L=1,3 M1=NODES(MOD(L, 3)+1,KHIGH) M2=NODES(MOD(L+1,3)+1,KHIGH) IF ((M2.EQ.N1).AND.(M1.EQ.N2)) THEN KELE=KHIGH RETURN ENDIF 120 CONTINUE ENDIF IF ((KLOW.GT.1).OR.(KHIGH.LT.NUMEL)) GO TO 100 RETURN END C C C REAL FUNCTION FLTLEN (PHI1,PHI2,RDETH,THETA1,THETA2) C C CALCULATES LENGTH OF GREAT CIRCLE SEGMENT BETWEEN C POINT (THETA1,PHI1) AND POINT (THETA2,PHI2), C IN PHYSICAL LENGTH UNITS (RADIANS*PLANET RADIUS). C DOUBLE PRECISION AB AB =SIN(THETA1)*SIN(THETA2)*COS(PHI1)*COS(PHI2)+ + SIN(THETA1)*SIN(THETA2)*SIN(PHI1)*SIN(PHI2)+ + COS(THETA1)*COS(THETA2) AB=ACOS(AB) FLTLEN=AB*RDETH RETURN END C C C SUBROUTINE READPM (INPUT,IUNIT7, IUNIT8, NAMES, NUMPLT, OFFMAX, + OUTPUT,ACREEP, ALPHAT, BCREEP, BIOT , + BYERLY, CCREEP, CFRIC , CONDUC, + DCREEP, ECREEP, EVERYP, FFRIC , GMEAN , + GRADIE, ICONVE, IPVREF, + MAXITR, OKDELV, OKTOQT, ONEKM, RADIO , + RDETH, REFSTR, RHOAST, RHOBAR, RHOH2O, + TADIAB, + TAUMAX, TEMLIM, TITLE3, TRHMAX, TSURF, + ZBASTH) C C READS STRATEGIC AND TACTICAL INPUT PARAMETERS FROM DEVICE IUNIT7, C AND ECHOES THEM ON DEVICE IUNIT8 WITH ANNOTATIONS. C CHARACTER*2 NAMES,PLTREF CHARACTER*80 TITLE3 LOGICAL EVERYP DIMENSION ACREEP(2),ALPHAT(2),BCREEP(2),CCREEP(2),CONDUC(2), + DCREEP(2),NAMES(NUMPLT),RADIO(2), + RHOBAR(2),TAUMAX(2),TEMLIM(2),VECTOR(2) C WRITE (IUNIT8,1) IUNIT7 1 FORMAT (/ /' ATTEMPTING TO READ PARAMETERS FROM UNIT',I3/) READ (IUNIT7,2) TITLE3 2 FORMAT (A80) WRITE (IUNIT8,3) TITLE3 3 FORMAT (/' TITLE OF PARAMETER SET ='/' ',A80) WRITE (IUNIT8,4) 4 FORMAT (/' **************************************************'/ + ' IT IS THE USERS RESPONSIBILITY TO INPUT ALL OF THE'/ + ' FOLLOWING NUMERICAL QUANTITIES IN CONSISTENT UNITS,'/ + ' SUCH AS SYSTEM-INTERNATIONAL (SI) OR CM-G-S (CGS).'/ + ' NOTE THAT TIME UNIT MUST BE THE SECOND (HARD-CODED).'/ + ' **************************************************'/ + /' ========== STRATEGIC PARAMETERS (DEFINE THE REAL', + '-EARTH PROBLEM) ======'/) READ (IUNIT7,*) FFRIC WRITE (IUNIT8,20) FFRIC 20 FORMAT (' ', F10.3,' COEFFICIENT OF FRICTION ON FAULTS') READ (IUNIT7,*) CFRIC WRITE (IUNIT8,30) CFRIC 30 FORMAT (' ', F10.3,' COEFFICIENT OF FRICTION WITHIN BLOCKS') READ (IUNIT7,*) BIOT BIOT = MAX(0.0,MIN(1.0,BIOT)) WRITE (IUNIT8,40) BIOT 40 FORMAT (' ',F10.4,' EFFECTIVE-PRESSURE (BIOT) COEFFICIENT,', + ' 0.0 TO 1.0') READ (IUNIT7,*) BYERLY BYERLY = MAX(0.0,MIN(0.99,BYERLY)) IF (OFFMAX.GT.0.) THEN WRITE (IUNIT8,41) BYERLY 41 FORMAT (' ',F10.4,' BYERLY COEFFICIENT (0. TO 0.99) ='/ + 11X,' FRACTIONAL FRICTION REDUCTION ON MASTER', + ' FAULT'/ + 11X,' (OTHER FAULTS HAVE LESS REDUCTION, IN', + ' PROPORTION TO'/ + 11X,' THEIR TOTAL PAST OFFSETS)') ELSE WRITE (IUNIT8,42) BYERLY 42 FORMAT (' ',F10.4,' BYERLY COEFFICIENT (NOT USED IN', + ' THIS RUN,'/ + 11X,' AS ALL OFFSETS ARE ZERO)') ENDIF CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,ACREEP) WRITE (IUNIT8,50) ACREEP(1),ACREEP(2) 50 FORMAT (' ',1P, E10.2,'/',E10.2,' A FOR CREEP = ', + 'PRE-EXPONENTIAL SHEAR', + ' STRESS CONSTANT FOR CREEP. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,BCREEP) WRITE (IUNIT8,60) BCREEP(1),BCREEP(2) 60 FORMAT (' ', F10.0,'/',F10.0,' B FOR CREEP =(ACTIVATION ', + 'ENERGY)/R/N', + ' IN K. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,CCREEP) WRITE (IUNIT8,70) CCREEP(1),CCREEP(2) 70 FORMAT (' ',1P, E10.2,'/',E10.2,' C FOR CREEP = DERIVATIVE OF B', + ' WITH RESPECT TO DEPTH. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,DCREEP) WRITE (IUNIT8,80) DCREEP(1),DCREEP(2) 80 FORMAT (' ',1P, E10.2,'/',E10.2,' D FOR CREEP = MAXIMUM SHEAR ', + 'STRESS UNDER ANY CONDITIONS. (CRUST/MANTLE)') READ (IUNIT7,*) ECREEP WRITE (IUNIT8,90) ECREEP 90 FORMAT (' ', F10.6,' E FOR CREEP = STRAIN-RATE EXPONENT FOR', + ' CREEP (1/N). (SAME FOR CRUST AND MANTLE!)') READ (IUNIT7,*) TADIAB, GRADIE WRITE (IUNIT8,92) TADIAB, GRADIE 92 FORMAT (' ',F10.0,1P,E10.2,' INTERCEPT AND GRADIENT OF UPPER' + ,' MANTLE ADIABAT (K, K/M)') READ (IUNIT7,*) ZBASTH WRITE (IUNIT8,94) ZBASTH 94 FORMAT (' ',1P,E10.2,' DEPTH OF BASE OF ASTHENOSPHERE') C READ (IUNIT7,952) PLTREF 952 FORMAT(A2) WRITE (IUNIT8,954) PLTREF 954 FORMAT(' ',A2,9X,'PLTREF: PLATE DEFINING VELOCITY ', + 'REFERENCE FRAME (EU, NA, AF, ...)') IPVREF=0 DO 956 I=1,NUMPLT IF (NAMES(I).EQ.PLTREF) IPVREF=I 956 CONTINUE IF (IPVREF.EQ.0) THEN WRITE (IUNIT8,958) (NAMES(I),I=1,NUMPLT) 958 FORMAT (/' *** ERR','OR in parameter input file: ***' + /' In line 13 (after ZBASTH, before ICONVE),' + /' in the first two columns of the line,' + /' define the velocity reference frame by' + /' entering one of the following plate names:' + /' ',26(A2,1X)) CALL PAUSE() STOP END IF C READ (IUNIT7,*) ICONVE WRITE (IUNIT8,96) ICONVE 96 FORMAT (' ',I10,' ICONVE: CODE FOR LOWER MANTLE FLOW ', + '(0=NONE)') READ (IUNIT7,*) TRHMAX WRITE (IUNIT8,101) TRHMAX 101 FORMAT (' ',1P,E10.2,' LIMIT ON HORIZONTAL TRACTIONS', + ' APPLIED TO BASE OF PLATE') C CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,VECTOR) TAUMAX(1)=VECTOR(1) TAUMAX(2)=VECTOR(2) C PROVIDE FOR OLD PARAMETER FILES WITH ONLY ONE TAUMAX: IF (TAUMAX(2).LE.0.0) TAUMAX(2)=TAUMAX(1) WRITE (IUNIT8,106) TAUMAX(1),TAUMAX(2) 106 FORMAT (' ',1P,E7.1,',',E7.1, + ' TAUMAX: SEA/LAND UPPER LIMIT', + ' ON INTEGRATED SUBDUCTION DRAG (N/m)') IF ((TAUMAX(1).LT.0.).OR.(TAUMAX(2).LT.0.)) THEN WRITE (IUNIT8,107) 107 FORMAT (/' *** ERR','OR in parameter input file: ***' + /' TAUMAX may not be negative.') CALL PAUSE() STOP END IF C READ (IUNIT7,*) RHOH2O WRITE (IUNIT8,110) RHOH2O 110 FORMAT (' ',1P,E10.3,' DENSITY OF GROUNDWATER, LAKES, & OCEANS') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,RHOBAR) WRITE (IUNIT8,120) RHOBAR(1),RHOBAR(2) 120 FORMAT (' ',1P,E10.3,'/',E10.3,' MEAN DENSITY,', + ' CORRECTED TO 0 DEGREES KELVIN. (CRUST/MANTLE)') READ (IUNIT7,*) RHOAST WRITE (IUNIT8,130) RHOAST 130 FORMAT (' ',1P,E10.3,' DENSITY OF ASTHENOSPHERE') READ (IUNIT7,*) GMEAN WRITE (IUNIT8,140) GMEAN 140 FORMAT (' ',1P,E10.3,' MEAN GRAVITATIONAL ACCELERATION', + ' (LENGTH/SEC**2)') READ (IUNIT7,*) ONEKM WRITE (IUNIT8,150) ONEKM 150 FORMAT (' ',1P,E10.3,' NUMBER OF LENGTH UNITS NEEDED TO', + ' MAKE 1 KILOMETER'/11X, + ' (E.G., 1000. IN SI, 1.E5 IN CGS)') READ (IUNIT7,*) RDETH WRITE (IUNIT8,155) RDETH 155 FORMAT (' ',1P,E10.3,' RADIUS OF THE EARTH') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,ALPHAT) WRITE (IUNIT8,160) ALPHAT(1),ALPHAT(2) 160 FORMAT (' ',1P,E10.2,'/',E10.2,' VOLUMETERIC THERMAL ', + 'EXPANSION OF CRUST', + ' (1/VOL)*(D.VOL/D.T). (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,CONDUC) WRITE (IUNIT8,170) CONDUC(1),CONDUC(2) 170 FORMAT (' ',1P,E10.2,'/',E10.2,' THERMAL CONDUCTIVITY, ENERGY/', + 'LENGTH/SEC/DEG. (CRUST/MANTLE)') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,RADIO) WRITE (IUNIT8,180) RADIO(1),RADIO(2) 180 FORMAT (' ',1P,E10.2,'/',E10.2,' RADIOACTIVE HEAT PRODUCTION', + ' ENERGY/VOLUME/SEC. (CRUST/MANTLE)') READ (IUNIT7,*) TSURF WRITE (IUNIT8,185) TSURF 185 FORMAT (' ', F10.0,' SURFACE TEMPERATURE, ON', + ' ABSOLUTE SCALE') CALL READN (INPUT,IUNIT7,IUNIT8,2, + OUTPUT,TEMLIM) WRITE (IUNIT8,190) TEMLIM(1),TEMLIM(2) 190 FORMAT (' ', F10.0,'/',F10.0,' CONVECTING TEMPERATURE (TMAX), ON' + ,' ABSOLUTE SCALE. (CRUST/MANTLE)') WRITE (IUNIT8,199) 199 FORMAT (/' ========== TACTICAL PARAMETERS (HOW TO REACH ', + 'THE SOLUTION) =========='/) READ (IUNIT7,*) MAXITR WRITE (IUNIT8,200) MAXITR 200 FORMAT (' ',I10,' MAXIMUM ITERATIONS WITHIN VELOCITY SOLUTION') READ (IUNIT7,*) OKTOQT WRITE (IUNIT8,210) OKTOQT 210 FORMAT (' ',F10.6,' ACCEPTABLE FRACTIONAL CHANGE IN VELOCITY ', + '(STOPS ITERATION EARLY)') READ (IUNIT7,*) REFSTR WRITE (IUNIT8,220) REFSTR 220 FORMAT (' ',1P,E10.2,' EXPECTED MEAN VALUE OF SHEAR STRESS IN', + ' CRUST'/' ',10X, + ' (USED TO INITIALIZE AND SET STIFFNESS LIMITS)') READ (IUNIT7,*) OKDELV WRITE (IUNIT8,230) OKDELV 230 FORMAT (' ',1P,E10.2,' MAGNITUDE OF VELOCITY ERR0RS ALLOWED', + ' DUE TO FINITE STIFFNESS'/11X, + '(SUCH ERR0RS MAY APPEAR IN SUCH FORMS AS:'/11X, + ' 1. FICTICIOUS BASAL SLIP OF CRUST OVER MANTLE'/11X, + ' 2. ERRONEOUS CONVERGENCE/DIVERGENCE AT VERTICAL FAULTS'/ + 11X, + ' 3. VELOCITY EFFECT OF FICTICIOUS VISCOUS COMPLIANCES'/11X, + ' HOWEVER, VALUES WHICH ARE TOO SMALL WILL CAUSE ILL-CONDITIONED' + /11X, + ' LINEAR SYSTEMS AND STRESS ERR0RS, ', + 'AND MAY PREVENT CONVERGENCE!)') READ (IUNIT7,*) EVERYP WRITE (IUNIT8,240) EVERYP 240 FORMAT (' ',L10,' SHOULD NODAL VELOCITIES BE OUTPUT EVERY STE', + 'P? (FOR CONVERGENCE STUDIES)') WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE FANGLS (INPUT, PHI,THETA, + OUTPUT, FANGLE) C C CALCULATE THE ARGUMENTS (ANGLES COUNTERCLOCKWISE FROM +THETA) C AT BOTH ENDS OF AN ARC OF A GREAT CIRCLE. C DOUBLE PRECISION FPOINT COMMON /SFAULT/ FPOINT DIMENSION FANGLE(2),FPOINT(7),PHI(2),THETA(2) DG180=3.141592654 A1=SIN(THETA(1))*COS(PHI(1)) A2=SIN(THETA(1))*SIN(PHI(1)) A3=COS(THETA(1)) B1=SIN(THETA(2))*COS(PHI(2)) B2=SIN(THETA(2))*SIN(PHI(2)) B3=COS(THETA(2)) S=0.99 XX=S*A1+(1.0-S)*B1 YY=S*A2+(1.0-S)*B2 ZZ=S*A3+(1.0-S)*B3 XVAL=SQRT(XX*XX+YY*YY+ZZ*ZZ) XX=XX/XVAL YY=YY/XVAL ZZ=ZZ/XVAL DX=XX-A1 DY=YY-A2 DZ=ZZ-A3 SITA=THETA(1) PHAI=PHI(1) S1=COS(SITA)*COS(PHAI) S2=COS(SITA)*SIN(PHAI) S3=-SIN(SITA) P1=-SIN(PHAI) P2=COS(PHAI) DXX=DX*S1+DY*S2+DZ*S3 DYY=DX*P1+DY*P2 FANGLE(1)=ATAN2F(DYY,DXX) S=0.01 XX=S*A1+(1.0-S)*B1 YY=S*A2+(1.0-S)*B2 ZZ=S*A3+(1.0-S)*B3 XVAL=SQRT(XX*XX+YY*YY+ZZ*ZZ) XX=XX/XVAL YY=YY/XVAL ZZ=ZZ/XVAL DX=B1-XX DY=B2-YY DZ=B3-ZZ SITA=ACOS(ZZ) PHAI=ATAN2F(YY,XX) IF(PHAI.LT.0.0) PHAI=2.0*DG180+PHAI S1=COS(SITA)*COS(PHAI) S2=COS(SITA)*SIN(PHAI) S3=-SIN(SITA) P1=-SIN(PHAI) P2=COS(PHAI) DXX=DX*S1+DY*S2+DZ*S3 DYY=DX*P1+DY*P2 FANGLE(2)=ATAN2F(DYY,DXX) RETURN END C C C SUBROUTINE FNODAL (INPUT,IUNITT,MXFEL,MXNODE,NFL,NODEF, + XNODE,YNODE, + OUTPUT,FPFLT) C C CALCULATES VECTOR NODAL FUNCTIONS AT INTERPOLATION POINT ON A C GREAT CIRCLE FAULT ELEMENT C DOUBLE PRECISION FPHI DIMENSION FPHI(4,7),FPFLT(2,2,2,7,MXFEL),FPP(2,2,2,7), + NODEF(4,MXFEL),PHI(2),THETA(2), + XNODE(MXNODE),YNODE(MXNODE) COMMON /FPHIS/ FPHI C SNCCOP: SIN(SITA)*COS(PHAI) AT INTERPOLATION POINT C SNCSNP SIN(SITA)*SIN(PHAI) AT INTERPOLATION POINT C DO 900 I=1,NFL N1=NODEF(1,I) N2=NODEF(2,I) THETA(1)=XNODE(N1) THETA(2)=XNODE(N2) PHI(1)=YNODE(N1) PHI(2)=YNODE(N2) CALL SNODAL (INPUT,PHI,THETA, + OUTPUT,FPP) DO 800 M=1,7 DO 500 J=1,2 DO 400 K=1,2 DO 300 L=1,2 FPFLT(L,K,J,M,I)=FPP(L,K,J,M) 300 CONTINUE 400 CONTINUE 500 CONTINUE SNCCOP=FPHI(1,M)*SIN(THETA(1))*COS(PHI(1))+ + FPHI(2,M)*SIN(THETA(2))*COS(PHI(2)) SNCSNP=FPHI(1,M)*SIN(THETA(1))*SIN(PHI(1))+ + FPHI(2,M)*SIN(THETA(2))*SIN(PHI(2)) COSM =FPHI(1,M)*COS(THETA(1))+FPHI(2,M)*COS(THETA(2)) PP=SQRT(SNCCOP*SNCCOP+SNCSNP*SNCSNP+COSM*COSM) COSM=COSM/PP SITA=ACOS(COSM) IF(SITA.LE.0.0.OR.SITA.GE.3.141592654) THEN SITA=SITA*57.29577951 WRITE(IUNITT,220) M,I,SITA 220 FORMAT('LATITUDE OF INTEGRATION POINT',I5, + ' OF FAULT ELEMENT', + I5,' IS OUT RANGE', + E14.4) ENDIF 800 CONTINUE 900 CONTINUE RETURN END C C C SUBROUTINE SNODAL (INPUT,PHI,THETA, + OUTPUT,FPP) C C CALCULATES VECTOR NODAL FUNCTION AT INTERPOLATION POINT ON A C GREAT CIRCLE OF A SIDE OF A FINITE ELEMENT C DOUBLE PRECISION FPHI,PP DIMENSION FPHI(4,7),FPP(2,2,2,7),PHI(2),THETA(2) COMMON /FPHIS/ FPHI C SNCCOP: SIN(SITA)*COS(PHAI) AT INTERPOLATION POINT C SNCSNP SIN(SITA)*SIN(PHAI) AT INTERPOLATION POINT C X1=SIN(THETA(1))*COS(PHI(1)) Y1=SIN(THETA(1))*SIN(PHI(1)) Z1=COS(THETA(1)) X2=SIN(THETA(2))*COS(PHI(2)) Y2=SIN(THETA(2))*SIN(PHI(2)) Z2=COS(THETA(2)) XN=X1+X2 YN=Y1+Y2 ZN=Z1+Z2 XYZN=SQRT(XN*XN+YN*YN+ZN*ZN) XN=XN/XYZN YN=YN/XYZN ZN=ZN/XYZN DD=X1*XN+Y1*YN+Z1*ZN DO 800 M=1,7 XX=FPHI(1,M)*X1+FPHI(2,M)*X2 YY=FPHI(1,M)*Y1+FPHI(2,M)*Y2 ZZ=FPHI(1,M)*Z1+FPHI(2,M)*Z2 PP=SQRT(XX*XX+YY*YY+ZZ*ZZ) XX=XX/PP YY=YY/PP ZZ=ZZ/PP SITAJ=ACOS(ZZ) PHAIJ=ATAN2F(YY,XX) RN=XX*XN+YY*YN+ZZ*ZN PPM=RN/DD CSCS=COS(SITAJ)*COS(PHAIJ) CSSN=COS(SITAJ)*SIN(PHAIJ) SNSN=SIN(SITAJ)*SIN(PHAIJ) SNC=SIN(SITAJ) SNP=SIN(PHAIJ) CSP=COS(PHAIJ) DO 500 J=1,2 FP=FPHI(J,M)*PPM FPP(1,1,J,M)=( COS(THETA(J))*COS(PHI(J))*CSCS + +COS(THETA(J))*SIN(PHI(J))*CSSN + +SIN(THETA(J))*SNC)*FP FPP(2,1,J,M)=(-SIN(PHI(J))*CSCS+COS(PHI(J))*CSSN)*FP FPP(1,2,J,M)=(-COS(THETA(J))*COS(PHI(J))*SNP + +COS(THETA(J))*SIN(PHI(J))*CSP)*FP FPP(2,2,J,M)=( SIN(PHI(J))*SNP + +COS(PHI(J))*CSP)*FP 500 CONTINUE 800 CONTINUE RETURN END C C C SUBROUTINE GETNET (INPUT,IUNIT7,IUNIT8, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + NFAKEN,NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000,OFFMAX,OFFSET, + TITLE1,XNODE,YNODE, + WORK,CHECKE,CHECKF,CHECKN) C C READ FINITE ELEMENT GRID FROM UNIT "IUNIT7". C ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8". C CHARACTER*80 TITLE1 LOGICAL ALLOK,BRIEF C C NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL CHECKE,CHECKF,CHECKN C DIMENSION CHECKE(MXEL),CHECKF(MXFEL),CHECKN(MXNODE), + DQDTDA(MXNODE),ELEDAT(3,MXEL),ELEV(MXNODE), + FDIP(2,MXFEL), NODEF(4,MXFEL), + NODES(3,MXEL),OFFSET(MXFEL), + XNODE(MXNODE),YNODE(MXNODE) DIMENSION DIPS(3),VECTOR(7) C WRITE (IUNIT8,1) IUNIT7 1 FORMAT (/ /' Attempting to read finite element grid from unit', + I3/) READ (IUNIT7,2) TITLE1 2 FORMAT (A80) WRITE (IUNIT8,3) TRIM(TITLE1) 3 FORMAT(/' Title of finite element grid ='/' ',A) C C READ NUMBER OF NODES. C INPUT NODAL LOCATIONS (X,Y), ELEVATIONS, & HEAT-FLOW (ONE RECORD C PER NODE). C (OPTION "BRIEF" SUPPRESSES MOST OUTPUT) C READ (IUNIT7,*) NUMNOD,NREALN,NFAKEN,N1000,BRIEF BRIEF=.TRUE. C IF (NUMNOD.NE.(NREALN+NFAKEN)) THEN WRITE (IUNIT8,4) NUMNOD,NREALN,NFAKEN 4 FORMAT (/' ERR0R: NUMNOD (',I6,') IS NOT EQUAL TO SUM' + /' OF NREALN (',I6,') AND NFAKEN (',I6,').') CALL PAUSE() STOP ENDIF C IF (NREALN.GT.N1000) THEN WRITE (IUNIT8,5) NREALN,N1000 5 FORMAT (/' ERR0R: NREALN (',I6,') IS GREATER THAN' + /' N1000 (',I6,').') CALL PAUSE() STOP ENDIF C IF (NUMNOD.GT.MXNODE) THEN WRITE (IUNIT8,10) NUMNOD 10 FORMAT(/' INCREASE PARAMETER MAXNOD TO BE AT LEAST' + /' THE NUMBER OF NODES (',I6,') AND RECOMPILE.') CALL PAUSE() STOP ENDIF C IF (BRIEF) THEN WRITE (IUNIT8,35) 35 FORMAT(/' (SINCE OPTION BRIEF=.TRUE., GRID WILL NOT BE ', + 'ECHOED HERE. BE CAREFUL!!!)') ELSE WRITE (IUNIT8,40) NUMNOD 40 FORMAT (/' THERE ARE',I5,' NODES IN THE GRID') WRITE (IUNIT8,50) 50 FORMAT (/ + ' NODE E-LONGITUDE N-LATITUDE', + ' THETA PHI ELEVATION', + ' HEAT-FLOW'/) ENDIF DO 90 K=1,NUMNOD CHECKN(K)=.FALSE. 90 CONTINUE DO 100 K=1,NUMNOD CALL READN (INPUT,IUNIT7,IUNIT8,5, + OUTPUT,VECTOR) INDEX=VECTOR(1)+0.5 IF (INDEX.GT.NREALN) THEN IF ((INDEX.LE.N1000).OR. + (INDEX.GT.(N1000+NFAKEN))) THEN WRITE (IUNIT8,91) INDEX 91 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER: ',I6) CALL PAUSE() STOP ENDIF ENDIF PLON=VECTOR(2) PLAT=VECTOR(3) IF (ABS(PLAT).GT.90.01) THEN WRITE (IUNIT8,92) INDEX 92 FORMAT (/' ERR0R: ABS(LATITUDE) > 90 AT NODE ',I6) CALL PAUSE() STOP ENDIF IF (ABS(PLAT).GT.89.99) THEN WRITE (IUNIT8,93) INDEX 93 FORMAT (/' ERR0R: NODE ',I6,' LIES ON A POLE.' + /' THIS IS A SINGULAR POINT OF THE' + ,' SPHERICAL COORDINATE SYSTEM.' + /' MOVE THIS NODE, AT LEAST SLIGHTLY.') CALL PAUSE() STOP ENDIF XI=(90.0-PLAT)*0.017453292 YI=PLON*0.017453292 ELEVI=VECTOR(4) QI=VECTOR(5) IF (INDEX.LE.NREALN) THEN I=INDEX ELSE I=NREALN+INDEX-N1000 ENDIF CHECKN(I)=.TRUE. XNODE(I)=XI YNODE(I)=YI ELEV(I)=ELEVI DQDTDA(I)=QI IF (QI.LT.0.) THEN WRITE (IUNIT8,96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL PAUSE() STOP ENDIF IF (.NOT.BRIEF) THEN WRITE (IUNIT8,99) INDEX,PLON,PLAT,XI,YI,ELEVI,QI 99 FORMAT (' ',I10,0P,2F12.3,2F11.5,1P,2E10.2) ENDIF 100 CONTINUE ALLOK=.TRUE. DO 101 I=1,NUMNOD ALLOK=ALLOK.AND.CHECKN(I) 101 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,102) 102 FORMAT(' THE FOLLOWING NODES WERE NEVER READ:') DO 104 I=1,NUMNOD IF (I.LE.NREALN) THEN INDEX=I ELSE INDEX=N1000+I-NREALN ENDIF IF (.NOT.CHECKN(I)) WRITE(IUNIT8,103)INDEX 103 FORMAT (' ',36X,I6) 104 CONTINUE CALL PAUSE() STOP ENDIF C C READ TRIANGULAR ELEMENTS C READ (IUNIT7,*) NUMEL IF (NUMEL.GT.MXEL) THEN WRITE (IUNIT8,108) NUMEL 108 FORMAT(/' INCREASE PARAMETER MAXEL TO BE AT LEAST EQUAL' + /' TO THE NUMBER OF ELEMENTS (',I6,') AND RECOMPILE.') CALL PAUSE() STOP ENDIF DO 109 K=1,NUMEL CHECKE(K)=.FALSE. 109 CONTINUE IF (.NOT.BRIEF) WRITE (IUNIT8,110) NUMEL 110 FORMAT(/' THERE ARE ',I6,' TRIANGULAR CONTINUUM ELEMENTS.'/ + ' (NODE NUMBERS FOR EACH ARE GIVEN CORNERS, COUNTER', + 'CLOCKWISE'/ / + ' ELEMENT C1 C2 C3') DO 200 K=1,NUMEL C (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) C (READN METHOD IS USED TO ALLOW FOR POSSIBLE ELEMENT DATA.) CALL READN (INPUT,IUNIT7,IUNIT8,7, + OUTPUT,VECTOR) I=NINT(VECTOR(1)) IF ((I.LT.1).OR.(I.GT.NUMEL)) THEN WRITE (IUNIT8,111) I 111 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I6) CALL PAUSE() STOP ENDIF NODES(1,I)=NINT(VECTOR(2)) NODES(2,I)=NINT(VECTOR(3)) NODES(3,I)=NINT(VECTOR(4)) ELEDAT(1,I)=VECTOR(5) ELEDAT(2,I)=VECTOR(6) ELEDAT(3,I)=VECTOR(7) CHECKE(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,120) I,(NODES(J,I),J=1,3), + (ELEDAT(J,I),J=1,3) 120 FORMAT (' ',I6,':',3I10,1P,E8.1,0P,F7.1,1P,E8.1) DO 130 J=1,3 N=NODES(J,I) IF (N.GT.NREALN) N=NREALN+(N-N1000) IF ((N.LE.0).OR.(N.GT.NUMNOD)) THEN WRITE (IUNIT8,125) NODES(J,I) 125 FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.') CALL PAUSE() STOP ENDIF NODES(J,I)=N 130 CONTINUE 200 CONTINUE ALLOK=.TRUE. DO 201 I=1,NUMEL ALLOK=ALLOK.AND.CHECKE(I) 201 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,202) 202 FORMAT (' THE FOLLOWING ELEMENTS WERE NEVER READ:') DO 204 I=1,NUMEL IF (.NOT.CHECKE(I)) WRITE(IUNIT8,203)I 203 FORMAT (' ',39X,I6) 204 CONTINUE CALL PAUSE() STOP ENDIF C C READ FAULT ELEMENTS C READ (IUNIT7,*) NFL IF (NFL.GT.MXFEL) THEN WRITE (IUNIT8,220)NFL 220 FORMAT (/' INCREASE PARAMETER MAXFEL TO BE AT LEAST EQUAL' + /' TO THE NUMBER OF FAULTS (',I6,') AND RECOMPILE.') CALL PAUSE() STOP ENDIF OFFMAX=0. DO 222 I=1,NFL CHECKF(I)=.FALSE. 222 CONTINUE IF (.NOT.BRIEF) WRITE(IUNIT8,230) NFL 230 FORMAT(/ /' THERE ARE ',I6,' GREAT CIRCLE FAULT ELEMENTS.') IF ((.NOT.BRIEF).AND.(NFL.GT.0)) WRITE(IUNIT8,231) 231 FORMAT (/' (THE 4 NODE NUMBERS DEFINING EACH ELEMENT MUST BE', + ' IN A COUNTERCLOCKWISE ORDER:'/ + ' N1, AND N2 ARE IN LEFT-TO-RIGHT SEQUENCE ON THE', + ' NEAR SIDE,'/ + ' THEN N3 IS OPPOSITE N2, N4 IS OPPOSITE N1 '/, + ' (FAULT DIPS ARE GIVEN AT N1, N2, ', + ' IN DEGREES FROM HORIZONTAL;'/ + ' POSITIVE DIPS ARE TOWARD N1, AND N2, RESPECTIVELY, '/ + ' WHILE NEGATIVE DIPS ARE TOWARD N4, AND N3.)'/ + ' (THE ARGUMENT OF THE FAULT TRACE IS GIVEN AT N1 AND N2,'/ + ' IN DEGREES COUNTERCLOCKWISE FROM THE THETA AXIS.)'/ + ' OFFSET IS THE TOTAL PAST SLIP OF THE FAULT.'/ / + ' ELEMENT N1 N2 N3 N4 DIP1 DIP2', + ' OFFSET'/) 240 FORMAT (' ',I6,':',4I5,1X,2F6.1,1X,F9.0) DO 300 K=1,NFL OFF=0. READ(IUNIT7,*) I,(NODEF(J,K),J=1,4),(DIPS(L),L=1,2),OFF IF ((I.LT.1).OR.(I.GT.NFL)) THEN WRITE (IUNIT8,241) I 241 FORMAT (/' ERR0R: ILLEGAL FAULT NUMBER: ',I6) CALL PAUSE() STOP ENDIF CHECKF(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,240) I,(NODEF(J,I),J=1,4), + (DIPS(L),L=1,2),OFF DO 250 J=1,4 N=NODEF(J,I) IF (N.GT.NREALN) N=NREALN+(N-N1000) IF ((N.LE.0).OR.(N.GT.NUMNOD)) THEN WRITE (IUNIT8,243) NODEF(J,I),I 243 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER (',I6, + ') IN FAULT ',I6) CALL PAUSE() STOP ENDIF NODEF(J,I)=N 250 CONTINUE DO 260 L=1,2 IF (ABS(DIPS(L)).GT.90.) THEN WRITE(IUNIT8,252) DIPS(L) 252 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN', + ' RANGE OF -90. TO +90. DEGREES.'/ + ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE', + ' HORIZONAL;'/ + ' A + PREFIX (OR NONE) INDICATES A DIP', + ' TOWARD THE N1-N2 SIDE;'/ + ' A - PREFIX INDICATES A DIP TOWARD', + ' THE N4-N3 SIDE.)') CALL PAUSE() STOP ENDIF IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L) FDIP(L,I)=DIPS(L)*0.017453293 260 CONTINUE IF (OFF.LT.0.) THEN WRITE (IUNIT8,280) OFF 280 FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2, + ' FOR FAULT ELEMENT',I6/ + ' OFFSETS MAY NOT BE NEGATIVE.') CALL PAUSE() STOP ENDIF OFFSET(I)=OFF OFFMAX=MAX(OFFMAX,OFF) 300 CONTINUE ALLOK=.TRUE. DO 301 I=1,NFL ALLOK=ALLOK.AND.CHECKF(I) 301 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,302) 302 FORMAT(' THE FOLLOWING FAULTS WERE NEVER READ:') DO 304 I=1,NFL IF (.NOT.CHECKF(I)) WRITE(IUNIT8,303)I 303 FORMAT(' ',36X,I6) 304 CONTINUE CALL PAUSE() STOP ELSE IF (OFFMAX.GT.0.) THEN WRITE (IUNIT8,400) OFFMAX 400 FORMAT (/' GREATEST FAULT OFFSET READ WAS ',1P,E10.2) ELSE WRITE (IUNIT8,401) 401 FORMAT (/' SINCE FAULT OFFSETS ARE ALL ZERO,', + ' INPUT PARAMETER BYERLY WILL HAVE NO EFFECT.') ENDIF ENDIF IF (.NOT. BRIEF) WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE PUTNET (INPUT,IUNITO, + BRIEF,DQDTDA,ELEDAT,ELEV,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE, + chemical_Delta_rho_list, + cooling_curvature_list) C C WRITES FINITE ELEMENT GRID TO UNIT "IUNITO". C CHARACTER*80 TITLE1 LOGICAL ANYEMU,BRIEF C DIMENSION chemical_Delta_rho_list(MXNODE), + cooling_curvature_list(MXNODE), + DQDTDA(MXNODE),ELEDAT(3,MXEL),ELEV(MXNODE), + FDIP(2,MXFEL),NP(4), + NODEF(4,MXFEL),NODES(3,MXEL), + OFFSET(MXFEL),TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(2) C WRITE(*,"(/ /' READY TO CREATE OUTPUT .FEG FILE ON UNIT ',I3)") WRITE(*,"(' Please enter a (new) filename with .feg ending: ')") WRITE (IUNITO,1) TITLE1 1 FORMAT (A80) C WRITE (IUNITO,2) NUMNOD, NREALN,NFAKEN,N1000,BRIEF 2 FORMAT(4I8,L8,' (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)') C DO 100 I=1,NUMNOD IF (I.LE.NREALN) THEN IPRINT=I ELSE IPRINT=N1000+(I-NREALN) ENDIF PLAT=90.-XNODE(I)*57.29577951 PLON=YNODE(I)*57.29577951 WRITE (IUNITO,91) I,PLON,PLAT,ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I),chemical_Delta_rho_list(I), + cooling_curvature_list(I) 91 FORMAT (I6,0P,2F11.5,1P,6E10.2) 100 CONTINUE C WRITE (IUNITO,110) NUMEL 110 FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') C ANYEMU=.FALSE. DO 140 I=1,NUMEL DO 130 K=1,3 IF (ELEDAT(K,I).NE.0.0) THEN ANYEMU=.TRUE. GO TO 141 END IF 130 CONTINUE 140 CONTINUE 141 DO 200 I=1,NUMEL DO 150 K=1,3 IF (NODES(K,I).LE.NREALN) THEN NP(K)=NODES(K,I) ELSE NP(K)=N1000+(NODES(K,I)-NREALN) ENDIF 150 CONTINUE IF (ANYEMU) THEN WRITE (IUNITO,155) I,(NP(K),K=1,3),(ELEDAT(K,I),K=1,3) 155 FORMAT (I5,3I6,1P,E8.1,0P,F7.1,1P,E8.1) ELSE WRITE (IUNITO,160) I,(NP(K),K=1,3) 160 FORMAT (I5,3I6) END IF 200 CONTINUE C WRITE (IUNITO,210) NFL 210 FORMAT (I10,' (NFL = NUMBER OF CURVILINEAR FAULT ELEMENTS)') C DO 300 I=1,NFL DO 220 K=1,4 IF (NODEF(K,I).LE.NREALN) THEN NP(K)=NODEF(K,I) ELSE NP(K)=N1000+(NODEF(K,I)-NREALN) ENDIF 220 CONTINUE DO 230 K=1,2 DIPS(K)=FDIP(K,I) DIPS(K)=DIPS(K)*57.2958 IF (DIPS(K).GT.90.01) DIPS(K)=DIPS(K)-180. 230 CONTINUE WRITE (IUNITO,250) I,(NP(K),K=1,4), + (DIPS(K),K=1,2),OFFSET(I) 250 FORMAT (I5,4I6,2F6.1,1P,E10.2) 300 CONTINUE C RETURN END C C C SUBROUTINE PAUSE() IMPLICIT NONE WRITE(*,"(' Press [Enter]..'\)") READ(*,*) RETURN END