! NOTE: All of these lines are in the NeoKinema.f90 source code! !========================================================================= ! NeoKinema.f90 ! ! A code which uses Neotectonic Kinematic data (geodesy, fault slip ! rates, stress directions) as additional constraints on a viscous-shell ! model of the lithosphere with plate-rotation boundary conditions, ! and computes present long-term-average horizontal velocities ! and distributed permanent strain rates. ! ! in the Fortran 90 language ! (for version and date, search below for "version =") ! ! by Peter Bird ! Professor Emeritus ! Department of Earth and Space Sciences ! University of California ! Los Angeles, CA 90095-1567 ! ! contacts for Peter Bird: phone (310) 825-1126 ! pbird@ess.ucla.edu ! ! (c) Copyright 1997, 2001, 2002, 2003, 2007, 2008, 2010, 2012, 2014 by ! Peter Bird and the Regents of the University of California ! !========================================================================= ! ! 1. GENERAL DESCRIPTION: ! ------------------------ ! This code uses Neotectonic Kinematic data (geodesy, fault slip ! rates, stress directions) as additional constraints on a viscous-shell ! model of the lithosphere with plate-rotation boundary conditions, ! and computes present long-term-average horizontal velocities ! and fault slip rates, and also permanent strain rates between faults. ! Although information on faults is a major part of the input, ! a continuum finite element approximation is used to describe ! the velocity field. In output graphics (usually prepared with ! companion program NeoKineMap) fault activity can be suggested by ! overprinted icons and/or by overprinted actual traces with ! information about estimated slip-rates. (The plotted slip rates can ! be either a-priori model inputs, or a-posteriori model outputs.) ! Also, NeoKineMap is capable of displaying maps of velocity magnitude ! which show all fault discontinuities, using some kludgy, ad-hoc ! post-processing code. If this is not desired, the distributed ! or "smeared" offsets across faults can be plotted to rigorously ! represent the internal solution mechanism in NeoKinema. ! !========================================================================= ! ! 2. ALGORITHM: !--------------- ! ! The algorithm of NeoKinema version 2.0~2.1 was described in ! Appendix S1 to Liu & Bird [2008; Geophys. J. Int., 172(2), ! 779-797, doi: 10.1111/j.1365-246X.2007.03640.x ]. ! ! The same file is available at URL of: ! http://peterbird.name/oldFTP/NeoKinema/ ! ! The method by which most-compressive horizontal principal stress ! directions are interpolated on the sphere was given by: ! Bird & Li [1996), J. Geophys. Res., v. 101, #B3, 5435-5443]. ! ! Model predictions of geodetic data are optimized under the criteria ! of minimizing (P - D) N (P - D), where: ! P is the vector of predicted horizontal velocity components at benchmarks ! D is the vector of actual* horizontal velocity components at benchmarks ! N is the normal matrix, or inverse of the covariance matrix of D. ! ! It is not necessary to pre-compute the normal matrix N; if needed (v.3+) ! NeoKinema will compute this from the covariance matrix of D. ! The covariance matrix of D can be represented by the block-diagonal ! error ellipses (for each benchmark in isolation) contained in the .gps file, ! or it can be supplemented by a full covariance matrix read from ! a .gp2 file. In either case, terms may be added to the covariance matrix ! to indicate that the velocity frame of reference is free-floating. ! Therefore, it is not necessary to use .gps data files in the same ! reference frame as the boundary condition--it is only necessary that all ! the lines in the .gps data file have the SAME velocity reference frame. ! ! *As the iteration of the solution proceeds, the geodetic data are corrected ! by adding velocities equal to the estimated long-term average of coseismic ! displacements caused by all faults within the model domain. In order ! to make this correction, estimated depth ranges of locking are read ! in two places: {1} default values from the parameter input file; and ! {2} local per-fault values, if available, from the right-hand columns ! of the f*.nki file (where negative numbers signal ignorance, in which ! case the default values {1} are used). ! !========================================================================== ! ! 3. FILE-NAMING CONVENTIONS ! ! NOTE: Adhering to these conventions will make use of ! the associated graphics package NeoKineMap MUCH EASIER! !----------------------------------------------------------------- ! Some pre-existing file-name conventions are retained: ! ! .dig indicates DIGitized polylines, from Digitise ! (or from coordinate-transformation utility Projector), ! such as digitized fault traces, or coastlines, ! or political boundaries. ! (such files can be read by NeoKineMap and FiniteMap as well). ! ! .feg indicates a Finite Element Grid, from OrbWin or OrbWeaver, ! plus post-processing utility program OrbNumber ! (such files can be read by NeoKineMap and FiniteMap as well). ! ! v*.out indicates a Velocity Output file, with 3 lines of titles, ! followed by (v_South, v_East) for each node, in m/s ! (such files can be read by NeoKineMap and FiniteMap as well). ! v_[token].out is the long-term-average velocity field ! (comparable to fault offset rates and plate tectonics); and ! v_interseismic-[token].out is the short-term interseismic field ! (comparable to GPS and other geodetic data). ! ! .gps indicates a geodetic-velocity file, of the sort ! that can be plotted by FiniteMap and NeoKineMap, ! or converted to a new velocity reference frame by ! ReframeGPS. ! ! .gp2 indicates an appendix to a .gps file, giving the ! covariance matrix for the geodetic velocity components. ! ! New file formats created for NeoKinema have one of two suffixes: ! ! .nki for NeoKinema Input (must be created before running NeoKinema); ! ! .nko for NeoKinema Output (created by NeoKinema). ! ! Specific file types are indicated by the first letter: ! ! .nki (Input): ! ! p*.nki Parameters, such as data weights and default fault locking depths, ! and names of other input files. ! ! b*.nki Boundary conditions ! ! f*.nki Fault offset-rates, with uncertainties (possibly huge!) ! ! s*.nki Stress directions ! ! .nko (Output): ! ! t*.nko Text record of the progress and success of the run ! ! f*.nko Offset-rates of faults estimated by NeoKinema ! ! h*.nko Heave-rates of faults estimated by NeoKinema ! (In version 1.3+, slip rates are also given to ! support PROGRAM Long_Term_Seismicity.) ! ! s*.nko Stress directions interpolated by NeoKinema ! ! e*.nko Continuum strain-rates, excluding faulting ! ! g*.nko Geodetic velocities of benchmarks, with velocity ! reference frame (optionally) redefined by NeoKinema, ! and with estimated mean coseismic velocities added ! to produce estimated long-term-average velocities. ! ! The first line of your p*.nki Parameter file lists a "name token" ! which will be built into all output file names. For example, ! ! _test01 results in the creation of t_test01.nko ! ! 2001-01 results in the creation of t2001-01.nko ! ! It is suggested (but not required) that you build the same token ! into the names of any input files that you prepare, IF they are ! created anew for this particular model. (If they are likely to be ! re-used in many model runs, this is not advised because it requires ! keeping redundant copies of the file.) !========================================================================== ! ! 4. OVERVIEW OF INPUT DATASETS ! ! !--------------------------------------------------------- ! p*.nki short file of parameter values, such as the weight factors ! for fault-slip and continuum constraints, the names ! of all other input files, etc. REQUIRED <1> ! ! f*.nki Fault names, slip-senses, and offset rates (or priors). OPTIONAL <2> ! As of NeoKinema v.3, "bracket" limits on offset rates may be added. ! ! f*.dig Digitized fault traces, with serial numbers to ! match f*.nki. REQUIRED if f.nki is used. <3> ! As of NeoKinema version 3, additional information on fault dip ! and depth of the "trace" (for blind thrusts) may be included. ! ! s*.nki Most-compressive horizontal principal stress azimuths ! (from centroid moment tensors of earthquakes, ! from dikes and veins, or from cluster analysis of faults ! with slickensides. Or, less reliably, from clusters of ! folds). OPTIONAL <9> ! ! *.gps Geodetic velocities at benchmarks, with error ! ellipses. OPTIONAL <10> ! ! *.gp2 Covariance matrix for geodetic velocity components. <10> ! (Note: This file is OPTIONAL and can be omitted, even if ! a .gps file was provided.) ! ! *.feg Finite element grid of connected spherical-triangle ! elements on the Earth's spherical surface. REQUIRED <11> ! Filename must have extension ".feg" or ".FEG". ! Typically created with interactive utility program ! ORBWEAVE for PCs with DOS or Windows (any version). ! The file begins with an ordered list of nodes ! with East longitude, North lattitude for each. ! Following is a list of triangular elements, each ! defined by the numbers of the 3 corner nodes. ! (It is important that these grids have been processed ! through utility program ORBNUMBR to reduce bandwidth.) ! The last line of the file contains "0" to show that ! there are no fault elements in an .feg file for NeoKinema. ! ! b*.nki A boundary-conditions file indicating which nodes (at ! least 2!) have specified velocities. <12> ! ! Notes: It is permissible (but not recommended) to include ! computer, drive, and path information in filenames ! (up to a maximum of 80 characters total for each). ! In practice, I ADVISE YOU to create a separate folder ! for this computation. Then, run NeoKinema from this folder, ! using short, simple filenames. ! This method will also make it easier to run the graphic program ! NeoKineMap, which assumes that input files are grouped in one folder. ! If any of the input datasets is null, enter 'none' ! for the filename. (However, you must always provide ! parameters, finite element grid, and boundary conditions.) ! !============================================================================= ! ! 4. DATA PREPARATION / INPUT FILE FORMATS ! *** General Note: Do not exceed 132 characters in any line of any file! ! -------------------------------------------------------------------------------------- ! ! -Build Parameter file p*.nki containing these lines; comments to right are OK. ! (V----beginning each in column 1): ! _test01 [name token for this run of NeoKinema] ! 1.00E3 L0 = length of fault trace whose offset rate gets unit weight (in m) ! 4.00E8 A0 = area of continuum whose stiffness & isotropy get unit weight (in m**2) ! 40 refinements (for nonlinearity) allowed in this solution (try 20-40) ! 1.0E-15 mu_ = scalar measure of typical anelastic strain rates in continuum (/s) ! 3.2E-17 xi_ = small strain-rate increment, in /s (Not TOO small! Try 3.2E-17 /s.) ! 20. sigma, in degrees, of angle between [heave vectors of dip-slip faults] & [trace-normal direction] ! 6371. radius of planet, in kilometers ! 1. 12. minimum and maximum default locking depths of intraplate faults, in km (may be overridden by f*.nki) ! 14. 40. minimum and maximum default locking depths of subduction zones, in km (may be overridden by f*.nki) ! FALSE switch: Do active faults give sigma_1h direction data? ! f_test01.nki filename of fault offset-rates (or, 'none') ! f_test01.dig filename of digitized fault traces (or, 'none') ! s_test01.nki filename of horizontal principal stress directions (or, 'none') ! WUSC002.gps filename of geodetic velocity data (or, 'none') ! WUSC002.gp2 filename of covariance matrix for geodetic velocity components (or, 'none') ! FALSE switch: Is velocity reference frame allowed to free-float? ! TRUE conservative_geodetic_adjustment? (uses geologic slip rates; not self-consistent) ! test01.feg filename of finite element grid (required) ! b_test01.nki filename of boundary conditions for finite element grid (required) ! NA plate defining velocity reference frame for type-4 boundary conditions ! FALSE dump_all_solutions? Creates velocity log file for studying convergence. ! ------------------------------------------------------------------------------------- ! ! -Build files f*.nki and s*.nki (both optional) as follows: ! A general feature of all of these files is that they are tables ! with one datum per line. ! To permit reading the table, there are two lines of headers before ! the actual data: ! -the Fortran 90 FORMAT needed to read the table, e.g., ! "(A6,1X,A50,F12.0,F12.0,L2,F6.1,F6.1,F6.1,F6.1)" for f*.nki ! (BE SURE TO INCLUDE formats for the 2 optional columns ! on the right, WHETHER ACTUALLY FILLED-IN OR NOT), ! or "(A8,1X,A8,1X,F8.3,F8.3,1X,I3,1X,A1)" in the case of s*.nki. ! [Special note on these 2 FORMATs; do *NOT* use format-repetition ! integers, as in "2F8.3"; instead, list each item separately, ! as in "F8.3,F8.3". This is necessary to avoid confusing my ! (rather primitive) code that will interpret the format: {1} in ! order to figure out whether your stress-quality indices are numeric ! or alphabetical; and {2} in order to change precision of output ! to 0.001 mm/a in f*.nko. ! -human-readable abbreviated column headings ! (which will be copied into the output files). ! NOTE: Create FORMATs carefully. Be sure that all input numbers ! have an included decimal point to guard against decimal point ! misalignment if the format is changed! If there is a chance ! that some input values do NOT have an included decimal, then ! be sure they are right-aligned, and use a corresponding format ! item like F12.0 or F6.0 (but not F12.3, for example). ! ! After these two header lines, the contents are different in ! each case: ! ! f*.nki = Fault offset-rate components ! For each fault of regional scale: ! 1. Identifier of fault and slip component*, consisting of 6 bytes: ! - first byte is always "F" (in column 1); ! - bytes #2-#5 are a 4-digit integer, used to locate ! the digitized fault trace within a .DIG file which contains ! many traces. Use leading 0's if number is less than 1000. ! - byte #6 determines the sense and component of slip ! whose rate is provided in this datum: ! R = Right-lateral (dextral) heave rate; ! L = Left-lateral (sinistral); ! D = Divergent or Detachment (extensional) heave rate; ! N = Normal (extensional) throw rate; ! P = thrust Plate or naPpe (convergent) heave rate; ! T = Thrust (convergent) throw rate; ! S = Subduction (convergent) heave rate. ! (*Note: More than one component can be entered for an oblique-slip ! fault, simply by using one line for the strike-slip component (R or L) ! and another line for the dip-slip component (D, N, P, T, or S). ! To reduce confusion, it is probably best to enter these two lines ! as adjacent data, but NeoKinema does not actually require this.) ! It is NOT permitted to enter more than one strike-slip component ! or more than one dip-slip component for any single fault trace. ! 2. Descriptive text with fault name, location, ... ! 3. Offset rate (slip, throw, or heave rate), in mm/year. ! Always a positive number, or zero. ! In the case of thrust (T) or normal (N) faulting, this number ! is the relative vertical offset (throw) rate ! of the two sides, which can usually be measured more ! accurately than slip rate. ! In all other cases (D, N, P, T, S) this number gives the ! relative horizontal velocity (heave rate) component ! perpendicular to the fault trace. ! 4. Standard error (sigma_; 68%-confidence ; half of 95%- ! confidence) of offset rate, in mm/year. May be large. ! However, may NOT be zero or negative. ! NOTE: It is perfectly reasonable to enter an offset rate of "0.0" ! and a large uncertainty (e.g., "49.9") for any fault for which ! you have no information. Then, this fault will be free to slip ! at any rate that will improve the fit of other data. (However, ! it will also be free to slip in the "wrong" direction! ! 5. Creeping? (T/F). If this fault is known to be creeping, at more ! than half of its long-term rate, then insert "T" in this column. ! This will indicate that nearby geodetic velocities should NOT ! be corrected for elastic strain accumulation along this fault. ! (This logical flag will also be copied to the h_token_nko ! output file, where it will eliminate any seismicity footprint ! for this fault in runs of PROGRAM Long_Term_Seismicity.) ! Otherwise, insert "F" for the normal case of a stick-slip fault. ! 6. Minimum/upper limit of seismogenic locked patch, in km of depth. ! Use a negative number if not known ("-1.0") and NeoKinema will ! substitute the appropriate value from the parameter input file ! p*.nki. {This entry is not used if Creeping? flag = T.} ! 7. Maximum/lower limit of seismogenic locked patch, in km of depth. ! Use a negative number if not known ("-1.0") and NeoKinema will ! substitute the appropriate value from the parameter input file ! p*.nki. {This entry is not used if Creeping? flag = T.} ! - - - - OPTIONAL additional columns added in NeoKinema version 3.01: - - - ! 8. Lower limit on offset rate in column 3, in mm/a. Must be lower. ! 9. Upper limit on offset rate in column 3, in mm/a. Must be higher. ! [Note: It is permissible to include a lower limit in column 8 but ! omit the upper limit. You may even wish to give a lower limit ! of "0.0" for all faults, to prevent them slipping in the ! wrong direction. However, if you wish to specify only ! an upper limit, you will need to enter a dummy place-holder ! lower limit, which can be very negative (e.g., "-999.").] ! ! s*.nki = most-compressive horizontal principal Stress directions ! For each assemblage of stress indicators: ! 1. Text-string #1, which could be a reference ! in short form (e.g., "Jones, 1991")? ! 2. Text-string #2, which might give location and state abbreviation? ! 3. East longitude, in decimal degrees (e.g., -106.92) ! 4. North latitude, in decimal degrees (e.g., 43.81) ! 5. Azimuth, measured clockwise from North, in degrees; ! may be either integer (I format) or real (F/E/D format). ! 6. Standard deviation (sigma_; +- for 68%-confidence; half of +- for ! 95%-confidence) of stress azimuth, in degrees, integer or real, ! OR: ! letter quality index (A, B, C, D, E) from World Stress Map ! (in which case the corresponding FORMAT item #7 is an A, not I or F). ! Note that stress-direction data should not be restricted to those ! within the area of the F-E grid (.feg file). Additional data up to ! 22 degrees arc distance from the grid edge can provide additional ! constraints on interpolated stress directions, and should be provided ! as input if available. (In this respect, NeoKinema differs from ! Restore, which is not able to utilize stress directions outside the ! .feg grid area.) ! ! -f*.dig, the file of digitised fault traces, must have some combination ! of the following formats (basic, or more complex): ! V---(column 1). [16 sample lines of f*.dig follow this line.] ! F0489 ! -1.11661E+02,+3.92411E+01 ! -1.11636E+02,+3.92543E+01 ! -1.11615E+02,+3.92747E+01 ! -1.11592E+02,+3.92896E+01 ! -1.11570E+02,+3.93063E+01 ! -1.11549E+02,+3.93249E+01 ! -1.11525E+02,+3.93436E+01 ! -1.11503E+02,+3.93622E+01 ! -1.11484E+02,+3.93824E+01 ! -1.11462E+02,+3.94046E+01 ! -1.11443E+02,+3.94248E+01 ! -1.11434E+02,+3.94274E+01 ! *** end of line segment *** ! F0532RT Doozy dextral transpressive fault, MI ! dig_degrees 72.1 ! -1.12405E+02,+3.91402E+01 ! -1.123795+02,+3.91117E+01 ! -1.12365E+02,+3.90894E+01 ! -1.12355E+02,+3.90844E+01 ! -1.12317E+02,+3.90801E+01 ! [... and so on...] ! Each fault trace must be introduced by a label line written ! by WRITE (nn,"('F',I4)") fltnum ! where the INTEGER :: fltnum is used to tie the trace to ! data in file f*.nki. ! (If you wish, you can follow the F1234 number with one or two ! bytes indicating the sense of slip on the fault, such as ! F1234R or F1234RT. These notations are IGNORED by NeoKinema, ! which takes the sense of each slip component from the f*.nki ! file. However, these notations are USED by NeoKineMap for ! plotting fault traces in the appropriate color(s).) ! You are also encouraged to include the fault name on this line. ! Beginning with NeoKinema v3.02, there may be an optional ! second header line giving the fault dip, in any of these formats: ! "dip_degrees 45" or "dip_degrees 45.7" or "dip_degrees45" or ! "dip_degess45.7". What you can NOT enter is: ! "dip_degrees45SomethingElse" because there must be white-space ! immediately to the right of the dip number. (White-space ! characters include the new-line bytes CR and/or LF.) ! The following lines must give (longitude, latitude) of each ! digitised point along the trace, in degrees, according to: ! FORMAT(1X,SP,1P,E12.5,',',E12.5) ! Notice that the first column must be blank, and the leading ! + sign on positive numbers is required (use SP in FORMAT). ! East longitude is positive; West of Greenwich is negative. ! North latitude is positive; Southern hemisphere is negative. ! Each fault trace is concluded by a record with ! "*** end of line segment ***" starting in column 1. ! The order in which faults are digitised is unimportant. ! The number of points in each segment is unimportant ! (but there must be at least 2). ! The order in which they are numbered is unimportant, except ! that the 6-byte (6-character) identifiers must tie to F.DAT. ! The easiest way to create these files is to use my program ! DIGITISE (which accepts output from a digitiser through the ! serial port COM1 or COM2) and one of my map-projection ! utility programs like PROJECTOR which converts ! from the (x,y) coordinates of the map projection to ! (longitude, latitude) coordinates. ! ! *.gps = Geodetic-velocity data file. ! Each .gps file will have 3 lines of headers: ! ------------------------------------------------------ ! File name and source(s) ! FORTRAN format for reading the data lines that follow the headers ! [ column header labels, in a standard order: ] ! E_lon_deg N_lat_deg v_E_mmpa v_N_mmpa v_E_sigma v_N_sigma correlation frame identifier(s) ! ------------------------------------------------------------------------------------------ ! with the (obvious) meanings: ! E_lon_deg = longitude, in degrees from Greenwich meridian, with East positive ! N_lat_deg = latitude, in degrees from equator, with North positive ! v_E_mmpa = velocity component to East, in millimeters per year ! v_N_mmpa = velocity component to North, in millimeters per year ! v_E_sigma = standard deviation (1-sigma) of v_E_mmpa, also in mm/a ! v_N_sigma = standard deviation (1-sigma) of v_N_mmpa, also in mm/a ! correlation = coefficient of correlation between v_E_mmpa and v_N_mmpa ! reference_frame = reference frame for velocity, left-justified, limited to 15 bytes ! identifier(s) = optional station name and/or source reference, if a compilation ! Following these headers there is one line of data per benchmark. ! ! *.gp2 Covariance matrix for velocity components listed in the ! preceding .gps file. Order of degrees of freedom is: ! East component before North component for each benchmark, ! and benchmarks in same order as in the .gps file. ! Note that the .gp2 file is symmetrical, so it is sufficient ! to enter the diagonal and upper (or lower) triangle. ! Zero values do not need to be entered. ! Values can be entered in any order. ! Each line in the .gp2 file has: ! irow jcolumn variance ! where: irow is an integer ! jcolumn is an integer ! variance is a real number in units of (mm/a)**2. ! Note that variance is required to be positive for all diagonal entries; ! zero is not acceptable on the diagonal. ! If "variance" does not agree with information from the .gps file, ! then information from the .gps file is over-written and replaced. ! -------------------------------------------------------- ! ! -Build the *.feg file with ORBWEAVE (using no fault elements), and ! renumber for minimum bandwidth with ORBNUMBR. ! -------------------------------------------------------- ! ! -Build the b*.nki file to go with *.feg and to specify which nodes ! have fixed velocities (at least 2!). This file begins with a title line; ! the bulk of the file is read with the Fortran list-directed ! input method, so column spacing and alignment are not important, ! but plate names (if used) must be inside quotation marks. ! Notice that my post-processing utility program OrbNumber, ! applied to any NeoKinema-type .FEG grid file, will produce ! an ordered list of boundary nodes; you can use this as the ! basis for your b*.nki file. ! For each node whose velocity you want to fix ! provide one line of b*.nki with EITHER: ! ! ordinal_integer node_number 2 velocity_in_m/s azimuth_in_degrees_clockwise_from_North, e.g.: ! 1 431 2 3.109E-09 90.0 ! ! --== OR ==-- ! ordinal_integer node_number 4 plate_ID_for_this_node, e.g.: ! 2 432 4 "PA" ! ! The leading ordinal_integer can be anything you like, but may not ! be omitted; I use it to count the boundary nodes. ! The following node number must refer to a boundary node in *.feg. ! The node numbers you use must be the NEW node numbers ! assigned by my renumbering utility OrbNumber; in order to see ! these, load the output *.feg file from OrbNumber back ! into OrbWin or OrbWeaver, select the Nodes command, and simply wave ! the mouse near the desired node to see its number displayed ! at the bottom center of the screen. ! The code "2" means that both components of velocity are fixed ! by the velocity magnitude and azimuth you have specified; ! velocities are in meters per second, in any reference frame you like! ! The azimuth of the velocity is measured in degrees, clockwise ! from local North. Negative (counterclockwise) azimuths are allowed. ! The code "4" means that velocities should be computed from the ! PB2002 global plate model of Bird [2003; Geochemistry, Geophysics, ! Geosystems]. (The reference frame for velocities was previously ! specified in the last line of the parameter input file.) ! If you wish, you may also list other boundary nodes with code "0" ! which will just leave them free. No other codes are supported now. ! If relative plate rotations are available for the ! surrounding plates, then it is certainly desirable to ! provide boundary conditions all the way around the grid! ! Otherwise, one stable interior side is fixed, and the ! rest left free to move as determined by the geologic data. ! A free edge may be appropriate to represent the edge of the ! overriding plate in a subduction zone. ! -------------------------------------------------------- !==================================================================== ! ! 5. OVERVIEW OF OUTPUT DATASETS: ! !---------------------------------------------- ! *Text dataset: t*.nko contains reports on progress. <21> ! The information is much the same as what you see on the screen, ! so do not be concerned if messages fly by too fast to read. ! This dataset is created in each run. Be careful- old t*.nko's ! will be overwritten if they are not moved or renamed! ! ! *Interpolated stresses: s*.nko <22> ! These are created to record the stress directions ! interpolated to each element center, for plotting by NeoKineMap. ! Produced only if stress directions were input, and/or if faults ! were designated as stress indicators. ! There are 4 columns in the file, in this order: ! ! element_number successful_interpolation? azimuth standard_deviation ! ! Element numbers refer to the .feg file used in the run of NeoKinema. ! Successful_interpolation is a logical value (T/F). ! Azimuth is in degrees clockwise from North. ! Standard deviation (sigma) is also in degrees. ! (If successful_interpolation = F, the following two values ! are meaningless, and should not be used in graphics.) ! ! *Model offset-rates for faults: f*.nko <23> ! Same format as f*.nki, but limits on locking depth are omitted. ! Instead two new columns are added at the right-hand side ! to show the model prediction, and prediction error in units of sigma. ! Rates are mean heave or throw rates depending on fault sense; all are in mm/a. ! This is the low-resolution form of fault offset-rate ! output; for the high-resolution heave rates, see h*.nko. ! Note that NeoKinema takes no special action to prevent faults ! from slipping in the wrong direction (although this becomes quite ! unlikely if the ratio of neotectonic rate to uncertainty exceeds ! 2-3). Slip on faults can reverse over time, and the sign of the ! neotectonic rate is not always known with certainty; one may wish ! to see what sense NeoKinema predicts. If "wrong-way" slipping ! on any fault is unacceptable, one can simply delete that fault's ! rate(s) from the f*.nki file, and run NeoKinema again. ! ! *Model heave-rates of faults: h*.nko <24> ! Each fault trace is divided into a number of short segments ! (a segment is the intersection of a fault trace with a triangular finite element), ! and the segments will appear in seemingly random order in the file. ! Each segment is described by a line in the file, whose initial (left-hand) part looks like: ! F9078T heave-rate = 1.878 mm/a in (-120.399, 34.460)-(-120.391, 34.460) = element 1887 ! This information is used by NeoKineMap to plot model fault heave rates in detail. ! The full output line contains additional information on the right, e.g.: ! F9078T heave-rate = 1.878 mm/a in (-120.399, 34.460)-(-120.391, 34.460) = element 1887 creeping?: F, slip-rate = 1.999 mm/a ! and this information on fault creep (T/F) and slip rate is used by ! PROGRAM Long_Term_Seismicity ! ! *Reframed, unlocked geodetic velocities: g*.nko <27>. ! If a *.gps file was provided as input, a g*.nko file will be written ! as output. It may be different from the input file in 3 ways: ! (1) Re-framing of the velocities (if allowed) into the velocity reference ! frame which allows them to best-fit the current NeoKinema solution. ! Note that errors will remain; the velocities in g*.nko will ! not generally be equal to the computed model velocities around ! them (although they should be close). ! (Note that this correction is optional, and only occurs if ! parameter "floating_frame" = .TRUE. in the parameter input file. ! (2) Deletion of any benchmarks which did not fall within the area ! of the finite-element grid in the current .feg file. ! (3) Addition of estimated long-term-average coseismic velocities ! to benchmarks, based on the current fault slip rates in the ! model, and the fault locking depths read from the parameter ! input file, or from the right-hand columns of the f*.nki file. ! (Note that this correction is also optional; if you set Creeps? = T ! for all faults in f*.nki, no such correction will be performed. ! Use this option where faults creep freely ! at all times, or where the authors of the *.gps file have already ! corrected the velocities to add estimated mean coseismic rates.) ! ! *Neotectonic (long-term-average) velocities: v*.out <25>. ! A dataset of velocities at the nodes of the finite-element ! grid "*.feg" is also produced. This permits ! plotting diagrams of the neotectonics using NeoKineMap (or FiniteMap). ! Maps produced from this output file are comparable to plate tectonics ! and long-term-average fault offset rates. ! ! *Interseismic (short-term) velocities: v_interseismic_*.out <25>. ! A dataset of velocities at the nodes of the finite-element ! grid "*.feg" is also produced. This permits ! plotting diagrams of the neotectonics using NeoKineMap (or FiniteMap). ! Maps produced from this output file are comparable to geodetic ! (e.g., GPS) velocities collected in years with no major earthquakes. ! ! *Distributed permanet strain-rates of the elements: e*.nko <26>. ! For every continuum finite element in x_feg, a logical value (T/F) ! indicates whether the strain-rate in that element contained a ! contribution from faulting. If so, T is shown, followed by the ! non-faulting or distributed part of the strain-rate tensor. ! If not, then F is shown, followed by the strain-rate tensor (that ! could just as well have been computed from v*.out). ! All strain-rate tensors are horizontal-plane only (2 x 2) ! and are in (theta = S, phi = E) coordinates, so that the ! 3 values given are: theta-theta or N-S, theta-phi or SE, phi-phi or E-W. ! Note: All strain-rates discussed here are permanent, or anelastic, ! strain-rates which are achieved by some combination of frictional ! plasticity at shallow depths and/or dislocation creep at greater ! depths. Thus, they are considered to be independent of time, ! unlike the elastic part of the strain-rate which reverses sense ! during each earthquake cycle, and tends to a long-term mean ! rate of zero. ! !========================================================================= version = 'NeoKinema: Version 3.1 of 6 May 2014' !GPBversion======================================================================= ! VERSION HISTORY ! Version 1.0, Summer 2002: released only to Zhen Liu (who used ! it for the preliminary modeling of the Persia-Tibet-Burma orogen). !---------------------------------------------------------------------------- ! *Eliminated redundant shadow strike-slip pseudo-datum for any fault ! trace which already has a user-specified strike-slip component. ! *Added damping during the updating of shadow strike-slip rate sigmas. ! *Sigma_perSecond, associated with enforcing no continuum shear ! straining on the stress principal axis directions ! was limited to be no less than xi_, in order to protect the ! condition number of the linear system and enhance reproducibility ! of iterated results (thus allowing for better convergence). ! *Input parameter dump_all_solutions added to control production of ! velocity-solution log file v_log.nko, which is intended for use ! with programs Analyze_Velocity_Evolution and FiniteMap. ! *Input parameter conservative_geodetic_adjustment added to provide ! an option in which the fault-unlocking velocity adjustment at ! geodetic benchmarks is based on the input/prior geologic slip rate, not ! the self-consistent current slip rate estimate. This may be necessary ! to preserve stability of the iteration of the solution, in cases of ! benchmarks sitting over two antithetic thrusts (or detachments) which ! dip together and intersect in the seismogenic depth range. !---------------------------------------------------------------------------- ! Version 1.1, 2003.10.06: also released to Zhen Liu (2nd-round Persia-Tibet-Burma). ! I also used it for reconnaissance (unpublished) Gorda-California-Nevada ! and SCEC sub-region models through September 2004. !---------------------------------------------------------------------------- ! *Fixed bug in code that was written to prevent ill-conditioning of ! linear system, by limiting highest eigenvalues: ! Previously, I subjected sigma_perSecond to a lower limit of xi_ ! to achieve this. However, when I later used the (inflated) ! sigma_perSecond to evaluate prediction error, the stress-direction ! errors were incorrectly biased downward. (The effect was largest ! for high values of xi_.) The fix is to leave sigma_perSecond ! unchanged, but introduce a new "effective_sigma_perSecond" to prevent ! ill-conditioning. The old, correct sigma_perSecond is used to assess ! the post-fit residuals (prediction errors) in continuum stress directions. !--------------------------------------------------------------------------------- ! Version 1.2, 2004.10.26: Tested fix for effectiveness, and found that overall ! change to best solutions was minor; kept on file. !--------------------------------------------------------------------------------- ! *Added label "!Write_h_token_nko" to help find code that writes this ! file. (It is not in a separate SUBROUTINE.) Altered this code to ! provide the additional REAL value slip_rate_mmpa for every dip-slip ! fault, and for every strike-slip fault in the input data file. ! However, no slip rate is provided for shadow strike-slip components ! on dip-slip faults; instead, this motion is merged with the dip-slip ! rate by Pythagorean theorom to give the total slip rates quoted. ! This information is needed by PROGRAM Long_Term_Seismicity. !--------------------------------------------------------------------------------- ! Version 1.3, 2004.10.26: used for 2004 SCEC Annual Report and 2005 SCEC proposal ! (N.B. These runs have erroneous labels of "NK v.2"!) !--------------------------------------------------------------------------------- ! *Altered comments to clarify that "D" and "N" are just 2 different ! ways of entering fault offset rates (as heave or throw components) for ! normal faults; and likewise, "P" and "T" are just 2 different ways ! of entering fault offset rates (as heave or throw components) for thrust ! faults, and these class distinctions do NOT imply different fault dips, ! or any other geometric or behavioral distinction. ! *Set compiled-in parameters of: ! REAL, PARAMETER :: normal_dip_degrees = 55.0 ! REAL, PARAMETER :: thrust_dip_degrees = 20.0 ! REAL, PARAMETER :: subduction_dip_degrees = 14.0 ! consistent with Bird & Kagan [2004] Table 5 values used in Long_Term_Seismicity, ! because these values are used in the Write_h_token_nko code section, ! and consistency between programs in a chain is highly desirable. !---------------------------------------------------------------------------- ! Version 1.4, 2004.11.22: used to test simple plate boundary cases with ! PROGRAM Long_Term_Seismicity for consistency with ! Bird & Kagan [2004, BSSA]. !---------------------------------------------------------------------------- ! *Eliminated input parameter "geodesy_weight" and set it internally to 1.0. ! The weight of one geodetic velocity component at one geodetic benchmark ! is now the standard unit in which other adjustable weights are measured. ! (N.B.: I do not consider a "weight" to include the associated factor ! of (sigma)**(-2), or any geometric factors specific to a datum site.) ! *Added new user-controlled dimensional weighting parameters: ! L0 = length of fault trace whose offset rate gets unit weight (in m) ! A0 = area of continuum whose stiffness & isotropy get unit weight (in m**2) ! and recoded to eliminate hidden dimensional factors of finite element size. !--------------------------------------------------------------------------------- ! Version 2.0, 2004.12.07: Used for final RELM/publication verion of the ! Gorda-California-Nevada orogen ! [Bird & Liu, 2007, Seismol. Res. Lett.]. ! Also used by Zhen Liu for most of the models ! (but not the final preferred model) ! in the Persia-Tibet-Burma orogen project ! [Liu & Bird, 2008, Geophys. J. Int.] !--------------------------------------------------------------------------------- ! *Modified Euler rotation-rate poles (all expressed relative to PA) ! for plates AM, AR, IN, SU: replaced the PB2002 [Bird, 2003, G^3] ! values with new values from last column of Table 1 in ! Liu & Bird [2008, Geophys. J. Int. = Persia-Tibet- ! Burma model]. Sources for this update (selected by Liu): ! -Sella et al. [2002] for AR, IN; ! -Kreemer et al. [2003] for AM, SU. !--------------------------------------------------------------------------------- ! Version 2.1, 2007.08.14 used for the final preferred Persia-Tibet-Burma model ! of Liu & Bird [2008, Geophys. J. Int.] !--------------------------------------------------------------------------------- ! *Added test of input stress-direction azimuths for range -360~+360, ! after having endless grief over undetected azimuths ! of "999" degrees in the World Stress Map dataset. ! *Subduction zone faults (S) are allowed to slip ! obliquely EVEN IF parameter "sigma_offnormal_degrees" ! permitting oblique slip on other dip-slip faults is ! set to a small value or zero. [2007.11.21] ! *Corrected FORMAT in subprogram Write_f_token_nki() to provide 0.001 mm/a precision ! in echoed datum offset rates and their standard deviations. [2008.01.10] ! *Modified input section to permit reading at most ONE strike-slip offset rate ! (sense R or L) and at most ONE dip-slip offset rate (sense D, N, P, S, or T) ! for any given fault trace (e.g., F4253). Note that allowing multiple rates ! would cause stacking of multiple virtual faults along the trace, rather than ! weighting of multiple opinions about the single fault, as the user might ! have intended! [2008.01.17] ! *Added two right-hand columns to f*.nki input file, after logical indicator ! for fault creep, with the upper/smaller seismogenic locking depth (in km), ! and the lower/deeper seismogenic locking depth (in km). ! These new columns are NOT optional; old input files must be upgraded by ! adding such columns! (However, negative values like "-1.0" can be used ! for unknown locking depths, which will then be replaced with the default ! values in the parameter input file p_*.nki. [2008.01.21] ! [Also, note that 2 MORE optional columns were later added; see v.3 below.] ! *Cosmetic changes in error-measure output: Old norms "L0, L1, L2" ! renamed to "N0, N1, N2" to avoid confusion with input "reference length" ! L0 used for weighting fault offset-rate data. ! "Fault" error measures are still reported, but now called "Offset-rate" errors. ! Also, printing of "Global" error measures is now suppressed, ! because with the addition of alternative "Potency" error measures ! to replace misleading "Fault" error measures, it is no longer ! clear how a "Global" error measure should be defined. ! *Added new "Potency-rate" error measures to final output. These are just ! the N0, N1, and N2 norms of [(model_offset_rate - datum_offset_rate)/ ! (datum_standard_error)] * [(model_potency_rate)/MEAN(model_potency_rates)]. ! The difference from the old "Faults" error measures (now called "Offset-rate" ! error measures) is the insertion of the second, potency-rate factor, ! in place of the simpler trace-length factor used in Offset-rate errors. ! [PB, 2008.01.30] !--------------------------------------------------------------------------------- ! Version 2.2, 2008.01.30: Used for: ! -Bird [2009, JGR] models of the western U.S. and adjacent offshore regions. ! -Howe & Bird [2010, Tectonophysics] models of the Alpine-Aegean orogen. !--------------------------------------------------------------------------------- ! *General revision of Euler poles used in type-4 velocity boundary conditions: ! -Adopted "Guadalupe" Euler pole for NA-PA from Gonzalez-Garcia et al. [2003] ! (endorsed by Bird [2009]). ! -Added CA-NA of DeMets et al. [2000] (endorsed by Mann et al. [2002]) to ! this revised NA-PA to get a better CA-PA pole. ! -Corrected BU-PA pole to be consistent with SU-PA of Kreemer [2003] ! (as in NeoKinema v. 2.0+) plus BU-SU of PB2002 [Bird, 2003]. ! {However, note that Robinson & Bird [2011?] later found this ! BU-SU rotation from PB2002 to be too fast, by about a factor of 2.} ! -Added capability to handle unexpected two-character plate abbreviations ! (such as "C?" or "C2" replacing "CA") in the boundary conditions ! b_*.nki file, in lines that specify type-4 boundary conditions. ! The user is now prompted to provide the Euler pole for any new ! (or revised) plate motions, relative to the current reference frame ! (variable "reference_plate_c2" of the input parameter p_*.nki file). ! This allows altering boundary velocities without recompiling ! NeoKinema. A record of user input is copied to the log file. ! -Added an automatically-generated table listing all Euler poles used in ! generating type-4 velocity boundary conditions in the current run. ! This table is also copied to the log file. ! *Added warning when any assigned uncertainty (sigma, standard deviation) of ! any fault offset rate exceeds 50 mm/a (which may cause ill-conditioning). ! User must manually override warning to continue, and a record is left ! in the log file. ! *Added warning(s) when any geodetic benchmark falls into the same finite ! element as a fast-moving fault (offset rate > 1 mm/a), because this ! is likely to cause solution errors of same order as the offset rate. ! There is both a prospective warning (based on input, prior offset rates) ! and a retrospective warning (based on posterior, output offset rates). ! After any prospective warning, user must choose whether to ignore warning. ! Both sets of warnings are copied to the log file. ! *Now allowing fault offset-senses and traces to be used as stress-direction ! pseudo-data (assuming user specifies faults_give_sigma_1h = .TRUE.) ! even when the target offset rate is zero. (Zero is often assigned ! where there is no slip-rate datum for a fault, but in fact ! the slip sense may be well-established.) ! *A minimum of 6 solution refinements by iteration is now enforced; ! this is a safety measure, because naieve users cannot always judge ! which types of problem require iteration for accurate solution. ! *Master SUBROUTINE Solve_for_vw now CALLs Prediction in every pass ! except the first. This solves a non-updating bug that appeared only ! in some artifically simple test problems. ! *Changed the computation of revised (long-term, fault-unlocked) geodetic ! velocity targets, so it is now based on fault locking from the surface ! down to the lower locking depth. (The upper locking depth is still ! read from the f_*.nki file and passed on to Long_Term_Seismicity ! for use in seismicity forecasts, but it is no longer used within ! NeoKinema.) The reason is that shallow strips of fault plane above ! the uppper (shallower) locking depth may creep a bit, and may not ! contribute much (or any) stress drop during earthquakes, but they creep ! at only a tiny fraction of the long-term relative microplate velocity. ! Therefore, for purposes of correcting geodetic velocities it is better ! to approximate these shallow strips as locked, rather than approximate ! them as creeping at the full long-term relative microplate velocity. !--------------------------------------------------------------------------------- ! Version 2.3, 2010.01.15: Used for: ! -modeling of the Ninety East-Sumatra orogen with Tom Robinson. ! -modeling of the Hispaniola orogen with Rafal Jankowski. ! -posted on web site and sent to Zhen Liu, GeoPentech, etc. !--------------------------------------------------------------------------------- ! *Added one more output file resulting from each run: ! v_interseismic_[token].out is now added to represent short-term ! interseismic velocities of nodes. This can be displayed in NeoKineMap, ! for comparison with GPS strain-rates, or with other codes that do not ! compute long-term velocities. It can be differentiated in ! NeoKineMap to display short-term interseismic strain-rates. ! Also, it can be used with Strainrate_exporter.f90/.exe to convert ! short-term interseismic strain-rates to the gridded-values format ! used by David Sandwell et al. for their comparison studies. ! (However, note that these strain-rates will be an undifferentiated ! mixture of elastic and permanent.) !--------------------------------------------------------------------------------- ! Version 2.4 of 10 March 2010: ! -used for GPS/UCERF3 comparison study run by Kaj Johnson, Wayne Thatcher, ! David Sandwell, and Liz Hearn in March-April 2010. ! -used for UCERF3 block-modeling exercises (under the same organizers) in ! Spring-Fall 2011. ! -used for the fault-based UCERF3064~UCERF3077 deformation models submitted to ! UCERF3 in 2011-2012. ! -used for exploratory block-models of the western US (through WUS3_003) ! prepared for the USGS NSHMP/WUS project in 2012.10. !--------------------------------------------------------------------------------- ! MAJOR REVISIONS FOR VERSION 3: ! -Large geodetic matrices (covariance, and normal) are only formed and used ! if required by user input: using_GPS_matrices = using_gp2_file.OR.floating_frame ! This permits MAJOR decreases in bandwidth and MAJOR increases in speed, and ! also frees up memory, allowing 32-bit computers to handle >>1,500 GPS sites. ! -Optional brackets (lower and upper limits) on each fault offset-rate can now ! be prescribed using 2 extra columns in the f*.nki input file. ! These are enforced by automatically increasing the weight on the geologic ! target offset rate (for nonconforming faults only) by successive factors ! of 2 in iterations during the latter 2/3rds of the model run. ! [N.B. Although the 2 new columns are optional in most lines of the file, ! there MUST be entries for these columns in the fixed-FORMAT at the top.] ! -Adjusted iteration strategy, and raised minimum iterations (refinements) to 12. ! During the first 1/3 of iterations, the program runs without any nonlinear ! constraints (which are irreversible). During the second 1/3 of iterations, ! it adds any brackets on fault offset rates. During the final 1/3, it boxes ! any continuum elements where eDot_1H would otherwise be off by 90 degrees, ! and continues to survey fault offset rates for possible bracket tightening. ! -Input file of digitized fault traces (f*.dig) is now scanned for an optional ! additional header line after the F1234 number, sense(s), and title of each ! fault. This new header line may contain the flag "dip_degrees" followed by a ! number (either integer or decimal), and the indicated dip will be used in ! all computations. Faults without this optional header will continue to dip ! 90 degrees for strike-slip, or according to these pre-coded values: ! REAL, PARAMETER :: normal_dip_degrees = 55.0 ! REAL, PARAMETER :: thrust_dip_degrees = 20.0 ! REAL, PARAMETER :: subduction_dip_degrees = 14.0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity. ! -Made dips of "shadow" strike-slip offset-rate components equal to the dips ! of their parent dip-slip offset-rate components, whether that dip came ! from a pre-coded value or from a "dip_degrees" flag. !--------------------------------------------------------------------------------- !--------------------------------------------------------------------------------- ! Version 3.0, 2012.12.06: Used for production runs of NSHM-WUS fault-based models ! prepared for USGS, to be one branch of the deformation- ! model tree in the Western U. S. (WUS) portion of the ! 2014 update of the National Seismic Hazard Map (NSHM). ! (These runs occurred in 2012.12..., although the map was ! not released until 2014.) !--------------------------------------------------------------------------------- !--------------------------------------------------------------------------------- ! *Corrected a logical bug that prevented the matrix of covariances of GPS ! velocity components from being read, even if the filename was given. ! *Changed the default behavior to always output the e[token].nko file, ! even when no faults are included in the model. While this file is ! not needed for plotting long-term permanent strain-rates with ! NeoKineMap, it might be useful for other purposes, such as seismic ! hazard assessment. It can always be deleted, if not wanted. !--------------------------------------------------------------------------------- ! Version 3.1, 2014.05.06: Provided to Michele Carafa for modeling of Europe. !---------------------------------------------------------------------------------