! Restore2.f90 ! ! A code which uses geologic data on fault offsets and strains in ! balanced cross sections, plus latitude anomalies and vertical- ! axis rotations from paleomagnetism, plus paleostress directions, ! to compute paleotectonic and palinspastic maps. It can also ! compute velocities and strain-rates (neotectonics) at the present ! or at an older epoch. ! ! in the Fortran 90 language ! (for version and date, see the 8th WRITE below) ! ! by Peter Bird ! Department of Earth and Space Sciences ! University of California ! Los Angeles, CA 90095-1567 ! ! contacts for Peter Bird: phone (310) 825-1126 ! Internet ! ! (c) Copyright 1997, 1999 by Peter Bird and the ! Regents of the University of California ! !========================================================================= ! ! 1. GENERAL DESCRIPTION: ! ------------------------ ! This program uses geologic and paleomagnetic data to compute ! paleotectonic flow and deformation patterns, and integrates ! these backward over time to create palinspastic restorations. ! Although information on faults is a major part of the input, ! a continuum finite element approximation is used to describe ! the model results; fault activity can be suggested by ! overprinted icons and/or by overprinted actual traces with ! information about estimated slip-rates. ! !========================================================================= ! ! 2. ALGORITHM: !--------------- ! P. Bird (1998), Kinematic history of the Laramide orogeny ! in latitudes 35-49N, western United States, ! Tectonics, vol. 17, #5, p. 780-801. ! !========================================================================== ! ! 3. OVERVIEW OF INPUT DATASETS ! (using suggested filenames for cross-reference in this ! documentation; most actual filenames may be different) ! ! !--------------------------------------------------------- !"PARAMETE[RS].DAT":A short file of strategic parameter values, such as ! the names of other input files, the allowable stain- !"paramete[rs].dat" rate for the a-priori stiffness, size of timesteps, ! number of iterations, et cetera. <1> !"Paramete[rs].dat" Note: This is the only filename which is fixed. ! The 8.3-character form can be used under DOS or MVS. ! Under Unix, almost any upper- or lower-case name ! will be found (except odd ones like "pARAMETERS.dAT"!). ! "F.DAT": Fault names, offsets, and age limits on movement. <2> ! "F.DIG": Digitized fault traces, with codes to match F.DAT. <3> ! "C.DAT": Present and restored lengths of balanced cross- ! sections. <4> ! "P.DAT": Paleomagnetic paleolatitude anomalies (from inclination ! anomalies remaining after structural corrections), and ! vertical-axis rotations (from declination ! anomalies remaining after structural and plate-tectonic ! corrections). <8> ! "S.DAT": Most-compressive horizontal principal stress azimuths ! (from dikes and veins, or from cluster analysis of faults ! with slickensides- or, less reliably, from clusters of ! folds). <9> ! "y.DIG": Digitised traces of any fiducial lines which are to be ! retrodeformed (state lines, coast lines, present ! parallels and meridians, etc.) <10> ! Typically created by utility programs DIGITISE ! for PCs with DOS, with coordinate transformations ! by program PROJECTOR. ! Limit the name 'y' to 4 characters or less, or it ! will be truncated. ! "x.FEG": Finite element grid of connected spherical-triangle ! elements on the Earth's spherical surface. <11> ! Filename must have extension ".feg" or ".FEG". ! Typically created with interactive utility program ! ORBWEAVE for PCs with DOS or Windows (3.x, 95, 98, NT). ! 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 OrbNumber to reduce bandwidth.) ! Limit the name 'x' to 4 characters or less, or it ! will be truncated. ! "x.BCS"; A boundary-conditions file indicating which nodes (at ! least 2!) have specified (usually zero) velocities. <12> ! Notes: Names "x" and "y" are chosen by the user ! and are entered in the "PARAMETE[RS].DAT" file, ! but will be truncated to 4 characters maximum ! when creating the names of output files. ! 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 directory ! for this computation, including any restarts or additional ! iterations. Then, run Restore from this directory, using ! short, simple filenames. Otherwise, the number of files ! may be overwhelming, and their locations will be scattered!) ! If any of the input datasets is null, enter 'none' ! for the filename. (However, you must provide parameters, ! finite element grid, and boundary conditions at least.) ! !============================================================================= ! ! 4. DATA PREPARATION / INPUT FILE FORMATS ! *** General Note: Do not exceed 132 characters in any line of any file! ! --------------------------------------------------------------------------- ! -Build a file "PARAMETE[RS].DAT" containing these lines; comments to right are OK. ! (V----beginning each in column 1): ! "0. geologic age at which this run starts, in Ma ( 0.? ) ! 85. greatest geologic age (end of run), in Ma ( >= 0. ) ! 5. timestep, in Ma ( >0., unless start_time == end_time ) ! 3 refinements (for nonlinearity) within each timestep ! 0 iterations of entire history PRIOR TO this computational run ! 10 iterations of entire history IN THIS computational run ! 100 number of iterations of history planned in ALL runs ! TRUE watch? (adds 1~4 rate files/iteration, for convergence studies) ! TRUE map_set? Do you want 2-4 output files/timestep in last iteration? ! 5.E-17 standard deviation of nominally zero strainrates (mu_), /s ! 1.0E-16 scale strain-rate of "rigid" blocks, /s (for early iterations only) ! 1.6E-12 scale slip-rate sigma of faults, m/s (for early iterations only) ! 1.6E-12 scale rate uncertainty of BCS's, m/s (for early iterations only) ! 3.2E-10 scale N-S drift uncertainty of PM, m/s (for early iterations only) ! 1.4E-16 scale spin uncertainty of PM, radians/s (for early iterations only) ! 1.E-17 small strain-rate increment (xi_), in /s ! 6371. radius of planet, in kilometers ! TRUE switch: Do new active faults give sigma_1h direction data? ! F78DEMO.DAT filename of fault offsets (or, 'none') ! F78DDEMO.DIG filename of digitized fault traces (or, 'none') ! C78.DAT filename of balanced cross-sections (or, 'none') ! P78.DAT filename of paleomagnetic data (or, 'none') ! S78.DAT filename of horizontal principal stress directions (or, 'none') ! STATES.DIG filename of fiducial lines to be restored (or, 'none') ! DEMO78A.FEG filename of finite element grid #1 (required) ! [ DEMO78B.FEG (*optionally, list additional finite element grids) ] ! DEMO78A.BCS filename of boundary conditions for finite element grid #1 (required) ! [ DEMO78B.BCS (*if there is more than 1 .FEG, each must have a .BCS file) ] ! ------------------------------------------------------------------------- ! *Optionally, you can list up to max_fegs = 20 finite element grid files. ! When the first grid becomes exessively distorted, the second will be ! used. When it is excessively distorted, the third will be used, and so ! on. When the last grid is excessively distorted, the program stops. ! You can list the same grid multiple times; this works well in the interior, ! but will not match the changing shape of the continental margin; for that, ! you need hand-edited custom grids based on untangling the topology of ! the previous deformed grid without changing its overall shape. Be sure to ! renumber each grid for minimum bandwidth, with OrbNumber. For each ! .FEG file listed, there must be a corresponding .BCS boundary-conditions ! file. This file must use the new (re-numbered) node identifiers produced ! by OrbNumber. ! ------------------------------------------------------------------------- ! -Build files F.DAT, C.DAT, P.DAT, and S.DAT (all optional) as follows. ! A general feature of all of these files is that they are tables ! with one datum per line (except that a paleomagnetic site has both ! a latitude anomaly and a vertical axis rotation on one 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,F11.2,F11.3,2F11.1)" in the case of F.DAT, and ! -human-readable abbreviated column headings ! (which will be copied into the output files). ! NOTE: Create FORMATs carefully, as they are also used to WRITE files ! that will be needed to restart a calculation. Be sure that there are ! enough significant digits provided in the "F11.2" type formats! ! Also, be sure that all input numbers have an included decimal point ! to guard against decimal point misalignment if the format is changed! ! After these two header lines, the contents are different in ! each case: ! F.DAT = Fault Offsets ! For each fault of regional scale: ! 1. Identifier of fault, 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 is the sense of offset (measurable component): ! R = right-lateral (dextral); L = left-lateral (sinistral); T = ! thrust (of ~25 degree dip); P = low-angle thrust plate or nappe; ! N = normal (high-dip, ~65 degrees); D = detachment (low-dip). ! Thus, a typical identifier might be "F0059T". ! 2. Descriptive text with fault name, location, ... ! 3. Amount of offset, in km. (Always a positive number. ! Letter above gives sense. In the case of thrust (T) or normal ! (N) faulting, give the relative vertical offset (throw) of ! the two sides, which can usually be measured more ! accurately than slip. ! In case of a low-angle thrust plate or nappe (P), give the ! amount of crustal shortening. ! In case of a detachment fault, listric normal fault, or domino- ! style set of rotating normal faults (D), give the net crustal ! extension. ! 4. Standard deviation (sigma_; 68%-confidence ; half of 95%- ! confidence ) of offset, in km ! 5. Maximum age of faulting, in Ma. ! 6. Minimum age of faulting, in Ma. If not known, enter 0. ! C.DAT = Strains Computed from Balanced Cross Sections ! For each balanced section: ! 1. Reference, in short form (e.g., "Jones, 1991"). ! 2. (present) East longitude of West end, in decimal degrees ! (typically negative, e.g., -106.92) ! 3. (present) North latitude of West end, in decimal degrees ! (e.g., 43.81) ! 4. (present) East longitude of East end, in decimal degrees ! (typically negative, e.g., -104.12) ! 5. (present) North latitude of East end, in decimal degrees ! (e.g., 44.03) ! 6. Map code (e.g., "C0001"). ! 7. Present length of section, in km. ! 8. Restored length of section, in km. ! 9. Standard deviation (sigma_; 68%-confidence; half of 95%- ! confidence) of restored length, in km ! 10. Age of restoration epoch, in Ma. ! 11. If available, geologic time by which all deformation was ! known to be over (e.g., from overlap assemblages). In Ma. ! If not available, enter 0; do not just omit! ! P.DAT = Paleomagnetic Latitude Anomalies and Vertical-Axis ! Rotations ! For each virtual geomagnetic pole (preferably from multiple ! samples): ! 1. Reference. If from IAGA Paleomagnetic database ! (McElhinny & Lock, 1995) then use prefix IAGA: "IAGA: Adams ! & Eve, 1901". If from my Notebook database, then use ! prefix NB: "NB: Jones, 1991". ! 2. (present) East longitude of sample(s), in decimal ! degrees (typically negative, e.g., -106.92) ! 3. (present) North latitude of sample(s), in decimal ! degrees (e.g., 43.81) ! 4. Paleolatitude anomaly, in degrees (paleolatitudes to the South ! of present are typical along the west coast; these are ! defined as negative paleolatitude anomalies). ! 5. Standard deviation (sigma_; 68%-confidence; half of 95%- ! confidence) of latitude anomaly, in degrees ! 6. Net vertical-axis rotation since magnetization, in ! degrees (counterclockwise rotations going from past to present ! are considered positive). ! 7. Standard deviation (sigma_; 68%-confidence; half of 95%- ! confidence) of vertical-axis rotation, in degrees ! 8. Maximum age of magnetization, in Ma ! 9. Minimum age of magnetization, in Ma ! 10. Longitude of geomagnetic North Pole (in a normal epoch) ! at the time of magnetization, in the reference frame used ! for velocity boundary conditions, in degrees (E = +). ! 11. Latitude of geomagnetic North Pole (in a normal epoch) ! at the time of magnetization, in the reference frame used ! for velocity boundary conditions, in degrees (N = +). ! S.DAT = Most-Compressive Horizontal Principal Stress Directions ! For each assemblage of stress indicators: ! 1. Reference, in short form (e.g., "Jones, 1991"). ! 2. Location and state abbreviation. ! 3. Map code (e.g., "S0001") ! 4. (present) East longitude, in decimal degrees (typically ! negative, e.g., -106.92) ! 5. (present) North latitude, in decimal degrees (e.g., ! 43.81) ! 6. (present) Azimuth, measured clockwise from North, in ! degrees ! 7. Standard deviation (sigma_; 68%-confidence; half of 95%- ! confidence) of stress azimuth, in degrees ! 8. Maximum age of dikes/veins/joints/faults, in Ma ! 9. Minimum age of dikes/veins/joints/faults, in Ma ! 10. Enter "Window" if the age of a single stress indicator is ! bracketed between (8.) and (9.); otherwise enter "Stage" if ! multiple stress indicators were used to show that the ! stress direction was constant from time (8.) to time (9.). ! -------------------------------------------------------- ! -Build the x.FEG file with ORBWEAVE (using no fault elements), and ! renumber for minimum bandwidth with OrbNumber. ! -------------------------------------------------------- ! -Build the x.BCS file to go with x.FEG and to specify which nodes ! have fixed velocities (at least 2!). This file is read with ! list-directed input, so format is not critical. For each ! node whose velocity you want to fix (and you can list these ! in any order), provide one line of x.BCS with: ! node_number Southward-velocity Eastward_velocity, e.g. ! " 129 -1.375E-09 +2.541E-10" ! Specify these velocities in meters/second. (However, in ! most cases, you will probably want to specify zero velocities: ! " 129 0. 0. " ! Please understand that it is NOT necessary or desirable to ! provide boundary conditions all the way around the grid! ! Usually, one stable interior side is fixed, and the ! rest left free to move as determined by the geologic data. ! This is appropriate because integrating strain is like solving ! a first-order differential equation for displacement, and ! first-order equations should only have a boundary condition on ! one side of the domain! (If you use conditions on both sides, ! you will get a solution which matches them, but the strain in ! unfaulted elements may be unreasonably large!) ! The node numbers you use must be the NEW node numbers ! assigned by renumbering utility OrbNumber; in order to see ! these, load the output x.feg file from OrbNumber back ! into ORBWEAVE, 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. ! -------------------------------------------------------- ! -F.DIG, the file of digitised fault traces, must have exactly ! the following format: ! V---(column 1). [16 sample lines of F.DIG follow this line.] ! F0489N ! -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 SEGMENT *** ! F0532T ! -1.12405E+02,+3.91902E+01 ! [... and so on...] ! Each fault trace must be introduced by a label line written ! by WRITE (nn,"('F',I4,A1)") fltnum, fltsns ! where the INTEGER :: fltnum is used to tie the trace to ! data in file F.DAT, and the CHARACTER*1 :: fltsns is either: ! 'T' for thrust, 'P' for low-angle thrust plates or nappes, ! 'N' for normal, 'D' for detachment or listric normal or domino-set, ! 'R' for dextral, or 'L' for sinstral. ! (Note: N faults have ~65 degree dip; D faults have varying dip or ! ~0 degrees dip at present, but may have dipped more steeply when ! active. T faults have ~25 degree dip; P faults have very low ! dips (except in ramps). ! See the paper listed under 2. Algorithm.) ! 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 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. ! Do not bother to digitise a lot of faults that will not fall ! within your finite element grids; these merely waste memory ! during the run(s), and are omitted from output files because ! they cannot be restored. ! -------------------------------------------------------- ! -y.DIG, the fiducial-lines file, has exactly the same format ! as F.DIG above, except that identifying text label(s) before each ! line segment are optional (if present, they will not be ! transferred to the output file). This file might contain ! present-day state lines, coastlines, meridians of longitude, ! parallels of latitude, etc.. They will be restored to their ! former posititions for use in graphics (which are your ! responsibility). This file is optional. If you wish to ! create one, use my programs DIGITISE and PROJECTOR, ! as described above under file F.DIG. ! Do not bother to digitise a lot of lines that will not fall ! within your finite element grids; these merely waste memory ! during the run(s), and are omitted from output files because ! they cannot be restored. ! !==================================================================== ! ! 5. OVERVIEW OF OUTPUT DATASETS: ! !---------------------------------------------- ! *Text dataset: REPORT.txt 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 REPORT.txt's ! will be overwritten if they are not moved or renamed! ! ! THE FOLLOWING ARE PRODUCED IN EVERY TIMESTEP OF THE LAST ITERATION ! -=IF=- the input parameter "map_set" is .TRUE. ("T"): ! ! *Palinspastic datasets: ! "Fmmnn.DIG": Digitised fault traces, restored. <22> ! (assuming that fault data were used in solution). ! "xmmnn.FEG": Finite-element grid, retrodeformed. <23> ! "ymmnn.DIG": Digitised fiducial lines, restored. <24> ! (assuming that fiducial points were read). ! These formats are unchanged. However, all positions (decimal ! degrees of East longitude and North latitude) are modified ! back to the time of the report. ! The fiducial time "nn" (a geologic age in Ma, rounded to the ! nearest integer) is inserted in each filename created. ! (Note: if the age in Ma is over 99, the tens digit will be a ! letter, allowing values from A0 = 100 to Z9 = 359 in two bytes.) ! The iteration number "mm" refers to the global iteration of the ! entire solution, and helps to keep very similar datasets ! catalogued without confusion. The same convention using values ! up to Z9 is used for iteration #. (The parameter old_iterations in ! the Parameters.dat input file is needed to keep these bytes accurate.) ! Points in "F.DIG", "S.DAT", and "y.DIG" which did not fall into ! the area of the finite element grid "x.FEG" will simply be ! deleted, since they cannot be restored. ! ! *Paleotectonic (or neotectonic) datatset: ! A dataset of velocities of the nodes of the finite-element ! grid "xmmnn.FEG" is also produced: xmmnn.VEL <25>. This permits ! plotting diagrams of the paleotectonics using RetroMap. ! Note that the velocities are NOT the averages over the timestep, ! but are the values at the computational-end (geologic-beginning) ! so that they tie to the node positions in "xnnn.FEG". ! Also, velocities point forward in time (toward the present), not ! back in time; they show what was happening at a previous epoch. ! ! *Strain-history datasets: ! See explanation of "mmnn" in filenames in "palinspastic datasets" above. ! These files have some lines identical to the corresponding input datasets: ! Pmmnn.DAT <28> is expanded from P.DAT <8> ! Smmnn.DAT <29> is expanded from S.DAT <9> ! but have other lines added to show what has happened in the run(s). ! The extra lines are distinguished by a special character (+, *, &, $) in column 1. ! + lines: These record the position of the datum nn Ma ago, ! in the reference frame which is used to define boundary conditions. ! The first number is East longitude in degrees, the second is ! North latitude in degrees. ! Fmmnn.DAT lacks these lines because the information is in Fmmnn.DIG. ! Cmmnn.DAT has two such lines per datum for the two ends of the section. ! Pmmnn.DAT has one such line per datum. ! Smmnn.DAT has one such line per datum. ! This information is needed to restart a history in the middle. ! * lines: These record the rates of offset/displacement/rotation ! during the most recent iteration of the history (left), ! followed in each case by a new goal (right) for future computations. ! The first number is the age in Ma at the young end of some timestep. ! The second number is the age in Ma the old end of the same timestep. ! The third number is the model rate during that step. ! The fourth number is the goal for future computations. ! The units are derived from the original data line above, ! with 1 m.y. used as the time unit to define rates: ! Fmmnn.DAT has fault offset rates in km/m.y.. Always positive. ! Cmmnn.DAT has cross-section length rates in km/m.y.. Extension ! from past to present has a positive sign. ! Pmmnn.DAT has latitude rates in degrees/m.y.. If the point moved ! toward the S paleopole in going from past to present, ! this rate is entered with a positive sign. ! Smmnn.DAT does not have * lines. ! This information is very important, as it is the only way that ! a computation restarted for more iterations can remember what ! it learned in previous iterations! ! & lines: These are like * lines except they record additional progress: ! Fmmnn.DAT does not have & lines. ! Cmmnn.DAT does not have & lines. ! Pmmnn.DAT has rotation rates in degrees/m.y.. ! Counterclockwise rotation from past to present is +. ! Smmnn.DAT does not have & lines. ! $ lines: These record a paleo-azimuth indicator nn Ma ago. ! The reference frame is that used to define boundary conditions. ! Fmmnn.DAT does not have $ lines. ! Cmmnn.DAT does not have $ lines. ! Pmmnn.DAT does not have $ lines. ! Smmnn.DAT has paleo-azimuths of sigma_1h. As in S.DAT, ! azimuth is in degrees clockwise from North. ! ! THE FOLLOWING ARE PRODUCED AT THE END OF THE LAST ITERATION !(or, at the end of EVERY iteration if parameter "watch" is .TRUE. ('T'): ! ! *Strain-history datasets: ! See explanation of "mmnn" in filenames in "palinspastic datasets" above. ! These files have some lines identical to the corresponding input datasets: ! Fmmnn.DAT <26> is expanded from F.DAT <2> ! Cmmnn.DAT <27> is expanded from C.DAT <4> ! Pmmnn.DAT <28> is expanded from P.DAT <8> ! Smmnn.DAT <29> is expanded from S.DAT <9> ! but have other lines added to show what has happened in the run(s). ! The extra lines are distinguished by a special character (+, *, &, $) in column 1. ! + lines: These record the position of the datum nn Ma ago, ! in the reference frame which is used to define boundary conditions. ! The first number is East longitude in degrees, the second is ! North latitude in degrees. ! Fmmnn.DAT lacks these lines because the information is in Fmmnn.DIG. ! Cmmnn.DAT has two such lines per datum for the two ends of the section. ! Pmmnn.DAT has one such line per datum. ! Smmnn.DAT has one such line per datum. ! This information is needed to restart a history in the middle. ! * lines: These record the rates of offset/displacement/rotation ! during the most recent iteration of the history (left), ! followed in each case by a new goal (right) for future computations. ! The first number is the age in Ma at the young end of some timestep. ! The second number is the age in Ma the old end of the same timestep. ! The third number is the model rate during that step. ! The fourth number is the goal for future computations. ! The units are derived from the original data line above, ! with 1 m.y. used as the time unit to define rates: ! Fmmnn.DAT has fault offset rates in km/m.y.. Always positive. ! Cmmnn.DAT has cross-section length rates in km/m.y.. Extension ! from past to present has a positive sign. ! Pmmnn.DAT has latitude rates in degrees/m.y.. If the point moved ! toward the S paleopole in going from past to present, ! this rate is entered with a positive sign. ! Smmnn.DAT does not have * lines. ! This information is very important, as it is the only way that ! a computation restarted for more iterations can remember what ! it learned in previous iterations! ! & lines: These are like * lines except they record additional progress: ! Fmmnn.DAT does not have & lines. ! Cmmnn.DAT does not have & lines. ! Pmmnn.DAT has rotation rates in degrees/m.y.. ! Counterclockwise rotation from past to present is +. ! Smmnn.DAT does not have & lines. ! $ lines: These record a paleo-azimuth indicator nn Ma ago. ! The reference frame is that used to define boundary conditions. ! Fmmnn.DAT does not have $ lines. ! Cmmnn.DAT does not have $ lines. ! Pmmnn.DAT does not have $ lines. ! Smmnn.DAT has paleo-azimuths of sigma_1h. As in S.DAT, ! azimuth is in degrees clockwise from North. ! ! THE FOLLOWING ARE PRODUCED IN THE EVENT THAT THE PROGRAM TERMINATES ! EARLY (because the grid becomes too strained, and no other grids are ! available to read in): ! ! *Palinspastic datasets (described above): ! "Fmmnn.DIG": Digitised fault traces, restored. <22> ! (assuming that fault data were used in solution). ! "xmmnn.FEG": Finite-element grid, retrodeformed. <23> ! "ymmnn.DIG": Digitised fiducial lines, restored. <24> ! (assuming that fiducial points were read, AND ! that input parameter "map_set" is .TRUE.) ! ! *Strain-history datasets (described above): ! Fmmnn.DAT <26> is expanded from F.DAT (if any) <2> ! Cmmnn.DAT <27> is expanded from C.DAT (if any) <4> ! Pmmnn.DAT <28> is expanded from P.DAT (if any) <8> ! Smmnn.DAT <29> is expanded from S.DAT (if any) <9> ! Note that in the case of early termination, the rates for ! completed timesteps are actual model predictions, but the rates ! for the remaining timesteps are goals. When Restore is restarted ! with a new grid, these goals will be re-used. !========================================================================= ! ! Fortran 90 free-form source code begins: ! ---------------------------------------------------------------- ! NOTE: Some kludgy coding intended to fix the fatal memory-leaks ! caused by the bugs in MicroSoft Fortran PowerStation 4.0 Pro: ! 1. No array segments are used as actual arguments. ! Instead, temporary vectors are used for input. ! (To locate all these kludges, search on the letters "tv".) ! 2. Array-valued functions were converted to SUBROUTINES (ouch!): ! (Cross, Del_Gjxy_del_thetaphi, Gjxy, Interpolate, Local_Phi, ! Local_Theta, Moved_by_vel, Moved_by_vw, Step_aside, Unitise, ! XYZ_from_lonlat). ! 3. TYPE(is123)-valued function "Inside" was replaced with ! SUBROUTINE Internal, which returns an integer argument ! and 3 real arguments (not a TYPE(is123) argument)! ! ----------------------------------------------------------------- !******************************************************************* PROGRAM Restore2 USE NUMERICAL_LIBRARIES ! DIGITAL version of International Mathematics Subroutine Library: ! LSLPB: solves systems of linear equations with real ! coefficients arranged in a special band storage mode. ! LSLSF: solves systems of linear equations with symmetric ! matrix of real coefficients in uncompressed form. ! EVCRG: eigenvectors and eigenvalues of a real matrix. IMPLICIT NONE ! All variable names must be declared. !-------------------------------------------------------------------- TYPE :: is123 ! element & internal coordinates of any point on surface INTEGER :: element REAL, DIMENSION(3) :: s END TYPE is123 TYPE :: crack ! Intersection of a fault offset rate datum with ! a fault segment. Some segments have > 1 crack, due to ! > 1 contiguous datum time windows in this timestep; ! many fault segments have no crack (inactive). INTEGER :: datum ! offset datum index INTEGER :: segment ! segment index REAL :: s_ ! sliprate, converted, rules according to 'sense' REAL :: sigma_ ! standard deviation of sliprate, ditto REAL, DIMENSION(3) :: H ! see (11)-(17) of Bird (1998?) CHARACTER(1) :: sense ! fault type: R, L, T, N, D END TYPE crack TYPE :: needle ! everything one needs to know about a stress datum REAL, DIMENSION(3) :: location ! Cartesian unit vector REAL :: azimuth ! clockwise from N, in radians REAL :: sigma ! uncertainty, radians REAL :: relevance ! 0.0 to 1. END TYPE needle !-------------------------------------------------------------------- !CONSTANTS: Fixed conversion factors, in alphabetical order: INTEGER, PARAMETER :: bytes_per_int = 4 ! descriptive, not prescriptive! INTEGER, PARAMETER :: bytes_per_real = 4 ! descriptive, not prescriptive! !Check documentation of your compiler to see if the above are correct. !If they are wrong, program runs the same, but array-size reporting will be off. INTEGER, PARAMETER :: bytes_per_is = bytes_per_int + 3 * bytes_per_real INTEGER, PARAMETER :: bytes_per_crack = 2 * bytes_per_int + 5 * bytes_per_real + 1 INTEGER, PARAMETER :: bytes_per_Mb = 1024 * 1024 REAL, PARAMETER :: cot_thrust_dip = 2.14451 ! 1./TAN(25.) REAL, PARAMETER :: cot_normal_dip = 0.46631 ! 1./TAN(65.) REAL, PARAMETER :: deg_per_rad = 180./3.141592654 REAL, PARAMETER :: m_per_km = 1000. INTEGER, PARAMETER :: max_fegs = 20 !(could be increased) REAL, PARAMETER :: s_per_Ma = 1000000.*365.25*24.*60.*60. REAL, PARAMETER :: Pi = 3.14159265 !-------------------------------------------------------------------- !VARIABLES: All global variables except arrays, in alphabetical order: REAL :: allowance ! grace period (s) for association of neotec solution with stage stress data LOGICAL :: any_action ! is x_active(j,i) TRUE for any j? INTEGER :: bcs_count ! # of boundary-condition nodes REAL :: ccw ! counterclockwise rotation, radians CHARACTER(80) :: c_dat ! filename, X-section data INTEGER :: c_dat_count ! number of X-section data CHARACTER(134) :: c_dat_format ! to read X-section data CHARACTER(134) :: c_dat_titles ! to write X-section data REAL :: c_scale ! scale stretch-rate uncertainty to start iteration, m/s CHARACTER(1) :: c, c1 !(temporary) CHARACTER(4) :: c4 !(temporary) CHARACTER(5) :: c5 !(temporary) CHARACTER(6) :: c6 !(temporary) CHARACTER(30) :: c30, c30a !(temporary) CHARACTER(47) :: c47 !(temporary) CHARACTER(50) :: c50 !(temporary) CHARACTER(80) :: c80 !(temporary) CHARACTER(134) :: c134 !(temporary) LOGICAL :: check_if ! any node lies "on" a fault trace INTEGER :: c_in_time_and_space ! number of data within time, domain windows CHARACTER(10) :: clock_time ! wall clock time INTEGER :: complete_timesteps = 0 INTEGER :: crack_count ! total number of cracks in this timestep INTEGER :: current_feg ! ordinal number of loaded .feg file CHARACTER(8) :: date ! date program is run INTEGER :: Delta_node ! greatest difference between connected node #s INTEGER :: Delta_node_feg ! lower limit of Delta_node, based on x_feg INTEGER :: Delta_node_last = 0 ! memory of previous Delta_node REAL :: Deltat_ ! timestep, in s REAL :: end_time ! greatest geologic age, in s LOGICAL :: eof ! SUBR sends signal to calling prog REAL :: exponent ! = n_ / last_iter LOGICAL :: faults_give_sigma_1h ! are faults stress indicators? CHARACTER(80) :: f_dig ! filename, digitized traces INTEGER :: f_dig_count ! # of points in f_dig CHARACTER(80) :: f_dat ! filename, offset data INTEGER :: f_dat_count ! # of data in f_dat CHARACTER(134) :: f_dat_format ! to read offset data CHARACTER(134) :: f_dat_titles ! to print offset data INTEGER :: f_highest ! max fault index INTEGER :: f_in_time_and_space ! number of data within time, domain windows REAL :: f_scale ! scale slip-rate uncertainty to start iteration, m/s REAL :: factor ! converts tabulated offset rate to internal REAL :: floor = 10. * EPSILON(floor) LOGICAL :: folding ! finite element grid has folded REAL :: gamma_ ! azimuth, clockwise from N, radians LOGICAL :: get_feg ! memo that new .feg file is needed LOGICAL :: got_index ! during reading of f_dig LOGICAL :: got_point ! during reading of f_dig REAL :: half_R2 ! R**2/2. INTEGER :: i,i1,i2,i3,i4,i5,i6 ! (temporary) INTEGER :: in_count ! # of fault offsets in L0, L1, L2 LOGICAL :: in_trace ! during reading of f_dig, y_dig INTEGER :: iteration ! index for repetitions of whole history INTEGER :: j, j1, j2, j3 ! (temporary) INTEGER :: k ! (temporary) INTEGER :: l_ ! finite element index INTEGER :: last_iter ! planned final iteration of a set INTEGER :: lda ! leading DIMENSION of ABCDEF, for LSLPB INTEGER :: line ! line number of any input file LOGICAL :: loading ! need to (re)read maps and datasets & do s1s2s3 INTEGER :: m ! (temporary) INTEGER :: max_iter ! # of iterations of history INTEGER :: memory ! bytes of memory in all arrays REAL :: misfits ! a prediction error, in sigmas CHARACTER(4) :: mmnn ! iteration (I2) and timestep (I2) LOGICAL :: map_set ! do you want complete graphical detail for every ! timestep of the last iteration? REAL :: mu_ ! Bird (1998?) REAL :: mu_scale ! scale strain-rate uncertainty to start iteration INTEGER :: n_ ! timestep number, absolute scale, present - >past INTEGER :: n_1 ! 1st timestep number, absolute scale INTEGER :: n_refine ! # refinements of each v solution INTEGER :: n_t_per_it ! timesteps/iteration; <= num_timesteps INTEGER :: ncoda ! # of codiagonals in ABCD (part of ABCDEF) INTEGER :: ndof ! # of degrees of freedom, =2*num_nod LOGICAL :: neotec ! do present velocities only INTEGER :: num_bad ! # of nodes lying on fault traces INTEGER :: num_ele ! # of finite elements INTEGER :: num_fegs ! # of .feg and .bcs files in Paramete[rs].dat INTEGER :: num_nod ! # of nodes in x_feg INTEGER :: num_timesteps ! # of timesteps in one iteration of history REAL :: offs ! fault offset, scalar form REAL :: overlap ! time common to two time windows, s REAL :: overlap_threshold ! minimum overlap to consider a "stage" stress datum relevant CHARACTER(80) :: p_dat ! filename, paleomagnetic data INTEGER :: p_dat_count ! number of paleomagnetic sites CHARACTER(134) :: p_dat_format ! to read paleomagnetic data CHARACTER(134) :: p_dat_titles ! to write paleomagnetic data REAL :: p_drift_scale ! scale N-S drift uncertainty to start iteration, m/s INTEGER :: p_in_time_and_space ! number of data within time, domain windows REAL :: p_spin_scale ! scale rotation uncertainty to start iteration, radians/s REAL :: P_weight ! weighting of paleomagnetic data LOGICAL :: paleotec ! do palinspastic reconstructions INTEGER :: past_iterations ! iterations of history in previous runs REAL :: R ! Bird (1998?) REAL :: r1, r2, r3, r4 !(temporary) INTEGER :: read_status ! did READ work? INTEGER :: s ! indicates stress datum CHARACTER(80) :: s_dat ! filename, stress directions INTEGER :: s_dat_count ! number of paleostress data from s_dat CHARACTER(134) :: s_dat_format ! to read paleostress data CHARACTER(134) :: s_dat_titles ! to write paleostress data INTEGER :: seg_count ! number of fault segments (element ^ trace) REAL :: south ! Sward motion of paleomag site, in m REAL :: start_time ! younest geologic age, in s (usually 0.) LOGICAL :: stress_ever ! are there stress data for any timestep? LOGICAL :: stress_now ! any stress constraints this timestep? REAL :: stretch ! actual extension of 1 cross-section REAL :: sum !(temporary) REAL :: t, t0, t1, t2, t3,t4 !(temporary) REAL :: time0, time1 ! age at young, old end of timestep, s INTEGER :: total_iterations ! = iteration + past_iterations (q.v.) REAL, DIMENSION(3):: tv,tvi,tvo ! temporary 3-vectors used in kludgy ! coding which avoids passing array ! sub-sections as actual arguments, ! because under MicroSoft Fortran PowerStation ! 4.0 Pro this causes a memory leak! LOGICAL :: watch ! outputs all rates files each iteration REAL :: W ! Bird (1998?) CHARACTER(80) :: y_dig ! filename of digitised fiducial lines INTEGER :: y_dig_count ! number of digitised points in y_dig CHARACTER(80) :: x_feg, x_vel ! current .feg and .vel names (before mmnn) REAL :: x1, x2, x3, x4 !(temporary) REAL :: xi_ ! Bird (1998?) CHARACTER(5) :: zone ! time zone !-------------------------------------------------------------------- !VARS !ARRAYS, in alphabetical order: REAL, DIMENSION(:),ALLOCATABLE :: a_ ! area of plane triangle element beneath surface, m**2 !(1:num_ele = element index l_) REAL, DIMENSION(:),ALLOCATABLE :: adjustments ! fraction of histories adjusted by equation (27) of Bird(1998?), % !(1:max_iter = iteration index) REAL, DIMENSION(:,:), ALLOCATABLE :: ABCDEF ! coefficient matrix of linear system used to solve for velocity ! components; see (4) and (5) of Bird (1998?). ! Stored in codiagonal band symmetric storage mode of Microsoft IMSL, ! in which an extra column is added on right to contain right-hand-side ! vector of the linear system. Also, extra rows are provided at top. ! Logically-diagonal entries are found in column 1, rows ncoda+1:ncoda+ndof=lda. ! The right-hand-side vector is in column ncoda+2, rows ncoda+1:ncoda+ndof=lda. !(1:ncoda+ndof=lda; 1:ncoda+2). INTEGER, DIMENSION(:),ALLOCATABLE :: boundary_node ! index number of node with specified velocities !(1:bcs_count = boundary-condition index) LOGICAL(1), DIMENSION(:),ALLOCATABLE :: boxed ! does this element require boxing for correct sigma_1h? !(1:num_ele = element index) REAL, DIMENSION(:,:),ALLOCATABLE :: center ! Cartesian unit vectors at center of finite elements !(1:3 = x,y,z; 1:num_ele = element index l_) LOGICAL(1), DIMENSION(:,:),ALLOCATABLE :: c_active ! .TRUE. indicates that the X-section datum is applicable to the timestep !(1:num_timesteps = timestep index; 1:c_dat_count = X-section index) CHARACTER(5), DIMENSION(:),ALLOCATABLE :: c_code ! short name of X-section used on master map ! (1:c_dat_count = X-section index) TYPE(is123), DIMENSION(:,:),ALLOCATABLE :: c_end_is ! locations of both ends of cross-section in internal coordinates !(1:2 = west,east end; 1:c_dat_count = X-section index) REAL, DIMENSION(:,:,:),ALLOCATABLE :: c_end_now ! current location of ends of each cross-section (integrated); ! both ends are Cartesian unit vectors from Earth center: ! (1:3 = x,y,z; 1:2 = west,east end; 1:c_dat_count = X-section index) REAL, DIMENSION(:,:,:),ALLOCATABLE :: c_end_0 ! present location of ends of each cross-section ! both ends are Cartesian unit vectors from Earth center; ! (1:3 = x,y,z; 1:2 = west,east end; 1:c_dat_count = X-section index) REAL, DIMENSION(:,:),ALLOCATABLE :: c_err ! 3 norms of cross-section error (each normalized by sigma before combining): !(0:2 = L0,L1,L2 norm; 1:max_iter = iteration index). REAL, DIMENSION(:,:),ALLOCATABLE :: c_goal ! target rates of extension of a X-section in each timestep, in m/s ! (1:num_timesteps = timestep index; 1:c_dat_count = X-section index) REAL, DIMENSION(:,:),ALLOCATABLE :: c_length ! present and past length of X-section, in m ! (1:2 = present, past; 1:c_dat_count = X-section index) CHARACTER(47), DIMENSION(:), ALLOCATABLE :: c_ref ! bibliographic reference for each X-section datum ! (1:c_dat_count = X-section index) REAL, DIMENSION(:,:),ALLOCATABLE :: c_rate ! rates of extension of a X-section in each timestep, in m/s ! (1:num_timesteps = timestep index; 1:c_dat_count = X-section index) REAL, DIMENSION(:), ALLOCATABLE :: c_rate_sigma_ ! uncertainty in rate of extension of cross-section, m/s !(1:c_dat_count = X-section index) REAL, DIMENSION(:), ALLOCATABLE :: c_sigma_ ! standard deviation of restored length of X-section (and of stretch) ! (1:c_dat_count = X-section index) REAL, DIMENSION(:), ALLOCATABLE :: c_stretch ! extension of X-section from past to present, in m ! (1:c_dat_count = X-section index) REAL, DIMENSION(:), ALLOCATABLE :: c_t_max ! restored age of X-section, in s ! (1:c_dat_count = X-section index) REAL, DIMENSION(:), ALLOCATABLE :: c_t_min ! age of overlap assemblage on restored X-section, in s ! (1:c_dat_count = X-section index) REAL, DIMENSION(:,:),ALLOCATABLE :: condition ! boundary velocities, in m/s !(1:2 = theta(South),phi(East) components; ! 1:bcs_count = boundary condition index) INTEGER, DIMENSION(:,:),ALLOCATABLE :: crack_index ! count and pointer to cracks active in each element; !(1:2 = number, position of 1st one in local_crack; ! 1:num_ele = element index) INTEGER, DIMENSION(8) :: datetimenumber ! for output from DATE_AND_TIME REAL, DIMENSION(:,:),ALLOCATABLE :: dot ! fiducial points read from y_dig, ! as Cartesian unit vectors !(1:3 = x,y,z; 1:y_dig_count = fiducial point index) LOGICAL(1), DIMENSION(:), ALLOCATABLE :: dot_last ! If T, given point is the last in a line segment. !(1:y_dig_count = fiducial point index) TYPE(is123), DIMENSION(:), ALLOCATABLE :: dot_is ! locations of digitised fiducial points in internal coordinates !(1:y_dig_count = fiducial point index) REAL, DIMENSION(:,:), ALLOCATABLE :: duplicate ! see ABCDEF above; this is an extra copy. REAL, DIMENSION(:), ALLOCATABLE :: ele_azim ! azimuth of most-compressive horizontal principal stress, at ! each element center, in radians clockwise from North. !(1:num_ele = element index l_) REAL, DIMENSION(:), ALLOCATABLE :: ele_q ! relevance of the stress datum (or interpolated) stress ! in each element for a particular timestep, 0.0 to 1.0 !(1:num_ele = element index l_) REAL, DIMENSION(:), ALLOCATABLE :: ele_sigma ! standard deviation of azimuth of most-compressive horizontal ! principal stress, at each element center, in radians. !(1:num_ele = element index l_) REAL, DIMENSION(:,:),ALLOCATABLE :: ele_strainrate ! continuum strain-rate tensor (not including fault strain) ! at center of each element, saved to be output for plotting ! by RetroMap; no other use by Restore. !(1:3 = theta-theta or N-S, theta-phi or SE, phi-phi or E-W; ! 1:num_ele = element index l_) LOGICAL(1), DIMENSION(:), ALLOCATABLE :: ele_stressed ! is there a useful interpolated value of ele_azim, with associated ! del_az_for_90pc <= 0.7854 radians (45 degrees), for this element? !(1:num_ele = element index l_) LOGICAL(1), DIMENSION(:,:),ALLOCATABLE :: f_active ! .TRUE. indicates that the offset datum is applicable to the timestep !(1:num_timesteps = timestep index; 1:f_dat_count = offset index) CHARACTER(1), DIMENSION(:), ALLOCATABLE :: f_dat_code ! N = Normal; P = Promoted (Quaternary rate used for whole first ! timestep); D = Demoted (rate for 1st timestep ignored since another ! rate for same fault trace was Promoted). !(1:f_dat_count = offset datum index) REAL, DIMENSION(:,:),ALLOCATABLE :: f_divide ! temporary storage for accumulating (over all active segments) ! the numerator and denominator that will be used to compute the ! mean slip rate corresponding to each fault offset datum !(1:2 = numerator/denominator; 1:f_dat_count = offset index) REAL, DIMENSION(:,:),ALLOCATABLE :: f_err ! 3 norms of fault offset error (each normalized by sigma before combining): !(0:2 = L0,L1,L2 norm; 1:max_iter = iteration index). REAL, DIMENSION(:,:),ALLOCATABLE :: f_goal ! target rates of fault offset in each timestep, in m/s; ! interpretation depends on sense. ! (1:num_timesteps = timestep index; 1:f_dat_count = offset datum index) LOGICAL(1), DIMENSION(:), ALLOCATABLE :: f_new ! was this the first stage of movement on the fault? ! (1:f_dat_count = offset datum index) REAL, DIMENSION(:,:),ALLOCATABLE :: f_rate ! rates of fault offset in each timestep, in m/s; ! interpretation depends on sense. ! (1:num_timesteps = timestep index; 1:f_dat_count = offset datum index) REAL, DIMENSION(:), ALLOCATABLE :: f_rate_sigma_ ! uncertainty in fault offset rate, m/s !(1:f_dat_count = offset datum index) REAL, DIMENSION(:), ALLOCATABLE :: f_t_max ! maximum age of fault movement, in s ! (1:f_dat_count = offset datum index) REAL, DIMENSION(:), ALLOCATABLE :: f_t_min ! minimum age of fault movement, in s ! (1:f_dat_count = offset datum index) LOGICAL(1), DIMENSION(:), ALLOCATABLE :: f_2_in ! does this fault have at least 2 points inside the current grid? ! Notes: 1. Once f_2_in(i) is .FALSE., it cannot become .TRUE. ! in any later timestep, even though a new grid is read. ! However, a .TRUE. value can always become .FALSE.. ! Thus, only faults continuously deformed are ever output. ! 2. Faults for which f_2_in is .FALSE. are not output in either ! fmmnn.dat or fmmnn.dig. !(1:f_highest = trace index) CHARACTER(50), DIMENSION(:), ALLOCATABLE :: fault_name ! (1:f_dat_count = offset datum index) CHARACTER(80), DIMENSION(:), ALLOCATABLE :: gridname_bcs ! filenames of finite element grid (.feg) files CHARACTER(80), DIMENSION(:), ALLOCATABLE :: gridname_feg ! filenames of boundary condition files REAL, DIMENSION(21,0:29) :: ln_rel_prob ! matrix of natural logarithms of relative probabilities of ! difference angles between two sigma_1h directions at spatially ! separated sites. Values are roughly +0.8 to -0.8. ! First (row) subscript identifies the distance annulus, ! out to 22 degrees arc (but annuli are not of uniform width). ! Second (column) subscript identifies the size of the angular ! discrepancy, in 3-degree steps: column 0 is for differences of ! 0-3 degrees, and column 29 is for differences of 87-90 deg. TYPE(crack), DIMENSION(:), ALLOCATABLE :: local_crack ! a compilation (using structure type crack) of all needed ! information to describe active fault segment in this timestep, ! sorted by element number. See crack_index for location. !(1:crack_count = crack index). REAL, DIMENSION(:), ALLOCATABLE :: mu_nod ! uncertainty (sigma_) of the nominally-zero strain-rate of stable ! areas; if 0. is read, mu_ is substituted for that node ! (1:num_nod = node index) INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor ! list of neighboring finite elements !(1:3 = side crossed; 1:num_ele = element index l_) INTEGER, DIMENSION(:,:),ALLOCATABLE :: node ! list of nodes defining each element, in counterclockwise order ! as seen from outside the planet. !(1:3 = 3 corners; 1:num_ele = element index l_) REAL, DIMENSION(:), ALLOCATABLE :: offset ! amount of fault offset, in m (always positive; see sense) ! (1:f_dat_count = offset datum index) REAL, DIMENSION(:), ALLOCATABLE :: offset_sigma_ ! sigma_ of fault offset, in m (always positive) ! (1:f_dat_count = offset datum index) LOGICAL(1), DIMENSION(:,:),ALLOCATABLE :: p_active ! .TRUE. indicates that the paleomagnetic datum is applicable to the timestep !(1:num_timesteps = timestep index; 1:p_dat_count = offset index) REAL, DIMENSION(:), ALLOCATABLE :: p_ccw ! counterclockwise rotation of paleomagnetic site, from past ! to present, in radians ! (1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:,:),ALLOCATABLE :: p_ccw_err ! 3 norms of rotation error (each normalized by sigma before combining): !(0:2 = L0,L1,L2 norm; 1:max_iter = iteration index). REAL, DIMENSION(:,:),ALLOCATABLE :: p_ccw_goal ! target counterclockwise rotation rate of a paleomagnetic site, ! going from past to present, in radians/s !(1:num_timesteps = timestep index; 1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:,:),ALLOCATABLE :: p_ccw_rate ! model counterclockwise rotation rate of a paleomagnetic site, ! going from past to present, in radians/s !(1:num_timesteps = timestep index; 1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_ccw_rate_sigma_ ! uncertainty in counterclockwise rotation rate, radian/s !(1:p_dat_count = paleomagentic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_ccw_sigma_ ! standard deviation of counterclockwise rotation of paleomagnetic ! site, from past to present, in radians ! (1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:,:),ALLOCATABLE :: p_pole ! Paleo-North-pole at time of magnetization of paleomagnetic site, ! in reference frame used for velocity boundary conditions, ! expressed as Cartesian unit vector !(1:3 = x,y,z; 1:p_dat_count = paleomagnetic site index) CHARACTER(50), DIMENSION(:), ALLOCATABLE :: p_ref ! bibliographic reference for each paleomagnetic site ! (1:p_dat_count = paleomagnetic site index) TYPE(is123), DIMENSION(:), ALLOCATABLE :: p_site_is ! locations of paleomagnetic sites in internal coordinates !(1:p_dat_count = paleomagentic site index) REAL, DIMENSION(:,:), ALLOCATABLE :: p_site_now ! current location of paleomagnetic site (integrated); ! Cartesian unit vector from Earth center: ! (1:3 = x,y,z; 1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:,:), ALLOCATABLE :: p_site_0 ! present location of paleomagnetic site; ! Cartesian unit vector from Earth center: ! (1:3 = x,y,z; 1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_south ! distance that paleomagnetic site has drifted South, in m ! (1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:,:),ALLOCATABLE :: p_south_err ! 3 norms of paleolatitude error (each normalized by sigma before combining): !(0:2 = L0,L1,L2 norm; 1:max_iter = iteration index). REAL, DIMENSION(:,:),ALLOCATABLE :: p_south_goal ! target velocities toward paleo-South of a paleomagnetic site, in m/s !(1:num_timesteps = timestep index; 1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:,:),ALLOCATABLE :: p_south_rate ! model velocities toward paleo-South of a paleomagnetic site, in m/s !(1:num_timesteps = timestep index; 1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_south_sigma_ ! standard deviation of distance that paleomagnetic site has ! drifted South, in m ! (1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_south_rate_sigma_ ! uncertainty in Southward drift rate, m/s !(1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_t_max ! mean age of magnetization, in seconds !(1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_t_min ! this is the age, in s, at which paleomagnetic sites were sampled; ! so all values are 0.; provided as a necessary actual parameter. !(1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_t1 ! maximum age of magnetization (averaged with p_t2 to give ! p_t_max), in seconds !(1:p_dat_count = paleomagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: p_t2 ! minimum age of magnetization (averaged with p_t1 to give ! p_t_max), in seconds !(1:p_dat_count = paleomagnetic site index CHARACTER(14), DIMENSION(18) :: param_name ! different possible names of the parameters.dat file. REAL, DIMENSION(3) :: pole ! temporary Cartesian unit vector REAL, DIMENSION(:,:,:), ALLOCATABLE :: rate_err ! 3 norms of rate error (each normalized by sigma before combining): !(0:2 = L0,L1,L2 norm; 0:num_timesteps = mean for iteration, each step; ! 1:max_iter = iteration index). REAL, DIMENSION(:,:),ALLOCATABLE :: s_activity ! relevance (q) of a stress datum in the timestep: 0.0 to 1.0 !(1:num_timesteps = timestep index; 1:s_dat_count = offset index) REAL, DIMENSION(:), ALLOCATABLE :: s_azim_now ! (integrated) azimuth of most compressive horizontal principal stress, ! in radians clockwise from North (in the reference frame ! used to define velocity boundary conditions) !(1:s_dat_count = paleostress site index) REAL, DIMENSION(:), ALLOCATABLE :: s_azim_0 ! present azimuth of most compressive horizontal principal stress, ! in radians clockwise from North !(1:s_dat_count = paleostress site index) CHARACTER(5), DIMENSION(:), ALLOCATABLE :: s_code ! master-map location memo for each paleostress datum ! (1:s_dat_count = paleostress index) CHARACTER(30), DIMENSION(:), ALLOCATABLE :: s_loc ! geographic location memo for each paleostress datum ! (1:s_dat_count = paleostress index) CHARACTER(30), DIMENSION(:), ALLOCATABLE :: s_ref ! bibliographic reference for each paleostress datum ! (1:s_dat_count = paleostress index) REAL, DIMENSION(:), ALLOCATABLE :: s_sigma_ ! standard deviation of azimuth of most compressive ! horizontal principal stress, in radians !(1:s_dat_count = paleostress site index) TYPE(is123), DIMENSION(:,:), ALLOCATABLE :: s_site_is ! locations of paleostress sites in internal coordinates !(1:2 = site,neighbor@azimuth_cw_from_N; ! 1:s_dat_count = paleostress site index) REAL, DIMENSION(:,:,:), ALLOCATABLE :: s_site_now ! current location of paleostress site (integrated); ! Cartesian unit vector from Earth center: ! (1:3 = x,y,z; 1:2 = site,neighbor@azimuth_cw_from_N; ! 1:s_dat_count = paleostress site index) REAL, DIMENSION(:,:), ALLOCATABLE :: s_site_0 ! present-day location of paleostress site; ! Cartesian unit vector from Earth center: ! (1:3 = x,y,z; 1:s_dat_count = paleostress site index) LOGICAL(1), DIMENSION(:), ALLOCATABLE :: s_stage ! .TRUE. indicates that paleostress is valid anytime ! from s_t_max to s_t_min (not just SOME time) !(1:s_dat_count = paleostress site index) REAL, DIMENSION(:), ALLOCATABLE :: s_t_max ! maximum age of paleostress, in s ! (1:s_dat_count = paleostress site index) REAL, DIMENSION(:), ALLOCATABLE :: s_t_min ! minimum age of paleostress, in s ! (1:s_dat_count = paleostress site index) INTEGER, DIMENSION(:,:),ALLOCATABLE :: seg_def ! defines a fault segment by its fault trace # and element # !(1:2 = trace, element; 1:seg_count = segment index) REAL, DIMENSION(:,:,:),ALLOCATABLE :: seg_end !Cartesian unit vectors at each end of fault segment !(a fault segment is the intersection of an element with a trace) !(1:3 = xyz, 1:2 = beginning/end, 1:seg_count = segment index) TYPE(is123), DIMENSION(:,:),ALLOCATABLE :: seg_end_is ! element number and internal coordinates at each end of fault segment !(a fault segment is the intersection of an element with a trace) !(1:2 = beginning/end, 1:seg_count = segment index) REAL, DIMENSION(:), ALLOCATABLE :: seg_eta_ ! eta_ (+1 or -1) for each fault segment, depending on whether ! isolated node u_ is to right or to left. !(1:seg_count = segment index) REAL, DIMENSION(:), ALLOCATABLE :: seg_kappa_ ! kappa_ (relative length) of each fault segment !(1:seg_count = segment index) INTEGER, DIMENSION(:), ALLOCATABLE :: seg_u_ ! u_ = 1, 2, or 3 to identify isolated node of segment !(1:seg_count = segment index) CHARACTER(1), DIMENSION(:), ALLOCATABLE :: sense ! T = thrust, N = normal, D = detachment, R = dextral, L = sinistral ! (1:f_dat_count = offset datum index) REAL, DIMENSION(:,:),ALLOCATABLE :: trace ! all fault traces in one long list; access through trace_loc. ! (1:3 = x,y,z components of unit vector; 1:f_dig_count = in order read) TYPE(is123), DIMENSION(:), ALLOCATABLE :: trace_is ! locations of fault traces in internal coordinates !(1:f_dig_count = in order read) INTEGER, DIMENSION(:,:),ALLOCATABLE :: trace_loc ! gives locations within traces where one trace is found, ! both in terms of digitised points (in traces) and in terms of ! fault segments (in seg_def, seg_end, seg_end_is); ! (1:4 = first, last entry in traces; first, last segment; ! 0:f_highest = trace index (ties f_dat to f_dig)) CHARACTER(1), DIMENSION(:), ALLOCATABLE :: trace_type ! records 1-character fault sense from f_deg, for writing !(1:f_highest = trace index) LOGICAL(1), DIMENSION(:), ALLOCATABLE :: twisted ! a paleomagnetic site is designated twisted if if shares a ! finite element with an active strike-slip fault at any time ! younger than its magnetization age. Twisted data are not ! used as data, and no corresponding model prediction is output. ! A site never becomes un-twisted once it is twisted ! (even though a different finite-element grid is read in). ! (1:p_dat_count = palemagnetic site index) REAL, DIMENSION(:), ALLOCATABLE :: u_flag ! indicator of singularity returned by LSLPB: 0. or 1. !(1:ndof) REAL, DIMENSION(3) :: vec1, vec2 ! temporary Cartesian unit vectors REAL, DIMENSION(:), ALLOCATABLE :: vw0 ! alternating theta (South) and phi (East) velocity components ! at finite element nodes, in m/s, at young end of timestep !(1:2:num_nod = position in solution vector of linear system) REAL, DIMENSION(:), ALLOCATABLE :: vw1 ! alternating theta (South) and phi (East) velocity components ! at finite element nodes, in m/s, at old end of timestep !(1:2:num_nod = position in solution vector of linear system) REAL, DIMENSION(:), ALLOCATABLE :: vw_add ! alternating theta (South) and phi (East) velocity component ! increments at finite lement nodes, in m/s, for the corrector !(1:2:num_nod = position in solution vector of linear system) REAL, DIMENSION(:), ALLOCATABLE :: vw_mean ! alternating theta (South) and phi (East) velocity components ! at finite element nodes, in m/s, averaging predictor & corrector !(1:2:num_nod = position in solution vector of linear system) INTEGER, DIMENSION(:), ALLOCATABLE :: which_trace ! which trace index goes with this fault offset datum? ! use value to read actual trace location from trace_loc ! (1:f_dat_count = fault offset datum index) REAL, DIMENSION(:,:),ALLOCATABLE :: xyz_nod ! positions of node of the finite-element grid (integrated), ! as Cartesian unit vectors from center of a unit sphere. !(1:3 = x,y,z; 1:num_nod = node index) !-------------------------------------------------------------------- ! ! Initialize array of relative probabilities of angular discrepancies, ! as a function of angular (arc) distance between stress indicators. ! (Bird & Li, 1996, J. Geophys. Res., v.101, #B3, 5435-5443) ! We formed 21 annuli for 0-22 deg. angular distance, using ! epsilon = 0.4 and 150 bins for the whole 180-degree range.. ! We used 30 3-degree sectors for beta, covering 0-90 degrees. ! Using all 6000 data of the World Stress Map (Zoback, 1992), we ! counted the frequency of certain differences (beta). ! In each annulus (row), we divided by the total and multiplied by ! 30 to get relative probabilities of order 1. ! Finally, we took the natural logarithm of all numbers, to speed ! formation of extended products. ln_rel_prob( 1,0:29) = & (/+0.790,+0.682,+0.639,+0.606,+0.553,+0.478,+0.449,+0.349,+0.318,+0.216,+0.179,+0.086,-0.006,-0.063,-0.192,& -0.261,-0.320,-0.401,-0.440,-0.546,-0.607,-0.644,-0.725,-0.727,-0.791,-0.791,-0.790,-0.781,-0.869,-0.843/) ln_rel_prob( 2,0:29) = & (/+0.522,+0.474,+0.449,+0.440,+0.378,+0.351,+0.317,+0.262,+0.252,+0.171,+0.110,+0.091,+0.015,-0.018,-0.088,& -0.134,-0.185,-0.244,-0.280,-0.309,-0.379,-0.408,-0.442,-0.394,-0.465,-0.447,-0.450,-0.456,-0.475,-0.462/) ln_rel_prob( 3,0:29) = & (/+0.416,+0.409,+0.379,+0.360,+0.357,+0.291,+0.276,+0.225,+0.191,+0.149,+0.123,+0.077,+0.021,-0.022,-0.045,& -0.109,-0.145,-0.180,-0.214,-0.256,-0.262,-0.295,-0.311,-0.331,-0.379,-0.381,-0.378,-0.387,-0.379,-0.433/) ln_rel_prob( 4,0:29) = & (/+0.361,+0.353,+0.312,+0.318,+0.296,+0.274,+0.256,+0.208,+0.191,+0.150,+0.112,+0.055,+0.031,-0.010,-0.046,& -0.067,-0.144,-0.143,-0.191,-0.262,-0.248,-0.311,-0.284,-0.342,-0.293,-0.270,-0.266,-0.309,-0.302,-0.347/) ln_rel_prob( 5,0:29) = & (/+0.273,+0.275,+0.267,+0.246,+0.239,+0.220,+0.210,+0.150,+0.126,+0.104,+0.094,+0.062,+0.051,-0.003,-0.028,& -0.054,-0.081,-0.109,-0.130,-0.158,-0.189,-0.217,-0.214,-0.263,-0.233,-0.225,-0.239,-0.254,-0.253,-0.210/) ln_rel_prob( 6,0:29) = & (/+0.343,+0.332,+0.346,+0.305,+0.270,+0.241,+0.213,+0.182,+0.186,+0.143,+0.117,+0.080,+0.038,+0.031,-0.025,& -0.056,-0.092,-0.127,-0.143,-0.200,-0.219,-0.289,-0.243,-0.303,-0.294,-0.315,-0.332,-0.343,-0.359,-0.358/) ln_rel_prob( 7,0:29) = & (/+0.327,+0.387,+0.314,+0.279,+0.294,+0.240,+0.236,+0.193,+0.173,+0.133,+0.103,+0.084,+0.060,+0.018,-0.047,& -0.075,-0.113,-0.123,-0.160,-0.223,-0.241,-0.295,-0.280,-0.289,-0.252,-0.313,-0.291,-0.341,-0.333,-0.329/) ln_rel_prob( 8,0:29) = & (/+0.290,+0.264,+0.282,+0.259,+0.279,+0.226,+0.203,+0.203,+0.163,+0.125,+0.092,+0.084,+0.062,+0.018,-0.009,& -0.020,-0.076,-0.118,-0.137,-0.143,-0.223,-0.212,-0.246,-0.267,-0.315,-0.276,-0.291,-0.277,-0.296,-0.333/) ln_rel_prob( 9,0:29) = & (/+0.259,+0.254,+0.252,+0.290,+0.223,+0.221,+0.215,+0.195,+0.172,+0.135,+0.098,+0.071,+0.053,-0.004,+0.007,& -0.062,-0.058,-0.078,-0.135,-0.157,-0.160,-0.205,-0.233,-0.259,-0.263,-0.263,-0.281,-0.265,-0.324,-0.323/) ln_rel_prob(10,0:29) = & (/+0.213,+0.208,+0.198,+0.190,+0.192,+0.195,+0.152,+0.158,+0.153,+0.121,+0.105,+0.106,+0.084,+0.048,+0.039,& +0.005,-0.040,-0.057,-0.109,-0.145,-0.158,-0.184,-0.215,-0.228,-0.237,-0.258,-0.230,-0.204,-0.266,-0.289/) ln_rel_prob(11,0:29) = & (/+0.179,+0.139,+0.164,+0.204,+0.183,+0.167,+0.164,+0.158,+0.157,+0.099,+0.071,+0.042,+0.086,+0.031,+0.000,& -0.001,-0.041,-0.069,-0.084,-0.118,-0.144,-0.170,-0.180,-0.192,-0.181,-0.181,-0.194,-0.193,-0.223,-0.201/) ln_rel_prob(12,0:29) = & (/+0.128,+0.192,+0.136,+0.158,+0.159,+0.132,+0.102,+0.118,+0.125,+0.139,+0.109,+0.038,+0.093,+0.033,-0.024,& +0.002,-0.056,-0.067,-0.091,-0.100,-0.178,-0.119,-0.147,-0.146,-0.164,-0.172,-0.156,-0.141,-0.196,-0.158/) ln_rel_prob(13,0:29) = & (/+0.078,+0.114,+0.085,+0.128,+0.090,+0.126,+0.089,+0.079,+0.046,+0.040,+0.037,+0.036,+0.010,-0.033,+0.019,& -0.072,-0.051,-0.071,-0.064,-0.087,-0.071,-0.065,-0.072,-0.077,-0.071,-0.067,-0.040,-0.056,-0.075,-0.085/) ln_rel_prob(14,0:29) = & (/+0.098,+0.069,+0.104,+0.106,+0.117,+0.099,+0.113,+0.130,+0.080,+0.047,+0.080,+0.042,+0.024,+0.021,-0.004,& -0.001,-0.049,-0.068,-0.032,-0.053,-0.060,-0.081,-0.110,-0.096,-0.080,-0.096,-0.115,-0.116,-0.139,-0.146/) ln_rel_prob(15,0:29) = & (/+0.187,+0.133,+0.181,+0.145,+0.138,+0.095,+0.140,+0.128,+0.083,+0.067,+0.085,+0.052,+0.051,+0.038,-0.010,& -0.066,+0.004,-0.021,-0.075,-0.098,-0.105,-0.150,-0.121,-0.106,-0.162,-0.149,-0.158,-0.174,-0.188,-0.160/) ln_rel_prob(16,0:29) = & (/+0.156,+0.172,+0.132,+0.139,+0.155,+0.117,+0.097,+0.087,+0.088,+0.107,+0.057,+0.056,+0.038,+0.026,+0.001,& -0.016,-0.038,-0.079,-0.077,-0.081,-0.105,-0.113,-0.077,-0.145,-0.143,-0.160,-0.157,-0.138,-0.139,-0.150/) ln_rel_prob(17,0:29) = & (/+0.130,+0.136,+0.131,+0.136,+0.107,+0.134,+0.119,+0.109,+0.123,+0.075,+0.027,+0.018,+0.046,+0.041,+0.029,& -0.036,-0.025,-0.028,-0.081,-0.054,-0.091,-0.070,-0.122,-0.113,-0.144,-0.117,-0.190,-0.162,-0.155,-0.145/) ln_rel_prob(18,0:29) = & (/+0.118,+0.097,+0.089,+0.092,+0.074,+0.089,+0.108,+0.122,+0.103,+0.104,+0.030,+0.053,+0.049,+0.055,-0.010,& +0.010,-0.051,-0.048,-0.056,-0.097,-0.077,-0.146,-0.086,-0.107,-0.098,-0.093,-0.107,-0.098,-0.122,-0.117/) ln_rel_prob(19,0:29) = & (/+0.135,+0.067,+0.070,+0.114,+0.110,+0.103,+0.084,+0.040,+0.102,+0.038,+0.011,+0.056,+0.013,+0.003,+0.005,& -0.024,+0.004,-0.059,-0.008,-0.022,-0.080,-0.064,-0.028,-0.144,-0.089,-0.083,-0.125,-0.123,-0.092,-0.108/) ln_rel_prob(20,0:29) = & (/+0.082,+0.047,+0.077,+0.054,+0.066,+0.051,+0.078,+0.049,+0.042,+0.010,+0.055,+0.015,+0.047,+0.049,+0.020,& +0.033,+0.020,-0.008,-0.002,-0.017,-0.044,-0.036,-0.081,-0.050,-0.070,-0.107,-0.091,-0.112,-0.098,-0.139/) ln_rel_prob(21,0:29) = & (/+0.051,+0.076,+0.037,+0.028,+0.037,+0.049,+0.027,+0.043,+0.031,+0.014,+0.022,+0.022,+0.023,+0.001,+0.005,& +0.020,-0.014,-0.001,-0.004,-0.006,-0.038,-0.036,-0.036,-0.067,-0.037,-0.034,-0.028,-0.042,-0.071,-0.096/) !-------------------------------------------------------------------- memory = 0 ! total of allocated arrays, in bytes ! write the header and initial time stamp PRINT "(' ')" ! All PRINTs have a space in 1st byte for "carriage-control"! PRINT "(' Starting Restore; for details see REPORT.txt')" OPEN (UNIT = 21, FILE = "REPORT.txt", STATUS = "REPLACE", & ACTION = "WRITE") WRITE (21, "('===========================================================')") WRITE (21, "('A record of one run of program Restore')") WRITE (21, "('(palinspastic restoration by integration of paleotectonics)')") WRITE (21, "('by Peter Bird')") WRITE (21, "(' Dept. of Earth and Space Sciences')") WRITE (21, "(' University of California')") WRITE (21, "(' Los Angeles, CA 90095-1567')") WRITE (21, "(' version 2.1, 5 June 1999')") ! see Relevance.doc WRITE (21, "('-----------------------------------------------------------')") CALL DATE_AND_TIME (date, clock_time, zone, datetimenumber) WRITE (21,"('Run began on ',I4,'.',I2,'.',I2,' at ',I2,':',I2,':',I2)") & datetimenumber(1), datetimenumber(2), datetimenumber(3), & datetimenumber(5), datetimenumber(6), datetimenumber(7) WRITE (21, "('-----------------------------------------------------------')") ! parameter section (with immediate conversions to SI units) PRINT "(' Reading parameters for this run')" WRITE (21, "('Begin reading parameters for this run')") WRITE (21, "(' Paramete[rs].dat = parameter file')") param_name = (/ 'PARAMETE.DAT ', 'PARAMETE.Dat ', 'PARAMETE.dat ', & & 'Paramete.DAT ', 'Paramete.Dat ', 'Paramete.dat ', & & 'paramete.DAT ', 'paramete.Dat ', 'paramete.dat ', & & 'PARAMETERS.DAT', 'PARAMETERS.Dat', 'PARAMETERS.dat', & & 'Parameters.DAT', 'Parameters.Dat', 'Parameters.dat', & & 'parameters.DAT', 'parameters.Dat', 'parameters.dat' /) DO i = 1, 18 OPEN (UNIT = 1, FILE = param_name(i), STATUS = "OLD", ACTION = "READ", & & PAD = "YES", IOSTAT = read_status) IF (read_status == 0) EXIT END DO IF (read_status /= 0) THEN WRITE (21, "('ERROR: Could not locate file PARAMETE[RS].DAT')") WRITE (21, "('(in either upper, lower, or mixed-case)')") STOP ENDIF line = 0 ! start time READ (1, *) t ; line = line + 1 IF (t < 0.) CALL Prevent ('negative age', line, "Parameters.dat") WRITE (21,"(F10.2,' geologic age at which this run starts, in Ma')") t start_time = t * s_per_Ma ! end time READ (1, *) t ; line = line + 1 IF (t < 0.) CALL Prevent ('negative age', line, "Parameters.dat") WRITE (21,"(F10.2,' geologic age of older end of history, in Ma')") t end_time = t * s_per_Ma IF (end_time < start_time) THEN PRINT "(' Error; age in line 2 must be .GE. ( >= ) age in line 1.')" WRITE (21,"('Error; age in line 2 must be .GE. ( >= ) age in line 1.')") STOP END IF neotec = (start_time == end_time) paleotec = .NOT. neotec READ (1, *) t ; line = line + 1 ! time step IF (paleotec) THEN IF (t <= 0.) CALL Prevent ('nonpositive timestep', line, "Parameters.dat") WRITE (21,"(F10.2,' time step, in Ma')") t ELSE WRITE (21,"(F10.2,' (this timestep will NOT be used)')") t t = 0. ENDIF Deltat_ = t * s_per_Ma ! number of refinements READ (1, *) n_refine ; line = line + 1 IF (n_refine < 0) CALL Prevent ('negative refinements', line, "Parameters.dat") WRITE (21,"(I10,' number of refinements of each velocity solution')") n_refine ! previous iterations READ (1, *) past_iterations ; line = line + 1 IF (past_iterations < 0) CALL Prevent ('negative past iterations', line, "Parameters.dat") WRITE (21,"(I10,' number of past iterations of whole history before this run')") past_iterations ! iterations in this run READ (1, *) max_iter ; line = line + 1 IF (max_iter <= 0) CALL Prevent ('nonpositive iteration limit', line, "Parameters.dat") IF (neotec .AND. (max_iter > 1)) THEN PRINT "(' Number of iterations set to 1')" WRITE (21,"('----------------------------------------------------------')") WRITE (21,"('Error: This run consumes no time (start_time == end_time),')") WRITE (21,"('and only computes velocities, so iteration is useless.')") WRITE (21,"('Your requested max_iter of ',I3,' has been set to 1.')") max_iter WRITE (21,"('If you want more precision, increase n_refine in line #4.')") WRITE (21,"('----------------------------------------------------------')") max_iter = 1 END IF IF ((start_time > 0.) .AND. (max_iter > 1)) THEN PRINT "(' Error: Inconsistent parameters; see REPORT.txt.')" WRITE (21,"('----------------------------------------------------------')") WRITE (21,"(' Error: Inconsistent parameters.')") WRITE (21,"('You have requested a starting age > 0. and more than 1 iteration.')") WRITE (21,"('These are inconsistent because they require different data files.')") WRITE (21,"('If your data files are present-day, set start_time to 0.')") WRITE (21,"('If your data files are paleo-data, then set max_iter to 1.')") WRITE (21,"('This run cannot continue.')") WRITE (21,"('----------------------------------------------------------')") STOP END IF ALLOCATE ( adjustments(max_iter) ) ! total number of iterations planned WRITE (21,"(I10,' number of iterations of whole history in this run')") max_iter READ (1, *) last_iter WRITE (21, "(I10,' number of iterations planned (total of all runs)')") last_iter ! watch convergence(?) of rate histories? READ (1, *) watch WRITE (21, "(L10,1X,'that 1-4 rate files will be written/iteration for convergence watch.')") watch ! output map information each timestep of last iteration? READ (1, *) map_set WRITE (21, "(L10,1X,'that 2-4 output files will be written/timestep in last iteration.')") map_set ! strain-rate uncertainty for rigid blocks READ (1, *) mu_ ; line = line + 1 IF (mu_ <= 0.) CALL Prevent ('nonpositive mu_', line, "Parameters.dat") IF (mu_ < SQRT(1.1 * TINY(mu_))) CALL Prevent ('mu_**2 will underflow!', line, "Parameters.dat") IF (mu_ > 1.E-10) CALL Prevent ('unreasonably large mu_', line, "Parameters.dat") WRITE (21,"(1P,E10.2,' standard deviation of nominally zero strainrates (mu_), /s')") mu_ ! set of scale rate uncertainties to start off the iteration READ (1, *) mu_scale ; line = line + 1 WRITE (21, "(1P,E10.2,' scale strain-rate of rigid blocks, /s (for early iterations only)')") mu_scale IF (paleotec.AND.(mu_scale <= 0.)) CALL Prevent ('nonpositive mu_scale', line, "Parameters.dat") READ (1, *) f_scale ; line = line + 1 WRITE (21, "(1P,E10.2,' scale slip-rate sigma of faults, m/s (for early iterations only)')") f_scale IF (paleotec.AND.(f_scale <= 0.)) CALL Prevent ('nonpositive f_scale', line, "Parameters.dat") READ (1, *) c_scale ; line = line + 1 WRITE (21, "(1P,E10.2,' scale rate uncertainty of BCSs, m/s (for early iterations only)')") c_scale IF (paleotec.AND.(c_scale <= 0.)) CALL Prevent ('nonpositive c_scale', line, "Parameters.dat") READ (1, *) p_drift_scale ; line = line + 1 WRITE (21, "(1P,E10.2,' scale N-S drift uncertainty of PM, m/s (for early iterations only)')") p_drift_scale IF (paleotec.AND.(p_drift_scale <= 0.)) CALL Prevent ('nonpositive p_drift_scale', line, "Parameters.dat") READ (1, *) p_spin_scale ; line = line + 1 WRITE (21, "(1P,E10.2,' scale spin uncertainty of PM, radians/s (for early iterations only)')") p_spin_scale IF (paleotec.AND.(p_spin_scale <= 0.)) CALL Prevent ('nonpositive p_spin_scale', line, "Parameters.dat") ! small strain-rate increment (xi_) for imposing stress-directions READ (1, *) xi_ ; line = line + 1 IF (xi_ <= 0.) CALL Prevent ('nonpositive xi_', line, "Parameters.dat") WRITE (21,"(1P,E10.2,' small strain-rate increment (xi_), /s')") xi_ ! radius of planet READ (1, *) t ; line = line + 1 IF (t <= 0.) CALL Prevent ('nonpositive R', line, "Parameters.dat") WRITE (21,"(1P,E10.2,' radius of the planet (R), in km')") t R = t * m_per_km half_R2 = (R**2)/2. ! do new active faults count as sigma_1h data? READ (1, *) faults_give_sigma_1h ; line = line + 1 WRITE (21,"(L10,' that active faults give stress directions')") & faults_give_sigma_1h ! names of additional input files (or, "none ") f_dat = Get_filename (unit = 1) ; line = line + 1 WRITE (21,"(' ',A)") TRIM(f_dat) WRITE (21,"(11X,'preceding line = filename of fault offsets')") IF (f_dat(1:5) == 'none ') THEN READ (1,*) ; line = line + 1 ! read and ignore f_dig = 'skipped' ELSE f_dig = Get_filename (unit = 1) ; line = line + 1 END IF WRITE (21,"(' ',A)") TRIM(f_dig) WRITE (21,"(11X,'preceding line = filename of digitised fault traces')") c_dat = Get_filename (unit = 1) ; line = line + 1 WRITE (21,"(' ',A)") TRIM(c_dat) WRITE (21,"(11X,'preceding line = filename of balanced cross-sections')") p_dat = Get_filename (unit = 1) ; line = line + 1 WRITE (21,"(' ',A)") TRIM(p_dat) WRITE (21,"(11X,'preceding line = filename of paleomagnetic data')") s_dat = Get_filename (unit = 1) ; line = line + 1 WRITE (21,"(' ',A)") TRIM(s_dat) WRITE (21,"(11X,'preceding line = filename of principal stress directions')") y_dig = Get_filename (unit = 1) ; line = line + 1 IF (neotec .OR. (.NOT.map_set)) y_dig = 'skipped' WRITE (21,"(' ',A)") TRIM(y_dig) WRITE (21,"(11X,'preceding line = filename of digitised fiducial lines')") CALL More_mem ('gridname_feg', max_fegs * 80) ALLOCATE ( gridname_feg (max_fegs) ) CALL More_mem ('gridname_bcs', max_fegs * 80) ALLOCATE ( gridname_bcs (max_fegs) ) gridname_feg(1) = Get_filename (unit = 1) ; line = line + 1 WRITE (21,"(' ',A)") TRIM(gridname_feg(1)) WRITE (21,"(11X,'preceding line = filename of finite element grid # 1')") CALL Test_file (name = gridname_feg(1), unit = 2) num_fegs = 1 DO i = 2, max_fegs c80 = Get_filename (unit = 1) ; line = line + 1 j1 = INDEX (c80, '.FEG ') j2 = INDEX (c80, '.feg ') IF ((j1 > 0) .OR. (j2 > 0)) THEN gridname_feg(i) = c80 num_fegs = num_fegs + 1 ELSE EXIT END IF END DO DO i = 2, num_fegs WRITE (21,"(' ',A)") TRIM(gridname_feg(i)) WRITE (21,"(11X,'preceding line = filename of finite element grid #',I2)") i CALL Test_file (name = gridname_feg(i), unit = 2) END DO gridname_bcs(1) = c80 DO i = 2, num_fegs gridname_bcs(i) = Get_filename (unit = 1) ; line = line + 1 END DO DO i = 1,num_fegs WRITE (21,"(' ',A)") TRIM(gridname_bcs(i)) WRITE (21,"(11X,'preceding line = filename of boundary conditions #',I2)") i CALL Test_file (name = gridname_bcs(i), unit = 2) END DO CLOSE (UNIT = 1) ! close PARAMETE[RS].DAT PRINT "(' Successfully read all run parameters')" WRITE (21, "('End Parameter Section')") WRITE (21, "('-----------------------------------------------------------')") ! Decide number of timesteps in advance, to preallocate rate arrays IF (paleotec) THEN num_timesteps = NINT(end_time / Deltat_) ELSE ! neotec num_timesteps = 1 END IF CALL More_mem ('rate_err', 3 * (1 + num_timesteps) * max_iter * bytes_per_real) ALLOCATE ( rate_err(0:2, 0:num_timesteps, 1:max_iter) ) rate_err = 0. ! Iterate the whole history!!! current_feg = 0 ! no grid loaded outer_loop: DO iteration = 1, max_iter IF (paleotec) THEN PRINT "(' Beginning iteration ',I3,' out of ',I3,' (this run), ',I3,' (total)')",& & iteration, max_iter, last_iter WRITE (21,"('Beginning iteration ',I3,' out of ',I3,' (this run), ',I3,' (total)')")& & iteration, max_iter, last_iter END IF IF (paleotec .AND. (start_time == 0.)) THEN total_iterations = iteration + past_iterations ELSE ! neotec, or else a restart total_iterations = MAX(1, past_iterations) END IF IF (paleotec) THEN IF (last_iter > 1) THEN exponent = MIN(1.000, REAL(total_iterations - 1) / REAL(last_iter - 1) ) ELSE exponent = 1.000 END IF PRINT "(' using exponent = ',F6.4)", exponent WRITE (21, "('using exponent = ',F6.4)") exponent ELSE exponent = 1.000 END IF loading = (iteration == 1) .OR. (current_feg /= 1) ! i.e., new .feg IF (loading) THEN ! Read input datasets (if loading; whenever there is a new .feg), ! with immediate conversion of quantities to SI units, except ! geographic positions to Cartesian unit vectors in a unit sphere. PRINT "(' ',4X,'Begin reading data input files')" WRITE (21,"(4X,'Begin reading data input files')") ! read f.dat IF (iteration == 1) THEN IF (f_dat(1:5) == 'none ') THEN f_dat_count = 0 f_highest = 0 ELSE PRINT "(' ',8X,'Reading fault offset data from ',A)", TRIM(f_dat) WRITE (21,"(8X,'Reading fault offset data from ',A)") TRIM(f_dat) OPEN (UNIT = 2, FILE = f_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") READ (2, "(A)") f_dat_format READ (2, "(A)") f_dat_titles ! Skim file and count number of data lines, highest fault index f_dat_count = 0 f_highest = 0 get_offset_lines: DO READ (2, "(A)", IOSTAT = read_status) c134 IF (read_status == 0) THEN ! read was successful IF ((c134(1:1) /= '+') .AND. & (c134(1:1) /= '*') .AND. & (c134(1:1) /= '&') .AND. & (c134(1:1) /= '$')) THEN f_dat_count = f_dat_count + 1 READ (c134, "(1X,I4,1X)") i f_highest = MAX (f_highest, i) END IF ELSE; EXIT get_offset_lines; END IF END DO get_offset_lines CLOSE (UNIT = 2) ! (will be re-read) ! allocate arrays CALL More_mem ('f_active', num_timesteps * f_dat_count * 1) ALLOCATE ( f_active(num_timesteps, f_dat_count) ) CALL More_mem ('which_trace', f_dat_count * bytes_per_int) ALLOCATE ( which_trace (f_dat_count) ) CALL More_mem ('sense', f_dat_count * 1) ALLOCATE ( sense (f_dat_count) ) CALL More_mem ('fault_name', f_dat_count * 50) ALLOCATE ( fault_name (f_dat_count) ) CALL More_mem ('offset', f_dat_count * bytes_per_real) ALLOCATE ( offset (f_dat_count) ) CALL More_mem ('offset_sigma_', f_dat_count * bytes_per_real) ALLOCATE ( offset_sigma_ (f_dat_count) ) CALL More_mem ('f_t_max', f_dat_count * bytes_per_real) ALLOCATE ( f_t_max (f_dat_count) ) CALL More_mem ('f_t_min', f_dat_count * bytes_per_real) ALLOCATE ( f_t_min (f_dat_count) ) CALL More_mem ('f_goal', f_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( f_goal(num_timesteps, f_dat_count) ) CALL More_mem ('f_rate', f_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( f_rate(num_timesteps, f_dat_count) ) f_rate = 0.0 CALL More_mem ('f_rate_sigma_', f_dat_count * bytes_per_real) ALLOCATE ( f_rate_sigma_(f_dat_count) ) CALL More_mem ('f_dat_code', f_dat_count * 1) ALLOCATE ( f_dat_code(f_dat_count) ) CALL More_mem ('f_divide', 2 * f_dat_count * bytes_per_real) ALLOCATE ( f_divide(2, f_dat_count) ) CALL More_mem ('f_err', 3 * max_iter * bytes_per_real) ALLOCATE ( f_err(0:2, max_iter) ) IF (faults_give_sigma_1h) THEN CALL More_mem ('f_new', f_dat_count) ALLOCATE ( f_new (f_dat_count) ) f_new = .TRUE. ! just initializing END IF ! fill arrays OPEN (UNIT = 2, FILE = f_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") ; line = 0 READ (2,*) ; line = line + 1 READ (2,*) ; line = line + 1 read_f_dat: DO i = 1, f_dat_count READ (2, f_dat_format) c6, c50, t1, t2, t3, t4 ; line = line + 1 READ (c6, "(1X,I4,1X)") i1 IF (i1 > f_highest) THEN PRINT "(' Illegally high trace index: ',I6)", i1 WRITE (21,"('Illegally high trace index: ',I6)") i1 STOP END IF which_trace(i) = i1 c = c6(6:6) IF (c == 't') c = 'T' IF (c == 'p') c = 'P' IF (c == 'n') c = 'N' IF (c == 'd') c = 'D' IF (c == 'r') c = 'R' IF (c == 'l') c = 'L' IF (.NOT.((c == 'T') .OR. (c == 'N') .OR. (c == 'R') .OR. (c == 'L') & .OR. (c == 'D') .OR. (c == 'P'))) THEN PRINT "(' Illegal slip sense: ',A1)", c WRITE (21,"('Illegal slip sense: ',A1)") c STOP END IF sense(i) = c fault_name(i) = c50 IF (t1 < 0.) CALL Prevent ('negative offset', line, f_dat) offset(i) = t1 * m_per_km IF (t2 <= 0.) CALL Prevent ('nonpositive sigma_', line, f_dat) offset_sigma_(i) = t2 * m_per_km IF (t3 <= 0.) CALL Prevent ('nonpositive maximum age', line, f_dat) f_t_max(i) = t3 * s_per_Ma IF (t4 < 0.) CALL Prevent ('negative minimum age', line, f_dat) f_t_min(i) = t4 * s_per_Ma IF (t3 <= t4) CALL Prevent ('null age span', line, f_dat) CALL Set_goal_A (index = i, total = offset, & & tmin = f_t_min, tmax = f_t_max, & & checkPD = .TRUE., & ! inputs & goal = f_goal, & ! output & active = f_active)! output f_rate_sigma_(i) = offset_sigma_(i) / (f_t_max(i) - f_t_min(i)) CALL Set_goal_B (index = i, & & checkPD = .TRUE., active = f_active, & & unit = 2, signal = '*', eof = eof, & & conversion = (m_per_km / s_per_Ma), & ! inputs & line = line, & ! modify & rate = f_rate, & ! modify & goal = f_goal) ! modify IF (eof) EXIT read_f_dat END DO read_f_dat PRINT "(' ',8X,I4,' fault-offset data were read')", f_dat_count WRITE (21,"(8X,I4,' fault-offset data were read')") f_dat_count CLOSE (UNIT = 2) ! close f_dat !scan for overlapping time windows; also set f_new? IF (f_dat_count > 1) THEN DO i = 1, f_dat_count - 1 DO j = i + 1, f_dat_count IF (which_trace(i) == which_trace(j)) THEN overlap = MIN (f_t_max(i) - f_t_min(i), & & f_t_max(j) - f_t_min(j), & & f_t_max(i) - f_t_min(j), & & f_t_max(j) - f_t_min(i)) IF (overlap > 0.) THEN WRITE (c6, "(I6)") which_trace(i) !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. DO k = 2, 5 IF (c6(k:k) == ' ') c6(k:k) = '0' END DO c6(1:1) = 'F' PRINT "(' Error: Two or more data concerning trace ',A/& & ' have overlapping time windows.'/& & ' Edit these data to make them contiguous.')", c6 WRITE (21,"('Error: Two or more data concerning trace ',A/& & ' have overlapping time windows.'/& & ' Edit these data to make them contiguous.')") c6 STOP ELSE ! no overlap; set f_new? IF (faults_give_sigma_1h) THEN ! set f_new of later offset FALSE IF (f_t_max(i) < f_t_max(j)) THEN f_new(i) = .FALSE. ELSE f_new(j) = .FALSE. END IF END IF END IF ! overlap / no overlap END IF ! >= 2 data on same trace END DO ! j = i+1, f_dat_count END DO ! i = 1, f_dat_count-1 END IF ! scan for overlapping times ! scan for Demotions (P and N already filled in) DO i = 1, f_dat_count ! may already be P DO j = 1, f_dat_count ! may need to be D IF (i /= j) THEN ! different datum lines IF (which_trace(i) == which_trace(j)) THEN !same trace IF (f_dat_code(i) == 'P') THEN ! one was promoted IF (f_active(1, j)) THEN ! conflict f_dat_code(j) = 'D' ! Demote. f_active(1, j) = .FALSE. f_goal(1, j) = 0. END IF END IF END IF END IF END DO END DO ! scan for Demotions WRITE (21,"(8X,'- - - - - - - - - - - - - - - - - - - - ')") END IF ! IF (f_dat /= 'none') END IF ! IF (iteration == 1) ! read f.dig IF ((f_dig(1:5) == 'none ') .OR. (f_dig(1:8) == 'skipped ')) THEN f_dig_count = 0 f_highest = 0 ELSE PRINT "(' ',8X,'Reading fault traces from ',A)", TRIM(f_dig) WRITE (21,"(8X,'Reading fault traces from ',A)") TRIM(f_dig) IF (iteration == 1) THEN OPEN (UNIT = 3, FILE = f_dig, STATUS = "OLD", ACTION = "READ", & PAD = "YES"); line = 0 ! Skim file and count number of data points f_dig_count = 0 loop_thru: DO READ (3, "(A)", IOSTAT = read_status) c50; line = line + 1 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN f_dig_count = f_dig_count + 1 ELSE IF ((c50(1:1) == 'F') .OR. (c50(1:1) == 'f')) THEN READ (c50,"(1X,I4)", IOSTAT = read_status) i if (read_status == 0) THEN IF (i < 0) CALL Prevent ('negative fault number', line, f_dig) f_highest = MAX (f_highest, i) END IF END IF ELSE; EXIT loop_thru; END IF END DO loop_thru CLOSE (UNIT = 3) ! (will be re-read) ! allocate arrays CALL More_mem ('trace', f_dig_count * 3 * bytes_per_real) ALLOCATE ( trace (3, f_dig_count) ) trace = 0. ! whole array ! CALL More_mem ('trace_is', f_dig_count * bytes_per_is) ALLOCATE ( trace_is(f_dig_count) ) CALL More_mem ('trace_loc', 4 * f_highest * bytes_per_int) ALLOCATE ( trace_loc (4, 0:f_highest) ) trace_loc = 0 ! whole array! CALL More_mem ('trace_type', f_highest * 1) ALLOCATE ( trace_type(f_highest) ) trace_type = ' ' CALL More_mem ('f_2_in', f_highest * 1) ALLOCATE ( f_2_in(f_highest) ) f_2_in = .TRUE. ! later statements can only change -> .FALSE. END IF ! IF (iteration == 1) ! (re)fill arrays OPEN (UNIT = 3, FILE = f_dig, STATUS = "OLD", ACTION = "READ", & PAD = "YES") ; line = 0 in_trace = .FALSE. i = 0 read_dig: DO READ (3, "(A)", IOSTAT = read_status) c50; line = line + 1 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN got_point = .TRUE. got_index = .FALSE. READ (c50,*) t1, t2 ! E longitude, N latitude IF (ABS(t2) > 90.001) THEN PRINT "(' Bad latitude ',F10.2,' in line ',I10,' of ',A)", & t2, line, TRIM(f_dig) WRITE (21,"('Bad latitude ',F10.2,' in line ',I10,' of ',A)") & t2, line, TRIM(f_dig) STOP END IF ELSE IF ((c50(1:1) == 'F') .OR. (c50(1:1) == 'f')) THEN got_point = .FALSE. READ (c50,"(1X,I4,A)") j2, c1 ! new trace number and sense got_index = .TRUE. trace_type(j2) = c1 ELSE ! most likely, '*** end of line segment ***' got_point = .FALSE. got_index = .FALSE. END IF ELSE; EXIT read_dig; END IF IF (in_trace) THEN IF (got_point) THEN i = i + 1 CALL Xyz_from_lonlat(t1, t2, tvo) trace(1:3,i) = tvo trace_loc(2,j1) = i ELSE ! *** end ... in_trace = .FALSE. ENDIF ELSE ! (not in_trace) IF (got_index) THEN j1 = j2 ! new index becomes current index ELSE IF (got_point) THEN i = i + 1 CALL Xyz_from_lonlat(t1, t2, tv) trace(1:3,i) = tv trace_loc(1,j1) = i in_trace = .TRUE. ELSE ! *** end ... END IF END IF ! (in_trace) END DO read_dig CLOSE (UNIT = 3) ! close f_dig PRINT "(' ',8X,I6,' fault-trace points were read')", f_dig_count WRITE (21,"(8X,I6,' fault-trace points were read')") f_dig_count WRITE (21,"(8X,'- - - - - - - - - - - - - - - - - - - - - ')") IF (iteration == 1) THEN ! Check that all necessary traces where actually read. j = 0 ! number of critical traces missing from f_dig DO i = 1, f_dat_count IF (trace_loc(2,which_trace(i)) == 0) j = j + 1 END DO IF (j > 0) THEN PRINT "(' Error: The following fault traces were missing:')" WRITE (21,"('Error: The following fault traces were missing:')") DO i = 1, f_dat_count IF (trace_loc(2,which_trace(i)) == 0) THEN WRITE (c4,"(I4)") which_trace(i) !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. DO j=1,4 IF (c4(j:j) == ' ') c4(j:j) = '0' END DO PRINT "(' F',A4,A1,' ',A)", & c4, sense(i), TRIM(fault_name(i)) WRITE(21,"(' F',A4,A1,' ',A)") & c4, sense(i), TRIM(fault_name(i)) END IF END DO STOP END IF ! IF (j > 0) END IF ! IF (iteration == 1) END IF ! IF (f_dig /= 'none'.OR. 'skipped') ! read c_dat IF (iteration == 1) THEN IF (c_dat(1:5) == 'none ') THEN c_dat_count = 0 ELSE PRINT "(' ',8X,'Reading cross-section data from ',A)", TRIM(c_dat) WRITE (21,"(8X,'Reading cross-section data from ',A)") TRIM(c_dat) OPEN (UNIT = 4, FILE = c_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") READ (4, "(A)") c_dat_format READ (4, "(A)") c_dat_titles ! Skim file and count number of data lines c_dat_count = 0 get_section_lines: DO READ (4, "(A)", IOSTAT = read_status) c134 IF (read_status == 0) THEN ! read was successful IF ((c134(1:1) /= '+') .AND. & (c134(1:1) /= '*') .AND. & (c134(1:1) /= '&') .AND. & (c134(1:1) /= '$')) THEN c_dat_count = c_dat_count + 1 END IF ELSE; EXIT get_section_lines; END IF END DO get_section_lines CLOSE (UNIT = 4) ! (will be re-read) ! allocate arrays CALL More_mem ('c_active', num_timesteps * c_dat_count * 1) ALLOCATE ( c_active(num_timesteps, c_dat_count) ) CALL More_mem ('c_ref', c_dat_count * 47) ALLOCATE ( c_ref(c_dat_count) ) CALL More_mem ('c_end_0', 3 * 2 * c_dat_count * bytes_per_real) ALLOCATE ( c_end_0(3, 2, c_dat_count) ) CALL More_mem ('c_end_is', 2 * c_dat_count * bytes_per_is) ALLOCATE ( c_end_is(2, c_dat_count) ) CALL More_mem ('c_end_now', 3 * 2 * c_dat_count * bytes_per_real) ALLOCATE ( c_end_now(3, 2, c_dat_count) ) CALL More_mem ('c_code', c_dat_count * 5) ALLOCATE ( c_code(c_dat_count) ) CALL More_mem ('c_length', 2 * c_dat_count * bytes_per_real) ALLOCATE ( c_length(2, c_dat_count) ) CALL More_mem ('c_stretch', c_dat_count * bytes_per_real) ALLOCATE ( c_stretch(c_dat_count) ) CALL More_mem ('c_sigma_', c_dat_count * bytes_per_real) ALLOCATE ( c_sigma_(c_dat_count) ) CALL More_mem ('c_t_max', c_dat_count * bytes_per_real) ALLOCATE ( c_t_max (c_dat_count) ) CALL More_mem ('c_t_min', c_dat_count * bytes_per_real) ALLOCATE ( c_t_min (c_dat_count) ) CALL More_mem ('c_goal', c_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( c_goal(num_timesteps, c_dat_count) ) CALL More_mem ('c_rate', c_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( c_rate(num_timesteps, c_dat_count) ) c_rate = 0.0 CALL More_mem ('c_rate_sigma_', c_dat_count * bytes_per_real) ALLOCATE ( c_rate_sigma_(c_dat_count) ) CALL More_mem ('c_err', 3 * max_iter * bytes_per_real) ALLOCATE ( c_err(0:2, 1:max_iter) ) ! fill arrays OPEN (UNIT = 4, FILE = c_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") ; line = 0 READ (4,*) ; line = line + 1 READ (4,*) ; line = line + 1 read_c_dat: DO i = 1, c_dat_count READ (4, c_dat_format) c47, x1, x2, x3, x4, c5, r1, r2, r3, t2, t1 ; line = line + 1 c_ref(i) = c47 CALL Xyz_from_lonlat(x1, x2, tv) c_end_0(1:3,1,i) = tv CALL Xyz_from_lonlat(x3, x4, tv) c_end_0(1:3,2,i) = tv c_end_now(1:3,1:2,i) = c_end_0(1:3,1:2,i) ! possibly overwritten below c_code(i) = c5 c_length(1,i) = r1 * m_per_km c_length(2,i) = r2 * m_per_km c_stretch(i) = c_length(1,i) - c_length(2,i) IF (r3 <= 0.) CALL Prevent ('nonpositive sigma_', line, c_dat) c_sigma_(i) = r3 * m_per_km IF (t2 <= t1) CALL Prevent ('null age range', line, c_dat) c_t_max(i) = t2 * s_per_Ma c_t_min(i) = t1 * s_per_Ma CALL Set_goal_A (index = i, total = c_stretch, & & tmin = c_t_min, tmax = c_t_max, checkPD = .FALSE., & ! inputs & goal = c_goal, active = c_active)! outputs c_rate_sigma_(i) = c_sigma_(i) / (c_t_max(i) - c_t_min(i)) !read any + lines, but use only if start_time > 0. READ (4,"(A)", IOSTAT = read_status) c134 IF (read_status /= 0) EXIT read_c_dat if (c134(1:1) == '+') THEN c134 = c134(2:134) // ' ' READ (c134,*) x1, x2 READ (4,"(A)") c134 c134 = c134(2:134) // ' ' READ (c134,*) x3, x4 IF (start_time > 0.) THEN CALL Xyz_from_lonlat(x1, x2, tv) c_end_now(1:3,1,i) = tv CALL Xyz_from_lonlat(x3, x4, tv) c_end_now(1:3,2,i) = tv END IF ELSE BACKSPACE(4) END IF !read any * lines CALL Set_goal_B (index = i, & & checkPD = .FALSE., active = c_active, & & unit = 4, signal = '*', eof = eof, & & conversion = (m_per_km / s_per_Ma), & ! inputs & line = line, & ! modify & rate = c_rate, & ! modify & goal = c_goal) ! modify IF (eof) EXIT read_c_dat END DO read_c_dat PRINT "(' ',8X,I4,' cross-section data were read')", c_dat_count WRITE (21,"(8X,I4,' cross-section data were read')") c_dat_count CLOSE (UNIT = 4) ! close c_dat WRITE (21,"(8X,'- - - - - - - - - - - - - - - - - - - - ')") END IF ! (c_dat /= 'none') END IF ! (iteration == 1) ! read p_dat IF (iteration == 1) THEN IF (p_dat(1:5) == 'none ') THEN p_dat_count = 0 ELSE PRINT "(' ',8X,'Reading paleomagnetic data from ',A)", TRIM(p_dat) WRITE (21,"(8X,'Reading paleomagnetic data from ',A)") TRIM(p_dat) OPEN (UNIT = 8, FILE = p_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") READ (8, "(A)") p_dat_format READ (8, "(A)") p_dat_titles ! Skim file and count number of data lines p_dat_count = 0 get_paleo_lines: DO READ (8, "(A)", IOSTAT = read_status) c134 IF (read_status == 0) THEN ! read was successful IF ((c134(1:1) /= '+') .AND. & (c134(1:1) /= '*') .AND. & (c134(1:1) /= '&') .AND. & (c134(1:1) /= '$')) THEN p_dat_count = p_dat_count + 1 END IF ELSE; EXIT get_paleo_lines; END IF END DO get_paleo_lines CLOSE (UNIT = 8) ! (will be re-read) ! allocate arrays CALL More_mem ('p_active', num_timesteps * p_dat_count * 1) ALLOCATE ( p_active(num_timesteps, p_dat_count) ) CALL More_mem ('p_ref', p_dat_count * 50) ALLOCATE ( p_ref(p_dat_count) ) CALL More_mem ('p_site_0', 3 * p_dat_count * bytes_per_real) ALLOCATE ( p_site_0(3, p_dat_count) ) CALL More_mem ('p_site_is', p_dat_count * bytes_per_is) ALLOCATE ( p_site_is(p_dat_count) ) CALL More_mem ('p_site_now', 3 * p_dat_count * bytes_per_real) ALLOCATE ( p_site_now(3, p_dat_count) ) CALL More_mem ('p_south', p_dat_count * bytes_per_real) ALLOCATE ( p_south(p_dat_count) ) CALL More_mem ('p_south_sigma_', p_dat_count * bytes_per_real) ALLOCATE ( p_south_sigma_(p_dat_count) ) CALL More_mem ('p_ccw', p_dat_count * bytes_per_real) ALLOCATE ( p_ccw(p_dat_count) ) CALL More_mem ('p_ccw_sigma_', p_dat_count * bytes_per_real) ALLOCATE ( p_ccw_sigma_(p_dat_count) ) CALL More_mem ('p_t2', p_dat_count * bytes_per_real) ALLOCATE ( p_t2(p_dat_count) ) CALL More_mem ('p_t1', p_dat_count * bytes_per_real) ALLOCATE ( p_t1(p_dat_count) ) CALL More_mem ('p_t_max', p_dat_count * bytes_per_real) ALLOCATE ( p_t_max(p_dat_count) ) CALL More_mem ('p_t_min', p_dat_count * bytes_per_real) ALLOCATE ( p_t_min(p_dat_count) ) p_t_min = 0. ! whole array; necessary as an actual parameter CALL More_mem ('p_pole', 3 * p_dat_count * bytes_per_real) ALLOCATE ( p_pole(3, p_dat_count) ) CALL More_mem ('p_south_goal', p_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( p_south_goal(num_timesteps, p_dat_count) ) CALL More_mem ('p_south_rate', p_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( p_south_rate(num_timesteps, p_dat_count) ) p_south_rate = 0.0 CALL More_mem ('p_south_rate_sigma_', p_dat_count * bytes_per_real) ALLOCATE ( p_south_rate_sigma_(p_dat_count) ) CALL More_mem ('p_ccw_goal', p_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( p_ccw_goal(num_timesteps, p_dat_count) ) CALL More_mem ('p_ccw_rate', p_dat_count * num_timesteps * bytes_per_real) ALLOCATE ( p_ccw_rate(num_timesteps, p_dat_count) ) p_ccw_rate = 0.0 CALL More_mem ('p_ccw_rate_sigma_', p_dat_count * bytes_per_real) ALLOCATE ( p_ccw_rate_sigma_(p_dat_count) ) CALL More_mem ('twisted', p_dat_count) ALLOCATE ( twisted(p_dat_count) ) CALL More_mem ('p_south_err', 3 * max_iter * bytes_per_real) ALLOCATE ( p_south_err(0:2,max_iter) ) CALL More_mem ('p_ccw_err', 3 * max_iter * bytes_per_real) ALLOCATE ( p_ccw_err(0:2,max_iter) ) twisted = .FALSE. ! whole array; later statements only -> .TRUE. ! fill arrays OPEN (UNIT = 8, FILE = p_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") ; line = 0 READ (8,*) ; line = line + 1 READ (8,*) ; line = line + 1 read_p_dat: DO i = 1, p_dat_count READ (8, p_dat_format) c50, x1, x2, r1, r2, r3, r4, t2, t1, x3, x4; line = line + 1 p_ref(i) = c50 CALL Xyz_from_lonlat(x1, x2, tv) p_site_0(1:3,i) = tv p_site_now(1:3,i) = p_site_0(1:3,i) ! possibly overwritten below p_south(i) = (r1 / deg_per_rad) * R IF (r2 <= 0.) CALL Prevent ('nonpositive sigma_ for latitude', line, p_dat) p_south_sigma_(i) = (r2 / deg_per_rad) * R p_ccw(i) = r3 / deg_per_rad IF (r4 <= 0.) CALL Prevent ('nonpositive sigma_ for rotation', line, p_dat) p_ccw_sigma_(i) = r4 / deg_per_rad IF (t2 < t1) THEN PRINT "(' Error in line ',I6,' of'/' ',A/' max.age ',F10.2,' is < min.age ',F10.2)", & line, TRIM(p_dat), t2, t1 WRITE (21,"('Error in line ',I6,' of'/A/'max.age ',F10.2,' is < min.age ',F10.2)") & line, TRIM(p_dat), t2, t1 STOP ENDIF p_t2(i) = t2 * s_per_Ma p_t1(i) = t1 * s_per_Ma p_t_max(i) = (p_t1(i) + p_t2(i))/2. ! mean age of magnetization p_t_min(i) = 0. ! rocks were sampled only yesterday CALL Xyz_from_lonlat(x3, x4, tv) p_pole(1:3,i) = tv CALL Set_goal_A (index = i, total = p_south, & & tmin = p_t_min, tmax = p_t_max, checkPD = .FALSE., & ! inputs & goal = p_south_goal, active = p_active)! output p_south_rate_sigma_(i) = p_south_sigma_(i) / (p_t_max(i) - p_t_min(i)) CALL Set_goal_A (index = i, total = p_ccw, & & tmin = p_t_min, tmax = p_t_max, checkPD = .FALSE., & ! inputs & goal = p_ccw_goal, active = p_active)! output p_ccw_rate_sigma_(i) = p_ccw_sigma_(i) / (p_t_max(i) - p_t_min(i)) !read any + lines, but use only if start_time > 0. READ (8,"(A)", IOSTAT = read_status) c134 IF (read_status /= 0) EXIT read_p_dat if (c134(1:1) == '+') THEN c134 = c134(2:134) // ' ' READ (c134,*) x1, x2 IF (start_time > 0.) THEN CALL Xyz_from_lonlat(x1, x2, tv) p_site_now(1:3,i) = tv END IF ELSE BACKSPACE(8) END IF !read any * and # lines and correct rates and goals CALL Set_goal_B (index = i, & & checkPD = .FALSE., active = p_active, & & unit = 8, signal = '*', eof = eof, & & conversion = (R / (deg_per_rad * s_per_Ma)), & ! inputs & line = line, & ! modify & rate = p_south_rate, & ! modify & goal = p_south_goal) ! modify IF (eof) EXIT read_p_dat CALL Set_goal_B (index = i, & & checkPD = .FALSE., active = p_active, & & unit = 8, signal = '&', eof = eof, & & conversion = (1. / (deg_per_rad * s_per_Ma)), & ! inputs & line = line, & ! modify & rate = p_ccw_rate, & ! modify & goal = p_ccw_goal) ! modify IF (eof) EXIT read_p_dat END DO read_p_dat PRINT "(' ',8X,I4,' paleomagnetic sites were read')", p_dat_count WRITE (21,"(8X,I4,' paleomagnetic sites were read')") p_dat_count CLOSE (UNIT = 8) ! close p_dat WRITE (21,"(8X,'- - - - - - - - - - - - - - - - - - - - ')") END IF ! IF (p_dat /= 'none') END IF ! IF (iteration == 1) ! read s_dat IF (iteration == 1) THEN IF (s_dat(1:5) == 'none ') THEN s_dat_count = 0 ELSE PRINT "(' ',8X,'Reading paleostress data from ',A)", TRIM(s_dat) WRITE (21,"(8X,'Reading paleostress data from ',A)") TRIM(s_dat) OPEN (UNIT = 9, FILE = s_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") READ (9, "(A)") s_dat_format READ (9, "(A)") s_dat_titles ! Skim file and count number of data lines s_dat_count = 0 get_stress_lines: DO READ (9, "(A)", IOSTAT = read_status) c134 IF (read_status == 0) THEN ! read was successful IF ((c134(1:1) /= '+') .AND. & (c134(1:1) /= '*') .AND. & (c134(1:1) /= '&') .AND. & (c134(1:1) /= '$')) THEN s_dat_count = s_dat_count + 1 END IF ELSE; EXIT get_stress_lines; END IF END DO get_stress_lines CLOSE (UNIT = 9) ! (will be re-read) ! allocate arrays CALL More_mem ('s_activity', num_timesteps * s_dat_count * bytes_per_real) ALLOCATE (s_activity(num_timesteps, s_dat_count) ) CALL More_mem ('s_ref', s_dat_count * 30) ALLOCATE ( s_ref(s_dat_count) ) CALL More_mem ('s_loc', s_dat_count * 30) ALLOCATE ( s_loc(s_dat_count) ) CALL More_mem ('s_code', s_dat_count * 5) ALLOCATE ( s_code(s_dat_count) ) CALL More_mem ('s_site_0', 3 * s_dat_count * bytes_per_real) ALLOCATE ( s_site_0(3, s_dat_count) ) CALL More_mem ('s_site_is', 2 * s_dat_count * bytes_per_is) ALLOCATE (s_site_is(2, s_dat_count) ) CALL More_mem ('s_site_now', 3 * 2 * s_dat_count * bytes_per_real) ALLOCATE ( s_site_now(3, 2, s_dat_count) ) CALL More_mem ('s_azim_0', s_dat_count * bytes_per_real) ALLOCATE ( s_azim_0(s_dat_count) ) CALL More_mem ('s_azim_now', s_dat_count * bytes_per_real) ALLOCATE ( s_azim_now(s_dat_count) ) CALL More_mem ('s_sigma_', s_dat_count * bytes_per_real) ALLOCATE ( s_sigma_(s_dat_count) ) CALL More_mem ('s_t_max', s_dat_count * bytes_per_real) ALLOCATE ( s_t_max(s_dat_count) ) CALL More_mem ('s_t_min', s_dat_count * bytes_per_real) ALLOCATE ( s_t_min(s_dat_count) ) CALL More_mem ('s_stage', s_dat_count * 1) ALLOCATE ( s_stage(s_dat_count) ) ! fill arrays OPEN (UNIT = 9, FILE = s_dat, STATUS = "OLD", ACTION = "READ", & PAD = "YES") ; line = 0 READ (9,*) ; line = line + 1 READ (9,*) ; line = line + 1 read_s_dat: DO i = 1, s_dat_count READ (9, s_dat_format) c30, c30a, c5, x1, x2, r1, r2, t2, t1, c6; line = line + 1 s_ref(i) = c30 s_loc(i) = c30a s_code(i) = c5 CALL Xyz_from_lonlat(x1, x2, tv) s_site_0(1:3,i) = tv s_azim_0(i) = r1 / deg_per_rad s_azim_now(i) = s_azim_0(i) ! possibly overwritten below s_site_now(1:3,1,i) = s_site_0(1:3,i) ! possibly overwritten below tvi = s_site_now(1:3, 1, i) CALL Step_aside(tvi, s_azim_0(i), tvo) s_site_now(1:3, 2, i) = tvo IF (r2 <= 0.) CALL Prevent ('nonpositive sigma_', line, s_dat) s_sigma_(i) = r2 / deg_per_rad IF (t2 <= t1) THEN PRINT "(' Error in line ',I6,' of'/' ',A/' ',F10.2,' is <= ',F10.2)", & line, TRIM(s_dat), t2, t1 WRITE (21,"('Error in line ',I6,' of'/A/F10.2,' is <= ',F10.2)") & line, TRIM(s_dat), t2, t1 STOP ENDIF s_t_max(i) = t2 * s_per_Ma s_t_min(i) = t1 * s_per_Ma s_stage(i) = (c6(1:1) == 'S') .OR. (c6(1:1) == 's') IF (paleotec) THEN DO j = 1, num_timesteps t1 = (j - 1) * Deltat_ t2 = j * Deltat_ overlap = MAX(0.0, MIN(t2 - t1, s_t_max(i) - s_t_min(i), & & t2 - s_t_min(i), s_t_max(i) - t1)) IF (s_stage(i)) THEN overlap_threshold = MIN(0.1 * s_per_Ma, 0.5 * Deltat_) IF (overlap >= overlap_threshold) THEN s_activity(j, i) = 1.0 ELSE s_activity(j, i) = 0.0 END IF ELSE ! "window" type datum s_activity(j, i) = overlap / (s_t_max(i) - s_t_min(i)) END IF END DO ELSE ! neotec IF (s_stage(i)) THEN allowance = 0.1 * s_per_Ma IF ((start_time > (s_t_min(i) - allowance)) .AND. & & (start_time < (s_t_max(i) + allowance))) THEN s_activity(1, i) = 1.0 ELSE s_activity(1, i) = 0.0 END IF ELSE ! "window" type datum s_activity(1, i) = 0.0 END IF ! stage(i) END IF ! paleotec or neotec !read any + line, but use only if start_time > 0. READ (9,"(A)", IOSTAT = read_status) c134 IF (read_status /= 0) EXIT read_s_dat if (c134(1:1) == '+') THEN c134 = c134(2:134) // ' ' READ (c134,*) x1, x2 IF (start_time > 0.) THEN CALL Xyz_from_lonlat(x1, x2, tv) s_site_now(1:3,1,i) = tv END IF ELSE BACKSPACE(9) END IF !read any $ line, but use only if start_time > 0. READ (9,"(A)", IOSTAT = read_status) c134 IF (read_status /= 0) EXIT read_s_dat if (c134(1:1) == '$') THEN c134 = c134(2:134) // ' ' READ (c134,*) t1 IF (start_time > 0.) THEN s_azim_now(i) = t1 / deg_per_rad tvi = s_site_now(1:3, 1, i) CALL Step_aside(tvi, s_azim_now(i), tvo) s_site_now(1:3,2,i) = tvo END IF ELSE BACKSPACE(9) END IF END DO read_s_dat PRINT "(' ',8X,I4,' paleostress sites were read')", s_dat_count WRITE (21,"(8X,I4,' paleostress sites were read')") s_dat_count CLOSE (UNIT = 9) ! close s_dat WRITE (21,"(8X,'- - - - - - - - - - - - - - - - - - - - ')") END IF ! IF (s_dat /= 'none') END IF ! IF (iteration == 1) ! read y_dig IF ((y_dig(1:5) == 'none ') .OR. (y_dig(1:7) == 'skipped')) THEN y_dig_count = 0 ELSE PRINT "(' ',8X,'Reading fiducial lines from ',A)", TRIM(y_dig) WRITE (21,"(8X,'Reading fiducial lines from ',A)") TRIM(y_dig) IF (iteration == 1) THEN OPEN (UNIT = 10, FILE = y_dig, STATUS = "OLD", ACTION = "READ", & PAD = "YES"); line = 0 ! Skim file and count number of data points y_dig_count = 0 y_loop_thru: DO READ (10, "(A)", IOSTAT = read_status) c50; line = line + 1 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN y_dig_count = y_dig_count + 1 END IF ELSE; EXIT y_loop_thru; END IF END DO y_loop_thru CLOSE (UNIT = 10) ! (will be re-read) ! allocate arrays CALL More_mem ('dot', 3 * y_dig_count * bytes_per_real) ALLOCATE ( dot(3, y_dig_count) ) CALL More_mem ('dot_is', y_dig_count * bytes_per_is) ALLOCATE ( dot_is(y_dig_count) ) CALL More_mem ('dot_last', (y_dig_count + 1) * 1) ALLOCATE ( dot_last(0:y_dig_count) ) END IF ! IF (iteration == 1) ! fill arrays OPEN (UNIT = 10, FILE = y_dig, STATUS = "OLD", ACTION = "READ", & PAD = "YES") ; line = 0 i = 0 ! will be used to count points again y_read_dig: DO READ (10, "(A)", IOSTAT = read_status) c50; line = line + 1 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN got_point = .TRUE. READ (c50,*) t1, t2 ! E longitude, N latitude IF (ABS(t2) > 90.001) THEN PRINT "(' Bad latitude ',F10.2,' in line ',I10,' of ',A)", & t2, line, TRIM(y_dig) WRITE (21,"('Bad latitude ',F10.2,' in line ',I10,' of ',A)") & t2, line, TRIM(y_dig) STOP END IF ELSE ! titles, or '*** end of line segment ***' got_point = .FALSE. END IF ELSE; EXIT y_read_dig; END IF IF (got_point) THEN i = i + 1 CALL Xyz_from_lonlat(t1, t2, tv) dot(1:3,i) = tv dot_last(i) = .FALSE. ELSE ! titles, or '*** end ...' dot_last(i) = .TRUE. ENDIF END DO y_read_dig CLOSE (UNIT = 10) ! close y_dig PRINT "(' ',8X,I6,' fiducial points were read')", y_dig_count WRITE (21,"(8X,I6,' fiducial points were read')") y_dig_count WRITE (21,"(8X,'- - - - - - - - - - - - - - - - - - - - - ')") END IF ! IF (y_dig /= 'none') PRINT "(' ',4X,'Successfully read all input datasets')" WRITE (21,"(4X,'Successfully read all input datasets')") stress_ever = (s_dat_count > 0) .OR. (faults_give_sigma_1h.AND.(f_dat_count > 0)) IF ((.NOT.stress_ever).AND.(n_refine > 0)) THEN n_refine = 0 PRINT "(' ',4X,'Lacking any stress data; refinement is useless: n_refine = 0.')" WRITE (21,"(4X,'Lacking any stress data; refinement is useless: n_refine = 0.')") END IF END IF ! (loading) ! Initialize integrated variables on 2nd and later iterations ! ( where we can assume that start_time = 0. ) IF (iteration > 1) THEN ! Note: The arrays "trace" (from f.dig) ! and "dot" (from y.dig) ! will be recomputed when we have ! re-read x.feg to rediscover original node locations. IF (c_dat_count > 0) c_end_now = c_end_0 ! whole array IF (p_dat_count > 0) p_site_now = p_site_0 ! whole array IF (s_dat_count > 0) THEN s_azim_now = s_azim_0 ! whole array DO i = 1, s_dat_count s_site_now(1:3, 1, i) = s_site_0(1:3, i) tvi = s_site_now(1:3, 1, i) CALL Step_aside(tvi, s_azim_0(i), tvo) s_site_now(1:3, 2,i) = tvo END DO END IF END IF ! (iteration > 1) ! [.feg and .bcs are read as needed within timestepping loop ] current_feg = 0 ! start at beginning of the list each iteration get_feg = .TRUE. ! better load one; if you've got one, it's deformed! IF (ALLOCATED (vw0)) vw0 = 0. ! so strain-rate is defined in Solve-for-vw IF (ALLOCATED (vw1)) vw1 = 0. ! so strain-rate is defined in Solve-for-vw IF (paleotec) THEN n_ = NINT(start_time / Deltat_) ELSE ! (neotec) n_ = 0 ! will be bumped to 1 END IF timestepping: DO IF (paleotec) THEN time0 = n_ * Deltat_ ELSE time0 = start_time END IF t0 = time0 / s_per_Ma ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ n_ = n_ + 1 ! not using DO index because we may repeat steps ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (paleotec) THEN time1 = n_ * Deltat_ t1 = time1 / s_per_Ma PRINT "(' ',4X,'Attempting timestep from ',F6.2,' to ',F6.2,' Ma')", t0, t1 WRITE (21,"(4X,'Attempting timestep from ',F6.2,' to ',F6.2,' Ma')") t0, t1 ELSE ! (neotec) time1 = time0 ! mostly for use by Write_x_vel END IF IF (get_feg) THEN get_feg = .FALSE. current_feg = current_feg + 1 ! progress through list IF (current_feg > 1) loading = .TRUE. IF (current_feg > num_fegs) THEN PRINT "(' Error: Needed but could not find another .FEG file, #',I2)", current_feg WRITE (21,"('Error: Needed but could not find another .FEG file, #',I2)") current_feg IF (current_feg == 1) THEN STOP ! (this shouldn't happen) ELSE ! there is an .feg in memory IF ((complete_timesteps == 0) .OR. (n_ == 1)) THEN PRINT "(' Provide valid .FEG and .BCS files and start from scratch.')" WRITE (21,"(' Provide valid .FEG and .BCS files and start from scratch.')") STOP ELSE PRINT "(' Writing output files necessary to restart:')" WRITE (21,"('Writing output files necessary to restart:')") mmnn = Mangle(total_iterations, time0) CALL Write_x_feg(shift = 0) IF (f_dig_count > 0) CALL Write_f_dig(shift = 0) IF (y_dig_count > 0) CALL Write_y_dig(shift = 0) IF (f_dat_count > 0) CALL Write_f_dat(shift = 0) IF (c_dat_count > 0) CALL Write_c_dat(shift = 0) IF (p_dat_count > 0) CALL Write_p_dat(shift = 0) IF (s_dat_count > 0) CALL Write_s_dat(shift = 0) PRINT "(' See further instructions in REPORT.txt.')" t1 = time0 / s_per_Ma WRITE (21,"('To resume this calculation:')") WRITE (21,"(' 1. Fix or replace deformed .FEG file.')") WRITE (21,"(' 2. Create corresponding .BCS file.')") WRITE (21,"(' 3. Correct PARAMETE[RS].DAT file:')") WRITE (21,"(' Use these, with the files just written, as input files.')") WRITE (21,"(' Set start_time to ',F7.3,' Ma in line 1.')") t1 WRITE (21,"(' Set old_iterations to ',I2,' in line 5.')") total_iterations WRITE (21,"(' Set max_iter to 1 in line 6.')") WRITE (21,"(' 4. Rerun Restore to complete this iteration.')") WRITE (21,"(' 5. Correct PARAMETE[RS].DAT file:')") WRITE (21,"(' Use present-day .DIG, .DAT files for input.')") WRITE (21,"(' Add any new .FEG, .BCS files to the corresponding lists.')") WRITE (21,"(' Set start_time to 0. in line 1.')") WRITE (21,"(' Set old_iterations to ',I2,' in line 5.')") total_iterations + 1 WRITE (21,"(' Set max_iter as high as desired in line 6.')") WRITE (21,"(' 6. Run Restore, and cross your fingers!')") STOP END IF ! (n_ > 1) END IF ELSE ! next name on list should be pre-tested as valid ! read .feg x_feg = gridname_feg(current_feg) PRINT "(' ',8X,'Reading ',A)", TRIM(x_feg) WRITE (21,"(8X,'Reading ',A)") TRIM(x_feg) ! create filename with .vel extension IF (INDEX (x_feg, '.FEG') > 0) THEN j = INDEX (x_feg, '.FEG') x_vel = x_feg(1:j) // 'VEL' ELSE IF (INDEX (x_feg, '.feg') > 0) THEN j = INDEX (x_feg, '.feg') x_vel = x_feg(1:j) // 'vel' ELSE PRINT "(' Error: No .FEG (or .feg) found in this filename'/' ',A)", TRIM(x_feg) WRITE (21,"('Error: No .FEG (or .feg) found in this filename'/A)") TRIM(x_feg) STOP END IF OPEN (UNIT = 11, FILE = x_feg, STATUS = 'OLD', ACTION = 'READ', PAD = 'YES') READ (11, *) ! skip title READ (11, *) num_nod ndof = 2 * num_nod line = 2 IF (loading) THEN ! only then do we want to thrash memory! IF (ALLOCATED (xyz_nod)) THEN DEALLOCATE (xyz_nod) DEALLOCATE (mu_nod) DEALLOCATE (vw0, vw1, vw_add, vw_mean) DEALLOCATE (u_flag) ELSE CALL More_mem ('xyz_nod', 3 * num_nod * bytes_per_real) CALL More_mem ('mu_nod', num_nod * bytes_per_real) CALL More_mem ('vw0', ndof * bytes_per_real) CALL More_mem ('vw1', ndof * bytes_per_real) CALL More_mem ('vw_add', ndof * bytes_per_real) CALL More_mem ('vw_mean', ndof * bytes_per_real) CALL More_mem ('u-flag', ndof * bytes_per_real) ENDIF ALLOCATE ( xyz_nod(3, num_nod) ) ALLOCATE ( mu_nod(num_nod) ) ALLOCATE ( vw0(ndof) ) vw0 = 0. ! so strain-rate is defined in first Solve-for-vw ALLOCATE ( vw1(ndof) ) vw1 = 0. ! so strain-rate is defined in first Solve-for-vw ALLOCATE ( vw_add(ndof) ) ALLOCATE ( vw_mean(ndof) ) ALLOCATE ( u_flag(ndof) ) END IF ! (loading) ! If there is any chance, check for nodes lying on fault traces! check_if = loading .AND. (f_dat_count > 0) IF (check_if) THEN PRINT "(' ',8X,'Checking that no node lies on any fault trace (slow!)')" WRITE (21,"(8X,'Checking that no node lies on any fault trace')") num_bad = 0 END IF DO i = 1, num_nod READ (11, *) i1, x1, x2, t; line = line + 1 IF ((i1 < 1).OR.(i1 > num_nod)) CALL Check_range('x_feg',line) CALL Xyz_from_lonlat(x1, x2, tv) xyz_nod(1:3, i1) = tv IF (check_if) THEN DO k = 1, f_highest j1 = trace_loc(1, k) j2 = trace_loc(2, k) IF (j2 > j1) THEN vec1 = trace(1:3, j1) DO j = j1 + 1, j2 vec2 = trace(1:3, j) CALL Cross(vec1, vec2, tv) CALL Unitise(tv, pole) tv = xyz_nod(1:3,i1) IF (ABS(Dot_3D(pole, tv)) < floor) THEN t1 = Arc_distance(tv, vec1) t2 = Arc_distance(tv, vec2) t3 = Arc_distance(vec1, vec2) IF ((t1 <= t3) .AND. (t2 <= t3)) THEN num_bad = num_bad + 1 PRINT "(' Error: Node ',I6,' lies on trace F',I4)", i1, k WRITE (21,"('Error: Node ',I6,' lies on trace F',I4)") i1, k END IF END IF vec1 = vec2 END DO END IF END DO END IF ! necessary to check for node on a fault trace IF (t == 0.) t = mu_ ! use default strain-rate uncertainty IF (t <= 0.) CALL Prevent ('nonpositive mu_', line, x_feg) IF (t < SQRT(1.1 * TINY(mu_))) CALL Prevent ('mu_**2 will underflow!', line, x_feg) IF (t > 1.E-10) CALL Prevent ('unreasonably large mu_', line, x_feg) mu_nod(i1) = t END DO ! i = 1, num_nod IF (check_if .AND. (num_bad > 0)) STOP IF (loading) THEN ! otherwise elements are the same READ (11, *) num_ele; line = line + 1 IF (p_dat_count > 0) THEN P_weight = SQRT(2. * num_ele) PRINT "(' ',8X,'Paleomagnetic weight factor P = ',1P,E9.2,' = SQRT(2*num_ele)')", P_weight WRITE (21,"(8X,'Paleomagnetic weight factor P = ',1P,E9.2,' = SQRT(2*num_ele)')") P_weight END IF IF (ALLOCATED (node)) THEN DEALLOCATE (node) DEALLOCATE (a_) DEALLOCATE (crack_index) IF (ALLOCATED (ele_stressed)) THEN DEALLOCATE (ele_stressed) DEALLOCATE (ele_azim) DEALLOCATE (ele_q) DEALLOCATE (ele_sigma) DEALLOCATE (ele_strainrate) END IF ELSE CALL More_mem ('node', 3* num_ele * bytes_per_int) CALL More_mem ('a_', num_ele * bytes_per_real) CALL More_mem ('crack_index', 2 * num_ele * bytes_per_int) IF (stress_ever) THEN CALL More_mem ('ele_azim', num_ele * bytes_per_real) CALL More_mem ('ele_q', num_ele * bytes_per_real) CALL More_mem ('ele_sigma', num_ele * bytes_per_real) CALL More_mem ('ele_stressed', num_ele) CALL More_mem ('boxed', num_ele) CALL More_mem ('ele_strainrate', 3 * num_ele * bytes_per_real) END IF ENDIF ALLOCATE ( node(3, num_ele) ) ALLOCATE ( a_(num_ele) ) ALLOCATE ( crack_index(2, num_ele) ) IF (stress_ever) THEN ALLOCATE ( ele_azim(num_ele) ) ALLOCATE ( ele_q(num_ele) ) ALLOCATE ( ele_sigma(num_ele) ) ALLOCATE ( ele_stressed(num_ele) ) ALLOCATE ( boxed(num_ele) ) ALLOCATE ( ele_strainrate(3, num_ele) ) END IF DO l_ = 1, num_ele READ (11, *) j, j1, j2, j3; line = line + 1 IF ((j < 1).OR.(j > num_ele).OR.(j1 < 1).OR.(j1 > num_nod) & & .OR.(j2 < 1).OR.(j2 > num_nod) & & .OR.(j3 < 1).OR.(j3 > num_nod))& & CALL Check_range('x_feg',line) node(1:3,j) = (/ j1, j2, j3 /) END DO END IF ! (loading) CLOSE (11) ! close x_feg ! read .bcs IF (loading) THEN PRINT "(' ',8X,'Reading ',A)", TRIM(gridname_bcs(current_feg)) WRITE (21,"(8X,'Reading ',A)") TRIM(gridname_bcs(current_feg)) OPEN (UNIT = 12, FILE = gridname_bcs(current_feg), STATUS = 'OLD',& & ACTION = 'READ', PAD = 'YES') bcs_count = 0 DO READ (12, *, IOSTAT = read_status) j, r1, r2 IF (read_status /= 0) EXIT bcs_count = bcs_count + 1 END DO IF (bcs_count < 2) THEN PRINT "(' Error: Provide at least 2 fixed boundary nodes in ',A)", & & TRIM(gridname_bcs(current_feg)) WRITE (21,"(' Error: Provide at least 2 fixed boundary nodes in ',A)") & & TRIM(gridname_bcs(current_feg)) STOP END IF IF (ALLOCATED (boundary_node)) THEN DEALLOCATE (boundary_node) DEALLOCATE (condition) ELSE CALL More_mem('boundary_node', bcs_count * bytes_per_int) CALL More_mem('condition', 2 * bcs_count * bytes_per_real) END IF ALLOCATE ( boundary_node(bcs_count) ) ALLOCATE ( condition(2, bcs_count) ) REWIND (12); line = 0 DO i = 1, bcs_count READ (12,*) j, r1, r2; line = line + 1 IF ((i < 1).OR.(i > num_nod)) CALL Check_range(gridname_bcs(current_feg),line) boundary_node(i) = j condition(1:2, i) = (/ r1, r2 /) END DO CLOSE (12) ! close x.bcs END IF CALL Plane_area (folding) !? IF (folding) THEN PRINT "(' ',8X,'This grid is ALREADY folded; starting timestep over')" WRITE (21,"(8X,'This grid is ALREADY folded; starting timestep over')") get_feg = .TRUE. n_ = n_ - 1 CYCLE timestepping ELSE IF (loading) THEN sum = 0. DO l_ = 1, num_ele sum = sum + a_(l_) END DO W = num_ele / sum PRINT "(' ',8X,'Continuum-stiffness weight factor W = ',1P,E9.2,' = 1/(mean element area)')", W WRITE (21,"(8X,'Continuum-stiffness weight factor W = ',1P,E9.2,' = 1/(mean element area)')") W END IF END IF !get internal coordinates for all integrated data points IF (loading) THEN IF ((f_dig_count + f_dat_count + c_dat_count + p_dat_count + s_dat_count + y_dig_count) > 0) THEN PRINT "(' ',8X,'Finding all data locations in grid coordinates')" WRITE (21,"(8X,'Finding all data locations in grid coordinates')") IF (ALLOCATED(neighbor)) THEN DEALLOCATE(neighbor) DEALLOCATE(center) ELSE CALL More_mem('neighbor', 3 * num_ele * bytes_per_int) CALL More_mem('center', 3 * num_ele * bytes_per_real) END IF ALLOCATE ( neighbor(3, num_ele) ) ALLOCATE ( center(3, num_ele) ) CALL Find_s1s2s3 ! count number of data actually in play f_in_time_and_space = 0 DO i = 1, f_dat_count IF (f_2_in(which_trace(i))) THEN any_action = .FALSE. DO j = 1, num_timesteps any_action = any_action .OR. f_active(j,i) END DO IF (any_action) f_in_time_and_space = f_in_time_and_space + 1 END IF END DO c_in_time_and_space = 0 DO i = 1, c_dat_count IF ((c_end_is(1,i)%element > 0).AND.(c_end_is(2,i)%element > 0)) THEN any_action = .FALSE. DO j = 1, num_timesteps any_action = any_action .OR. c_active(j,i) END DO IF (any_action) c_in_time_and_space = c_in_time_and_space + 1 END IF END DO p_in_time_and_space = 0 DO i = 1, p_dat_count IF (p_site_is(i)%element > 0) THEN any_action = .FALSE. DO j = 1, num_timesteps any_action = any_action .OR. p_active(j,i) END DO IF (any_action) p_in_time_and_space = p_in_time_and_space + 1 END IF END DO !Survey and then store fault segments !Note: Uses arrays "center" and "neighbor" created by Find-s1s2s3. !Also uses "trace_is" from Find-s1s2s3. IF (f_dat_count > 0) THEN PRINT "(' ',12X,'Counting fault segments')" WRITE (21,"(12X,'Counting fault segments')") CALL Def_seg (seg_count = seg_count, savem = .FALSE.) IF (seg_count > 0) THEN IF (ALLOCATED(seg_end)) THEN DEALLOCATE (seg_def) DEALLOCATE (seg_end) DEALLOCATE (seg_end_is) DEALLOCATE (seg_eta_) DEALLOCATE (seg_kappa_) DEALLOCATE (seg_u_) ELSE CALL More_mem('seg_def', 2 * seg_count * bytes_per_int) CALL More_mem('seg_end', 3 * 2 * seg_count * bytes_per_real) CALL More_mem('seg_end_is', 2 * seg_count * bytes_per_is) CALL More_mem('seg_eta_', seg_count * bytes_per_real) CALL More_mem('seg_kappa_', seg_count * bytes_per_real) CALL More_mem('seg_u_', seg_count * bytes_per_int) END IF ALLOCATE( seg_def(2, seg_count) ) ALLOCATE( seg_end(1:3, 2, seg_count) ) ALLOCATE( seg_end_is(2, seg_count) ) ALLOCATE( seg_eta_(seg_count) ) ALLOCATE( seg_kappa_(seg_count) ) ALLOCATE( seg_u_(seg_count) ) PRINT "(' ',12X,'Recording fault segments')" WRITE (21,"(12X,'Recording fault segments')") CALL Def_seg (seg_count = seg_count, savem = .TRUE.) END IF ! seg_count > 0 END IF END IF ! there is any data needing internal coordinates ! Check for paleomag. data sharing elements w/ active s-s faults. IF (p_dat_count > 0) THEN IF (seg_count > 0) THEN PRINT "(' ',8X,'Checking for paleomag data in same element as strike-slip fault(s)')" WRITE (21,"(8X,'Checking for paleomag data in same element as strike-slip fault(s)')") num_bad = 0 check_ps: DO i = 1, p_dat_count IF (time0 < p_t_max(i)) THEN ! sharing an element w/ s-s would twist datum DO j = 1, seg_count k = seg_def(1,j) ! fault trace # l_ = seg_def(2,j) ! element # IF (p_site_is(i)%element == l_) THEN DO m = 1, f_dat_count IF (which_trace(m) == k) THEN ! segment has rate info IF ((sense(m) == 'R').OR.(sense(m) == 'L')) THEN IF ((f_t_min(m) < p_t_max(i)).AND.(f_t_max(m) > time0)) THEN num_bad = num_bad + 1 twisted(i) = .TRUE. tv = p_site_0(1:3,i) CALL Lonlat_from_xyz (tv, x1, x2) WRITE (c4,"(I4)") k DO j2 = 1,3; IF (c4(j2:j2) == ' ') c4(j2:j2) = '0'; END DO PRINT "(' WARNING: Paleomag site ',A)", TRIM(p_ref(i)) PRINT "(' at ',F7.2,'E, ',F6.2,'N is twisted by F',A4,A1)", & & x1, x2, c4, sense(m) WRITE (21,"('WARNING: Paleomag site ',A)") TRIM(p_ref(i)) WRITE (21,"(' at ',F7.2,'E, ',F6.2,'N is twisted by F',A4,A1)")& & x1, x2, c4, sense(m) CYCLE check_ps END IF ! windows overlap END IF ! strike-slip END IF ! segment has a rate END DO ! m = 1, f_dat_count END IF ! segment and site are in same element END DO ! j = 1, seg_count END IF ! time0 < p_t_max(i) END DO check_ps ! i = 1, p_dat_count IF (num_bad > 0) THEN PRINT "(' ','Vertical-axis rotations at these sites will not be used.')" WRITE (21,"('Vertical-axis rotations at these sites will not be used.')") END IF ! some were twisted END IF ! seg_count > 0 END IF ! p_dat_count > 0 ! find Delta_node_feg = half-bandwidth, expressed in nodes (2 dof each) Delta_node_feg = 0 DO i = 1, num_ele i1 = node(1, i) i2 = node(2, i) i3 = node(3, i) j1 = MIN(i1, i2, i3) j3 = MAX(i1, i2, i3) Delta_node_feg = MAX(Delta_node_feg, j3 - j1) END DO ELSE ! (.NOT.loading) ! In this case, the iteration number is > 1, ! AND the first .feg was used for all timesteps. ! Therefore, the first .feg has been re-read for ! its node locations ONLY. We now have to restore ! positions back the the beginning, based on ! on internal coordinates, which are still valid. IF (f_dig_count > 0) THEN DO i = 1, f_dig_count IF (trace_is(i)%element > 0) THEN CALL Interpolate(trace_is(i), tv) trace(1:3,i) = tv END IF END DO END IF IF (y_dig_count > 0) THEN DO i = 1, y_dig_count IF (dot_is(i)%element > 0) THEN CALL Interpolate(dot_is(i), tv) dot(1:3,i) = tv END IF END DO END IF IF (seg_count > 0) THEN DO i = 1, seg_count CALL Interpolate(seg_end_is(1,i), tv) seg_end(1:3,1,i) = tv CALL Interpolate(seg_end_is(2,i), tv) seg_end(1:3,2,i) = tv END DO END IF END IF ! (loading / .NOT.loading) END IF ! next .feg on list is valid ELSE PRINT "(' ',8X,'Using old finite-element grid again')" WRITE (21,"(8X,'Using old finite-element grid again')") END IF ! (get_feg) ! Count cracks (~active fault segments), both total and by element. ! Note: A crack is not exactly an active segment; the same segment ! may be used for more than one crack when two or more offset data have ! contiguous time windows both falling into this timestep. crack_count = 0 crack_index = 0 ! whole array IF (f_dat_count > 0) THEN DO i = 1, f_dat_count IF (f_active(n_, i)) THEN ! datum is active k = which_trace(i) ! trace index j1 = trace_loc(3, k) ! 1st segment of trace j2 = trace_loc(4, k) ! last segment of trace IF (j1 > 0) THEN ! at least one segment exists DO j = j1, j2 IF ((seg_end(1,1,j) /= seg_end(1,2,j)) .OR. & & (seg_end(2,1,j) /= seg_end(2,2,j)) .OR. & & (seg_end(3,1,j) /= seg_end(3,2,j))) THEN ! ignore any zero-length segments(?) crack_count = crack_count + 1 l_ = seg_def(2, j) ! element number crack_index(1, l_) = crack_index(1, l_) + 1 END IF ! not a null segment END DO ! segments in trace END IF ! at least one segment exists END IF ! active END DO ! i = 1, f_dat_count IF (crack_count > 0) THEN ! Decide initial storage positions for local cracks in each element crack_index(2, 1) = 1 DO l_ = 2, num_ele crack_index(2, l_) = crack_index(2, l_ - 1) + crack_index(1, l_ - 1) END DO DO l_ = 1, num_ele IF (crack_index(1, l_) == 0) crack_index(2, l_) = 0 ! clear warning END DO ! allocate local crack storage IF (ALLOCATED(local_crack)) THEN DEALLOCATE (local_crack) ELSE CALL More_mem('local_crack', crack_count * bytes_per_crack) END IF ALLOCATE ( local_crack(crack_count) ) !PRINT "(' ',8X,'Recording active cracks in this timestep')" !WRITE (21,"(8X,'Recording active cracks in this timestep')") ! reset counts to 0 in each element so they can be used to keep track of open spots DO l_ = 1, num_ele crack_index(1, l_) = 0 END DO DO i = 1, f_dat_count ! offset datum index IF (f_active(n_, i)) THEN ! datum is active k = which_trace(i) ! trace index j1 = trace_loc(3, k) ! 1st segment of trace j2 = trace_loc(4, k) ! last segment of trace IF (j1 > 0) THEN ! at least one segment exists DO j = j1, j2 ! segment index IF ((seg_end(1,1,j) /= seg_end(1,2,j)) .OR. & & (seg_end(2,1,j) /= seg_end(2,2,j)) .OR. & & (seg_end(3,1,j) /= seg_end(3,2,j))) THEN ! ignore any zero-length null segments (?) l_ = seg_def(2, j) ! element number k = crack_index(2, l_) + crack_index(1, l_) !storage location crack_index(1, l_) = crack_index(1, l_) + 1 !bump up count local_crack(k)%datum = i ! to refer back when we have p_ local_crack(k)%segment = j ! to access changing length, azimuth, ! eta_, kappa_, u_, ... local_crack(k)%sense = sense(i) ! T, P, N, R, L, D IF (sense(i) == 'T') THEN factor = -cot_thrust_dip ELSE IF (sense(i) == 'P') THEN factor = -1. ELSE IF (sense(i) == 'N') THEN factor = cot_normal_dip ELSE IF (sense(i) == 'D') THEN factor = 1. ELSE IF (sense(i) == 'R') THEN factor = 1. ELSE IF (sense(i) == 'L') THEN factor = -1. ELSE ! should not happen! CALL Prevent ('illegal slip sense', i, 'array "sense"') ENDIF local_crack(k)%s_ = factor * f_goal(n_, i) local_crack(k)%sigma_ = ABS(factor) * f_rate_sigma_(i) END IF ! non-null segment END DO ! segments in trace END IF ! at least one segment exists END IF ! goal is not zero END DO ! i = 1, f_dat_count END IF ! crack_count > 0 END IF ! f_dat_count > 0 ! determine how active X-sections might have changed bandwidth, ! and reallocate linear system arrays. Delta_node = Delta_node_feg IF (c_dat_count > 0) THEN scan_sections: DO i = 1, c_dat_count j1 = c_end_is(1, i)%element j2 = c_end_is(2, i)%element IF ((j1 * j2) == 0) CYCLE scan_sections IF ((c_t_max(i) < time0) .OR. (c_t_min(i) > time1)) CYCLE scan_sections i1 = node(1, j1) i2 = node(2, j1) i3 = node(3, j1) i4 = node(1, j2) i5 = node(2, j2) i6 = node(3, j2) Delta_node = MAX(Delta_node, MAX(i1,i2,i3,i4,i5,i6) - MIN(i1,i2,i3,i4,i5,i6)) END DO scan_sections END IF ! c_dat_count > 0 IF (Delta_node /= Delta_node_last) THEN Delta_node_last = Delta_node PRINT "(' ',8X,'Delta_node is ',I4)", Delta_node WRITE (21,"(8X,'Delta_node is ',I4)") Delta_node ncoda = 2 * Delta_node + 1 lda = ndof + ncoda IF (ALLOCATED (ABCDEF)) THEN DEALLOCATE (ABCDEF) ELSE CALL More_mem ('ABCDEF matrix', lda * (ncoda + 2) * bytes_per_real) END IF ALLOCATE ( ABCDEF(lda, ncoda+2) ) IF (n_refine > 0) THEN IF (ALLOCATED (duplicate)) THEN DEALLOCATE (duplicate) ELSE CALL More_mem ('duplicate matrix', lda * (ncoda + 2) * bytes_per_real) END IF ALLOCATE ( duplicate(lda, ncoda+2) ) END IF END IF IF (paleotec) THEN PRINT "(' ',8X,'Solving for velocities at young end of timestep')" WRITE (21,"(8X,'Solving for velocities at young end of timestep')") ELSE ! neotec PRINT "(' Solving for velocities')" WRITE (21,"('Solving for velocities')") END IF ! Is there any kind of stress-direction constraint in this timestep? stress_now = .FALSE. ! (initialization only) IF (s_dat_count > 0) THEN all_stresses: DO s = 1, s_dat_count IF (s_activity(n_,s) > 0.) THEN stress_now = .TRUE. EXIT all_stresses END IF END DO all_stresses END IF ! s_dat_count > 0 IF (.NOT. stress_now) THEN IF (faults_give_sigma_1h) THEN IF (f_dat_count > 0) THEN all_faults: DO s = 1, f_dat_count IF (f_active(n_,s)) THEN stress_now = .TRUE. EXIT all_faults END IF END DO all_faults END IF END IF END IF vw0 = vw1 ! initial estimate, used in Solve-for-vw IF (stress_now) THEN boxed = .FALSE. ! Solve-for-vw can only turn it to TRUE CALL Solve_for_vw (passes = (1 + n_refine), vw = vw0) !* < * PREDICTOR * < * * < * * ELSE CALL Solve_for_vw (passes = 1, vw = vw0) !* < * PREDICTOR * < * * < * * END IF IF (paleotec) THEN ! move nodes in predictor and check for grid folding CALL Move_feg (vw = vw0) CALL Plane_area (folding) !? IF (folding) THEN PRINT "(' ',8X,'This grid would fold in predictor; starting timestep over')" WRITE (21,"(8X,'This grid would fold in predictor; starting timestep over')") get_feg = .TRUE. n_ = n_ - 1 CYCLE timestepping END IF ! move all integrated points to older positions: explicit part PRINT "(' ',8X,'Moving all points back in time: predictor')" WRITE (21,"(' ',8X,'Moving all points back in time: predictor')") CALL Move_data PRINT "(' ',8X,'Solving for velocities at older end of timestep')" WRITE (21,"(8X,'Solving for velocities at older end of timestep')") vw1 = vw0 ! initial estimate, used by Solve-for-vw CALL Solve_for_vw (passes = 1, vw = vw1) !* < * * < * CORRECTOR * < * * < * * < * * !correction velocities are half of changes DO i = 1, ndof vw_add (i) = (vw1(i) - vw0(i))/2. vw_mean(i) = (vw1(i) + vw0(i))/2. END DO ! move nodes in corrector and check for grid folding CALL Move_feg (vw = vw_add) CALL Plane_area (folding) !? IF (folding) THEN PRINT "(' ',8X,'This grid would fold in corrector; abandoning correction')" WRITE (21,"(8X,'This grid would fold in corrector; abandoning correction')") get_feg = .TRUE. vw_mean = vw0 !whole array ELSE PRINT "(' ',8X,'Adjusting for accelerations: corrector')" WRITE (21,"(' ',8X,'Adjusting for accelerations: corrector')") CALL Move_data END IF ELSE ! (neotec) vw_mean = vw0 ! whole array END IF ! Compute model predicted rates (p_) of this timestep CALL Prediction (vw = vw_mean, L0 = rate_err(0,n_,iteration), & & L1 = rate_err(1,n_,iteration), L2 = rate_err(2,n_,iteration)) ! provide complete basis for graphics if desired IF (neotec) THEN mmnn = Mangle(total_iterations, time0) mmnn(1:2) = "NT" CALL Write_x_vel(shift = 0, vw = vw0) complete_timesteps = 1 ELSE ! (paleotec) IF ((iteration == max_iter) .AND. map_set) THEN mmnn = Mangle(total_iterations, time1) CALL Write_x_feg(shift = 8) CALL Write_x_vel(shift = 8, vw = vw1) IF (f_dig_count > 0) CALL Write_f_dig(shift = 8) IF (y_dig_count > 0) CALL Write_y_dig(shift = 8) ! output positions and orientations of stress and paleomag. data, ! strictly for plotting by RetroMap: IF (p_dat_count > 0) CALL Write_p_dat(shift = 8) IF (s_dat_count > 0) CALL Write_s_dat(shift = 8) END IF t0 = time0 / s_per_Ma t1 = time1 / s_per_Ma PRINT "(' ',4X,'Completed timestep from ',F6.2,' to ',F6.2,' Ma')", t0, t1 WRITE (21,"(4X,'Completed timestep from ',F6.2,' to ',F6.2,' Ma')") t0, t1 complete_timesteps = complete_timesteps + 1 END IF ! (paleotec) IF (neotec .OR. (n_ >= num_timesteps)) EXIT timestepping END DO timestepping ! with index n_ ! An iteration of the whole history has been completed! IF (paleotec) THEN ! Describe errors for this iteration ! All rate errors (both data goals and a-priori zero-strain goals) rate_err(0, 0 , iteration) = 0. rate_err(1, 0 , iteration) = 0. rate_err(2, 0 , iteration) = 0. n_1 = NINT(start_time / Deltat_) + 1 n_t_per_it = num_timesteps - n_1 + 1 DO i = n_1, num_timesteps rate_err(0,0,iteration) = rate_err(0,0,iteration) + rate_err(0,i,iteration) rate_err(1,0,iteration) = rate_err(1,0,iteration) + rate_err(1,i,iteration) rate_err(2,0,iteration) = rate_err(2,0,iteration) + rate_err(2,i,iteration) END DO rate_err(0, 0, iteration) = rate_err(0, 0, iteration) / n_t_per_it rate_err(1, 0, iteration) = rate_err(1, 0, iteration) / n_t_per_it rate_err(2, 0, iteration) = rate_err(2, 0, iteration) / n_t_per_it PRINT "(' ',4X,'Mean rate error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)", & & rate_err(0,0,iteration),rate_err(1,0,iteration),rate_err(2,0,iteration) WRITE (21, "(4X,'Mean rate error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)") & & rate_err(0,0,iteration),rate_err(1,0,iteration),rate_err(2,0,iteration) ! Fault offset errors IF (f_in_time_and_space > 0) THEN f_err(0:2, iteration) = 0.0 ! initialize sums in_count = 0 DO i = 1, f_dat_count IF (f_2_in(which_trace(i)) .AND. (f_t_max(i) <= end_time)) THEN in_count = in_count + 1 IF (f_dat_code(i) == 'P') THEN ! Promoted type offs = f_rate(1, i) * f_t_max(i) ELSE ! Normal or Demoted types offs = 0. DO n_ = 1, num_timesteps IF (f_active(n_, i)) THEN offs = offs + f_rate(n_, i) * Deltat_ END IF END DO END IF ! Promoted or not IF (ABS((offs - offset(i))) > 2. * offset_sigma_(i)) & & f_err(0,iteration) = f_err(0,iteration) + 1. f_err(1,iteration) = f_err(1,iteration) + ABS(offs - offset(i)) / offset_sigma_(i) f_err(2,iteration) = f_err(2,iteration) + ((offs - offset(i)) / offset_sigma_(i))**2 END IF ! fault is inside time and space windows END DO ! on all fault offset data IF (in_count > 0) THEN f_err(0,iteration) = f_err(0, iteration) / in_count f_err(1,iteration) = f_err(1, iteration) / in_count f_err(2,iteration) = SQRT(f_err(2, iteration) / in_count) PRINT "(' ',4X,'Fault offset error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)", & & f_err(0,iteration),f_err(1,iteration),f_err(2,iteration) WRITE (21, "(4X,'Fault offset error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)") & & f_err(0,iteration),f_err(1,iteration),f_err(2,iteration) END IF END IF ! Cross-section extension errors IF (c_in_time_and_space > 0) THEN c_err(0:2, iteration) = 0.0 ! initialize sums in_count = 0 DO i = 1, c_dat_count IF ( (c_end_is(1,i)%element > 0) .AND. & & (c_end_is(2,i)%element > 0) .AND. & & (c_t_max(i) <= end_time)) THEN in_count = in_count + 1 stretch = 0. DO n_ = 1, num_timesteps IF (c_active(n_, i)) THEN stretch = stretch + c_rate(n_, i) * Deltat_ END IF END DO misfits = ABS(stretch - c_stretch(i)) / c_sigma_(i) IF (misfits > 2.) c_err(0,iteration) = c_err(0,iteration) + 1. c_err(1,iteration) = c_err(1,iteration) + misfits c_err(2,iteration) = c_err(2,iteration) + misfits**2 END IF ! section is inside time and space windows END DO ! on all cross-section data IF (in_count > 0) THEN c_err(0,iteration) = c_err(0, iteration) / in_count c_err(1,iteration) = c_err(1, iteration) / in_count c_err(2,iteration) = SQRT(c_err(2, iteration) / in_count) PRINT "(' ',4X,'Cross-section error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)", & & c_err(0,iteration),c_err(1,iteration),c_err(2,iteration) WRITE (21, "(4X,'Cross-section error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)") & & c_err(0,iteration),c_err(1,iteration),c_err(2,iteration) END IF END IF ! c_dat_count > 0 ! Paleolatitude errors IF (p_in_time_and_space > 0) THEN p_south_err(0:2, iteration) = 0.0 ! initialize sums in_count = 0 DO i = 1, p_dat_count IF ( (p_site_is(i)%element > 0) .AND. & & (p_t_max(i) <= end_time)) THEN in_count = in_count + 1 south = 0. DO n_ = 1, num_timesteps IF (p_active(n_, i)) & & south = south + p_south_rate(n_, i) * Deltat_ END DO misfits = ABS(south - p_south(i)) / p_south_sigma_(i) IF (misfits > 2.) p_south_err(0,iteration) = p_south_err(0,iteration) + 1. p_south_err(1,iteration) = p_south_err(1,iteration) + misfits p_south_err(2,iteration) = p_south_err(2,iteration) + misfits**2 END IF ! site is inside time and space windows END DO ! on all paleomagnetic sites IF (in_count > 0) THEN p_south_err(0,iteration) = p_south_err(0, iteration) / in_count p_south_err(1,iteration) = p_south_err(1, iteration) / in_count p_south_err(2,iteration) = SQRT(p_south_err(2, iteration) / in_count) PRINT "(' ',4X,'Paleolatitude error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)", & & p_south_err(0,iteration),p_south_err(1,iteration),p_south_err(2,iteration) WRITE (21,"(4X,'Paleolatitude error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)") & & p_south_err(0,iteration),p_south_err(1,iteration),p_south_err(2,iteration) END IF END IF ! Vertical-axis rotation errors IF (p_dat_count > 0) THEN p_ccw_err(0:2, iteration) = 0.0 ! initialize sums in_count = 0 DO i = 1, p_dat_count IF ( (p_site_is(i)%element > 0) .AND. & & (p_t_max(i) <= end_time)) THEN IF (.NOT.twisted(i)) THEN in_count = in_count + 1 ccw = 0. DO n_ = 1, num_timesteps IF (p_active(n_, i)) & & ccw = ccw + p_ccw_rate(n_, i) * Deltat_ END DO misfits = ABS(ccw - p_ccw(i)) / p_ccw_sigma_(i) IF (misfits > 2.) p_ccw_err(0,iteration) = p_ccw_err(0,iteration) + 1. p_ccw_err(1,iteration) = p_ccw_err(1,iteration) + misfits p_ccw_err(2,iteration) = p_ccw_err(2,iteration) + misfits**2 END IF ! .NOT.twisted END IF ! site is inside time and space windows END DO ! on all paleomagnetic sites IF (in_count > 0) THEN p_ccw_err(0,iteration) = p_ccw_err(0, iteration) / in_count p_ccw_err(1,iteration) = p_ccw_err(1, iteration) / in_count p_ccw_err(2,iteration) = SQRT(p_ccw_err(2, iteration) / in_count) PRINT "(' ',4X,'Vertical-axis rotation error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)", & & p_ccw_err(0,iteration),p_ccw_err(1,iteration),p_ccw_err(2,iteration) WRITE (21,"(4X,'Vertical-axis rotation error: L0 = ',F5.3,', L1 = ',F7.3,', L2 = ',F7.3)") & & p_ccw_err(0,iteration),p_ccw_err(1,iteration),p_ccw_err(2,iteration) END IF END IF PRINT "(' ','End of iteration ',I3,' out of ',I3)", iteration, max_iter WRITE (21,"('End of iteration ',I3,' out of ',I3)") iteration, max_iter END IF IF (paleotec) THEN ! Adjust goals for next iteration based on actual rates, which ! may come from iteration just finished or from restart datasets. IF ((f_dat_count + c_dat_count + p_dat_count) > 0) THEN PRINT "(' Assigning revised goals for the next iteration...')" j1 = 0; j2 = 0 ! numerator and denominator of %-adjusted END IF IF (f_dat_count > 0) THEN CALL New_goal (count = f_dat_count, total = offset, & & tmin = f_t_min, tmax = f_t_max, checkPD = .TRUE., rate = f_rate, & ! inputs & goal = f_goal, active = f_active, n_adjusted = j3) ! outputs j1 = j1 + j3 j2 = j2 + f_in_time_and_space END IF IF (c_dat_count > 0) THEN CALL New_goal (count = c_dat_count, total = c_stretch, & & tmin = c_t_min, tmax = c_t_max, checkPD = .FALSE., rate = c_rate, & ! inputs & goal = c_goal, active = c_active, n_adjusted = j3) ! outputs j1 = j1 + j3 j2 = j2 + c_in_time_and_space END IF IF (p_dat_count > 0) THEN CALL New_goal (count = p_dat_count, total = p_south, & & tmin = p_t_min, tmax = p_t_max, checkPD = .FALSE., rate = p_south_rate, & ! inputs & goal = p_south_goal, active = p_active, n_adjusted = j3) ! outputs j1 = j1 + j3 j2 = j2 + p_in_time_and_space CALL New_goal (count = p_dat_count, total = p_ccw, & & tmin = p_t_min, tmax = p_t_max, checkPD = .FALSE., rate = p_ccw_rate, & ! inputs & goal = p_ccw_goal, active = p_active, n_adjusted = j3) ! outputs j1 = j1 + j3 j2 = j2 + p_in_time_and_space END IF IF ((f_dat_count + c_dat_count + p_dat_count) > 0) THEN t = 100. * REAL(j1) / REAL(j2) PRINT "('+','Assigning revised goals for the next iteration: ',F6.2,'% adjusted')", t WRITE (21,"('Assigning revised goals for the next iteration: ',F6.2,'% adjusted')") t adjustments(iteration) = t END IF END IF ! paleotec ! output histories of strain, to permit more iterations, or scoring of model, ! or studies of model convergence. IF (watch .OR. neotec .OR. (iteration == max_iter)) THEN mmnn = Mangle(total_iterations, time1) IF (neotec) mmnn(1:2) = "NT" IF (f_dat_count > 0) CALL Write_f_dat(shift = 0) IF (c_dat_count > 0) CALL Write_c_dat(shift = 0) IF (p_dat_count > 0) CALL Write_p_dat(shift = 0) IF (paleotec .AND. (s_dat_count > 0)) CALL Write_s_dat(shift = 0) END IF END DO outer_loop ! iterations of the whole history ! Review error histories, for convergence studies IF (paleotec) THEN ! review history of goal adjustments PRINT "(' ','Review of goal adjustments in each iteration:')" WRITE (21,"('Review of goal adjustments in each iteration:')") PRINT "(' ','Iteration %-adjusted')" WRITE (21,"('Iteration %-adjusted')") DO i = 1, max_iter j = i + past_iterations PRINT "(' ',I5,F13.2,'%')", j, adjustments(i) WRITE (21,"(I5,F13.2,'%')") j, adjustments(i) END DO ! all rate errors (combines all data with a-priori zero-strains) PRINT "(' ','Review of mean rate errors in each iteration:')" WRITE (21,"('Review of mean rate errors in each iteration:')") PRINT "(' ','(these will not always converge, as rate goals shift)')" WRITE (21,"('(these will not always converge, as rate goals shift)')") PRINT "(' ','Iteration L0 L1 L2')" WRITE (21,"('Iteration L0 L1 L2')") DO i = 1, max_iter j = i + past_iterations PRINT "(' ',I5,F8.3,F7.3,F7.3)", j, rate_err(0,0,i),rate_err(1,0,i),rate_err(2,0,i) WRITE (21,"(I5,F8.3,F7.3,F7.3)") j, rate_err(0,0,i),rate_err(1,0,i),rate_err(2,0,i) END DO ! fault offset errors IF (f_dat_count > 0) THEN PRINT "(' ','Review of fault offset errors in each iteration:')" WRITE (21,"('Review of fault offset errors in each iteration:')") PRINT "(' ','Iteration L0 L1 L2')" WRITE (21,"('Iteration L0 L1 L2')") DO i = 1, max_iter j = i + past_iterations PRINT "(' ',I5,F8.3,F7.3,F7.3)", j, f_err(0,i),f_err(1,i),f_err(2,i) WRITE (21,"(I5,F8.3,F7.3,F7.3)") j, f_err(0,i),f_err(1,i),f_err(2,i) END DO END IF ! cross-section errors IF (c_dat_count > 0) THEN PRINT "(' ','Review of cross-section errors in each iteration:')" WRITE (21,"('Review of cross-section errors in each iteration:')") PRINT "(' ','Iteration L0 L1 L2')" WRITE (21,"('Iteration L0 L1 L2')") DO i = 1, max_iter j = i + past_iterations PRINT "(' ',I5,F8.3,F7.3,F7.3)", j, c_err(0,i),c_err(1,i),c_err(2,i) WRITE (21,"(I5,F8.3,F7.3,F7.3)") j, c_err(0,i),c_err(1,i),c_err(2,i) END DO END IF ! paleolatitude errors IF (p_dat_count > 0) THEN PRINT "(' ','Review of paleolatitude errors in each iteration:')" WRITE (21,"('Review of paleolatitude errors in each iteration:')") PRINT "(' ','Iteration L0 L1 L2')" WRITE (21,"('Iteration L0 L1 L2')") DO i = 1, max_iter j = i + past_iterations PRINT "(' ',I5,F8.3,F7.3,F7.3)", j, p_south_err(0,i),p_south_err(1,i),p_south_err(2,i) WRITE (21,"(I5,F8.3,F7.3,F7.3)") j, p_south_err(0,i),p_south_err(1,i),p_south_err(2,i) END DO ! vertical-axis rotation errors PRINT "(' ','Review of vertical-axis rotation errors in each iteration:')" WRITE (21,"('Review of vertical-axis rotation errors in each iteration:')") PRINT "(' ','Iteration L0 L1 L2')" WRITE (21,"('Iteration L0 L1 L2')") DO i = 1, max_iter j = i + past_iterations PRINT "(' ',I5,F8.3,F7.3,F7.3)", j, p_ccw_err(0,i),p_ccw_err(1,i),p_ccw_err(2,i) WRITE (21,"(I5,F8.3,F7.3,F7.3)") j, p_ccw_err(0,i),p_ccw_err(1,i),p_ccw_err(2,i) END DO END IF ! p_dat_count > 0 END IF ! paleotec; need iteration history ! final time stamp, and close CALL DATE_AND_TIME (date, clock_time, zone, datetimenumber) WRITE (21,"('Run ended on ',I4,'.',I2,'.',I2,' at ',I2,':',I2,':',I2)") & datetimenumber(1), datetimenumber(2), datetimenumber(3), & datetimenumber(5), datetimenumber(6), datetimenumber(7) WRITE (21, "('===========================================================')") CLOSE (UNIT = 21) ! close REPORT.txt PRINT "(' Successful termination of Restore; see REPORT.txt.')" PRINT "(' ')" !============================================================================== CONTAINS ! Any code that would otherwise be typed repeatedly SUBROUTINE Add_datum(prefix,f_,g_,goal,A,C,D,E,F) ! Adds a datum to the linear system of one element. ! Note that only the diagonal and lower triangle are filled; ! the upper triangle is symmetric. REAL, INTENT(IN) :: prefix, goal REAL, DIMENSION(3), INTENT(IN) :: f_,g_ REAL, DIMENSION(3,3), INTENT(INOUT) :: A,C,D REAL, DIMENSION(3), INTENT(INOUT) :: E,F INTEGER i, j DO i = 1, 3 DO j = 1, 3 IF (j <= i) THEN ! diagonal and lower triangle only A(i,j) = A(i,j) + prefix * f_(i) * f_(j) D(i,j) = D(i,j) + prefix * g_(i) * g_(j) END IF ! All of B lies in the upper triangle; transpose of C ! B(i,j) = B(i,j) + prefix * f_(i) *g_(j) ! All of C lies in lower triangle C(i,j) = C(i,j) + prefix * g_(i) * f_(j) END DO ! on j = 1, 6 E(i) = E(i) + prefix * f_(i) * goal F(i) = F(i) + prefix * g_(i) * goal END DO ! on i = 1, 6 END SUBROUTINE Add_datum REAL FUNCTION Arc_distance (a, b) ! distance in radians along a great circle from a to b, ! which are both Cartesian unit vectors REAL, DIMENSION(3), INTENT(IN) :: a, b REAL, DIMENSION(3) :: tv REAL :: cosa, sina cosa = Dot_3D( a, b ) CALL Cross(a, b, tv) sina = Magnitude ( tv ) Arc_distance = ATAN2(sina, cosa) END FUNCTION Arc_distance REAL FUNCTION ATAN2F(y, x) ! works like ATAN2 but corrects for case of (0, 0). ! returns inverse tangent in radians. REAL, INTENT(IN) :: y, x IF (y == 0.) THEN IF (x == 0.) THEN ATAN2F = 0. ELSE ATAN2F = ATAN2(y, x) END IF ELSE ATAN2F = ATAN2(y, x) END IF END FUNCTION ATAN2F SUBROUTINE Check_range (filename, line) CHARACTER(*):: filename INTEGER :: line PRINT "(' Error: Integer out of range in line ',I5,' of '/' ',A)", line, TRIM(filename) WRITE (21,"('Error: Integer out of range in line ',I5,' of '/A)") line, TRIM(filename) STOP END SUBROUTINE Check_range SUBROUTINE Components(l_,G,vw,v_,w_) ! Given element l_ and local nodal functions G, ! computes S and E components of velocity (v_,w_) ! from the long-vector form vw. INTEGER, INTENT(IN) :: l_ REAL, DIMENSION(3,2,2), INTENT(IN) :: G REAL, DIMENSION(:), INTENT(IN) :: vw REAL, INTENT(OUT) :: v_, w_ INTEGER :: iv, iw, j v_ = 0. w_ = 0. DO j = 1, 3 iv = 2 * node(j, l_) - 1 iw = iv + 1 v_ = v_ + G(j,1,1) * vw(iv) + G(j,2,1) * vw(iw) w_ = w_ + G(j,1,2) * vw(iv) + G(j,2,2) * vw(iw) END DO END SUBROUTINE Components SUBROUTINE Cross (a, b, c) REAL, DIMENSION(3) :: a, b, c ! vector cross product: c = a x b c(1) = a(2)*b(3) - a(3)*b(2) c(2) = a(3)*b(1) - a(1)*b(3) c(3) = a(1)*b(2) - a(2)*b(1) END SUBROUTINE Cross SUBROUTINE Def_seg (seg_count, savem) ! Defines fault segments !(a fault segment is the intersection of a trace with an element) ! IF (savem) THEN segments are recorded; otherwise just counted. !NOTE: Logic does not consider possible segments until one ! digitised point is actually in an element! That is, if ! a fault trace is defined by two points very far apart, ! spanning the model domain, then it could be missed. ! However, once first segment on fault is found, we try hard ! not to miss any small internal ones; it is OK if there is ! not a digitised point in each element crossed! INTEGER, INTENT(OUT) :: seg_count LOGICAL, INTENT(IN) :: savem LOGICAL :: no_cut, punt REAL :: s1, s2, s3, smin INTEGER :: element, i, j, jin, jout, j1, j2, k, kside, lastele, limit, newele REAL, DIMENSION(3) :: tv, tv1, tv2, vector, last_vec TYPE(is123) :: ist, lastis LOGICAL cutin, cutout ! does segment touch element sides? LOGICAL busy ! indicates segment has been opened already !decide on a plausibility limit for the number of segments limit = 0 DO i = 1, f_highest IF (trace_loc(1, i) > 0) limit = limit + 1 END DO limit = 2 * limit * num_ele punt = .FALSE. seg_count = 0 DO i = 1, f_highest trace_loc(3, i) = 0 ! in case no segments trace_loc(4, i) = 0 ! are found in grid j1 = trace_loc(1, i) j2 = trace_loc(2, i) IF ((j1 > 0) .AND. (j2 > j1)) THEN busy = .FALSE. DO j = j1, j2 ! scan along trace i 1000 element = trace_is(j)%element IF (busy) THEN ! looking for end of segment IF (element == lastele) THEN ! still in same element IF (j == j2) THEN ! end of the line ! busy = .FALSE. will be reset above IF (savem) THEN cutout = .FALSE. trace_loc(4, i) = seg_count seg_end(1:3, 2, seg_count) = trace(1:3, j) seg_end_is(2, seg_count) = trace_is(j) CALL Finish_seg (seg_count, cutin, jin, cutout, jout, no_cut) IF (no_cut) THEN seg_count = seg_count - 1 IF (trace_loc(4,i) > trace_loc(3,i)) THEN trace_loc(4,i) = trace_loc(4,i) - 1 ELSE trace_loc(3,i) = 0 trace_loc(4,i) = 0 END IF END IF END IF ! savem ELSE last_vec = trace(1:3, j) lastis = trace_is(j) END IF ! end/not end of the line ELSE ! we have crossed out of lastele busy = .FALSE. tv = trace(1:3, j) CALL Find_out (element = lastele, & & inside = last_vec, & & outside = tv, & & border = vector, & & coords = ist) IF (savem) THEN cutout = .TRUE. seg_end(1:3, 2, seg_count) = vector seg_end_is(2, seg_count) = ist jout = Which_zero(ist) CALL Finish_seg (seg_count, cutin, jin, cutout, jout, no_cut) IF (no_cut) THEN seg_count = seg_count - 1 IF (trace_loc(4,i) > trace_loc(3,i)) THEN trace_loc(4,i) = trace_loc(4,i) - 1 ELSE trace_loc(3,i) = 0 trace_loc(4,i) = 0 END IF END IF END IF ! savem ! Now, have we crossed into another ele? smin = 999. DO k = 1, 3 IF (ABS(ist%s(k)) < smin) THEN kside = k smin = ABS(ist%s(k)) END IF END DO newele = neighbor(kside, lastele) IF (newele == 0) THEN ! have left domain; wait for re-entry(?) ELSE ! still in domain; start new element CALL Dumb_s123 (newele, vector, s1, s2, s3) ist%element = newele ist%s(1) = s1 ist%s(2) = s2 ist%s(3) = s3 busy = .TRUE. seg_count = seg_count + 1 IF (seg_count > limit) CALL Dump_seg(limit, savem, punt) IF (punt) RETURN IF (savem) THEN trace_loc(4, i) = trace_loc(4, i) + 1 cutin = .TRUE. seg_def(1, seg_count) = i ! trace index seg_def(2, seg_count) = newele seg_end(1:3, 1, seg_count) = vector seg_end_is(1, seg_count) = ist jin = Which_zero(ist) END IF ! savem lastele = newele last_vec = vector lastis = ist GOTO 1000 END IF ! still in the domain; new element END IF ! we have crossed out of lastele ELSE ! looking for beginning (initial point, or entry) IF (element > 0) THEN !first point has been found busy = .TRUE. seg_count = seg_count + 1 IF (seg_count > limit) CALL Dump_seg(limit, savem, punt) IF (punt) RETURN lastele = element last_vec = trace(1:3, j) lastis = trace_is(j) IF (savem) THEN trace_loc(3, i) = seg_count ! 1st seg for trace trace_loc(4, i) = seg_count ! last = 1st seg_def(1, seg_count) = i ! trace index seg_def(2, seg_count) = element END IF ! savem IF (j == j1) THEN ! 1st point already inside IF (savem) THEN cutin = .FALSE. seg_end(1:3, 1, seg_count) = last_vec seg_end_is(1, seg_count) = lastis END IF ! savem ELSE ! 1st point is border crossing tv1 = trace(1:3, j) tv2 = trace(1:3, j - 1) CALL Find_out (element = element, & & inside = tv1, & & outside = tv2, & & border = vector, & & coords = ist) IF (savem) THEN cutin = .TRUE. seg_end(1:3, 1, seg_count) = vector seg_end_is(1, seg_count) = ist jin = Which_zero(ist) END IF ! savem END IF ! 1st point is j1 / border crossing IF (j == j2) THEN ! this is also the end ! busy = .FALSE. will be reset IF (savem) THEN cutout = .FALSE. trace_loc(4, i) = seg_count seg_end(1:3, 2, seg_count) = trace(1:3, j) seg_end_is(2, seg_count) = trace_is(j) Call Finish_seg (seg_count, cutin, jin, cutout, jout, no_cut) IF (no_cut) THEN seg_count = seg_count - 1 IF (trace_loc(4,i) > trace_loc(3,i)) THEN trace_loc(4,i) = trace_loc(4,i) - 1 ELSE trace_loc(3,i) = 0 trace_loc(4,i) = 0 END IF END IF END IF ! savem END IF ! 1st digitised point is also the end END IF ! beginning found END IF ! busy / not busy END DO ! on j = j1, j2; points of trace trace_loc(4, i) = seg_count END IF ! this trace is in memory END DO ! on i = 1, f_highest (all possible trace indeces) END SUBROUTINE Def_seg SUBROUTINE Del_Gjxy_del_thetaphi (l_, r_, dG) INTEGER, INTENT(IN) :: l_ ! element number REAL, DIMENSION(3), INTENT(IN) :: r_ ! position vector REAL, DIMENSION (3,2,2,2), INTENT(OUT) :: dG ! computes array of 2 derivitives of each of the 2 components of ! each of the 6 nodal functions for element l_ at ! position r_ (Cartesian unit vector). ! Results are in 1./radian (dimensionless), NOT 1./m or 1./degree ! ! It is user's responsibility that element l_ contains r_. SAVE ! allows fast re-entry when l_ is unchanged. INTEGER :: l_last = 0 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2)? INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function? INTEGER :: m ! 1:2 = theta (S-ward) or phi (E-ward) derivitive? REAL, DIMENSION(3,2) :: del_r_ ! theta- and phi-derivitives of r_ (in 3-D) REAL, DIMENSION(3,2) :: local ! local Theta, Phi unit vectors at r_ (xyz, SE) REAL, DIMENSION(3,2,2) :: del_local ! theta-, phi- derivitives of local REAL, DIMENSION(3,3) :: corner ! positions vector of corner nodes (xyz, 123) REAL, DIMENSION(3,3,2) :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, tv3, vfa, vfb ! temporary vector factors REAL :: cos_phi, cos_theta, phi, sin_phi, sin_theta INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ DO j = 1, 3 corner(1:3, j) = xyz_nod(1:3, node(j, l_)) tvi = corner(1:3, j) CALL Local_Theta(tvi, tvo) post(1:3, j, 1) = tvo CALL Local_Phi (tvi, tvo) post(1:3, j, 2) = tvo END DO END IF ! begin calculations which depend on r_ CALL Local_Theta(r_, tv) local(1:3,1) = tv CALL Local_Phi(r_, tv) local(1:3,2) = tv ! Note: these functions will catch polar points; don't test again phi = ATAN2(r_(2), r_(1)) cos_phi = COS(phi) sin_phi = SIN(phi) cos_theta = r_(3) sin_theta = SQRT(r_(1)**2 + r_(2)**2) del_r_(1:3,1) = local(1:3,1) ! d.r_/d.theta = Theta del_r_(1:3,2) = local(1:3,2) * sin_theta ! d.r_/d.phi = Phi * SIN(theta) del_local(1:3,1,1) = - r_(1:3) ! d.Theta/d.theta = - r_ del_local(1:3,1,2) = local(1:3,2) * cos_theta ! d.Theta/d.phi = Phi * COS(theta) del_local(1:3,2,1) = (/ 0., 0., 0. /) ! d.Phi/d.theta = 0 del_local(1:3,2,2) = (/ -cos_phi, -sin_phi, 0. /) ! d.Phi/d.phi = (-COS(phi),-SIN(phi,0) DO j = 1, 3 ! 3 corner nodes of element i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1 = corner(1:3, i1) tv2 = corner(1:3, i2) tv3 = corner(1:3, i3) CALL Cross(tv2, tv3, vfa) vfb = vfa / Dot_3D (tv1, vfa) DO x = 1, 2 ! unit velocity at node is S or E DO y = 1, 2 ! S- or E- component of nodal function tv1 = post(1:3, j, x) tvi = local(1:3, y) DO m = 1, 2 ! theta- or phi-derivitive tv = del_r_(1:3, m) tvo = del_local(1:3, y, m) dG(j, x, y, m) = & & (Dot_3D(tv,vfb)*Dot_3D(tv1,tvi)) + & & (Dot_3D(r_,vfb)*Dot_3D(tv1,tvo)) END DO END DO END DO END DO END SUBROUTINE Del_Gjxy_del_thetaphi REAL FUNCTION Dot_3D (a_, b_) ! Dot product of 3-component vectors. REAL, DIMENSION(3) :: a_, b_ Dot_3D = a_(1)*b_(1) + a_(2)*b_(2) + a_(3)*b_(3) END FUNCTION Dot_3D SUBROUTINE Dumb_s123 (element, vector, s1, s2, s3) ! Finds s1, s2, s3 coordinates of position vector "in" element ! (whether the point is actually in the element or NOT). INTEGER, INTENT(IN) :: element REAL, DIMENSION(3), INTENT(IN) :: vector REAL, INTENT(OUT) :: s1, s2, s3 INTEGER :: i1, i2, i3 REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL :: d1, dc, t IF (element == 0) THEN CALL Prevent('element = 0', 1, 'Dumb_s123') END IF i1 = node(1, element) i2 = node(2, element) i3 = node(3, element) !shorten(?) vector to just touch plane element -> v1 tv1 = center(1:3, element) dc = Dot_3D(vector, tv1) IF (dc <= 0.) CALL Prevent('"Internal" vector >= 90 deg. from element', 1, 'Dumb_s123') tv2 = xyz_nod(1:3, i1) d1 = Dot_3D(tv2, tv1) t = d1 / dc v1 = t * vector tvi = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) tvo = v1(1:3) - xyz_nod(1:3,i3) CALL Cross(tvi, tvo, tv) s1 = Dot_3D(tv1, tv) * half_R2 / a_(element) tvi = xyz_nod(1:3,i1) - xyz_nod(1:3,i3) tvo = v1(1:3) - xyz_nod(1:3,i1) CALL Cross(tvi, tvo, tv) s2 = Dot_3D(tv1, tv) * half_R2 / a_(element) s3 = 1.00 - s1 - s2 END SUBROUTINE Dumb_s123 SUBROUTINE Dump_seg(limit, savem, punt) ! Called for debugging information: dump of segments list, ! when number becomes excessive. INTEGER, INTENT(IN) :: limit LOGICAL, INTENT(IN) :: savem LOGICAL, INTENT(OUT) :: punt INTEGER :: i IF (savem) THEN PRINT "(' Error: Number of fault segments exceeds plausible limit of ',I8)", limit PRINT "(' Probable infinite loop in Def_seg.')" PRINT "(' See file REPORT.txt for dump of segments.')" WRITE (21,"('Error: Number of fault segments exceeds plausible limit of ',I8)") limit WRITE (21,"('(2 * number of elements * number of traces)')") WRITE (21,"('Probable infinite loop in Def_seg.')") WRITE (21,"('Dump of segments follows')") WRITE (21,"(/'Trace Element I1 S1 S2 S3 I2 S1 S2 S3')") DO i = 1, limit WRITE (21, "(I5,I7,I5,3F6.4,I5,3F6.4)") & & seg_def(1, i), seg_def(2,i), & & seg_end_is(1,i)%element, seg_end_is(1,i)%s(1:3), & & seg_end_is(2,i)%element, seg_end_is(2,i)%s(1:3) END DO STOP ELSE punt = .TRUE. PRINT "(' Aborting search for fault segments at count of ',I8)", limit ENDIF END SUBROUTINE Dump_seg SUBROUTINE E_rate(l_, G, dG, theta_, vw, eps_dot) ! evaluate strain-rate INTEGER, INTENT(IN) :: l_ ! element number REAL, DIMENSION(3,2,2) :: G ! nodal functions @ selected point REAL, DIMENSION(3,2,2,2):: dG ! derivitives of nodal functions REAL, INTENT(IN) :: theta_ ! colatitude, radians REAL, DIMENSION(:), INTENT(IN) :: vw REAL, DIMENSION(3), INTENT(OUT) :: eps_dot INTEGER :: iv, iw, j REAL :: cott, csct, prefix eps_dot = 0. ! (1..3) cott = 1. / TAN(theta_) csct = 1. / SIN(theta_) prefix = 1./R ! global variable; planet radius DO j = 1, 3 iv = 2 * node(j, l_) - 1 ! global index array iw = iv + 1 ! epsilon_dot_sub_theta_theta eps_dot(1) = eps_dot(1) + & & vw(iv) * prefix * dG(j,1,1,1) + & & vw(iw) * prefix * dG(j,2,1,1) ! epsilon_dot_sub_theta_phi eps_dot(2) = eps_dot(2) + & & vw(iv) * prefix * 0.5 * (csct * dG(j,1,1,2) + dG(j,1,2,1) - cott * G(j,1,2)) + & & vw(iw) * prefix * 0.5 * (csct * dG(j,2,1,2) + dG(j,2,2,1) - cott * G(j,2,2)) ! epsilon_dot_sub_phi_phi eps_dot(3) = eps_dot(3) + & & vw(iv) * prefix * (csct * dG(j,1,2,2) + cott * G(j,1,1)) + & & vw(iw) * prefix * (csct * dG(j,2,2,2) + cott * G(j,2,1)) END DO ! 3 local nodes END SUBROUTINE E_rate SUBROUTINE Find_out (element, inside, outside, border, coords) ! Given a directed arc of a great circle from point "inside" ! (which is inside, or on the border of "element") to point ! "outside" (both positions Cartesian unit vectors), finds ! the last point in contact with "element" and reports its ! Cartesian vector as "border" and its internal coordinates as "coords". INTEGER, INTENT(IN) :: element REAL, DIMENSION(3), INTENT(IN) :: inside, outside REAL, DIMENSION(3), INTENT(OUT):: border TYPE(is123), INTENT(OUT) :: coords INTEGER :: i, side, sidea, sideb INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL, DIMENSION(3) :: s, far, frac ! NOT vectors REAL :: distance, slope coords%element = element CALL Dumb_s123 (element, inside, s(1), s(2), s(3)) CALL Pull_in(s) CALL Dumb_s123 (element, outside, far(1), far(2), far(3)) DO i = 1, 3 slope = far(i) - s(i) IF (slope < 0.) THEN frac(i) = -s(i) / slope ELSE frac(i) = 999. ENDIF END DO distance = MINVAL(frac) array = MINLOC(frac) side = array(1) sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(side) = 0. s(sidea) = s(sidea) + distance * (far(sidea) - s(sidea)) s(sideb) = s(sideb) + distance * (far(sideb) - s(sideb)) CALL Pull_in(s) coords%s(1:3) = s CALL Interpolate(coords, border) END SUBROUTINE Find_out SUBROUTINE Find_s1s2s3 !determines internal coordinates of all integrated positions; ! for convenience, also turns off c_active, p_active if outside grid. REAL, DIMENSION(3) :: tv, v1 INTEGER :: a, b, back1, back2, element, lastel INTEGER :: i, j, jold, j1, j2, jt, k, l_, m, n REAL :: s1, s2, s3 CHARACTER(61) :: bar_graph LOGICAL :: debug = .FALSE. ! find element center points DO l_ = 1, num_ele v1 = xyz_nod(1:3,node(1,l_)) + xyz_nod(1:3,node(2,l_)) + xyz_nod(1:3,node(3,l_)) CALL Unitise(v1, tv) center(1:3, l_) = tv END DO ! find neighboring elements on all 3 sides of each neighbor = 0 homes: DO l_ = 1, num_ele sides: DO j = 1, 3 ! 3 sides k = 1 + MOD (j, 3) a = node(k, l_) ! 1st node along side b = node(1 + MOD (k, 3), l_) ! 2nd node along side strangers: DO m = 1, num_ele IF ((node(1, m) == b) .AND. (node(2, m) == a)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF ((node(2, m) == b) .AND. (node(3, m) == a)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF ((node(3, m) == b) .AND. (node(1, m) == a)) THEN neighbor(j, l_) = m EXIT strangers END IF END DO strangers END DO sides END DO homes IF (f_dig_count > 0) THEN bar_graph(1:41) = ' Locating fault traces (slow) ' IF (debug) WRITE (21, "('trace_is before correction:')") IF (debug) WRITE (21, "(' point element s1 s2 s3')") DO i = 42, 61 bar_graph(i:i) = CHAR(176) END DO PRINT "(' ',A)", bar_graph PRINT "('+',A)", bar_graph(1:41) jold = 0 DO i = 1, f_dig_count tv = trace(1:3, i) CALL Internal(tv, l_, s1, s2, s3) trace_is(i)%element = l_ trace_is(i)%s(1) = s1 trace_is(i)%s(2) = s2 trace_is(i)%s(3) = s3 IF (debug) WRITE (21, "(I6,I8,3F6.2)") i, trace_is(i)%element, & & (trace_is(i)%s(j3),j3=1,3) j = (20 * i) / f_dig_count IF (j > jold) THEN bar_graph(j+41:j+41) = CHAR(219) PRINT "('+',A)", bar_graph(1:j+41) jold = j END IF END DO ! Scan for cases where a trace wanders from element a briefly ! into element b and then back into a. Change element assignments ! of b -> a. This does not move the trace, but is does prevent ! serious problems with fault segmentation later. IF (debug) WRITE (21, "('Corrections to trace_is:')") IF (debug) WRITE (21, "('trace point old_ele old_s1 old_s2 old_s3 new_ele new_s1 new_s2 new_s3')") DO i = 1, f_highest j1 = trace_loc(1, i) j2 = trace_loc(2, i) IF ((j1 > 0) .AND. (j2 > j1+1)) THEN lastel = -1 ! deliberately invalid initial values back1 = -1 back2 = -1 DO j = j1, j2 element = trace_is(j)%element IF (element /= lastel) THEN back2 = back1 ! cascade back1 = lastel IF ((element == back2).AND.(element > 0)) THEN ! Problem ! jt = j - 1 fix_back: DO IF (trace_is(jt)%element == back1) THEN trace_is(jt)%element = back2 IF (back2 > 0) THEN tv = trace(1:3, jt) CALL Dumb_s123 (element = back2, vector = tv, s1 = s1, s2 = s2, s3 = s3) IF (debug) WRITE (21, 9921) i, jt, back1, trace_is(jt)%s(1), trace_is(jt)%s(2), & & trace_is(jt)%s(3), back2, s1, s2, s3 9921 FORMAT (I5,I6,I8,3F7.2,I8,3F7.2) trace_is(jt)%s(1) = s1 trace_is(jt)%s(2) = s2 trace_is(jt)%s(3) = s3 ELSE trace_is(jt)%s = (/ 0.0, 0.0, 0.0 /) END IF ELSE EXIT fix_back END IF jt = jt - 1 END DO fix_back END IF END IF ! prepare to loop lastel = element END DO END IF END DO ENDIF IF (f_highest > 0) THEN DO i = 1, f_highest a = trace_loc(1, i) b = trace_loc(2, i) IF ((a > 0) .AND. (b > a)) THEN n = 0 ! count points of trace inside the grid DO j = a, b IF (trace_is(j)%element > 0) n = n + 1 END DO IF (n < 2) f_2_in(i) = .FALSE. ELSE f_2_in(i) = .FALSE. END IF END DO END IF IF (c_dat_count > 0) THEN PRINT "(' ',12X,'Locating cross-section ends')" DO i = 1, c_dat_count tv = c_end_now(1:3, 1, i) CALL Internal (tv, l_, s1, s2, s3) c_end_is(1, i)%element = l_ c_end_is(1, i)%s(1) = s1 c_end_is(1, i)%s(2) = s2 c_end_is(1, i)%s(3) = s3 tv = c_end_now(1:3, 2, i) CALL Internal (tv, l_, s1, s2, s3) c_end_is(2, i)%element = l_ c_end_is(2, i)%s(1) = s1 c_end_is(2, i)%s(2) = s2 c_end_is(2, i)%s(3) = s3 IF ((l_ == 0) .OR. (c_end_is(1, i)%element == 0)) THEN DO j = 1, num_timesteps c_active(j, i) = .FALSE. END DO END IF END DO ENDIF IF (p_dat_count > 0) THEN PRINT "(' ',12X,'Locating paleomagnetic sites')" DO i = 1, p_dat_count tv = p_site_now(1:3, i) CALL Internal (tv, l_, s1, s2, s3) p_site_is(i)%element = l_ p_site_is(i)%s(1) = s1 p_site_is(i)%s(2) = s2 p_site_is(i)%s(3) = s3 IF (l_ == 0) THEN DO j = 1, num_timesteps p_active(j, i) = .FALSE. END DO END IF END DO ENDIF IF (s_dat_count > 0) THEN PRINT "(' ',12X,'Locating paleostress sites')" DO i = 1, s_dat_count tv = s_site_now(1:3, 1, i) CALL Internal (tv, l_, s1, s2, s3) s_site_is(1, i)%element = l_ s_site_is(1, i)%s(1) = s1 s_site_is(1, i)%s(2) = s2 s_site_is(1, i)%s(3) = s3 tv = s_site_now(1:3, 2, i) CALL Internal (tv, l_, s1, s2, s3) s_site_is(2, i)%element = l_ s_site_is(2, i)%s(1) = s1 s_site_is(2, i)%s(2) = s2 s_site_is(2, i)%s(3) = s3 END DO ENDIF IF (y_dig_count > 0) THEN bar_graph(1:44) = ' Locating fiducial points (slow) ' DO i = 45, 61 bar_graph(i:i) = CHAR(176) END DO PRINT "(' ',A)", bar_graph PRINT "('+',A)", bar_graph(1:44) jold = 0 DO i = 1, y_dig_count tv = dot(1:3, i) CALL Internal (tv, l_, s1, s2, s3) dot_is(i)%element = l_ dot_is(i)%s(1) = s1 dot_is(i)%s(2) = s2 dot_is(i)%s(3) = s3 j = (17 * i) / y_dig_count IF (j > jold) THEN bar_graph(j+44:j+44) = CHAR(219) PRINT "('+',A)", bar_graph(1:j+44) jold = j END IF END DO ENDIF END SUBROUTINE Find_s1s2s3 SUBROUTINE Finish_seg (segment, cutin, jin, cutout, jout, no_cut) ! completes segment description by deciding u_, eta_, kappa_ ! IF (cutin), uses jin; else, computes it. ! IF (cutout), uses jout; else, computes it. ! IF (jin == jout) sets flag no_cut = T; segment should be junked. ! Refers to global arrays seg_... LOGICAL, INTENT(IN) :: cutin, cutout INTEGER, INTENT(INOUT) :: jin, jout INTEGER, INTENT(IN) :: segment LOGICAL, INTENT(OUT) :: no_cut INTEGER :: element, n REAL, DIMENSION(3) :: c, in, out, tv, tv1, tv2, vn, vs TYPE(is123) :: ist LOGICAL :: debug = .FALSE. ! IF (debug), table of segments goes to REPORT.txt. REAL :: lat1, lat2, lon1, lon2 element = seg_def(2, segment) IF (cutin) THEN in = seg_end(1:3, 1, segment) ELSE ! project backward tv2 = seg_end(1:3, 2, segment) tv1 = seg_end(1:3, 1, segment) CALL Find_out (element = element, & & inside = tv2, & & outside = tv1, & & border = in, & & coords = ist) jin = Which_zero(ist) END IF IF (cutout) THEN out = seg_end(1:3, 2, segment) ELSE ! project forward tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) CALL Find_out (element = element, & & inside = tv1, & & outside = tv2, & & border = out, & & coords = ist) jout = Which_zero(ist) END IF IF (jin == jout) THEN no_cut = .TRUE. ELSE no_cut = .FALSE. n = 6 - jin - jout ! using local node numbering 1-3 IF ((n >=1) .AND. (n <= 3)) THEN seg_u_(segment) = n ELSE PRINT "(' Error: In Finish_seg, jin =',I2,', jout =',I2)", jin, jout PRINT "(' cutin = ',L1,', cutout = ',L1)", cutin, cutout PRINT "(' element =',I5,', segment =',I6)", element, segment WRITE (21,"(' Error: In Finish_seg, jin =',I2,', jout =',I2)") jin, jout WRITE (21,"(' cutin = ',L1,', cutout = ',L1)") cutin, cutout WRITE (21,"(' element =',I5,', segment =',I6)") element, segment STOP ENDIF IF (cutin .AND. cutout) THEN seg_kappa_(segment) = 1. ELSE tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) seg_kappa_(segment) = Arc_distance(tv1, tv2) / & & Arc_distance(in, out) ENDIF n = node(n, element) ! changing n to absolute node number vn = xyz_nod(1:3, n) - seg_end(1:3, 1, segment) vs = seg_end(1:3, 2, segment) - seg_end(1:3, 1, segment) CALL Cross (vn, vs, c) tv = center(1:3, element) IF (Dot_3D(c, tv) > 0.) THEN seg_eta_(segment) = +1. ELSE seg_eta_(segment) = -1. END IF END IF IF (debug) THEN ! A table is printed in file REPORT.txt, listing all segments. tv = seg_end(1:3, 1, segment) CALL Lonlat_from_xyz (tv, lon1, lat1) tv = seg_end(1:3, 2, segment) CALL Lonlat_from_xyz (tv, lon2, lat2) WRITE (21, & &"(' seg& & def(1) def(2)& & tr_lo(3)& & tr_lo(4)& & seg_end_is(1)& & seg_end_is(2)& & seg_end(1)& & seg_end(2)& & cin jin cout jout no_cut& & u_ eta_ kappa_')") WRITE (21, 9921) & &segment,& &seg_def(1,segment),seg_def(2,segment),& &trace_loc(3,seg_def(1,segment)),& &trace_loc(4,seg_def(1,segment)),& &seg_end_is(1,segment)%element,(seg_end_is(1,segment)%s(j3),j3=1,3),& &seg_end_is(2,segment)%element,(seg_end_is(2,segment)%s(j3),j3=1,3),& &lon1,lat1,& &lon2,lat2,& &cutin,jin,cutout,jout,no_cut,& &seg_u_(segment),seg_eta_(segment),seg_kappa_(segment) 9921 FORMAT (I4,& &I7,I7,& &I9,& &I9,& &I4,3F5.2,& &I4,3F5.2,& &F7.2,F6.2,& &F7.2,F6.2,& &L4,I4,L5,I5,L7,& &I4,F5.0,F7.2) END IF ! debug END SUBROUTINE Finish_seg REAL FUNCTION Get_azimuth (v1, v2) REAL, DIMENSION(3), INTENT(IN) :: v1, v2 REAL, DIMENSION(3) :: Phi, step, Theta REAL :: del_phi_, del_theta_, lon, lat ! v1 and v2 are Cartesian unit vectors. ! Result is azimuth from v1 to v2, measured at v1, in radians, ! and measured clockwise from North. IF (v1(1) == v2(1)) THEN IF (v1(2) == v2(2)) THEN IF (v1(3) == v2(3)) THEN CALL Lonlat_from_xyz(v1, lon, lat) PRINT "(' Error: Get_azimuth of coincident points at ',F8.3,'E, ',F7.3,'N')", lon, lat WRITE (21,"('Error: Get_azimuth of coincident points at ',F8.3,'E, ',F7.3,'N')") lon, lat STOP END IF END IF END IF step = v2 - v1 CALL Local_Phi (v1, Phi) CALL Local_Theta(v1, Theta) del_phi_ = Dot_3D(step, Phi) del_theta_ = Dot_3D(step, Theta) Get_azimuth = ATAN2(del_phi_, -del_theta_) END FUNCTION Get_azimuth CHARACTER(80) FUNCTION Get_filename (unit) ! Obtains a filename from the beginning of the next line; ! name may be padded with blanks on both sides, and ! may be followed by comments. Note that ! "unit" should have been opened with PAD = "YES". CHARACTER(132) :: buffer INTEGER :: i LOGICAL :: past INTEGER, INTENT (IN) :: unit ! Fortran device number READ (unit,"(A)") buffer buffer = ADJUSTL(buffer) ! left-justify past = .FALSE. ! will be T when past end of filename blank_right: DO i=2,132 IF ((buffer(i:i) == ' ') .OR. & (buffer(i:i) == ',') .OR. & (buffer(i:i) == '=') .OR. & (buffer(i:i) == ':')) past = .TRUE. if (past) buffer(i:i) = ' ' END DO blank_right IF (((buffer(1:1) == 'N') .OR. (buffer(1:1) == 'n')) .AND. & ((buffer(2:2) == 'O') .OR. (buffer(2:2) == 'o')) .AND. & ((buffer(3:3) == 'N') .OR. (buffer(3:3) == 'n')) .AND. & ((buffer(4:4) == 'E') .OR. (buffer(4:4) == 'e')) .AND. & (buffer(5:5) == ' ')) buffer = 'none' Get_filename = buffer(1:80) END FUNCTION Get_filename SUBROUTINE Gjxy (l_, r_, G) INTEGER, INTENT(IN) :: l_ ! element number REAL, DIMENSION(3), INTENT(IN) :: r_ ! position vector REAL, DIMENSION (3,2,2), INTENT(OUT) :: G ! computes matrix of 6 vector nodal functions for element l_ at ! position r_ (Cartesian unit vector). ! It is user's responsibility that element l_ contains r_. SAVE ! allows fast re-entry when l_ is unchanged. INTEGER :: l_last = 0 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2) INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function REAL, DIMENSION(3,2) :: local ! local unit vectors at r_ (xyz, SE) REAL, DIMENSION(3,3) :: corner ! positions vector of corner nodes (xyz, 123) REAL, DIMENSION(3,3,2) :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) REAL, DIMENSION(3) :: tvi, tvo, tv1, tv2, tv3, vf ! temporary vector factor REAL :: f_sup_j ! as in Kong and Bird (1995) [ j == k ] INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ DO j = 1, 3 corner(1:3, j) = xyz_nod(1:3, node(j, l_)) tvi = corner(1:3, j) CALL Local_Theta(tvi, tvo) post(1:3, j, 1) = tvo CALL Local_Phi (tvi, tvo) post(1:3, j, 2) = tvo END DO END IF ! begin computations which depend on r_ CALL Local_Theta(r_, tvo) local(1:3,1) = tvo CALL Local_Phi(r_, tvo) local(1:3,2) = tvo DO j = 1, 3 i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1 = corner(1:3, i1) tv2 = corner(1:3, i2) tv3 = corner(1:3, i3) CALL Cross(tv2, tv3, vf) f_sup_j = Dot_3D(r_, vf) / Dot_3D (tv1, vf) DO x = 1, 2 tv1 = post(1:3, j, x) DO y = 1, 2 tv2 = local(1:3, y) G(j, x, y) = f_sup_j * Dot_3D(tv1, tv2) END DO END DO END DO END SUBROUTINE Gjxy CHARACTER(80) FUNCTION Insert (filename, mmnn) ! truncates left part of filename to 4 or less bytes, and adds 'mmnn'; ! should be able to handle complicated names like: ! "../project/restore/NA5Ablah.FEG " ! (in either DOS or Unix) CHARACTER(80) :: filename CHARACTER(4) :: mmnn INTEGER :: left_frame, right_frame, old_stub, new_stub left_frame = MAX (INDEX (filename,'\',.TRUE.), INDEX (filename,'/',.TRUE.)) right_frame = INDEX (filename,'.',.TRUE.) IF ((right_frame == 0) .OR. (right_frame < left_frame)) THEN right_frame = INDEX (filename,' ') END IF old_stub = (right_frame - left_frame) - 1 new_stub = MIN (old_stub, 4) Insert = filename(1 : (left_frame + new_stub)) // mmnn // & TRIM(filename(right_frame : 80)) END FUNCTION Insert SUBROUTINE Internal (b_, iele, s1, s2, s3) REAL, DIMENSION(3), INTENT(IN) :: b_ INTEGER, INTENT(OUT) :: iele REAL, INTENT(OUT) :: s1, s2, s3 ! Input is a Cartesian unit vector in the unit sphere. ! Output is element number and s1, s2, s3 internal coordinates of the ! plane triangle finite element. INTEGER :: back1, back2, i, iet, l_ REAL :: r2, r2min, s1t, s2t, s3t REAL, DIMENSION(3) :: s_temp, tv ! establish defaults (not found) in case of quick exit iele = 0 ! default s1 = 0.0; s2 = 0.0; s3 = 0.0 ! default !find closest element center to initialize search r2min = 4.01 DO l_ = 1, num_ele r2 = (b_(1) - center(1,l_))**2 +(b_(2) - center(2,l_))**2 +(b_(3) - center(3,l_))**2 IF (r2 < r2min) THEN r2min = r2 iet = l_ END IF END DO ! If closest element center is more than 1 radian away, give up. tv = center(1:3, iet) IF (Dot_3D(b_, tv) < 0.540) RETURN ! initialize search memory (with impossible numbers) back1 = -1 back2 = -2 is_it_here: DO ! first, check for infinite loop between 2 elements! IF (iet == back2) THEN ! in loop; force location in one or the other! CALL Dumb_s123 (iet, b_, s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! normal operation CALL Dumb_s123 (iet, b_, s1t, s2t, s3t) IF ((s1t < s2t) .AND. (s1t < s3t)) THEN ! s1 is most negative; most critical IF (s1t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t < s1t) .AND. (s2t < s3t)) THEN ! s2 is most negative; most critical IF (s2t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3 is most negative; most critical IF (s3t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF END IF END IF ! in/not in a loop END DO is_it_here ! successful completion iele = iet s1 = s1t s2 = s2t s3 = s3t END SUBROUTINE Internal SUBROUTINE Interpolate (coordinates, v) ! Input is an internal coordinate set. ! Output is v = unit Cartesian (xyz) location vector TYPE(is123) :: coordinates REAL, DIMENSION(3), INTENT(OUT) :: v INTEGER :: i1, i2, i3, iele REAL :: s1, s2, s3 REAL, DIMENSION(3) :: vt iele = coordinates%element i1 = node(1, iele) i2 = node(2, iele) i3 = node(3, iele) s1 = coordinates%s(1) s2 = coordinates%s(2) s3 = coordinates%s(3) vt = s1 * xyz_nod(1:3,i1) + s2 * xyz_nod(1:3,i2) + s3 * xyz_nod(1:3,i3) CALL Unitise(vt, v) END SUBROUTINE Interpolate SUBROUTINE Local_Phi (b_, Phi) ! returns local East-pointing unit vector in Cartesian coordinates ! for location b_; not intended to work at the poles! REAL, DIMENSION(3), INTENT(IN) :: b_ REAL, DIMENSION(3), INTENT(OUT) :: Phi REAL, DIMENSION(3) :: temp IF (b_(1) == 0.) THEN IF (b_(2) == 0.) THEN PRINT "(' Error: Local_Phi was requested for N or S pole.')" WRITE (21,"('Error: Local_Phi was requested for N or S pole.')") STOP END IF END IF temp(1) = - b_(2) temp(2) = b_(1) temp(3) = 0. CALL Unitise(temp, Phi) END SUBROUTINE Local_Phi SUBROUTINE Local_Theta (b_, Theta) ! returns local South-pointing unit vector in Cartesian coordinates ! for location b_; not intended to work at the poles! REAL, DIMENSION(3), INTENT(IN) :: b_ REAL, DIMENSION(3), INTENT(OUT) :: Theta REAL, DIMENSION(3) :: temp REAL :: equat, new_equat equat = SQRT(b_(1)**2 + b_(2)**2) !equatorial component IF (equat == 0.) THEN PRINT "(' Error: Local_Theta was requested for N or S pole.')" WRITE (21,"('Error: Local_Theta was requested for N or S pole.')") STOP END IF new_equat = b_(3) ! swap components in a meridional plane temp(3) = - equat ! " temp(1) = new_equat * b_(1) / equat ! partition new equatorial component temp(2) = new_equat * b_(2) / equat ! " CALL Unitise(temp, Theta) END SUBROUTINE Local_Theta SUBROUTINE Lonlat_from_xyz (b_, lon, lat) REAL, DIMENSION(3), INTENT(IN) :: b_ ! Cartesian unit vector from center of planet REAL, INTENT(OUT) :: lon, lat REAL :: equat equat = b_(1)**2 + b_(2)**2 IF (equat == 0.) THEN ! N or S pole lon = 0. ! arbitrary convention IF (b_(3) > 0.) THEN lat = 90. ELSE lat = -90. END IF ELSE lat = ATAN2 (b_(3), SQRT(equat)) lon = ATAN2 (b_(2), b_(1)) lon = lon * deg_per_rad lat = lat * deg_per_rad END IF END SUBROUTINE Lonlat_from_xyz REAL FUNCTION Magnitude (b_) REAL, DIMENSION(3), INTENT(IN) :: b_ Magnitude = SQRT(b_(1)**2 +b_(2)**2 + b_(3)**2) END FUNCTION Magnitude CHARACTER(4) FUNCTION Mangle (iteration, time) INTEGER :: iteration, it REAL :: time, tMa CHARACTER(2) :: c2l, c2r tMa = time / s_per_Ma it = NINT(tMa) IF ((iteration < 360) .AND. (it < 360)) THEN IF (iteration < 100) THEN ! use decimal integers WRITE (c2l,"(I2)") iteration !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. ELSE ! Use A .. Z in the tens digit c2l(1:1)=CHAR(64 + INT(iteration/10.) - 9) ! A..Z ! Use 0 .. 9 in the ones digit c2l(2:2)=CHAR(48 + MOD(iteration,10)) ! 0..9 ENDIF IF (it < 100) THEN ! use decimal integers WRITE (c2r,"(I2)") it !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. ELSE ! Use A .. Z in the tens digit c2r(1:1)=CHAR(64 + INT(it/10.) - 9) ! A..Z ! Use 0 .. 9 in the ones digit c2r(2:2)=CHAR(48 + MOD(it,10)) ! 0..9 ENDIF ELSE PRINT "(' Error in Mangle: cannot handle ',2I10)", iteration, it WRITE (21,"(' Error in Mangle: cannot handle ',2I10)") iteration, it STOP END IF IF (c2l(1:1) == ' ') c2l(1:1) = '0' IF (c2r(1:1) == ' ') c2r(1:1) = '0' Mangle = c2l // c2r END FUNCTION Mangle SUBROUTINE More_mem (array_name, bytes_added) ! Keeps track of total array allocation CHARACTER(*) :: array_name ! literal text INTEGER :: bytes_added REAL :: Mb_added, Mb_total CHARACTER(80) :: buffer memory = memory + bytes_added Mb_added = REAL(bytes_added) / REAL(bytes_per_Mb) Mb_total = REAL(memory) / REAL(bytes_per_Mb) WRITE (buffer,"('Allocated ',A,' =',F8.3,' Mb, total',F8.3,' Mb')") & array_name, Mb_added, Mb_total !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. buffer = ADJUSTR(buffer) PRINT "(A)", buffer WRITE (21,"(A)") buffer END SUBROUTINE More_mem SUBROUTINE Move_data INTEGER :: i TYPE(is123) :: m1 REAL, DIMENSION(3) :: tvo, tv1, tv2 ! relocates all integrated positions and orientations (except xyz_nod). IF (ALLOCATED(center)) THEN m1%s(1:3) = 1./3. DO l_ = 1, num_ele m1%element = l_ CALL Interpolate(m1, tvo) center(1:3, l_) = tvo END DO END IF IF (f_dig_count > 0) THEN DO i = 1, f_dig_count IF (trace_is(i)%element > 0) THEN CALL Interpolate(trace_is(i), tvo) trace(1:3, i) = tvo END IF END DO END IF IF (seg_count > 0) THEN DO i = 1, seg_count IF (seg_end_is(1, i)%element > 0) THEN CALL Interpolate(seg_end_is(1, i), tvo) seg_end(1:3, 1, i) = tvo END IF IF (seg_end_is(2, i)%element > 0) THEN CALL Interpolate(seg_end_is(2, i), tvo) seg_end(1:3, 2, i) = tvo END IF END DO END IF IF (c_dat_count > 0) THEN DO i = 1, c_dat_count IF (c_end_is(1,i)%element > 0) THEN CALL Interpolate(c_end_is(1, i), tvo) c_end_now(1:3, 1, i) = tvo END IF IF (c_end_is(2,i)%element > 0) THEN CALL Interpolate(c_end_is(2, i), tvo) c_end_now(1:3, 2, i) = tvo END IF END DO END IF IF (p_dat_count > 0) THEN DO i = 1, p_dat_count IF (p_site_is(i)%element > 0) THEN CALL Interpolate(p_site_is(i), tvo) p_site_now(1:3, i) = tvo END IF END DO END IF IF (s_dat_count > 0) THEN DO i = 1, s_dat_count IF ((s_site_is(1,i)%element > 0).AND.(s_site_is(2,i)%element > 0)) THEN CALL Interpolate(s_site_is(1, i), tvo) s_site_now(1:3, 1, i) = tvo CALL Interpolate(s_site_is(2, i), tvo) s_site_now(1:3, 2, i) = tvo tv1 = s_site_now(1:3, 1, i) tv2 = s_site_now(1:3, 2, i) s_azim_now(i) = Get_azimuth(tv1, tv2) END IF END DO END IF IF (y_dig_count > 0) THEN DO i = 1, y_dig_count IF (dot_is(i)%element > 0) THEN CALL Interpolate(dot_is(i), tvo) dot(1:3, i) = tvo END IF END DO END IF END SUBROUTINE Move_data SUBROUTINE Move_feg (vw) REAL, DIMENSION(ndof), INTENT(IN) :: vw INTEGER :: i REAL, DIMENSION(3) :: tvi, tvo DO i = 1, num_nod tvi = xyz_nod(1:3, i) CALL Moved_by_vw (tvi, vw(2*i-1), vw(2*i), .FALSE., tvo) xyz_nod(1:3,i) = tvo END DO END SUBROUTINE Move_feg SUBROUTINE Moved_by_vel (old_vec, velocity, forward, new_vec) ! Computes new position (Cartesian unit vector) ! after a great-circle finite rotation with fixed angular ! velocity which begins at position old_vec ! with initial velocity vector 'velocity' (Cartesian, m/s) ! and continues for time interval Deltat_. ! IF (.NOT.forward), rotation is in opposite direction. REAL, DIMENSION(3), INTENT(IN) :: old_vec, velocity LOGICAL, INTENT(IN) :: forward REAL, DIMENSION(3), INTENT(OUT):: new_vec REAL, DIMENSION(3) :: temp, v_unit REAL :: angle, cos_a, sin_a, vel_mag vel_mag = SQRT(velocity(1)**2 + velocity(2)**2 + velocity(3)**2) angle = vel_mag * Deltat_ / R ! Deltat_, R are global. IF (angle == 0.) THEN new_vec = old_vec ELSE IF (.NOT. forward) angle = - angle CALL Unitise(velocity, v_unit) sin_a = SIN(angle) cos_a = COS(angle) temp = cos_a * old_vec + sin_a * v_unit CALL Unitise(temp, new_vec) ! just for insurance END IF END SUBROUTINE Moved_by_vel SUBROUTINE Moved_by_vw (old_vec, v_, w_, forward, new_vec) ! Computes new position (Cartesian unit vector) ! after a great-circle finite rotation with fixed angular ! velocity which begins at position old_vec ! with initial velocity components v_ (theta, South) and ! w_ (phi, East), and continues for time interval Deltat_. ! IF (.NOT.forward), rotation is in opposite direction. REAL, DIMENSION(3), INTENT(IN) :: old_vec REAL, INTENT(IN) :: v_, w_ LOGICAL, INTENT(IN) :: forward REAL, DIMENSION(3), INTENT(OUT) :: new_vec REAL, DIMENSION(3) :: Phi, Theta, veloc REAL :: angle, vel_mag vel_mag = SQRT(v_**2 + w_**2) angle = vel_mag * Deltat_ / R ! Deltat_, R are global. IF (angle == 0.) THEN new_vec = old_vec ELSE CALL Local_Theta (old_vec, Theta) CALL Local_Phi (old_vec, Phi) veloc = v_ * Theta + w_ * Phi CALL Moved_by_vel (old_vec, veloc, forward, new_vec) END IF END SUBROUTINE Moved_by_vw SUBROUTINE New_goal (count, total, tmin, tmax, checkPD, rate, & ! inputs & goal, active, n_adjusted) ! outputs ! refer to equation (27) of Bird (1997?) and adjacent paragraphs INTEGER, INTENT(IN) :: count ! number of data in arrays REAL,DIMENSION(:), INTENT(IN) :: total ! offset, displacement, angle REAL,DIMENSION(:), INTENT(IN) :: tmin ! ending age REAL,DIMENSION(:), INTENT(IN) :: tmax ! starting age LOGICAL, INTENT(IN) :: checkPD ! check for P, D (fault offsets only) REAL,DIMENSION(:,:),INTENT(IN) :: rate ! results of previous iteration REAL,DIMENSION(:,:),INTENT(OUT) :: goal ! targets for next iteration LOGICAL(1),DIMENSION(:,:),INTENT(OUT):: active ! (may be reset by Set-goal-A) INTEGER, INTENT(OUT) :: n_adjusted ! number of data with rescaled (not reset) goals INTEGER :: i, j REAL :: factor, sum, sum_goal LOGICAL :: adjust n_adjusted = 0 DO i = 1, count IF (checkPD .AND. (f_dat_code(i) == 'P')) THEN ! reassign the same Promoted goal as in first iteration CALL Set_goal_A (i, total, tmin, tmax, checkPD, goal, active) ELSE ! usual case ... sum_goal = total(i) / Deltat_ IF (sum_goal > 0.) THEN ! we hope most rates are + sum = 0. adjust = .FALSE. DO j = 1, num_timesteps sum = sum + MAX(0., rate(j, i)) !count only positive rates adjust = adjust .OR. (rate(j,i) > goal(j,i)) END DO IF (sum > 0.) THEN IF (adjust) THEN n_adjusted = n_adjusted + 1 factor = sum_goal / sum IF ((tmax(i) <= end_time) .OR. (factor < 1.)) THEN ! normal case; scale goals from rates DO j = 1, num_timesteps goal(j,i) = MAX(0., factor * rate(j,i)) END DO ELSE ! hang-over datum AND factor > 1; don't scale; let drift. DO j = 1, num_timesteps goal(j,i) = MAX(0., rate(j,i)) END DO END IF END IF ! adjust ELSE ! wrong way; must reassign goals from scratch CALL Set_goal_A (i, total, tmin, tmax, checkPD, goal, active) END IF ELSE IF (sum_goal == 0.) THEN DO j = 1, num_timesteps goal(j,i) = 0. END DO ELSE ! goal is -; we hope most rates are - sum = 0. adjust = .FALSE. DO j = 1, num_timesteps sum = sum + MIN(0., rate(j, i)) !count only negative rates adjust = adjust .OR. (rate(j,i) < goal(j,i)) END DO IF (sum < 0.) THEN IF (adjust) THEN n_adjusted = n_adjusted + 1 factor = sum_goal / sum ! factor is still + = - / - IF ((tmax(i) <= end_time) .OR. (factor < 1.)) THEN ! normal case; scale goals from rates DO j = 1, num_timesteps goal(j,i) = MIN(0., factor * rate(j,i)) END DO ELSE ! hang-over datum AND factor > 1; don't scale; let drift. DO j = 1, num_timesteps goal(j,i) = MIN(0., rate(j,i)) END DO END IF END IF ! adjust ELSE ! wrong way; must reassign goals from scratch CALL Set_goal_A (i, total, tmin, tmax, checkPD, goal, active) END IF END IF ! checkPD and 'P' END IF ! goal is - IF (checkPD .AND. (f_dat_code(i) == 'D')) goal(1, i) = 0. END DO ! on i = datum number END SUBROUTINE New_goal SUBROUTINE Plane_area (folding) LOGICAL, INTENT(OUT) :: folding !puts areas of plane triangles (below surface) into array a_; !if any is zero or negative, reports folding INTEGER :: i1, i2, i3, l_ REAL, DIMENSION(3) :: a, b, c, t folding = .FALSE. DO l_ = 1, num_ele ! global, like most variables i1 = node(1,l_) i2 = node(2,l_) i3 = node(3,l_) a = xyz_nod(1:3,i2) - xyz_nod(1:3,i1) b = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) CALL Cross (a, b, c) t = xyz_nod(1:3,i1) + xyz_nod(1:3,i2) + xyz_nod(1:3,i3) IF (Dot_3D(t, c) > 0.) THEN a_(l_) = Magnitude(c) * half_R2 ELSE folding = .TRUE. RETURN END IF END DO END SUBROUTINE Plane_area SUBROUTINE Plug_in_33 (l_,A,B,C,D,E,F) ! Completes upper triangle of element (3-node x 3-node) matrix, then ! adds element matrix and element forcing vector to global system ! augmented-matrix ABCDEF (which has a proprietary storage layout). INTEGER, INTENT(IN) :: l_ ! element number, to access nodes REAL, DIMENSION(3,3):: A, B, C, D ! submatrices of element matrix REAL, DIMENSION(3) :: E, F ! subvectors of element forcing vector INTEGER :: i, it, j, jt, m, n INTEGER :: ABCDrow, ABCDcol, EFrow, EFcol ! statement functions ABCDrow(j) = j + ncoda ! statement function; ncoda is global ABCDcol(i, j) = (j - i) + 1 ! statement function EFrow(i) = i + ncoda ! statement function; ncoda is global EFcol = ncoda + 2 ! These statement functions are for ! codiagonal band symmetric storage mode of matrices ABCD and EF, ! per Microsoft version of IMSL. Note that element (row #i, col #j) ! of the idealized square matrix square_ABCD, or ! square_ABCD(i, j) is stored as ABCDEF(ABCDrow(j),ABCDcol(i, j)); ! and that only elements with j >= i (upper right) can be stored. ! Element (row i) of the linear vector EF is stored in ! ABCDEF(EFrow(i), EFcol). ! upper-triangle values filled in by symmetry A(1,2) = A(2,1); A(1,3) = A(3,1); A(2,3) = A(3,2) B = TRANSPOSE(C) D(1,2) = D(2,1); D(1,3) = D(3,1); D(2,3) = D(3,2) ! add element matrix to upper right part of global stiffness matrix: DO m = 1, 3 ! row within A i = 2 * node(m, l_) - 1 ! logical row in square coefficient matrix DO n = 1, 3 ! column of A j = 2 * node(n, l_) - 1 ! logical column in square coefficient matrix IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + A(m, n) END IF END DO END DO DO m = 1, 3 ! row within B i = 2 * node(m, l_) - 1 DO n = 1, 3 ! column of B j = 2 * node(n, l_) IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + B(m, n) END IF END DO END DO DO m = 1, 3 ! row within C i = 2 * node(m, l_) DO n = 1, 3 ! column of C j = 2 * node(n, l_) - 1 IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + C(m, n) END IF END DO END DO DO m = 1, 3 ! row within D i = 2 * node(m, l_) DO n = 1, 3 ! column of D j = 2 * node(n, l_) IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + D(m, n) END IF END DO END DO ! add element's forcing vectors to global system: DO m = 1, 3 i = 2 * node(m, l_) - 1 ! logical row ABCDEF(EFrow(i), EFcol) = ABCDEF(EFrow(i), EFcol) + E(m) i = i + 1 ABCDEF(EFrow(i), EFcol) = ABCDEF(EFrow(i), EFcol) + F(m) END DO END SUBROUTINE Plug_in_33 SUBROUTINE Plug_in_66 (l_1,l_2,A,B,C,D,E,F) ! Completes upper triangle of a super-element ! (6-node x 6-node) matrix, then ! adds element matrix and element forcing vector to global system ! augmented-matrix ABCDEF (which has a proprietary storage layout). INTEGER, INTENT(IN) :: l_1, l_2 ! element numbers, to access nodes REAL, DIMENSION(6,6):: A, B, C, D ! submatrices of super-element matrix REAL, DIMENSION(6) :: E, F ! subvectors of super-element forcing vector INTEGER :: i, it, j, jt, m, n INTEGER :: ABCDrow, ABCDcol, EFrow, EFcol ! statement functions ABCDrow(j) = j + ncoda ! statement function; ncoda is global ABCDcol(i, j) = (j - i) + 1 ! statement function EFrow(i) = i + ncoda ! statement function; ncoda is global EFcol = ncoda + 2 ! These statement functions are for ! codiagonal band symmetric storage mode of matrices ABCD and EF, ! per Microsoft version of IMSL. Note that element (row #i, col #j) ! of the idealized square matrix square_ABCD, or ! square_ABCD(i, j) is stored as ABCDEF(ABCDrow(j),ABCDcol(i, j)); ! and that only elements with j >= i (upper right) can be stored. ! Element (row i) of the linear vector EF is stored in ! ABCDEF(EFrow(i), EFcol). ! upper-triangle values filled in by symmetry A(1,2)=A(2,1); A(1,3)=A(3,1); A(1,4)=A(4,1); A(1,5)=A(5,1); A(1,6)=A(6,1) A(2,3)=A(3,2); A(2,4)=A(4,2); A(2,5)=A(5,2); A(2,6)=A(6,2) A(3,4)=A(4,3); A(3,5)=A(5,3); A(3,6)=A(6,3) A(4,5)=A(5,4); A(4,6)=A(6,4) A(5,6)=A(6,5) B = TRANSPOSE(C) D(1,2)=D(2,1); D(1,3)=D(3,1); D(1,4)=D(4,1); D(1,5)=D(5,1); D(1,6)=D(6,1) D(2,3)=D(3,2); D(2,4)=D(4,2); D(2,5)=D(5,2); D(2,6)=D(6,2) D(3,4)=D(4,3); D(3,5)=D(5,3); D(3,6)=D(6,3) D(4,5)=D(5,4); D(4,6)=D(6,4) D(5,6)=D(6,5) ! add element matrix to upper right part of global stiffness matrix: DO m = 1, 6 ! row within A IF (m <= 3) THEN i = 2 * node(m, l_1) - 1 ! logical row in square coefficient matrix ELSE i = 2 * node(m-3, l_2) - 1 ! logical row in square coefficient matrix END IF DO n = 1, 6 ! column of A IF (n <= 3) THEN j = 2 * node(n, l_1) - 1 ! logical column in square coefficient matrix ELSE j = 2 * node(n-3, l_2) - 1 ! logical column in square coefficient matrix END IF IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + A(m, n) END IF END DO END DO DO m = 1, 6 ! row within B IF (m <= 3) THEN i = 2 * node(m, l_1) - 1 ELSE i = 2 * node(m-3, l_2) - 1 END IF DO n = 1, 6 ! column of B IF (n <= 3) THEN j = 2 * node(n, l_1) ELSE j = 2 * node(n-3, l_2) END IF IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + B(m, n) END IF END DO END DO DO m = 1, 6 ! row within C IF (m <= 3) THEN i = 2 * node(m, l_1) ELSE i = 2 * node(m-3, l_2) END IF DO n = 1, 6 ! column of C IF (n <= 3) THEN j = 2 * node(n, l_1) - 1 ELSE j = 2 * node(n-3, l_2) - 1 END IF IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + C(m, n) END IF END DO END DO DO m = 1, 6 ! row within D IF (m <= 3) THEN i = 2 * node(m, l_1) ELSE i = 2 * node(m-3, l_2) END IF DO n = 1, 6 ! column of D IF (n <= 3) THEN j = 2 * node(n, l_1) ELSE j = 2 * node(n-3, l_2) END IF IF (j >= i) THEN it = ABCDrow(j) ! actual row (sic; i not involved) jt = ABCDcol(i, j) ! actual column ABCDEF(it, jt) = ABCDEF(it, jt) + D(m, n) END IF END DO END DO ! add element's forcing vectors to global system: DO m = 1, 6 IF (m <= 3) THEN i = 2 * node(m, l_1) - 1 ! logical row ELSE i = 2 * node(m-3, l_2) - 1 ! logical row END IF ABCDEF(EFrow(i), EFcol) = ABCDEF(EFrow(i), EFcol) + E(m) i = i + 1 ABCDEF(EFrow(i), EFcol) = ABCDEF(EFrow(i), EFcol) + F(m) END DO END SUBROUTINE Plug_in_66 SUBROUTINE Prediction (vw, L0, L1, L2) ! Computes actual model predictions of rates (p_). ! Also computes 3 norms (L0, L1, L2) of rate errors (all data ! types and a-priori merged) in non-dimensional sigma units. REAL, DIMENSION(:), INTENT(IN) :: vw REAL, INTENT(OUT):: L0, L1, L2 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: B, X !NOTE: Double precision is to prevent underflow/overflow ! during solution of linear system, not for great accuracy. REAL :: cott REAL :: csct INTEGER :: datum REAL, DIMENSION(3,2,2,2) :: dG REAL :: equat REAL, DIMENSION(3) :: eps_dot, eps_dot_c REAL :: error REAL :: error_count REAL, DIMENSION(3) :: f_, g_ REAL, DIMENSION(3,2,2) :: G REAL :: gamma_b, gamma_d INTEGER :: iv, iw, i1, i2, i3 INTEGER :: l_, LDA, loc REAL :: L0_sum, L1_sum, L2_sum REAL :: misfits REAL :: mu_of_r_ INTEGER :: N REAL, DIMENSION(3) :: outward INTEGER :: p REAL :: p_ REAL, DIMENSION(3) :: Phi DOUBLE PRECISION :: prefix REAL, DIMENSION(3) :: r_ REAL :: rho_ INTEGER :: segment REAL :: theta_ REAL, DIMENSION(3) :: Theta REAL, DIMENSION(3) :: tv, tv1, tv2 REAL :: sint REAL :: tant, t_sigma REAL :: v_, w_ REAL, DIMENSION(3) :: veloc INTEGER :: z_, Z L0_sum = 0. L1_sum = 0. L2_sum = 0. error_count = 0. ! Continuum strain-rate and active fault segments (cracks) IF (f_dat_count > 0) f_divide = 0. ! all numerators and denominators DO l_ = 1, num_ele i1 = node(1, l_) i2 = node(2, l_) i3 = node(3, l_) ! evaluate nodal function and derivitives at center of element tv = center(1:3, l_) CALL Gjxy(l_, tv, G) CALL Del_Gjxy_del_thetaphi(l_, tv, dG) equat = SQRT(tv(1)**2 + tv(2)**2) IF (equat == 0.) THEN PRINT "(' Error: center of element ',I5,' is N or S pole.')", l_ WRITE (21,"('Error: center of element ',I5,' is N or S pole.')") l_ STOP END IF theta_ = ATAN2(equat, tv(3)) sint = SIN(theta_) csct = 1. / sint tant = TAN(theta_) cott = 1. / tant CALL E_rate(l_, G, dG, theta_, vw, eps_dot) ! evaluate mu_ at center mu_of_r_ = (mu_nod(i1) + mu_nod(i2) + mu_nod(i3)) / 3. IF (f_dat_count > 0) THEN Z = crack_index(1, l_) ELSE Z = 0 END IF IF (Z > 0) THEN ! element has active cracks ! Allocate and build linear system N = Z + 3 + 3 ALLOCATE ( A(N,N) ) ALLOCATE ( B(N) ) ALLOCATE ( X(N) ) LDA = N A = 0. ! whole matrix B = 0. ! whole vector DO z_ = 1, Z loc = crack_index(2, l_) + z_ - 1 ! storage location in local_crack t_sigma = f_scale * ((local_crack(loc)%sigma_/f_scale)**exponent) A(z_, z_) = 2. / (t_sigma**2) A(z_, Z + 4) = local_crack(loc)%H(1) A(z_, Z + 5) = local_crack(loc)%H(2) A(z_, Z + 6) = local_crack(loc)%H(3) B(z_) = 2. * local_crack(loc)%s_ / (t_sigma**2) END DO ! z_ = 1, Z t_sigma = mu_scale * ((mu_of_r_ / mu_scale)**exponent) prefix = 2. / (t_sigma**2) A(Z + 1, Z + 1) = prefix A(Z + 2, Z + 2) = prefix A(Z + 3, Z + 3) = prefix A(Z + 1, Z + 3) = 0.5 * prefix A(Z + 1, Z + 4) = 1. A(Z + 2, Z + 5) = 1. A(Z + 3, Z + 6) = 1. B(Z + 4) = eps_dot(1) B(Z + 5) = eps_dot(2) B(Z + 6) = eps_dot(3) ! Solve a real symmetric system of linear equations without iterative refinement. ! Usage CALL DLSLSF (N, A, LDA, B, X) ! NOTE: Double-precision version will reduce underflow/overflow problems! ! Arguments ! N = Number of equations. (Input) ! A = N by N matrix containing the coefficient matrix of the symmetric linear system. (Input) ! Only the upper triangle of A is referenced. ! LDA = Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input) ! B = Vector of length N containing the right-hand side of the linear system. (Input) ! X = Vector of length N containing the solution to the linear system. (Output) DO z_ = 1, Z p_ = X(z_) ! PREDICTED SLIP-RATE OF THIS CRACK loc = crack_index(2, l_) + z_ - 1 ! storage location in local_crack segment = local_crack(loc)%segment tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) rho_ = R * Arc_distance(tv1, tv2) datum = local_crack(loc)%datum f_divide(1, datum) = f_divide(1, datum) + rho_ * p_ f_divide(2, datum) = f_divide(2, datum) + rho_ error = (p_ - local_crack(loc)%s_) / local_crack(loc)%sigma_ IF (ABS(error) > 2.) L0_sum = L0_sum + 1. L1_sum = L1_sum + ABS(error) L2_sum = L2_sum + error**2 error_count = error_count + 1. END DO ! z_ = 1, Z eps_dot_c(1) = X(Z + 1) eps_dot_c(2) = X(Z + 2) eps_dot_c(3) = X(Z + 3) IF (stress_ever) ele_strainrate(1:3,l_) = eps_dot_c(1:3) ! saved for Write_x_feg error = SQRT( eps_dot_c(1)**2 + & & eps_dot_c(1) * eps_dot_c(3) + & & eps_dot_c(3)**2 + & & eps_dot_c(2)**2 ) / mu_of_r_ IF (ABS(error) > 2.) L0_sum = L0_sum + 1.5 L1_sum = L1_sum + ABS(error) * 1.5 L2_sum = L2_sum + error**2 * 1.5 error_count = error_count + 1.5 DEALLOCATE ( A, B, X ) ELSE ! no active cracks in this element IF (stress_ever) ele_strainrate(1:3,l_) = eps_dot(1:3) ! saved for Write_x_feg error = SQRT( eps_dot(1)**2 + & & eps_dot(1) * eps_dot(3) + & & eps_dot(3)**2 + & & eps_dot(2)**2 ) / mu_of_r_ IF (ABS(error) > 2.) L0_sum = L0_sum + W * a_(l_) L1_sum = L1_sum + ABS(error) * W * a_(l_) L2_sum = L2_sum + error**2 * W * a_(l_) error_count = error_count + W * a_(l_) END IF ! active cracks / no active cracks END DO ! l_ = 1, num_ele ! nominal rate for each active fault is average over segments IF (f_dat_count > 0) THEN DO i = 1, f_dat_count IF (f_divide(2, i) > 0.) THEN IF (sense(i) == 'T') THEN factor = -cot_thrust_dip ELSE IF (sense(i) == 'P') THEN factor = -1. ELSE IF (sense(i) == 'N') THEN factor = cot_normal_dip ELSE IF (sense(i) == 'D') THEN factor = 1. ELSE IF (sense(i) == 'R') THEN factor = 1. ELSE IF (sense(i) == 'L') THEN factor = -1. ENDIF f_rate(n_, i) = (f_divide(1, i) / f_divide(2, i)) / factor END IF END DO ! i = 1, f_dat_count END IF ! f_dat_count > 0 ! Active cross-sections, if entirely within domain IF (c_dat_count > 0) THEN DO k = 1, c_dat_count IF (c_active(n_,k)) THEN IF ((c_end_is(1,k)%element > 0).AND. & &(c_end_is(2,k)%element > 0)) THEN tv1 = c_end_now(1:3, 1, k) tv2 = c_end_now(1:3, 2, k) gamma_b = Get_azimuth(tv1, tv2) gamma_d = Pi + Get_azimuth(tv2, tv1) ! Note: These both point from b to d, but are evaluated at ! different locations b and d respectively. ! WORK ON WEST END: l_ = c_end_is(1,k)%element CALL Gjxy(l_, tv1, G) CALL Components(l_,G,vw,v_,w_) CALL Local_Theta(tv1, Theta) CALL Local_Phi (tv1, Phi) veloc = v_ * Theta + w_ * Phi outward = COS(gamma_b) * Theta - SIN(gamma_b) * Phi p_ = Dot_3D(veloc, outward) ! WORK ON EAST END: l_ = c_end_is(2,k)%element CALL Gjxy(l_, tv2, G) CALL Components(l_,G,vw,v_,w_) CALL Local_Theta(tv2, Theta) CALL Local_Phi (tv2, Phi) veloc = v_ * Theta + w_ * Phi outward = -COS(gamma_d) * Theta + SIN(gamma_d) * Phi p_ = p_ + Dot_3D(veloc, outward) misfits = ABS(p_ - c_goal(n_,k)) / c_rate_sigma_(k) error_count = error_count + 1. IF (misfits > 2.) L0_sum = L0_sum + 1. L1_sum = L1_sum + misfits L2_sum = L2_sum + misfits**2 c_rate(n_,k) = p_ END IF ! this section is entirely in domain END IF ! this section is active in this timestep END DO ! k = 1, c_dat_count END IF ! c_dat_count > 0 ! Paleomagnetic sites IF (p_dat_count > 0) THEN DO p = 1, p_dat_count IF (p_active(n_,p)) THEN l_ = p_site_is(p)%element IF (l_ > 0) THEN r_ = p_site_now(1:3,p) CALL Gjxy(l_, r_, G) ! paleolatitude part tv = p_pole(1:3, p) gamma_ = Get_azimuth(r_, tv) CALL Components(l_,G,vw,v_,w_) p_ = v_ * COS(gamma_) - w_ * SIN(gamma_) misfits = ABS(p_ - p_south_goal(n_,p)) / p_south_rate_sigma_(p) error_count = error_count + 1. IF (misfits > 2.) L0_sum = L0_sum + 1. L1_sum = L1_sum + misfits L2_sum = L2_sum + misfits**2 p_south_rate(n_,p) = p_ ! vertical-axis rotation part IF (.NOT.twisted(p)) THEN CALL Del_Gjxy_del_thetaphi(l_, r_, dG) equat = SQRT(r_(1)**2 + r_(2)**2) theta_ = ATAN2(equat, r_(3)) sint = SIN(theta_) csct = 1. / sint tant = TAN(theta_) cott = 1. / tant prefix = 1. / (2. * R) p_ = 0. DO j = 1, 3 f_(j) = prefix * (G(j,1,2)*cott + dG(j,1,2,1) - csct*dG(j,1,1,2)) g_(j) = prefix * (G(j,2,2)*cott + dG(j,2,2,1) - csct*dG(j,2,1,2)) iv = 2 * node(j,l_) - 1 iw = iv + 1 p_ = p_ + vw(iv) * f_(j) + vw(iw) * g_(j) END DO misfits = ABS(p_ - p_ccw_goal(n_,p)) / p_ccw_rate_sigma_(p) error_count = error_count + 1. IF (misfits > 2.) L0_sum = L0_sum + 1. L1_sum = L1_sum + misfits L2_sum = L2_sum + misfits**2 p_ccw_rate(n_,p) = p_ END IF ! .NOT.twisted END IF ! l_ > 0 END IF ! p_active(n_,p) END DO ! p = 1, p_dat_count END IF ! p_dat_count > 0 ! Finish and report error measures: IF (error_count > 0.) THEN L0 = L0_sum / error_count L1 = L1_sum / error_count L2 = SQRT( L2_sum / error_count ) PRINT "(' ',8X,'Rate errors: L0 = ',F4.3,', L1 = ',F6.3,', L2 = ',F6.3)", L0, L1, L2 WRITE (21,"(8X,'Rate errors: L0 = ',F4.3,', L1 = ',F6.3,', L2 = ',F6.3)") L0, L1, L2 ELSE L0 = 0. ; L1 = 0. ; L2 = 0. END IF END SUBROUTINE Prediction SUBROUTINE Prevent (bad_thing, line, filename) INTEGER, INTENT(IN) :: line CHARACTER(*), INTENT(IN) :: bad_thing, filename PRINT "(' Error: ',A,' is illegal in line ',I6/' of ',A)", & TRIM(bad_thing), line, TRIM(filename) WRITE (21,"('Error: ',A,' is illegal in line ',I6/' of ',A)") & TRIM(bad_thing), line, TRIM(filename) STOP END SUBROUTINE Prevent SUBROUTINE Pull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. REAL, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.) THEN highest = MAXVAL(s) medium = 1.00 - lowest - highest IF (medium > 0.) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0. sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00 / (1.00 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0. s(sideb) = 0. END IF END IF END SUBROUTINE Pull_in SUBROUTINE Set_goal_A (index, total, tmin, tmax, checkPD, & !inputs & goal, active) ! outputs ! Initialize goals as uniform within tmin:tmax, ! except for discretization effect of chopping into timesteps. ! Acts only on datum #index, so typically called from within loop. INTEGER, INTENT(IN) :: index ! which datum? REAL,DIMENSION(:), INTENT(IN) :: total, tmin, tmax ! for all data LOGICAL, INTENT(IN) :: checkPD ! promotions and demotions? REAL,DIMENSION(:,:),INTENT(OUT) :: goal LOGICAL(1), DIMENSION(:,:), INTENT(OUT) :: active REAL :: epsilon, overlap, t0, t1 IF (checkPD) THEN ! applies only to faults IF ((tmin(index) == 0.).AND.(tmax(index) < Deltat_)) THEN f_dat_code(index) = 'P' ! Promoted ELSE f_dat_code(index) = 'N' ! Normal; Demotions done later. END IF END IF IF (neotec) THEN active(1, index) = (start_time >= tmin(index)) .AND. (start_time < tmax(index)) ! if model time is exactly on boundary, then only windows extending to older times ! will match; this is based on idea that neotec models actually span a tiny bit of time. IF (active(1, index)) THEN goal(1, index) = total(index) / (tmax(index) - tmin(index)) ELSE goal(1, index) = 0. END IF ELSE ! paleotec DO j = 1, num_timesteps t0 = (j - 1) * Deltat_ t1 = j * Deltat_ overlap = MIN(t1 - t0, tmax(index) - tmin(index), & & tmax(index) - t0, t1 - tmin(index)) active(j, index) = (overlap > 0.) IF (active(j, index)) THEN epsilon = total(index) * overlap / (tmax(index) - tmin(index)) goal(j, index) = epsilon / Deltat_ ! special kludge: Quaternary fault slip rates apply to whole first timestep IF (checkPD .AND. (j == 1) .AND. (f_dat_code(index) == 'P')) THEN goal(j, index) = total(index) / (tmax(index) - tmin(index)) END IF ELSE goal(j, index) = 0. END IF END DO END IF ! paleotec END SUBROUTINE Set_goal_A SUBROUTINE Set_goal_B (index, checkPD, active, unit, signal, eof, conversion, & !inputs & line, rate, goal) ! modify ! Overwrite the default rates and goals ! (which are constant-rate goals and 0 rates) ! if * or & lines are available in the input file. ! Acts only on datum #index, so typically called from within loop. INTEGER, INTENT(IN) :: index, unit LOGICAL, INTENT(IN) :: checkPD ! promotion/demotion of Q rates LOGICAL(1), DIMENSION(:,:),INTENT(IN) :: active CHARACTER(1), INTENT(IN) :: signal LOGICAL, INTENT(OUT) :: eof REAL, INTENT(IN) :: conversion INTEGER, INTENT(INOUT) :: line REAL,DIMENSION(:,:),INTENT(INOUT) :: goal, rate INTEGER :: j, j0, j1, read_status REAL :: overlap, r0, r1, r2, r3, t0, t1 CHARACTER(134) :: c134 LOGICAL :: new, process, synch seek_stars: DO READ (unit, "(A)", IOSTAT = read_status) c134; line = line + 1 eof = .FALSE. IF (read_status /= 0) THEN eof = .TRUE. RETURN ELSE IF (c134(1:1) == signal) THEN c134 = c134(2:134) // ' ' READ (c134, *) r0, r1, r2, r3 process = .NOT. (checkPD .AND. (f_dat_code(index) == 'P')) ! never change Promoted Quaternary goals IF (process) THEN r0 = r0 * s_per_Ma ! global r1 = r1 * s_per_Ma r2 = r2 * conversion r3 = r3 * conversion IF (paleotec) THEN j0 = NINT(r0 / Deltat_) ! global only if paleotec j1 = NINT(r1 / Deltat_) synch = ((j1 - j0) == 1) .AND. & & (ABS((j0 * Deltat_ - r0) / Deltat_) < 0.01) .AND. & & (ABS((j1 * Deltat_ - r1) / Deltat_) < 0.01) IF (synch) THEN ! 1-to-1 correspondance of input lines with timesteps rate(j1, index) = r2 goal(j1, index) = r3 ELSE ! not synch'ed ! more complex logic when steps are staggered or nested DO j = 1, num_timesteps ! global; = 1 IF (neotec) IF (active(j, index)) THEN t0 = (j - 1) * Deltat_ t1 = j * Deltat_ overlap = MIN (Deltat_, r1 - r0, t1 - r0, r1 - t0) IF (overlap > 0.) THEN ! does [interval) just read include t0 of this timestep ? new = (r0 <= t0) .AND. (r1 > t0) IF (new) THEN !overwrite; later records may add rate(j, index) = r2 * overlap / Deltat_ goal(j, index) = r3 * overlap / Deltat_ ELSE ! add rate(j, index) = rate(j, index) + r2 * overlap / Deltat_ goal(j, index) = goal(j, index) + r3 * overlap / Deltat_ END IF ! new / NOT new END IF ! overlap > 0. END IF ! active(j, index) END DO ! j = 1, num_timesteps END IF ! synch or asynchronous ELSE ! not paleotec, so neotec IF ((r0 <= start_time).AND.(r1 >= start_time)) THEN rate(1, index) = r2 goal(2, index) = r3 END IF END IF ! paleotec / neotec END IF ! process ELSE ! signal not found BACKSPACE (unit) RETURN END IF END DO seek_stars END SUBROUTINE Set_goal_B SUBROUTINE Sigma_1h(stress_count,needles,x_,azimuth,del_az_for_90pc) ! interpolates most-compressive horizontal principal stress ! direction "azimuth" (in radians clockwise from North), ! and reports the (one-sided, or half-) width of the sector ! that contains it with 90%-confidence, as "del_az_for_90pc" (radians). ! The method is the simple "independent-data" method of ! Bird & Li (1996), J. Geophys. Res., v. 101, #B3, 5435-5443. ! The information about the present global stress field is ! in global array ln_rel_prob. INTEGER, INTENT(IN) :: stress_count ! count of data TYPE(needle),DIMENSION(:),INTENT(IN):: needles ! data table REAL, DIMENSION(3),INTENT(IN) :: x_ ! Cartesian unit vector ! from center of Earth ! to interpolation point REAL, INTENT(OUT) :: azimuth, del_az_for_90pc INTEGER :: j, jl, jlo, jr, jro, left, n, neighbors, right INTEGER, DIMENSION(1) :: peak REAL, DIMENSION(3) :: a_ REAL, DIMENSION(0:59) :: probability REAL :: fraction, gamma_a_, gamma_x_, oldsum, prediction, q_, sum, theta_ neighbors = 0 probability = 0. ! initialize array DO i = 1, stress_count a_ = needles(i)%location theta_ = Arc_distance(x_, a_) n = 1.00001 + 150. * (0.5 - 0.5 * COS(theta_))**0.6 IF (n <= 21) THEN neighbors = neighbors + 1 gamma_x_ = Get_azimuth(x_, a_) gamma_a_ = Get_azimuth(a_, x_) + 3.141593 prediction = MOD((needles(i)%azimuth - gamma_a_ + gamma_x_ + 18.84955), 3.141593) q_ = needles(i)%relevance left = 19.0985 * prediction + 0.5 left = MOD(left,60) right = MOD(left + 1, 60) IF (q_ > 0.99) THEN ! omit factor for more speed DO j = 0, 29 jl = MOD((left - j + 60), 60) jr = MOD((right + j), 60) probability(jl) = probability(jl) + ln_rel_prob(n,j) probability(jr) = probability(jr) + ln_rel_prob(n,j) END DO ELSE IF (q_ > 0.0) THEN DO j = 0, 29 jl = MOD((left - j + 60), 60) jr = MOD((right + j), 60) probability(jl) = probability(jl) + q_ * ln_rel_prob(n,j) probability(jr) = probability(jr) + q_ * ln_rel_prob(n,j) END DO END IF ! q_ = or /= 1. END IF ! inside correlation horizon at 22 deg. END DO ! i = 1, stress_count IF (neighbors > 0) THEN sum = 0.0 DO j = 0, 59 probability(j) = EXP(probability(j)) sum = sum + probability(j) END DO probability = probability / sum peak = MAXLOC(probability) - 1 ! to compensate for (0:59), not (60)! azimuth = (peak(1) + 0.5) * 0.05236 jlo = peak(1) jro = jlo oldsum = 0. DO j = 1, 29 jl = MOD(peak(1) - j + 60, 60) jr = MOD(peak(1) + j, 60) sum = oldsum + 0.5 *(probability(jl) + probability(jlo) + & & probability(jr) + probability(jro)) IF (sum >= 0.90) THEN fraction = (0.90 - oldsum) / (sum - oldsum) del_az_for_90pc = 0.05236 * (j - 1. + fraction) RETURN END IF oldsum = sum jlo = jl jro = jr END DO ELSE azimuth = 0.0 END IF del_az_for_90pc = 1.41372 END SUBROUTINE Sigma_1h SUBROUTINE Solve_for_vw (passes, vw) INTEGER, INTENT(IN) :: passes REAL, DIMENSION(ndof), INTENT(INOUT) :: vw ! NOTE: vw is used BEFORE it is computed in the first pass ! through the stress section, IF (stress_now), and also to ! track convergence(?) of velocity solution during refinement. ! An array of zeros is OK, but undefined values are not! !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - TYPE(needle), DIMENSION(:), ALLOCATABLE :: needles !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REAL, DIMENSION(3,3) :: A, B, C, D REAL, DIMENSION(3) :: E, F REAL, DIMENSION(6,6) :: A6,B6,C6,D6 REAL, DIMENSION(6) :: E6,F6 REAL :: allowance, azimuth REAL :: big_one REAL :: boxed_frac REAL :: cosg, cos2g REAL :: cott, csct REAL, DIMENSION(3,2,2,2) :: dG REAL :: del_az_for_90pc ! one-sided, in radians REAL :: del_s_, difference, divisor REAL :: dV_frac COMPLEX, DIMENSION(3) :: eigenvalues COMPLEX, DIMENSION(3,3) :: eigenvectors REAL, DIMENSION(3) :: epsilon_dot REAL :: equat REAL :: error REAL :: eta_ REAL, DIMENSION(6) :: f_, g_ ! sometimes, only 1..3 are used REAL :: factor REAL :: floor = 2. * TINY(floor) REAL, DIMENSION(3,2,2) :: G REAL :: gamma_ TYPE(is123), DIMENSION(7) :: Gauss_point REAL, DIMENSION(7) :: Gauss_weight = (/ 0.225, & & 0.13239415, 0.13239415, 0.13239415, & & 0.12593918, 0.12593918, 0.12593918 /) REAL :: goal INTEGER :: h INTEGER :: i, i1, i2, i3, ita, ite INTEGER, PARAMETER :: ijob = 4 ! mode setting for LSLPB INTEGER :: j, jta, jte, k REAL :: kappa_ INTEGER :: l_, l_1, l_2 REAL, DIMENSION(3) :: lambda_ REAL, DIMENSION(3,3) :: Lambda INTEGER :: loc REAL :: lat, lon INTEGER :: m REAL :: mu_of_r_, mu_2 REAL :: new_V, num_boxed, num_stressed REAL :: overlap_threshold INTEGER :: p, pass LOGICAL :: plot_stress_at_nodes = .FALSE. REAL :: prefix REAL, DIMENSION(3) :: r_ REAL :: radius INTEGER :: s REAL :: s_ INTEGER :: segment CHARACTER(1) :: sense REAL :: sigma_ REAL :: sing, sin2g REAL :: sint INTEGER :: stress_count REAL :: stressed_frac REAL :: sum_base_n, sum_base_o, sum_diff, sum_err REAL :: t_sigma, tant, theta REAL :: top_value INTEGER :: twoi, u_ REAL, DIMENSION(3) :: tv, tv1, tv2 REAL :: vn, vo, wn, wo REAL, DIMENSION(3,3) :: V DOUBLE PRECISION, DIMENSION(3,3) :: V8 REAL :: W_over_R2 INTEGER :: Z, z_ INTEGER :: ABCDrow, ABCDcol, EFrow, EFcol ! statement functions ABCDrow(j) = j + ncoda ! statement function; ncoda is global ABCDcol(i, j) = (j - i) + 1 ! statement function EFrow(i) = i + ncoda ! statement function; ncoda is global EFcol = ncoda + 2 ! These statement functions are for ! codiagonal band symmetric storage mode of matrices ABCD and EF, ! per Microsoft version of IMSL. Note that element (row #i, col #j) ! of the idealized square matrix square_ABCD, or ! square_ABCD(i, j) is stored as ABCDEF(ABCDrow(j),ABCDcol(i, j)); ! and that only elements with j >= i (upper right) can be stored. ! Element (row i) of the linear vector EF is stored in ! ABCDEF(EFrow(i), EFcol). Gauss_point(1)%s(1:3) = (/ 0.33333333, 0.33333333, 0.33333333 /) Gauss_point(2)%s(1:3) = (/ 0.05971587, 0.47014206, 0.47014206 /) Gauss_point(3)%s(1:3) = (/ 0.47014206, 0.05971587, 0.47014206 /) Gauss_point(4)%s(1:3) = (/ 0.47014206, 0.47014206, 0.05971587 /) Gauss_point(5)%s(1:3) = (/ 0.79742698, 0.10128650, 0.10128650 /) Gauss_point(6)%s(1:3) = (/ 0.10128650, 0.79742698, 0.10128650 /) Gauss_point(7)%s(1:3) = (/ 0.10128650, 0.10128650, 0.79742698 /) W_over_R2 = W / R**2 ABCDEF = 0. ! global coefficient matrix AND forcing vector DO l_ = 1, num_ele i1 = node(1, l_) i2 = node(2, l_) i3 = node(3, l_) A = 0.; B = 0.; C = 0.; D = 0. ! element submatrices; see (4) ! through (9) of Bird (1997?). E = 0.; F = 0. ! forcing subvectors for element Z = crack_index(1, l_) IF (Z > 0) THEN ! active faulting! mu_of_r_ = mu_nod(i1) * Gauss_point(1)%s(1) + & & mu_nod(i2) * Gauss_point(1)%s(2) + & & mu_nod(i3) * Gauss_point(1)%s(3) t_sigma = mu_scale * ((mu_of_r_ / mu_scale)**exponent) mu_2 = t_sigma**2 V8(1,1) = mu_2 * 4./3. ! begin element covariance matrix V8(1,2) = 0.D0 V8(1,3) = mu_2 * (-2./3.) V8(2,1) = 0.D0 V8(2,2) = mu_2 V8(2,3) = 0.D0 V8(3,1) = V8(1,3) V8(3,2) = 0.D0 V8(3,3) = V8(1,1) epsilon_dot = 0. ! zero target strainrate vector ! evaluate nodal function and derivitives at center of element tv = center(1:3, l_) CALL Gjxy(l_, tv, G) CALL Del_Gjxy_del_thetaphi(l_, tv, dG) equat = SQRT(tv(1)**2 + tv(2)**2) IF (equat == 0.) THEN PRINT "(' Error: center of element ',I5,' is N or S pole.')", l_ WRITE (21,"('Error: center of element ',I5,' is N or S pole.')") l_ STOP END IF theta = ATAN2(equat, tv(3)) sint = SIN(theta) csct = 1. / sint tant = TAN(theta) cott = 1. / tant DO z_ = 1, Z ! loop on active traces in this element ! segment description loc = crack_index(2, l_) + z_ - 1 ! storage location in local_crack segment = local_crack(loc)%segment tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) gamma_ = Get_azimuth (tv1, tv2) sense = local_crack(loc)%sense s_ =local_crack(loc)%s_ del_s_ = f_scale * ((local_crack(loc)%sigma_ / f_scale)**exponent) u_ = seg_u_(segment) ! 1, 2, or 3? eta_ = seg_eta_(segment) kappa_ = seg_kappa_(segment) ! H vector of this trace prefix = eta_ * kappa_ / R sing = SIN(gamma_) cosg = COS(gamma_) IF ((sense == 'L') .OR. (sense == 'R')) THEN local_crack(loc)%H(1) = prefix * (dG(u_,1,1,1) * cosg - dG(u_,2,1,1) * sing) local_crack(loc)%H(2) = prefix * 0.5 * (dG(u_,1,1,2) * cosg/sint - dG(u_,2,1,2) * sing/sint + & & dG(u_,1,2,1) * cosg - dG(u_,2,2,1) * sing - & & cott * (G(u_,1,2) * cosg - G(u_,2,2) * sing)) local_crack(loc)%H(3) = prefix * (dG(u_,1,2,2) * cosg/sint - dG(u_,2,2,2) * sing/sint + & & cott * (G(u_,1,1) * cosg - G(u_,2,1) * sing)) ELSE IF ((sense == 'T') .OR. (sense == 'P') .OR. (sense == 'N') .OR. (sense == 'D')) THEN local_crack(loc)%H(1) = prefix * (dG(u_,1,1,1) * sing + dG(u_,2,1,1) * cosg) local_crack(loc)%H(2) = prefix * 0.5 * (dG(u_,1,1,2) * sing/sint + dG(u_,2,1,2) * cosg/sint + & & dG(u_,1,2,1) * sing + dG(u_,2,2,1) * cosg - & & cott * (G(u_,1,2) * sing + G(u_,2,2) * cosg)) local_crack(loc)%H(3) = prefix * (dG(u_,1,2,2) * sing/sint + dG(u_,2,2,2) * cosg/sint + & & cott * (G(u_,1,1) * sing + G(u_,2,1) * cosg)) ELSE ! should not happen CALL Prevent ('bad slip sense', 65, 'Solve for w') END IF ! increment strain-rate goals and variances DO i = 1, 3 epsilon_dot(i) = epsilon_dot(i) + s_ * local_crack(loc)%H(i) DO j = 1, 3 V8(i,j) = V8(i,j) + del_s_**2 * local_crack(loc)%H(i) * & & local_crack(loc)%H(j) END DO END DO END DO ! z_ = 1, Z ! scale V matrix to prevent underflows V(1:3,1:3) = REAL(V8(1:3,1:3) * 1.D30) ! eigen-analysis of V: CALL EVCRG (3, V, 3, eigenvalues, eigenvectors, 3) ! Compute all of the eigenvalues and eigenvectors of a real matrix. ! Usage: ! CALL EVCRG (N, A, LDA, EVAL, EVEC, LDEVEC) ! Arguments: ! N = Order of the matrix. (Input) ! A = Floating-point array containing the matrix. (Input) ! LDA = Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input) ! EVAL = Complex array of size N containing the eigenvalues of A in decreasing order of magnitude. (Output) ! EVEC = Complex array containing the matrix of eigenvectors. (Output) ! The #h eigenvector, corresponding to EVAL(h), is stored in the #h column. Each vector is normalized to have ! Euclidean length equal to the value one. ! LDEVEC = Leading dimension of EVEC exactly as specified in the dimension statement of the calling program. (Input) top_value = 0. DO h = 1, 3 ! the eigenvectors each define a new scalar datum ! undo the scaling applied previously to V lambda_(h) = 1.E-30 * REAL(eigenvalues(h)) top_value = MAX(top_value, lambda_(h)) END DO DO h = 1,3 IF (lambda_(h) <= 0.) THEN ! prevent needless halts due to loss-of-precision IF (ABS(lambda_(h)) <= (0.01 * top_value)) THEN lambda_(h) = mu_2 ELSE CALL Prevent('nonpositive eigenvalue',185,'Solve for vw') END IF ! large or small negative eigenvalue END IF ! nonpositive eigenvalue (shouldn't happen) sigma_ = SQRT(lambda_(h)) prefix = 1. / sigma_**2 goal = 0. DO m = 1, 3 Lambda(h, m) = REAL(eigenvectors(m, h)) ! sic transpose; see above goal = goal + Lambda(h, m) * epsilon_dot(m) END DO DO j = 1, 3 ! local node numbers f_(j) = (1./R) * (dG(j,1,1,1) * Lambda(h, 1) + & & 0.5 * (csct * dG(j,1,1,2) + dG(j,1,2,1) - cott * G(j,1,2)) * Lambda(h, 2) + & & (csct * dG(j,1,2,2) + cott * G(j,1,1)) * Lambda(h, 3)) g_(j) = (1./R) * (dG(j,2,1,1) * Lambda(h, 1) + & & 0.5 * (csct * dG(j,2,1,2) + dG(j,2,2,1) - cott * G(j,2,2)) * Lambda(h, 2) + & & (csct * dG(j,2,2,2) + cott * G(j,2,1)) * Lambda(h, 3)) END DO ! j = 1, 2, 3 (local node numbers) CALL Add_datum(prefix,f_,g_,goal,A,C,D,E,F) END DO ! eigenvalues h = 1, 2, 3 :: 3 more data for linear system ELSE ! use a-priori constraint ! Basic stiffness matrix of nonfaulting elements DO m = 1, 7 mu_of_r_ = mu_nod(i1) * Gauss_point(m)%s(1) + & & mu_nod(i2) * Gauss_point(m)%s(2) + & & mu_nod(i3) * Gauss_point(m)%s(3) t_sigma = mu_scale * ((mu_of_r_ / mu_scale)**exponent) prefix = W_over_R2 * a_(l_) * Gauss_weight(m) / (t_sigma**2) Gauss_point(m)%element = l_ CALL Interpolate(Gauss_point(m), r_) CALL Gjxy(l_, r_, G) CALL Del_Gjxy_del_thetaphi(l_, r_, dG) equat = SQRT(r_(1)**2 + r_(2)**2) IF (equat == 0.) THEN PRINT "(' Error: integration point ',I1,' of element ',I5,' is N or S pole.')", m, l_ WRITE (21,"('Error: integration point ',I1,' of element ',I5,' is N or S pole.')") m, l_ STOP END IF theta = ATAN2(equat, r_(3)) csct = 1. / SIN(theta) tant = TAN(theta) cott = 1. / tant DO i = 1, 3 DO j = 1, 3 IF (j <= i) THEN ! diagonal and lower triangle only A(i,j) = A(i,j) + prefix * & & ( 2.*dG(i,1,1,1)*dG(j,1,1,1)+ & & csct*(dG(i,1,1,1)*dG(j,1,2,2)+dG(i,1,2,2)*dG(j,1,1,1))+ & & cott*(dG(i,1,1,1)*G(j,1,1)+G(i,1,1)*dG(j,1,1,1))+ & & 2.*(csct*dG(i,1,2,2)+G(i,1,1)*cott)*(csct*dG(j,1,2,2)+G(j,1,1)*cott)+ & & 0.5*(csct*dG(i,1,1,2)+dG(i,1,2,1)-G(i,1,2)*cott)*(csct*dG(j,1,1,2)+dG(j,1,2,1)-G(j,1,2)*cott) ) D(i,j) = D(i,j) + prefix * & & ( 2.*dG(i,2,1,1)*dG(j,2,1,1)+ & & csct*(dG(i,2,1,1)*dG(j,2,2,2)+dG(i,2,2,2)*dG(j,2,1,1))+ & & cott*(dG(i,2,1,1)*G(j,2,1)+G(i,2,1)*dG(j,2,1,1))+ & & 2.*(csct*dG(i,2,2,2)+G(i,2,1)*cott)*(csct*dG(j,2,2,2)+G(j,2,1)*cott)+ & & 0.5*(csct*dG(i,2,1,2)+dG(i,2,2,1)-G(i,2,2)*cott)*(csct*dG(j,2,1,2)+dG(j,2,2,1)-G(j,2,2)*cott) ) END IF ! All of B lies in the upper triangle; transpose of C ! B(i,j) = B(i,j) + prefix * & ! & ( 2.*dG(i,1,1,1)*dG(j,2,1,1)+ & ! & csct*(dG(i,1,1,1)*dG(j,2,2,2)+dG(i,1,2,2)*dG(j,2,1,1))+ & ! & cott*(dG(i,1,1,1)*G(j,2,1)+G(i,1,1)*dG(j,2,1,1))+ & ! & 2.*(csct*dG(i,1,2,2)+G(i,1,1)*cott)*(csct*dG(j,2,2,2)+G(j,2,1)*cott)+ & ! & 0.5*(csct*dG(i,1,1,2)+dG(i,1,2,1)-G(i,1,2)*cott)*(csct*dG(j,2,1,2)+dG(j,2,2,1)-G(j,2,2)*cott) ) ! All of C lies in lower triangle C(i,j) = C(i,j) + prefix * & & ( 2.*dG(i,2,1,1)*dG(j,1,1,1)+ & & csct*(dG(i,2,1,1)*dG(j,1,2,2)+dG(i,2,2,2)*dG(j,1,1,1))+ & & cott*(dG(i,2,1,1)*G(j,1,1)+G(i,2,1)*dG(j,1,1,1))+ & & 2.*(csct*dG(i,2,2,2)+G(i,2,1)*cott)*(csct*dG(j,1,2,2)+G(j,1,1)*cott)+ & & 0.5*(csct*dG(i,2,1,2)+dG(i,2,2,1)-G(i,2,2)*cott)*(csct*dG(j,1,1,2)+dG(j,1,2,1)-G(j,1,2)*cott) ) END DO ! on j = 1, 3 END DO ! on i = 1, 3 END DO ! on 7 Gauss integration points END IF ! faults, OR a-priori stiffness in this element ! Add any paleomagnetic data in this element IF (p_dat_count > 0) THEN DO p = 1, p_dat_count IF (p_active(n_,p)) THEN IF (p_site_is(p)%element == l_) THEN r_ = p_site_now(1:3,p) CALL Gjxy(l_, r_, G) ! paleolatitude anomaly part tv = p_pole(1:3, p) gamma_ = Get_azimuth(r_,tv) sing = SIN(gamma_) cosg = COS(gamma_) DO j = 1, 3 f_(j) = G(j,1,1) * cosg - G(j,1,2) * sing g_(j) = G(j,2,1) * cosg - G(j,2,2) * sing END DO t_sigma = p_drift_scale * ((p_south_rate_sigma_(p) / p_drift_scale)**exponent) prefix = P_weight / (t_sigma**2) !GPB CALL Add_datum(prefix,f_,g_,p_south_goal(n_,p),A,C,D,E,F) ! vertical-axis rotation part IF (.NOT.twisted(p)) THEN CALL Del_Gjxy_del_thetaphi(l_, r_, dG) equat = SQRT(r_(1)**2 + r_(2)**2) theta = ATAN2(equat, r_(3)) sint = SIN(theta) csct = 1. / sint tant = TAN(theta) cott = 1. / tant prefix = 1. / (2. * R) DO j = 1, 3 f_(j) = prefix * (G(j,1,2)*cott + dG(j,1,2,1) - csct*dG(j,1,1,2)) g_(j) = prefix * (G(j,2,2)*cott + dG(j,2,2,1) - csct*dG(j,2,1,2)) END DO t_sigma = p_spin_scale * ((p_ccw_rate_sigma_(p) / p_spin_scale)**exponent) prefix = P_weight / (t_sigma**2) !GPB CALL Add_datum(prefix,f_,g_,p_ccw_goal(n_,p),A,C,D,E,F) END IF ! .NOT.twisted END IF ! paleomag site is in this element END IF ! p_active(n_,p) END DO ! p = 1, p_dat_count END IF ! p_dat_count > 0 CALL Plug_in_33 (l_,A,B,C,D,E,F) ! add element matrix,vector to global END DO ! l_ = 1,num_ele, doing a-priori stiffness & faults & paleomag ! Add any active cross-sections with both ends in model IF (c_dat_count > 0) THEN DO k = 1, c_dat_count IF (c_active(n_,k)) THEN IF ((c_end_is(1,k)%element > 0).AND. & &(c_end_is(2,k)%element > 0)) THEN A6 = 0.; B6 = 0.; C6 = 0.; D6 = 0. E6 = 0.; F6 = 0. l_1 = c_end_is(1,k)%element ! West end is in this element tv1 = c_end_now(1:3, 1, k) tv2 = c_end_now(1:3, 2, k) gamma_ = Get_azimuth(tv1, tv2) sing = SIN(gamma_) cosg = COS(gamma_) CALL Gjxy(l_1, tv1, G) DO j = 1,3 f_(j) = G(j,1,1) * cosg - G(j,1,2) * sing g_(j) = G(j,2,1) * cosg - G(j,2,2) * sing END DO l_2 = c_end_is(2,k)%element ! East end is in this element gamma_ = Pi + Get_azimuth(tv2, tv1) ! gamma_ is almost same as above, but evaluated at East end now sing = SIN(gamma_) cosg = COS(gamma_) CALL Gjxy(l_2, tv2, G) DO j = 1,3 f_(j+3) = -G(j,1,1) * cosg + G(j,1,2) * sing g_(j+3) = -G(j,2,1) * cosg + G(j,2,2) * sing END DO t_sigma = c_scale * ((c_rate_sigma_(k) / c_scale)**exponent) prefix = 1./(t_sigma**2) DO i = 1, 6 DO j = 1, 6 IF (j <= i) THEN ! diagonal and lower triangle only A6(i,j) = A6(i,j) + prefix * f_(i) * f_(j) D6(i,j) = D6(i,j) + prefix * g_(i) * g_(j) END IF ! All of B lies in the upper triangle; transpose of C ! B(i,j) = B(i,j) + prefix * f_(i) *g_(j) ! All of C lies in lower triangle C6(i,j) = C6(i,j) + prefix * g_(i) * f_(j) END DO ! on j = 1, 6 E6(i) = E6(i) + prefix * f_(i) * c_goal(n_, k) F6(i) = F6(i) + prefix * g_(i) * c_goal(n_, k) END DO ! on i = 1, 6 CALL Plug_in_66 (l_1,l_2,A6,B6,C6,D6,E6,F6) END IF ! this cross-section is within domain END IF ! this cross-section is active now END DO ! k = 1, c_dat_count END IF ! c_dat_count > 0 IF (passes > 0) THEN ! selection criterion could be > 1 (0R, > 0 to always print) ! Write headers for convergence report IF (stress_now) THEN PRINT "(' ',18X,'Refinement dV/V RMS(V) Stressed Boxed Mean_error(deg.)')" WRITE (21,"(19X,'Refinement dV/V RMS(V) Stressed Boxed Mean_error(deg.)')") ELSE PRINT "(' ',18X,'Refinement dV/V RMS(V)')" WRITE (21,"(19X,'Refinement dV/V RMS(V)')") END IF END IF many_passes: DO pass = 1, passes ! Shuffle copies of matrix, if needed IF (passes > 1) THEN IF (pass == 1) THEN duplicate = ABCDEF ! save constant part ELSE ABCDEF = duplicate ! retrieve constant part END IF END IF num_boxed = 0.; num_stressed = 0. IF (stress_now) THEN ! Initialize stress in each element IF (pass == 1) THEN !PRINT "(' ',12X,'Interpolating stress directions (slow!)')" ! Count the stress data (including cracks?) stress_count = 0 DO s = 1, s_dat_count IF (s_site_is(1,s)%element > 0) THEN IF (s_activity(n_,s) > 0.) stress_count = stress_count + 1 END IF END DO IF (faults_give_sigma_1h) THEN DO s = 1, crack_count i = local_crack(s)%datum IF (f_goal(n_,i) /= 0.) THEN IF (f_new(i)) stress_count = stress_count + 1 END IF END DO END IF ! Allocate the array to hold the data ALLOCATE ( needles(stress_count) ) ! Fill the array with data stress_count = 0 DO s = 1, s_dat_count IF (s_site_is(1,s)%element > 0) THEN IF (s_activity(n_,s) > 0.) THEN stress_count = stress_count + 1 needles(stress_count)%location = s_site_now(1:3,1,s) needles(stress_count)%azimuth = s_azim_now(s) needles(stress_count)%sigma = s_sigma_(s) needles(stress_count)%relevance = s_activity(n_,s) END IF END IF END DO IF (faults_give_sigma_1h) THEN DO s = 1, crack_count i = local_crack(s)%datum IF (f_new(i)) THEN IF (f_goal(n_,i) /= 0.) THEN stress_count = stress_count + 1 segment = local_crack(s)%segment tv = 0.5 * (seg_end(1:3,1,segment) + seg_end(1:3,2,segment)) CALL Unitise(tv, r_) needles(stress_count)%location = r_ tv1 = seg_end(1:3, 1, segment) tv2 = seg_end(1:3, 2, segment) gamma_ = Get_azimuth (tv1, tv2) sense = local_crack(s)%sense ! T, N, R, L, D IF ((sense == 'T') .OR. (sense == 'P')) THEN gamma_ = gamma_ + 1.5708 ! Pi / 2 ELSE IF (sense == 'R') THEN gamma_ = gamma_ + 0.7854 ! Pi / 4 ELSE IF (sense == 'L') THEN gamma_ = gamma_ - 0.7854 ! Pi / 4 END IF needles(stress_count)%azimuth = gamma_ needles(stress_count)%sigma = 0.3928 ! Pi / 8 !Treat as stage data: IF (paleotec) THEN overlap = MAX(0.0, MIN(time1 - time0, f_t_max(i) - f_t_min(i), & & time1 - f_t_min(i), f_t_max(i) - time0)) overlap_threshold = MIN(0.1 * s_per_Ma, 0.5 * Deltat_) IF (overlap >= overlap_threshold) THEN needles(stress_count)%relevance = 1.0 ELSE needles(stress_count)%relevance = 0.0 END IF ELSE ! neotec allowance = 0.1 * s_per_Ma IF ((start_time > (f_t_min(i) - allowance)).AND. & & (start_time < (f_t_max(i) + allowance))) THEN needles(stress_count)%relevance = 1.0 ELSE needles(stress_count)%relevance = 0.0 END IF END IF ! paleotec / neotec END IF ! f_goal(n_,i) /= 0. END IF ! new fault END DO ! all local cracks END IF ! faults_give_sigma_1h ! Interpolate stress to centers of elements with no (real) data ele_stressed = .TRUE. ! (just initializing) ele_q = 0. ! (ditto) ! No interpolation if element contains a datum; use best datum DO s = 1, s_dat_count l_ = s_site_is(1,s)%element IF (l_ > 0) THEN IF (s_activity(n_,s) > ele_q(l_)) THEN IF (ele_stressed(l_)) THEN ele_azim(l_) = s_azim_now(s) ele_q(l_) = s_activity(n_,s) ele_sigma(l_) = s_sigma_(s) END IF END IF END IF END DO IF ((plot_stress_at_nodes).AND.(passes > 1)) THEN PRINT "(' Writing Sigma_1h.feg to show stress at nodes; read with OrbWeave')" WRITE (21, "('Writing Sigma_1h.feg to show stress at nodes; read with OrbWeave')") OPEN (1, FILE = 'Sigma_1h.feg') WRITE (1, "('Use with OrbWeave to see stress at node locations')") WRITE (1,"(I5, I5,' 0 30000 T')") num_nod, num_nod DO i = 1, num_nod r_ = xyz_nod(1:3, i) CALL Sigma_1h(stress_count,needles,r_,azimuth,del_az_for_90pc) CALL Lonlat_from_xyz (r_, lon, lat) WRITE (1, "(I5,2F10.4,2F10.2,' 0.0 0.0')") i, lon, lat, azimuth*deg_per_rad, del_az_for_90pc*deg_per_rad END DO WRITE (1, "(I5)") num_ele DO i = 1, num_ele WRITE (1, "(4I5)") i, node(1,i),node(2,i),node(3,i) END DO WRITE (1, "(' 0')") CLOSE (1) END IF ! Interpolate to centers of remaining elements DO l_ = 1, num_ele IF (ele_q(l_) == 0.) THEN ! has no datum r_ = center(1:3,l_) CALL Sigma_1h(stress_count,needles,r_,azimuth,del_az_for_90pc) ele_azim(l_) = azimuth ele_sigma(l_) = del_az_for_90pc * 0.6079 ! (assumes Gaussian shape!) IF (del_az_for_90pc <= 0.7854) THEN ! +- 45 degrees @ 90%-confidence is cutoff ele_q(l_) = 1.00 ! (because low q was translated to increased del_az_for_90pc ! during the interpolation process ELSE ! too uncertain to bother ele_stressed(l_) = .FALSE. END IF END IF ! no datum in this element END DO ! l_ = 1, num_ele END IF ! pass == 1 (stress needs interpolating) sum_err = 0. ! total of angular errors, in radians ! Insert stress constraints (different in each refinement) DO l_ = 1, num_ele IF (ele_stressed(l_)) THEN num_stressed = num_stressed + ele_q(l_) A = 0.; B = 0.; C = 0.; D = 0. E = 0.; F = 0. r_ = center(1:3, l_) CALL Gjxy(l_, r_, G) CALL Del_Gjxy_del_thetaphi(l_, r_, dG) equat = SQRT(r_(1)**2 + r_(2)**2) theta = ATAN2(equat, r_(3)) csct = 1. / SIN(theta) tant = TAN(theta) cott = 1. / tant cos2g = COS(2. * ele_azim(l_)) sin2g = SIN(2. * ele_azim(l_)) prefix = 1. / (2. * R) DO j = 1, 3 f_(j) = prefix * ( (csct*dG(j,1,1,2) + dG(j,1,2,1) - cott*G(j,1,2))*cos2g + & & (dG(j,1,1,1) - csct*dG(j,1,2,2) - cott*G(j,1,1))*sin2g ) g_(j) = prefix * ( (csct*dG(j,2,1,2) + dG(j,2,2,1) - cott*G(j,2,2))*cos2g + & & (dG(j,2,1,1) - csct*dG(j,2,2,2) - cott*G(j,2,1))*sin2g ) END DO CALL E_rate(l_, G, dG, theta, vw, epsilon_dot) ! epsilon_alpha_beta = 0 (only if element not faulting!) IF (crack_index(1, l_) == 0) THEN radius = SQRT(epsilon_dot(2)**2 + (0.5*(epsilon_dot(1)-epsilon_dot(3)))**2) radius = MAX(radius, xi_) ! global parameter sigma_ = 2. * ele_sigma(l_) * radius sigma_ = MAX(sigma_, floor) prefix = ele_q(l_) / (sigma_**2) CALL Add_datum(prefix,f_,g_,0.,A,C,D,E,F) END IF ! element is not faulting ! assess angular errors in sigma_1h azimuth = 0.5 * ATAN2F(epsilon_dot(2),0.5*(epsilon_dot(3)-epsilon_dot(1))) error = azimuth - ele_azim(l_) error = MIN(MOD(error + 6.2832, 3.1416), MOD(-error + 6.2832, 3.1416)) sum_err = sum_err + error * ele_q(l_) ! epsilon_alpha_alpha < epsilon_beta_beta IF (((epsilon_dot(3) - epsilon_dot(1))*cos2g + 2. * epsilon_dot(2) * sin2g) < 0.) boxed(l_) = .TRUE. IF (boxed(l_)) THEN ! Needs constraint to enforce correct sense num_boxed = num_boxed + ele_q(l_) prefix = 1. / R DO j = 1, 3 f_(j) = prefix * ((csct*dG(j,1,2,2)+G(j,1,1)*cott-dG(j,1,1,1)) * cos2g + & & (csct*dG(j,1,1,2)+dG(j,1,2,1)-G(j,1,2)*cott) * sin2g) g_(j) = prefix * ((csct*dG(j,2,2,2)+G(j,2,1)*cott-dG(j,2,1,1)) * cos2g + & & (csct*dG(j,2,1,2)+dG(j,2,2,1)-G(j,2,2)*cott) * sin2g) END DO prefix = ele_q(l_) / (0.83 * xi_)**2 CALL Add_datum(prefix,f_,g_,xi_,A,C,D,E,F) END IF ! boxed(l_): sense constraint needed CALL Plug_in_33 (l_,A,B,C,D,E,F) ! add element matrix,vector to global END IF ! ele_stressed(l_) END DO ! l_ = 1, num_ele END IF ! stress_now ! Scale linear system (to prevent overflow/underflow) big_one = MAXVAL(ABCDEF) ! biggest should be in ABCD part, not EF part factor = 1./big_one ABCDEF = factor * ABCDEF ! scales both coefficients and rhs ! Impose boundary conditions, maintaining symmetry, but only ! reading or changing upper right (j >= i). DO m = 1, bcs_count ! South (theta) component of velocity k = 2 * boundary_node(m) - 1 ! logical row/col to be replaced j = k ! work on the column DO i = MAX(1, k - ncoda), k-1 ! upper part of column ! clear column, changing rhs ita = ABCDrow(j) jta = ABCDcol(i, j) ite = EFrow(i) jte = EFcol ABCDEF(ite, jte) = ABCDEF(ite, jte) - condition(1, m) * ABCDEF(ita, jta) ABCDEF(ita, jta) = 0. END DO i = k ! work on the row DO j = k+1, MIN(ndof, k + ncoda) ! right part of row ! clear row, changing rhs ita = ABCDrow(j) jta = ABCDcol(i, j) ite = EFrow(j) jte = EFcol ABCDEF(ite, jte) = ABCDEF(ite, jte) - condition(1, m) * ABCDEF(ita, jta) ABCDEF(ita, jta) = 0. END DO ABCDEF(ABCDrow(k), ABCDcol(k, k)) = 1. ! put 1 on diagonal ABCDEF(EFrow(k), EFcol) = condition(1, m) ! put BC velocity component in rhs ! East (phi) component of velocity k = 2 * boundary_node(m) ! logical row/col to be replaced j = k ! work on the column DO i = MAX(1, k - ncoda), k-1 ! upper part of column ! clear column, changing rhs ita = ABCDrow(j) jta = ABCDcol(i, j) ite = EFrow(i) jte = EFcol ABCDEF(ite, jte) = ABCDEF(ite, jte) - condition(2, m) * ABCDEF(ita, jta) ABCDEF(ita, jta) = 0. END DO i = k ! work on the row DO j = k+1, MIN(ndof, k + ncoda) ! right part of row ! clear row, changing rhs ita = ABCDrow(j) jta = ABCDcol(i, j) ite = EFrow(j) jte = EFcol ABCDEF(ite, jte) = ABCDEF(ite, jte) - condition(2, m) * ABCDEF(ita, jta) ABCDEF(ita, jta) = 0. END DO ABCDEF(ABCDrow(k), ABCDcol(k, k)) = 1. ! put 1 on diagonal ABCDEF(EFrow(k), EFcol) = condition(2, m) ! put BC velocity component in rhs END DO ! on boundary node index m CALL LSLPB (ndof, ABCDEF, lda, ncoda, ijob, u_flag) ! Usage: ! CALL LSLPB (N, A, LDA, NCODA, IJOB, U) ! Arguments: ! N = Order of the matrix. (Input) ! Must satisfy N > 0. ! A = Array containing the N by N positive definite band coefficient ! matrix and right hand side in MS-IMSL's ! codiagonal band symmetric storage mode. (Input/Output) ! The number of array columns must be at least NCODA + 2. ! The number of columns is not an input to this subprogram. ! LDA = Leading dimension of A exactly as specified in the ! dimension statement of the calling program. (Input) ! Must satisfy LDA >= N + NCODA. ! NCODA = Number of upper codiagonals of matrix A. (Input) ! Must satisfy NCODA >= 0 and NCODA < N. ! IJOB = Flag to direct the desired factorization or solving step. (Input) ! IJOB Meaning: ! 1 factor the matrix A and solve the system Ax = b, where b is stored in column NCODA + 2 of array A. ! The vector x overwrites b in storage. ! 2 solve step only. Use b as column NCODA + 2 of A. (The factorization step has already been done.) ! The vector x overwrites b in storage. ! 3 factor the matrix A but do not solve a system. ! 4,5,6 same meaning as with the value IJOB - 3. For efficiency, no ! error checking is done on values LDA, N, NCODA, and U(*). ! U = Array of flags that indicate any singularities of A, namely loss of positive-definiteness of a leading minor. (Output) ! A value U(I) = 0. means that the leading minor of dimension I is not positive-definite. Otherwise, U(I) = 1. ! Comments: ! Automatic workspace usage is: NCODA real numbers. IF (passes > 0) THEN ! Compare new with old solution sum_base_o = 0. sum_base_n = 0. sum_diff = 0. DO i = 1, num_nod twoi = 2 * i wo = vw(twoi) wn = ABCDEF(EFrow(twoi), EFcol) twoi = twoi - 1 vo = vw(twoi) vn = ABCDEF(EFrow(twoi), EFcol) sum_base_o = sum_base_o + vo**2 +wo**2 sum_base_n = sum_base_n + vn**2 +wn**2 sum_diff = sum_diff + (vn - vo)**2 + (wn - wo)**2 END DO new_V = SQRT(sum_base_n/num_nod) divisor = SQRT(MAX(sum_base_o,sum_base_n)/num_nod) difference = SQRT(sum_diff/num_nod) IF (divisor > 0.) THEN dV_frac = difference / divisor ELSE dV_frac = 0. END IF IF (stress_now) THEN stressed_frac = (100. * num_stressed) / num_ele boxed_frac = (100. * num_boxed) / num_ele IF (num_stressed > 0.) THEN error = deg_per_rad * sum_err / num_stressed ELSE error = 0. END IF PRINT "(' ',18X,I10,F8.5,1P,E9.2,0P,F8.2,'%',F7.2,'%',F11.2)", & & pass-1, dV_frac, new_V, stressed_frac, boxed_frac, error WRITE (21,"(19X,I10,F8.5,1P,E9.2,0P,F8.2,'%',F7.2,'%',F11.2)") & & pass-1, dV_frac, new_V, stressed_frac, boxed_frac, error ELSE PRINT "(' ',18X,I10,F8.5,1P,E9.2)", & & pass-1, dV_frac, new_V WRITE (21,"(19X,I10,F8.5,1P,E9.2)") & & pass-1, dV_frac, new_V END IF ! stress_now END IF ! passes > 1 (or 0?) ! Transfer solution to velocity supervector DO i = 1, ndof vw(i) = ABCDEF(EFrow(i), EFcol) END DO IF ((passes > 1).AND.(dV_frac < 0.00001)) EXIT many_passes END DO many_passes ! pass = 1, passes IF (stress_now) DEALLOCATE ( needles ) END SUBROUTINE Solve_for_vw SUBROUTINE Step_aside (b_, gamma_, new_vec) REAL, DIMENSION(3), INTENT(IN) :: b_ REAL, INTENT(IN) :: gamma_ REAL, DIMENSION(3), INTENT(OUT) :: new_vec ! gives position (Cartesian unit vector) close to b_, ! but displaced toward gamma_ (in radians, clockwise from N) REAL, DIMENSION(3) :: offset, Phi, Theta, v1 REAL :: radians = 0.002 ! offset ~13 km on Earth CALL Local_Phi (b_, Phi) CALL Local_Theta(b_, Theta) offset = Phi * SIN(gamma_) - Theta * COS(gamma_) v1 = b_ + (offset * radians) CALL Unitise(v1, new_vec) END SUBROUTINE Step_aside SUBROUTINE Test_file (name, unit) ! Tests existing(?) file for validity by briefly opening and then closing CHARACTER(*) :: name INTEGER :: unit, open_status OPEN (unit, file = name, STATUS = 'OLD', IOSTAT = open_status) IF (open_status == 0) THEN CLOSE (unit) RETURN ELSE PRINT "(' Error: Following filename is invalid:'/' ',A)", TRIM(name) WRITE (21,"('Error: Following filename is invalid:'/A)") TRIM(name) STOP END IF END SUBROUTINE Test_file SUBROUTINE Unitise (b_, unit_vec) REAL, DIMENSION(3), INTENT(IN) :: b_ REAL, DIMENSION(3), INTENT(OUT) :: unit_vec REAL :: length length = Magnitude (b_) IF (length /= 0.) THEN unit_vec(1) = b_(1)/length unit_vec(2) = b_(2)/length unit_vec(3) = b_(3)/length ELSE unit_vec(1) = 1. unit_vec(2) = 0. unit_vec(3) = 0. END IF END SUBROUTINE Unitise INTEGER FUNCTION Which_zero(is) ! Returns 1, 2, or 3 to identify which si is closest to zero. TYPE(is123), INTENT(IN) :: is INTEGER, DIMENSION(1) :: j REAL, DIMENSION(3) :: sabs sabs(1:3) = ABS(is%s(1:3)) j = MINLOC(sabs) Which_zero = j(1) END FUNCTION Which_zero SUBROUTINE Write_c_dat (shift) ! Writes cmmnn.dat output file with +, * lines. ! Note that X-sections (even partly) outside the grid are dropped. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt CHARACTER(80) :: filename INTEGER :: unit = 27 ! see comment lines at top of file INTEGER :: i REAL, DIMENSION(3) :: tv1, tv2 filename = Insert (c_dat, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') WRITE (unit, "(A)") c_dat_format WRITE (unit, "(A)") c_dat_titles DO i = 1, c_dat_count IF ((c_end_is(1, i)%element > 0).AND.(c_end_is(2, i)%element > 0)) THEN c47 = c_ref(i) tv1 = c_end_0(1:3, 1, i) CALL Lonlat_from_xyz (tv1, x1, x2) tv2 = c_end_0(1:3, 2, i) CALL Lonlat_from_xyz (tv2, x3, x4) c5 = c_code(i) r1 = c_length(1,i) / m_per_km r2 = c_length(2,i) / m_per_km r3 = c_sigma_(i) / m_per_km t2 = c_t_max(i) / s_per_Ma t1 = c_t_min(i) / s_per_Ma WRITE (unit, c_dat_format) c47, x1, x2, x3, x4, c5, r1, r2, r3, t2, t1 IF (paleotec .OR. (start_time > 0.)) THEN tv1 = c_end_now(1:3, 1, i) CALL Lonlat_from_xyz (tv1, x1, x2) WRITE (unit, "('+',1X,F8.3,1X,F7.3)") x1, x2 tv2 = c_end_now(1:3, 2, i) CALL Lonlat_from_xyz (tv2, x1, x2) WRITE (unit, "('+',1X,F8.3,1X,F7.3)") x1, x2 END IF CALL Write_rates_and_goals & & (index = i, unit = unit, signal = '*', & & conversion = (s_per_Ma / m_per_km), & & goal = c_goal, rate = c_rate, active = c_active) END IF ! section is in domain END DO ! on all sections CLOSE (unit) END SUBROUTINE Write_c_dat SUBROUTINE Write_f_dat (shift) ! Writes fmmnn.dat output file with * lines. ! Note that faults which do not have at least 2 points within the grid ! are dropped. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt CHARACTER(80) :: filename CHARACTER(6) :: c6 INTEGER :: unit = 26 ! see comment lines at top of file INTEGER :: i, j filename = Insert (f_dat, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') WRITE (unit, "(A)") f_dat_format WRITE (unit, "(A)") f_dat_titles DO i = 1, f_dat_count IF (f_2_in(which_trace(i))) THEN WRITE (c6, "(1X,I4,1X)") which_trace(i) !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. DO j = 2, 5 IF (c6(j:j) == ' ') c6(j:j) = '0' END DO c6(1:1) = 'F' c6(6:6) = sense(i) c50 = fault_name(i) t1 = offset(i) / m_per_km t2 = offset_sigma_(i) / m_per_km t3 = f_t_max(i) / s_per_Ma t4 = f_t_min(i) / s_per_Ma WRITE (unit, f_dat_format) c6, c50, t1, t2, t3, t4 CALL Write_rates_and_goals & & (index = i, unit = unit, signal = '*', & & conversion = (s_per_Ma / m_per_km), & & goal = f_goal, rate = f_rate, active = f_active) END IF END DO CLOSE (unit) END SUBROUTINE Write_f_dat SUBROUTINE Write_f_dig (shift) ! Writes fmmnn.dig output file. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt CHARACTER(80) :: filename CHARACTER(4) :: c4 INTEGER :: unit = 22 ! see comment lines at top of file INTEGER :: i, j, j1, j2, n REAL :: lon, lat REAL, DIMENSION(3) :: tv filename = Insert (f_dig, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') DO i = 1, f_highest j1 = trace_loc(1,i) j2 = trace_loc(2,i) IF (j2 > j1) THEN n = 0 ! count of trace points inside the grid DO j = j1, j2 IF (trace_is(j)%element > 0) n = n + 1 END DO IF (n >= 2) THEN WRITE (c4,"(I4)") i !BUG: Formatted internal WRITE causes memory leak ! under Microsoft Fortran Powerstation 4.0, ! but it will be unimportant in this case. IF (c4(1:1) == ' ') c4(1:1) = '0' IF (c4(2:2) == ' ') c4(2:2) = '0' IF (c4(3:3) == ' ') c4(3:3) = '0' IF (c4(4:4) == ' ') c4(4:4) = '0' WRITE (unit,"('F',A,A)") c4, trace_type(i) DO j = j1, j2 IF (trace_is(j)%element > 0) THEN tv = trace(1:3, j) CALL Lonlat_from_xyz (tv , lon, lat) WRITE (unit, "(1X,SP,1P,E12.5,',',E12.5)") lon, lat END IF END DO WRITE (unit,"('*** END OF SEGMENT ***')") END IF END IF END DO CLOSE (unit) END SUBROUTINE Write_f_dig SUBROUTINE Write_p_dat (shift) ! Writes pmmnn.dat output file with +, *, & lines. ! Note that points outside the grid are dropped, ! and twisted points have no vertical-axis rotation reported. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt CHARACTER(80) :: filename INTEGER :: unit = 28 ! see comment lines at top of file INTEGER :: i REAL, DIMENSION(3) :: tv filename = Insert (p_dat, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') WRITE (unit, "(A)") p_dat_format WRITE (unit, "(A)") p_dat_titles DO i = 1, p_dat_count IF (p_site_is(i)%element > 0) THEN c50 = p_ref(i) tv = p_site_0(1:3, i) CALL Lonlat_from_xyz (tv, x1, x2) r1 = (p_south(i) / R) * deg_per_rad r2 = (p_south_sigma_(i) / R) * deg_per_rad r3 = p_ccw(i) * deg_per_rad r4 = p_ccw_sigma_(i) * deg_per_rad t2 = p_t2(i) / s_per_Ma t1 = p_t1(i) / s_per_Ma tv = p_pole(1:3, i) CALL Lonlat_from_xyz (tv, x3, x4) WRITE (unit, p_dat_format) c50, x1, x2, r1, r2, r3, r4, t2, t1, x3, x4 IF (paleotec .OR. (start_time > 0.)) THEN tv = p_site_now(1:3, i) CALL Lonlat_from_xyz (tv, x1, x2) WRITE (unit, "('+',1X,F8.3,1X,F7.3)") x1, x2 END IF CALL Write_rates_and_goals & & (index = i, unit = unit, signal = '*', & & conversion = ((deg_per_rad * s_per_Ma) / R), & & goal = p_south_goal, rate = p_south_rate, active = p_active) IF (.NOT. twisted(i)) & & CALL Write_rates_and_goals & & (index = i, unit = unit, signal = '&', & & conversion = (deg_per_rad * s_per_Ma), & & goal = p_ccw_goal, rate = p_ccw_rate, active = p_active) END IF END DO CLOSE (unit) END SUBROUTINE Write_p_dat SUBROUTINE Write_rates_and_goals & & (index, unit, signal, conversion, goal, rate, active) ! Output histories of rates (left) and goals for next iteration (right). INTEGER, INTENT(IN) :: index, unit CHARACTER(1), INTENT(IN) :: signal REAL, INTENT(IN) :: conversion REAL,DIMENSION(:,:),INTENT(IN) :: goal, rate LOGICAL(1),DIMENSION(:,:),INTENT(IN):: active INTEGER :: j REAL :: g, r, t0, t1 DO j = 1, num_timesteps ! global; = 1 IF (neotec) IF (active(j, index)) THEN r = rate(j, index) * conversion g = goal(j, index) * conversion IF (paleotec) THEN t0 = ((j - 1) * Deltat_) / s_per_Ma t1 = (j * Deltat_) / s_per_Ma ELSE ! (neotec) t0 = start_time / s_per_Ma t1 = t0 END IF WRITE (unit, "(A,1X,F7.3,1X,F7.3,1X,F9.4,1X,F9.4)") signal, t0, t1, r, g END IF END DO END SUBROUTINE Write_rates_and_goals SUBROUTINE Write_s_dat (shift) ! Writes smmnn.dat output file with +, $ lines. ! Note that points outside the grid are dropped. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt CHARACTER(80) :: filename INTEGER :: unit = 29 ! see comment lines at top of file INTEGER :: i REAL, DIMENSION(3) :: tv filename = Insert (s_dat, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') WRITE (unit, "(A)") s_dat_format WRITE (unit, "(A)") s_dat_titles DO i = 1, s_dat_count IF ((s_site_is(1,i)%element > 0).AND.(s_site_is(2,i)%element > 0)) THEN c30 = s_ref(i) c30a = s_loc(i) c5 = s_code(i) tv = s_site_0(1:3, i) CALL Lonlat_from_xyz (tv, x1, x2) r1 = s_azim_0(i) * deg_per_rad r2 = s_sigma_(i) * deg_per_rad t2 = s_t_max(i) / s_per_Ma t1 = s_t_min(i) / s_per_Ma IF (s_stage(i)) THEN c6 = "Stage " ELSE c6 = "Window" END IF WRITE (unit, s_dat_format) c30, c30a, c5, x1, x2, r1, r2, t2, t1, c6 tv = s_site_now(1:3, 1, i) CALL Lonlat_from_xyz (tv, x1, x2) WRITE (unit, "('+',1X,F8.3,1X,F7.3)") x1, x2 r1 = s_azim_now(i) * deg_per_rad r1 = MOD( (r1 + 360.), 180.) WRITE (unit, "('$',1X,F5.1)") r1 END IF END DO CLOSE (unit) END SUBROUTINE Write_s_dat SUBROUTINE Write_x_feg (shift) ! Writes xmmnn.feg output file. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt CHARACTER(80) :: filename INTEGER :: unit = 23 ! see comment lines at top of file INTEGER :: i, l_ LOGICAL :: faulting REAL :: lon, lat, s1h_azimuth, s1h_sigma REAL, DIMENSION(3) :: tv filename = Insert (x_feg, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') WRITE (unit, "(A)") TRIM(filename)//', restored from '//TRIM(x_feg) WRITE (unit, "(I7,4X,I7,' 0 30000 T')") num_nod, num_nod DO i = 1, num_nod tv = xyz_nod(1:3, i) CALL Lonlat_from_xyz (tv, lon, lat) WRITE (unit, "(I7,1X,F9.4,1X,F8.4,1X,1P,E9.2,' 0. 0. 0.')") i, lon, lat, mu_nod(i) END DO WRITE (unit, "(I7)") num_ele DO l_ = 1, num_ele IF (stress_ever) THEN s1h_sigma = ele_sigma(l_) * deg_per_rad s1h_azimuth = ele_azim(l_) * deg_per_rad IF (s1h_azimuth < 0.) s1h_azimuth = s1h_azimuth + 180. IF (s1h_azimuth < 0.) s1h_azimuth = s1h_azimuth + 180. IF (s1h_azimuth < 0.) s1h_azimuth = s1h_azimuth + 180. IF (s1h_azimuth >= 180.) s1h_azimuth = s1h_azimuth - 180. IF (s1h_azimuth >= 180.) s1h_azimuth = s1h_azimuth - 180. IF (s1h_azimuth >= 180.) s1h_azimuth = s1h_azimuth - 180. faulting = (crack_index(1, l_) /= 0) WRITE (unit, "(4(I7,1X),L1,1X,0P,F4.0,1X,F4.1,1P,3E10.2,1X,L1)") & & l_, node(1,l_), node(2,l_), node(3,l_), & & ele_stressed(l_), s1h_azimuth, s1h_sigma, & & ele_strainrate(1:3,l_), faulting ELSE WRITE (unit, "(4(I7,1X))") & & l_, node(1,l_), node(2,l_), node(3,l_) END IF END DO WRITE (unit, "(' 0')") CLOSE (unit) END SUBROUTINE Write_x_feg SUBROUTINE Write_x_vel (shift, vw) ! Writes xmmnn.vel output file. INTEGER :: shift ! controls placement of announcement on screen, in REPORT.txt REAL, DIMENSION(:), INTENT(IN) :: vw CHARACTER(80) :: filename INTEGER :: unit = 25 ! see comment lines at top of file INTEGER :: i REAL :: tMa filename = Insert (x_vel, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'WRITE', FILE = filename, STATUS = 'REPLACE') WRITE (unit, "(A,', for use with')") TRIM(filename) IF (neotec) THEN WRITE (unit, "(A)") TRIM(x_feg) ELSE WRITE (unit, "('a deformed version of ',A)") TRIM(x_feg) END IF tMa = time1 / s_per_Ma CALL DATE_AND_TIME (date, clock_time, zone, datetimenumber) WRITE (unit,"(F6.2,' Ma; iteration ',I3,'; computed ',I4,'.',I2,'.',I2,' at ',I2,':',I2,':',I2)") & & tMa, total_iterations, & & datetimenumber(1), datetimenumber(2), datetimenumber(3), & & datetimenumber(5), datetimenumber(6), datetimenumber(7) DO i = 1, num_nod WRITE (unit, "(1P,E13.6,1X,E13.6)") vw(2 * i - 1), vw(2 * i) END DO CLOSE (unit) END SUBROUTINE Write_x_vel SUBROUTINE Write_y_dig (shift) ! Writes ymmnn.dig output file. ! Note that points outside the grid are dropped; ! segments which are reduced to less than 2 points are also dropped. INTEGER, INTENT(IN) :: shift ! controls placement of announcement on screen, in REPORT.txt INTEGER :: unit = 24 ! see comment lines at top of file INTEGER :: i, n_in_seg REAL :: lon, lat LOGICAL :: seg_open CHARACTER(80) :: filename REAL, DIMENSION(3) :: tv filename = Insert (y_dig, mmnn) ! both arguments global variables; typical IF (shift == 8) THEN PRINT "(' ',8X,'Writing ',A)", TRIM(filename) WRITE (21,"(8X,'Writing ',A)") TRIM(filename) ELSE PRINT "(' Writing ',A)", TRIM(filename) WRITE (21,"('Writing ',A)") TRIM(filename) END IF OPEN (unit, ACTION = 'READWRITE', FILE = filename, STATUS = 'REPLACE') ! note that READWRITE is needed to permit BACKSPACE, below seg_open = .FALSE. ; n_in_seg = 0 DO i = 1, y_dig_count ! global IF (dot_is(i)%element > 0) THEN ! point lies within grid; location known tv = dot(1:3, i) CALL Lonlat_from_xyz (tv, lon, lat) WRITE (unit, "(1X,SP,1P,E12.5,',',E12.5)") lon, lat seg_open = .TRUE. ; n_in_seg = n_in_seg +1 ELSE ! this point fell outside the grid; location undefined IF (seg_open) THEN IF (n_in_seg > 1) THEN WRITE (unit,"('*** END OF SEGMENT ***')") ELSE BACKSPACE (unit) END IF seg_open = .FALSE. ; n_in_seg = 0 END IF END IF IF (dot_last(i)) THEN IF (seg_open) THEN IF (n_in_seg > 1) THEN WRITE (unit,"('*** END OF SEGMENT ***')") ELSE BACKSPACE (unit) END IF seg_open = .FALSE. ; n_in_seg = 0 END IF END IF END DO CLOSE (unit) END SUBROUTINE Write_y_dig SUBROUTINE Xyz_from_lonlat (lon, lat, vector) REAL, INTENT(IN) :: lon, lat REAL, DIMENSION(3), INTENT(OUT) :: vector ! assuming longitude is East longitude in degrees, ! and latitude is North latitude in degrees, ! computes 3 components of Cartesian unit vector ! from center of planet to surface point (on unit sphere). REAL :: theta_, phi_, equat theta_ = (90. - lat) / deg_per_rad phi_ = lon / deg_per_rad equat = SIN(theta_) vector(1) = equat * COS(phi_) vector(2) = equat * SIN(phi_) vector(3) = COS(theta_) END SUBROUTINE Xyz_from_lonlat END PROGRAM Restore2