24 June 1998 TO: Fellow Earth Scientists and Hackers FROM: Peter Bird Department of Earth and Space Sciences University of California Los Angeles, CA 90045 (310) 825-1126 pbird@ess.ucla.edu RE: Introduction to software of the "LARAMY" group and associated utility programs THIS FILE CONTAINS A GENERAL INTRODUCTION TO THE SOFTWARE I AM OFFERING, WITH INSTRUCTIONS ON HOW TO OBTAIN RELATED PROGRAMS. THIS IS FOLLOWED BY DETAILED NOTES ON THE INDIVIDUAL FILES OF THE "LARAMY" SUBDIRECTORY. I currently offer four thin-plate finite-element programs: LARAMY finite-strain no faults flat-Earth crust/mantle FAULTS neotectonics with faults flat-Earth crust only PLATES neotectonics with faults flat-Earth crust/mantle SHELLS neotectonics with faults round-Earth crust/mantle The finite element program LARAMY may be copied and used without charge or restriction, except that (1) any user who obtains any publishable results is expected to acknowledge the source of the program, and (2) no charge may be made for copies of the software which are passed on to third parties. It was developed with support from the National Science Foundation and the University of California, which places it in the public domain. Also, I would be very pleased to see my 6 years of programming effort bring some benefit to others besides myself. Finally, I have no desire to be personally responsible for running simulations of every interesting problem for the rest of my career! This FORTRAN program was created to permit simulations of continental deformation in which a crust of either one or two rheologic layers (either one of which may enclose the brittle-ductile transition) overlies a mantle lithosphere and/or fluid asthenosphere and/or subducting slabs. While it has been customized for simulations of the Tertiary history of North America, the problem-specific data are mostly enclosed in DATA statements in two BLOCK DATA subprograms which can be replaced with relative ease. Its only fundamental limitations are that it is (1) a continuum program, with no fault discontinuities, but only an equivalent frictional plasticity; (2) rigorously isostatic, and neglects all bending stresses; and (3) written for a flat Earth. Before attempting to work with Laramy, it is essential to gain an overview of the program's structure, abilities, and limitations. Readers are STRONGLY ADVISED to skim the following two articles before proceding further: Bird, P. (1988) Formation of the Rocky Mountains, western United States: a continuum computer model, Science, v. 239, p. 1501-1507. Bird, P. (1989) New finite element techniques for modeling deformation histories of continents with stratified temperature-dependent rheologies, J. Geophys. Res., v. 94, p. 3967- 3990. The essential hardware required to run LARAMY is only about 5 megabytes of fast memory. The software needed to run LARAMY is just a FORTRAN 77 compiler (optimizing, if possible) and a couple of matrix-manipulation routines: one for solving large, banded linear systems, and one for finding eigenvalues and eigenvectors of small matrices. If you have IBM's ESSL package at your institution, you can use the same matrix routines that I did, without any change to the programs. A Fortran 90 compiler should also work, since they are supposed to be backward- compatible. (You might have to change some of my simple "END" statements to "END SUBROUTINE WHATEVER" to satisfy Fortran 90.) Here are some actual run times and costs for LARAMY: *IBM 3090 with Vector Facility, running MVS operating system: 3.5 hours (costing about $560 @ intramural rate; May 1989). *IBM RS/6000 model 590 ("phi"), running AIX operating sytem: 1.2 hours (costing about $48 @ intramural rate; June 1998). (Details of the job: Each layer [crust, mantle lithosphere] had 609 nodes arranged in 29 rows of 21 columns; each layer had 280 triangular elements. Duration of 85 Ma was modeled with 1 Ma time steps, giving detailed output every 5 Ma. Up to 10 iterations of the velocity solution were used in each timestep to converge to 2% precision. Up to 1000 diffusive smoothing steps were used per timestep to model lateral diffusion of lower crust.) To keep I/O simple, program LARAMY reads all input from a single (small) text file, and puts all output into two text files: (1) a user-friendly file formatted for line- or laser-printer (with 133- character lines), including "printer plots"; and (2) a packed-numbers file that is used for graphics post-processing and/or restarting a simulation. Thus, the program itself is NOT interactive or graphical, and does not require any particular hardware or software for its interfaces. The following accessory programs are available in other subdirectories of this archive. (Note: numbers 0 - 5 tie to Read_me.n filenames.) 0. Laramy2AI (added 1997) is a graphics program that takes output from LARAMY and produces single or multiple maps of various scalar, vector, or tensor variables-- as .AI files readable by Adobe Illustrator for Windows. This allows the maps to be previewed and edited(!) before printing, and also allows the printing to take place on a wide variety of hardware (including hardware run by service bureaus)! The .AI files produced can be read by Adobe Illustrator 4 for Windows 3.1, and by Adobe Illustrator 7 for Windows 95 or Windows NT. They probably can be read by Adobe Illustrator on Macs, too. The processing program Laramy2AI.exe runs under Windows 95 or Windows NT on Intel (or compatible) processors. Source code is also provided in Fortran 90 so that you can move this utility to Mac, Unix, or other systems. (If you compile it with XL Fortran for AIX, use the "-qnoescape" switch!) 2. GRAFIC is a two-part graphics package that post-processes finite element results to produce colored contour maps of all variables, with overlays to show vector and tensor directions and/or state lines, etc. The number-crunching part is run as a batch job, and produces graphics metafiles. The second part is run interactively from a graphics terminal to view the images at various scales, with or without the overlays, etc. These images can be photographed from the screen to make slides. This (older) graphics package requires an IBM utility package called Graphical Data Display Manager. (To my knowledge, this is only for mainframes, not PCs.) Its capabilities are roughly the same as Laramy2AI, but it is designed for a mainframe environment where you might keep libraries of 20-100 plots (10-50 Mb) per model run. 3. VERSATEC is comparable to GRAFIC, but produces either B/W or color output on a Versatec electrostatic plotter (for publications). It runs as a single batch job and is not interactive at all. Naturally, your institution must own own a Versatec plotter to use this option. If you bought the plotter, then I presume you also got the necessary library of graphics subroutines to drive it. 4. FIXER is a fully interactive package for correcting small parts of a finite element grid which become excessively strained, so that a run can continue. It can also be used to fine-tune the files representing initial conditions. (However, it is not a bells-and-whistles grid-generator like some expensive engineering programs offer; just an editor. You must have already created and input a grid with the same topology as the final product you want.) A color graphics monitor attached to your mainframe is required, but not a mouse. This also requires IBM's Graphical Data Display Manager. 5. MAPPER is a fully interactive package for making plate reconstructions of oceanic plates which are no longer on the Earth's surface (based on marine magnetic anomaly data and finite rotations) and saving the results in the form of DATA statements that can then be compiled into LARAMY in order to represent complex boundary conditions on the bottom of continents. This program also interactively displays the boundary conditions as they would be applied to a representative finite-element grid, so that unexpected or erroneous BC's can be corrected. A color graphics monitor attached to your mainframe is required, but not a mouse. This also requires IBM's Graphical Data Display Manager. The hitch with accessory programs #2-#5 is that, while they have been written in pure FORTRAN77, they make thousands of calls to proprietary graphics packages which may not be available to you. My packages GRAFIC, FIXER, and MAPPER make calls to subroutines of the IBM Graphical Data Display Manager software (sometimes known as "that G*D-D*M* software"); my VERSATEC program calls the proprietary package that is supplied with that brand of plotter. My advice to you is that if you only want to plot the output, you should use Laramy2AI. If you need to use FIXER or MAPPER (and don't have GDDM), it will probably be easiest and cheapest to hire a programmer to provide dummy interface subroutines that accept these calls and translate them into further calls to your local graphics software (rather than programming these tools from scratch). Or, If you have finite-element post-processors that you are satisfied with, another simple alternative is to write an "output translator" which reads the packed-numbers (device-9) output from LARAMY and puts it into the format required by your existing package. Please understand that I will expect you to communicate with me if you find significant bugs or opportunities for program improvement, and that if you use these programs to obtain publishable results I will expect the source of the program to be acknowledged in your manuscript. I also hope that you will remember me when you have reprints to distribute. =========================================================== NOTES ON THE FILES OF THE "LARAMY" SUBDIRECTORY OLD_LARAMY.for = Fortran77 source code. (Note: "OLD_" refers to the original version of the rheologic code in SUBROUTINE VISCOS and SUBROUTINE DIAMND. Someday I will rewrite these to parallel my improved versions of SHELLS, PLATES, and FAULTS.) compile.jcl = IBM Job Control Language used to compile LARAMY (on UCLA's MVS system) and create a load module. (Note: on an AIX system with XL Fortran you only have to enter "f77com -lessl OLD_LARAMY.f" .) runner.jcl = IBM Job Control Language used to run an MVS load module. OLD_LARAMY.aix = AIX or UNIX commands to run load module OLD_LARAMY.exe, demonstrating the proper method of connecting input and output files. IN9828 = sample parameter input file, containing values that will give a decent simulation of the Tertiary history of North America. (Note that this file is also required to run every one of the graphics packages. Do not write over your only copy.) CORD11N = sample initial-condition file, often (but not necessarily) used to begin a run. This one represents North America 85,000,000 years ago. (Note that this file can be used to test out the graphics packages as well, as it is formatted just like output from LARAMY.) P9828A00 = sample text (device-6) output from OLD_LARAMY when run with input files IN9828 and CORD11N on an IBM RS/6000 under AIX. (This run took 1.2 hours.) Contains crude printer-plots of the variables every 5 Ma. T9828A00 = sample packed-arrays (device-9) output from OLD_LARAMY when run with input files IN9828 and CORD11N on an IBM RS/6000 under AIX. Can be used with any of the graphics programs described above to better visualize the output. ----------------------------------------------------- LARAMY I suggest that you print this file of about 14,500 lines on a high-speed line or laser printer attached to a mainframe, as it will amount to about 300 pages. You will need this printout as a reference. You will also need to make an index for yourself to stay sane. The first thing to notice is that the program contains a good deal of internal documentation. The first several pages explain what kinds of input and output files are associated with the Fortran device numbers 1, 6, 8, and 9. They identify the source and functions of the few external routines that are used. They also contain some PARAMETER statements that you may want to adjust before compiling, because these determine the size of the load module (about 5 megabytes, on an IBM 3090) and the dimensions of the largest grid allowable. (If you try to exceed these maxima, the subroutine SETDIM will detect the problem and print a report.) The next part of the main program is an extensive set of comment (C) lines containing a glossary of the array names, defining the variables AND their subscripts (don't you wish everyone would do this?). This should be your constant reference if you need to modify or (Heaven forbid!) debug the program. I suggest that you do NOT modify the program at all until AFTER you have successfully run a timestep or two using the input files as supplied, and inspected (and perhaps plotted) the output. Two aspects of the output on device 6 (line printer? remote terminal?) deserves explanation. First, it has line lengths up to 133 characters, with the first character reserved for printer-control characters (i.e., 0=double-space, 1=new-page, +=overstrike). This means you can print it readily on most high-speed/low-quality chain printers, or sideways (in landscape mode) on most laser printers. If you use a word-processing program to inspect the output, you will want to look at it with a word processor that does windowing instead of line-wrapping. Second, I should explain the format of the printer- plots. It was necessary to display maps of scalars (e.g., heat flow, crustal thickness, elevation), vectors (e.g., velocity), and second- rank tensors (e.g., strain-rate, finite strain, vertically-integrated stresses), by using only text characters. My solution to this was to use the numerals 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, * (* = 10) to represent magnitudes, according to a scale printed below the map. This means that each magnitude is uncertain by plus or minus 5% of the range of the map. Vectors will also have a direction indicated by a letter. These letters correspond to numbers on a clockface: A = 1 o'clock = 30 degrees clockwise from +y B = 2 o'clock = 60 degrees clockwise from +y C = 3 o'clock = 90 degrees clockwise from +y (+x axis) D = 4 o'clock = 120 degrees clockwise from +y E = 5 o'clock = 150 degrees clockwise " " F = 6 o'clock = 180 degrees clockwise " " G = 7 o'clock = 210 degrees clockwise H = 8 o'clock = 240 degrees clockwise I = 9 o'clock = 270 degrees clockwise J = 10 o'clock = 300 degrees clockwise K = 11 o'clock = 330 degrees clockwise L = 12 o'clock = 0 degrees clockwise so that the directions of the vectors are uncertain by plus or minus 15 degrees in this display. It is not really possible to display a tensor in this way, even a reduced tensor like the horizontal-plane strain rate with its three independent components (extension rate along X, extension rate along Y, and shear in the XY plane). My solution was to isolate the principal axis in the horizontal plane which has the greatest absolute value, and to plot it as if it were a vector. Of course, it is a funny kind of vector, because it could equally well be plotted in the opposite direction (180 degrees opposite). Therefore, we do not really need all the direction letters A-L, but only half of them. The convention I have adopted is that if the principal stress or strain- rate or finite strain is compressive, then I choose the end of the principal axis with a direction in the A-F range to plot, and if it is tensile then I choose the end in the G-L range to plot. This takes some time to get used to. These printer plots also each contain a single character #, which marks the origin of the x/y Cartesian coordinates. It will usually be somewhere in the middle of the problem. At this point, and at this point only, +y ("up" on the plot) corresponds to North, and +x ("right" on the plot) corresponds to East. I quite agree that these plots are not really satisfactory, but they give a "quick look" at solutions, even over complex networks, modem links, etc. which can only handle plain ASCII text. To obtain more precise and beautiful graphics output, see the other directories mentioned above, or write your own program. Most users of LARAMY will not be content to run North American cases. Therefore, they will need to slightly modify the source code in file LARAMY. Here is some essential background information to help you understand what you are doing. 1. FINITE ELEMENTS. Program Laramy and its associated utility programs work with only one kind of element: the 6-node, isoparametric plane triangle. This is a very general element, whose sides can be straight or curved. (If s is an internal variable that varies 0 to 1 along a side, then the side is located at: x = quadratic_function_#1(s); y = quadratic_ function_#2(s)). Do not even THINK about introducing other element types. I MIGHT be able to do it, but I haven't got two years to spare. Within each triangular element, the nodes are assigned local numbers 1-6, which are temporary aliases for their global node numbers. (The user does not have to know this, but the programmer might.) They are arranged generally counterclockwise like this: 1--------6--------3 2 | / /| | / / | | / / | | / 5 4 | / / | | / / | | / / | 4 5 / | | / 3---6----1 | / | / | / | / | / | / | / 2 Points within an element are often located by their 3 internal coordinates S1, S2, and S3. Regardless of the size or shape of the element, S1 varies from 0 (along the side 2-5-3) to 1 (at node #1). S2 varies from 0 (along side 3-6-1) to 1 (at node #2). S3 varies from 0 (along side 1-4-2) to 1 (at node #3). It follows that S1+S2+S3 = 1.00 identically for all points. A negative value of any coordinate means that you are outside the element. (Note: Mineralogists and petrologists will recognize their old friend, the "ternary diagram" in distorted form). If the element has bent sides, then the contours of S1 (and S2 and S3) will also be bent. 2. GRID TOPOLOGY. Another restriction in that Laramy and its associated utilities only work with a topologically-rectangular grid. Imagine two triangular elements placed together (hypotenuse to hypotenuse) to form a quadrilateral "cell". A CELL: *----*-----* | /| | / | | / | | / | * * * | / | | / | | / | | / | *-----*----* The grids in Laramy always consist of integer numbers of rows and columns of these quadrilateral cells. Every row is the same length (measured by number of cells). Every column is the same length (but not necessarily as long as the rows). In every cell of this rectangular grid, the hypotenuse slants the same way. If there are R rows of quadrilateral cells and C columns of quadrilateral cells, then the nodes are also in a topologically-rectangular array with (r = 2R+1) rows and (c = 2C+1) columns. What makes this a finite element program instead of a finite difference program is that this rectangular "net" of elements can be extremely distorted, as long as it is not cut, overlapped (even locally!), or overturned. (If you get any messages about "negative area" it means that the grid is locally folded; technically, the mapping from (x,y) to (S1,S2,S3) is no longer functional because it is not single-valued. This can happen even in a grid of 9 nodes and 2 elements if the strain in the elements is large enough and inhomogeneous enough.) It is important to understand the numbering of nodes and of elements in this rectangular array. One side is designated as the "top". Now, the other three sides will be referred to as "right", "bottom", and "left", in clockwise order. If you want to have subduction beneath the base of your model, it MUST come from the "left", so this will define which side is which. If the paper is rotated to place the "top" side on top, all the hypotenuses inside the cells MUST slant from lower left to upper right. The nodes are assigned global node numbers in the same order as text in a book: "left" to "right" within rows, with rows progressing down from "top" to "bottom". Thus, the #1 node is at the "top left" and the highest-numbered node (#rc) is at the "bottom right". Notice that the largest difference between the numbers of two nodes involved in any element (the bandwidth) is 2c = 4C+2, and that this is constant throughout the grid. The elements are numbered in a similar way, progressing "left" to "right" in each row of cells, and then "down" to the next row. Within each quadrilateral cell, the element toward the "top left" has the lower (odd) number, and the element toward the "bottom right" has the higher (even) number. You may be tempted to try to generalise Laramy for non- rectangular grids, but I advise against it. First, you would have to seek out all the subroutines in which the program assumes a rectangular grid or "assumes it knows" which elements are next to which, and which ones are on the boundary. The list would include subroutines READIN, SETDIM, GRIDDR, HOWFAR, SETEDT, ORTHO, LINKER, LOOKUP, BUILDF, BCS, SLIPBC, UNFOLD, SMOOTH, STENCL, DRAGIT, PANCAK, FETOFD, EBCS, and NET at least, and perhaps others I have forgotten. A more fundamental reason is that the algorithms of two very tricky subroutines absolutely require topological regularity. LINKER (and the associated LOOKUP) which find the element and internal coordinates matching an arbitrary (x,y) point use a search algorithm that only works on such grids. And PANCAK (with its several associated subroutines), which handles lateral extrusion of lower crust, uses a very-fast distorted-grid finite-difference program that also cannot work except with such grids. You would have to invent completely new algorithms for these functions, and you would have to make them FAST. Perhaps after you have drawn a few of your own grids, you might be worried that this restriction requires you to cram a lot of small elements together someplace where they are not really needed, and thus wastes money. Well, relax. The primary cost is usually that of solving the linear equations for velocity, and this cost depends more strongly on the bandwidth of the grid (bandwidth-squared) than it does on the total number of nodes (linear in this number). Since we are already paying for a certain bandwidth, it might as well be full of nodes and elements, eh? (More degrees of freedom can never degrade the accuracy of a solution.) 3. CARTESIAN COORDINATE SYTEM. We have already covered two coordinate systems: the local internal coordinates of one element (S1,S2,S3), and the topological (integer) coordinates of the node and element numbers and row numbers and column numbers of the grid. However, the primary coordinate system for posing and solving the differential equations is Cartesian: (x,y,z). This means that Laramy is a flat-Earth program, and cannot be used for global modeling without extensive revisions. Directions x and y are horizontal, with y 90 degrees counterclockwise from x. At the origin, +y points North, and +x points East. All plots produced by Laramy (and by the graphics programs in the other directories) will orient their output so that x is to the right any y is up. This cannot be modified except by recoding. Therefore, you should give some thought to the orientation of coordinates when beginning a new project. The restriction to a flat Earth is sometimes an actual convenience, because it means that grids can be drawn on a flat map (perhaps a conic projection of part of the round Earth), digitised, and used directly in Laramy without any transformations. Of course, the geometric distortion implicit in the base maps becomes a (small) source of error in the solutions. Although I previously referred to the topology of the grid as having a "top" and "bottom", there is NO requirement that the "top" be toward +y, or that the "right" be towards +x. You will rotate the topological super-rectangle of the grid so that subduction (if any) occurs along the "left" side. Remember, the sides of the grid, and the internal rows and columns, may be bent as much as necessary. (It should be noted that it is dangerous to think of +x as being East or +y as being North. Because of the curvature of the Earth, these identities can only hold for a single value of longitude, and not for a whole region.) Program Laramy makes use of round-Earth (latitude, longitude) coordinates in only one place, which you are free to ignore if you wish. This is inside BELOWY, a subroutine which tells Laramy what oceanic slabs are touching the continent from below, and how they are moving. BELOWY makes use of digitised maps of present-day positions of reconstructed seafloor, which are stored in (latitude, longitude) coordinates. Furthermore, to determine the paleo-positions of points on the subducting plates, it uses round-Earth finite rotations for accuracy. That means it must convert from (x,y) to (lat,lon), do the rotation, and then convert back to (x',y'). For this reason, subroutine READIN will ask for the latitude and longitude of your x/y origin (and the tangent latitude of your base-map projection) and pass this information to BELOWY; it is not used anywhere else. =========== CHECKLIST FOR APPLYING LARAMY TO NEW PROBLEMS ========= A) Select a base map. If it is a conic projection, note the tangent latitude to place in the parameter file. (If it is a more complex projection, you can just input a reasonable mean latitude for the map, and this will be almost the same thing.) B) Select a center of interest for your problem, and make this the x/y origin. This point will be "right side up" (North = +y, East = +x) in all plots. C) Draw the +y axis straight Northward from the origin. D) Draw the +x axis as a straight line 90 degrees clockwise from +y. This begins as East, but then diverges (unless you have a Mercator basemap, which is similar to a conic projection with a tangent at the equator). E) Decide whether you will work in cgs units or SI units. (For your first runs, you will probably use cgs, because my sample parameter files are in cgs, and it is not trivial to convert rheological units.) If you are using cgs, mark off distances along x and y in centimeters (e.g., 0, 2.E7, 4.E7,...). If you are using SI, mark off distances along x and y in meters (e.g., 0, 2.E5, 4. E5,...). F) Draw in a topologically-rectangular grid of cells, elements, and nodes on your base map, following the rules above. The "top" of the grid does not have to be toward the top of the base map. However, there is a rule that the "left" side of the grid MUST be the side with the trench (convergent margin) if you wish to have a side and/or basal boundary condition that represents subducting plates during any part of the model history! Note that the diagonals (hypotenuses) of the quadrilateral cells must all slant from "top right" to "bottom left". The grid may be as distorted as you like, but be wary of extremely distorted elements. (This is because elements will become more distorted during the run, and the program will crash if they "fold"; you will then have to hand-edit the grid with the program in the FIXER subdirectory.) I suggest that you use straight lines for all element sides (except those needed to represent trench boundaries) at first. Until you have some experience with costs, limit yourself to no more than 14 rows and 10 columns of quadrilateral cells. Cost rises as R * C**3 ! This means that in terms of rows and columns, the grid should be tall (large R) rather than fat (large C). Label the side with the first row as "top", the side with the last row as "bottom", the side with the first column (the trench side) as "left", and the side with the last columnn as "right". G) Make a copy of the parameter-input file (which will be read from device 1 by subroutine READIN) from the sample given (file IN9003), and begin to modify it for your problem. For example, change the title on line 2, which identifies a particular model in all the resulting plots and files. It is VERY IMPORTANT that two of the integer parameters should be correct. In line 45, parameter NELROW must be equal to the number of rows of quadrilateral cells in your grid (= R). Also, parameter NELCOL (line 46) must be equal to the number of columns of quadrilateral cells in your grid. Otherwise, the program will crash almost immediately, because it will find that the input files are the "wrong length". (Check this carefully if you get error messages saying that letters were read when numbers were expected.) On the other hand, there are 3 parameters used for automatic grid generation (height, width, and angle) that are irrelevant if you create your own grid. These parameters are currently found toward the end of IN9828. You should not delete these lines, but it is not necessary to understand them or modify them either! Just leave them alone. H) Digitise the node locations in the order previously explained (text order; from "left" to "right" in each row, and progressing from the "top" row first to the "bottom" row last). If necesary, scale the results for the scale of the map. Then, write a tiny program in any language to reformat these numbers as LARAMY expects them. The x values will become array XNODC, and the y values will make up array YNODC. Look at how these are formatted in file CORD11N, and write a small program to reformat your numbers in the same way: (1P,6E13.6). Don't forget to include the title, time, and array-name lines in the new file you are building; use CORD11N as a model. I) If there is no subduction beneath your model, your initial grid for the mantle lithosphere will probably be a duplicate of the initial grid for the crust. If so, just copy your two arrays and name the copies XNODM and YNODM, and arrange the arrays with their accompanying text labels in the (alphabetical) order that you find in file CORD9 (and in subroutines TAPE and GOON). However, if you want to have subduction, you may wish to have an initial offset between the crust and mantle-lithosphere grids, perhaps 200 km or more. In this case, you will need to draw a separate mantle grid and digitise it as well. Then merge the two sets of node locations into one file, as in CORD11N. If you DON'T want the crust and mantle layers of the continental lithosphere to cover the same area in the initial condition, then remember that EVERY mantle node and element must be UNDER a crustal element. (In order to guarantee that this rule is still obeyed as the models run, it is usually a good idea to use impermeable or immobile side boundary conditions on the side(s) toward which any basal subducted slabs are traveling.) It's OK if crustal elements do not have mantle elements beneath; they will be assumed to have asthenosphere and/or subducting slabs beneath. The cooling of the asthenosphere as it ages in contact with the crust will be computed and stored in array GEOTHA, so that the Moho temperature can decline over time; however, no new elements will be generated to represent this new mantle lithosphere, so its mechanical strength witll be neglected. Remember that it is REQUIRED that the crust and mantle grids have the same topology and number of nodes. If you really don't want a mantle boundary layer for some reason, you will have to make it very small and cram it against the inland side of the grid, I suppose. (I've never tried this.) J) Now you have the task of creating all the tedious layer thickness values, geotherms, initial finite strains, etc. that make up the rest of the initial-conditions file (see CORD11N). One neat way to do this is to set logical parameter OLDGRD to T in the parameter file (e.g., IN9828), and set RESTRT to F. Then the program will read ONLY your new title, time, and digitised node locations from device 8, and will generate the rest of the arrays based on your input parameters (read from device 5). Notice that you can have an initial cordillera built for you along the "left" (1st-column) side of the grid if you wish. If your problem involves subduction, you should definately use this option, which automatically simulates an active margin. (You can make the height of the cordillera as little as 1 centimeter, if you wish; but keep the width reasonable.) Set the program for a very short run; you might specify a duration (difference of beginning and ending times) of 0 seconds and the maximum number of iterations per timestep as 1. The program will soon stop and output the completed initial-conditions file on device 9. (If you only allowed one iteration, do not be concerned that the solution is "bad" and "unconverged". This only referes to the horizontal velocities, strain-rates, and stresses. The other variables are accurate.) K) If you want more details (e.g., present-day topography and heat flow) in your initial conditions, you can obtain them by editing the file that was just created on device 9 in the previous step. The programs in the FIXER package, in the subdirectory of that name, will allow to do interactive editing on colored contour maps on a graphics terminal, which is very sexy but slow. If you want to input a lot of geographic variation, you should write your own program to modify the initial-condition file, based on your own datasets. L) Modify the code inside Laramy as needed to prescribe your lateral (side) boundary conditions. This code is found in two places: In subroutine BUILDF, a pressure equal to the local lithostatic pressure is applied to all sides of the grid. (Lithostatic pressures were previously determined by subroutine SQUEZE, and BUILDF just extracts these values from array TAUZZ.) The reason for this is that in the absence of velocity boundary conditions, and in the absence of variations in topography, we want to obtain a "null solution" in which the rocks just sit there and there are no velocities or strain rates or shear stresses.) Do NOT be concerned that this pressure boundary condition conflicts with some velocity boundary conditions that you might want to impose; in finite-elements, velocity boundary conditions always overwrite and replace pressure boundary conditions! Therefore, you will not need to modify this routine unless you want a pressure boundary condition that is NOT locally lithostatic. To accomplish this, find the block of code in BUILDF delimited by IF (LEFTY.OR.RIGHTY) THEN.....ENDIF and the following block delimited by IF (UPPER.OR.LOWER) THEN......ENDIF. Within each block, add the desired pressure anomalies in the 3 statements that look like: TAUZZ1=TAUZZ(6,I) + DELTAU |<-- old code ->| |<-- you add-->| You can make DELTAU vary from place to place as you wish. Be aware that the quantity to add (DELTAU) should be the negative of the vertical integral of the pressure anomaly through the layer (crust or mantle). That is, to add 1 kilobar of anomalous compression on the edge of a 40 km crust, the cgs value of DELTAU is -1.E9*4.E6= -4.E15. This subprogram is called once for the crustal layer and once more for the mantle layer; to tell which layer is being processed, your added code can use the value of the logical variable MANTLE. For example, IF (MANTLE) THEN DELTAU = 0.0 ELSE DELTAU = -4.E15 ENDIF The other place where lateral (side) boundary conditions are set is in subprogram BCS, which dictates which nodes are restricted and/or forced to move. This is currently written to set velocities to zero on the "top" (first-row) side, "right" (last- column; usually inland) side, and "bottom" (last-row) side of the grid, implicitly leaving a pressure boundary conditions on the "left" (trench?) side. If you will read the code of BCS, you will see that it is trivial to change. Subprogram FIXNOD, which does the actual work of imposing the velocity boundary conditions at any single node, is very general and will probably not need to be changed. The only cases where FIXNOD would not be adequate would involve mixed boundary conditions (i.e., stress is a function of velocity), or free-slip impermeable boundaries that are NOT parallel to x or y; such complexities are "left as an exercise for the student" ! (Hint; study subprogram SLIPBC, which implements a nonlinear mixed boundary condition in an oblique direction, for Neogene models of North America only.) M) Decide what to do about the bottom boundary conditions. If you do not need any subducting slabs to contact the base of your model, then just set parameter IBELOW = 0 in the parameter- input file (e.g., IN9828). If you are modeling North America during any time in the last 85 million years, you can set parameter IBELOW = 1 or 2 and use the program with no changes. Notice that this will also activate routine SLIPBC to create a proto-San Andreas transform margin in Neogene times. If you wish to see exactly what this subprogram will do, copy the MAPPER group of utility programs and run them on an MVS/TSO system (or patch them into your own graphics software). The MAPPER programs will also allow you to modify the plate model, extend it earlier in time, etc. If you are brave, you can use it to create an entirely new plate model for another part of the world, beginning from a digitised map of present-day marine magnetic anomalies! A third alternative is to replace BELOW with a subprogram of your own creation that fills arrays SZZB, TOUCH, and VSLAB by an entirely different method. (See the definitions of the input and output arrays in the Glossary section at the beginning of the Laramy main program for help.) If you do create a new version of BELOW, remember to install this same version in the post- processing graphics programs GDDMCOMP (in the GRAFIC subdirectory) and VERSCOMP (in the VERSATEC subdirectory), and MAPCOMP (in the MAPPER subdirectory) so that these post-processors will be using the same plate model. Otherwise, your post-production plots of isostatic elevation, basal temperature, and basal shear stress will be all wrong! N) Consider rewriting SUBROUTINE SETEDT to provide an initialization of the strain-rates that is appropriate for you problem, or at least is isotropic. This initialization determines how the rheology is linearized for the first iteration of the first timestep. In some problems (basically, orogenies) this turns out not to be very important; the "memory" of this initial condition is quickly lost. However, problems of taphrogeny (horizontal extension) are inherently unstable, leading to a rift. Even the slight lateral heterogeneity introduced by the present SETEDT (which is customized for the Laramide orogeny) can affect the localization of strain in the entire time history! If you want to replace this SUBROUTINE with a simpler version that gives a laterally-homogeneous and isotropic initialization, then I suggest the following: SUBROUTINE SETEDT (NUMNOD,VM,VDECOL, + NUMEL,ERATEC,ERATEM,ESUMC,ESUMM, + XIPC,XIPM,YIPC,YIPM, + XNODC,XNODM,YNODC,YNODM) C C INITIALIZES: C * ESUMC AND ESUMM (TOTAL STRAIN MATRICES, CRUST AND MANTLE) C * ERATEC AND ERATEM (STRAIN-RATE TENSORS, CRUST AND MANTLE) C * VM (VELOCITY OF MANTLE LITHOSPHERE). C STRAIN-RATE AND MANTLE VELOCITY MAY NOT BE ZERO, OR C DIVIDE CHECKS WILL OCCUR IN THE ATTEMPT TO LINEARIZE C THE RHEOLOGY FOR THE FIRST ITERATION. C THE EXACT VALUES USED ARE NOT IMPORTANT -=IF=- SOLUTIONS C ARE ALLOWED TO CONVERGE FULLY; HOWEVER, IN A POORLY C CONVERGED SOLUTION, ANY LATERAL HETEROGENEITY IMPRINTED C BY THESE INITIAL CONDITIONS MAY REMAIN. C THEREFORE, LATERALLY-UNIFORM VALUES ARE RECOMMENDED. C DIMENSION ERATEC(4,7,NUMEL),ERATEM(4,7,NUMEL), + ESUMC(2,2,7,NUMEL),ESUMM(2,2,7,NUMEL), + VM(2,NUMNOD), + XIPC(7,NUMEL),XIPM(7,NUMEL), + XNODC(NUMNOD),XNODM(NUMNOD), + YIPC(7,NUMEL),YIPM(7,NUMEL), + YNODC(NUMNOD),YNODM(NUMNOD) DATA EDOT /1.0E-16/ DO 50 I=1,NUMNOD VM(1,I)=VDECOL VM(2,I)=0. 50 CONTINUE DO 100 M=1,7 DO 90 I=1,NUMEL ERATEC(1,M,I)= -EDOT ERATEC(2,M,I)= -EDOT ERATEC(3,M,I)=0.0 ERATEM(1,M,I)= -EDOT ERATEM(2,M,I)= -EDOT ERATEM(3,M,I)=0.0 ESUMC(1,1,M,I)=1. ESUMC(1,2,M,I)=0. ESUMC(2,1,M,I)=0. ESUMC(2,2,M,I)=1. ESUMM(1,1,M,I)=1. ESUMM(1,2,M,I)=0. ESUMM(2,1,M,I)=0. ESUMM(2,2,M,I)=1. 90 CONTINUE 100 CONTINUE RETURN END O) Review all of the parameters in the input file; check the rheology; etc. Now you should be ready! For the first "real" run, it is recommended that you set the timestep to 0.0 (and end-time equal to start-time) so you are only doing an instantaneous velocity solution and not integrating finite strains (yet). Allow a large number of iterations in the velocity solution (parameter MAXITR in line #54 of the parameter file >=20, and parameter OKTOQT in line #55 of the parameter file no more than 0.001). This should allow the velocities, strain- rates, and stresses to converge. Laramy will make printer-plots of the results for you to inspect. Are your boundary conditions actually working? Is the solution reasonable? P) Consider the values of the two artificial viscosity limits, VISMAX (line #56 of the parameter input file) and ETAMAX (line #57 of the parameter input file). If these are set too low, then artificial compliance will make the whole lithosphere artificially weak (low VISMAX) or will make the crust/mantle coupling too weak (low ETAMAX). On the other hand, values that are too high can cause convergence problems. A high VISMAX determines a "noise level" below which the solution cannot converge. An excessive ETAMAX can actually lead to the misleading appearance of convergence when it has not been achieved (in cases where the crust and mantle are basically not detached). What can happen is that each layer becomes so strongly bonded to the other that its solution cannot change very much from the last iteration, and so we appear to have a converged solution. I suggest that you read what I have written about this issue in Bird (1989; J. Geophys. Res., 94, 3967-3990) and Bird and Kong (1994; Geol. Soc. Am. Bull., 106, 159-174). Next, I suggest that you do a series of short runs to test the effect of each parameter. These runs can be only a single timestep each, and the timestep length can be zero. (Be sure to set MAXITR in line #54 to at least 20, and OKTOQT in line #55 to 0.001 or less, so that you are comparing solutions which are at least approximately converged.) In theory, what you should find is that as VISMAX or ETAMAX is increased, it initially changes the solution (by removing artificial viscous compliance), then there is a plateau where it doesn't affect the solution very much (this is where you want to be), and finally there is a region where a converged solution is hard to get. Q) After settling on your preferred VISMAX and ETAMAX (see above), do some additional tests to find out what number of iterations and convergence level are needed. That is, experiment with the parameters MAXITR (line #54 of the parameter file) and OKTOQT (line #55 of the parameter file). I used 10 and 0.02 (both dimensionless) in the Laramide simulation that is provided as an example, but that was ONLY BECAUSE I COULD NOT AFFORD BETTER CONVERGENCE on the slow computers of that time. Today, I would want 20 (or more) and 0.001 (or less). You should make some runs of one timestep each (and the timestep duration can be zero) to see how much the solution is affected by these values. Remember, OKTOQT is only the non-dimensional change in the solution from one iteration to the next, which is NOT the same as absolute error. If convergence is proceeding monotonically by small steps, then the residual absolute error can easily be 10 times greater. See Bird (1989, J. Geophys. Res., 94, 3967-3990) for extended discussion of these convergence issues. R) Finally you are ready to run an actual history!!!!!!!!!!!!!!!!!! Set the timestep to something like 3.15576E13 seconds (1 million years). Be sure to set the parameter file for a reasonable number of intermediate reports so that you can restart the calculation in the middle if anything goes wrong, and so that you can make graphics from these intermediate timesteps if everything goes right! Good luck! ----------------------------------------------------- compile.jcl (or) compile.aix This are short files of IBM Job Control Language (.jcl) or AIX (IBM's version of Unix; .aix) statements used to compile LARAMY. In compile.jcl, I store the load module as LARMOD(GO), which is a "partioned dataset" LARMOD with the single "member" GO. If you are on an IBM MVS system, you can run this with only minor changes (IRAAGPB and EFF9GPB are two of my account numbers; change them to yours everywhere(!) they occur). If your system is is another type, you may wish to ignore this file completely and start from scratch. However, you should still read the notes below about matrix-manipulation routines. In compile.aix, you can choose either the "debugging" or the "production" version of the compile statement. (Select a statement by removing the ### comment marks before it, and replacing ### comment marks before all others.) The debugging version is not optimized, so it works with the dbx debugger. The production version is 3-5x faster! Note that the compilation makes use of some routines from something called "ESSL". This is IBM's Engineering and Scientific Subroutine Library. If you don't have ESSL, your work will be a litte harder. Read the first page of the LARAMY code for the names of the matrix- manipulation routines that it calls. Then substitute appropriate calls to your local libraries. The most important routines are DGBF and DGBS called from SOLVER to solve the linear systems for the velocities of the nodes; this must be fast and must work with real, banded, assymetric coefficient matrices. If your computer has 4-byte precision as the default, then use double-precision routines as I have. If your computer (Cray?) has 8-byte numbers as the default, double precision is not needed and will be very costly. Unfortunately, different solvers assume different storage schemes for the banded coefficient matrix! To preserve portability, I have always referred to the coefficient matrix as a "vector" array, with a single subscript. The index to this vector (i.e., position in memory) is computed from the two logical (row, column) indeces by a statement function INDEXK which appears in many subroutines. Thus, if you change the banded-equation solver called by SOLVER, then you should also write a new statement function INDEXK and copy it into each subprogram that currently has such a statement (the name is always the same). [Incidentally, the incorrect compilation of LARAMY by some earlier version of VS Fortran seemed to be connected with improper vectorization of loops which invoked this statement function for indirect adressing. ] Less critical are the routines SGEF and SGES used to solve real linear systems with full matrices, and SGEEV used to find eigenfunctions and eigenvalues of a general real matrix. In LARAMY, these problems are never larger than 5 x 5, so it is not critical to use particularly fast or storage-efficient routines. ----------------------------------------------------- runner.jcl (or) OLD_LARAMY.aix This are short files of IBM Job Control Language (.jcl) or AIX (IBM's version of Unix; .aix) used to run the compiled finite-element program, which was produced by the previous job in file compile. -Input file IN9828 is to be associated with Fortran device 1, and is provided. -Input file CORD11N is to be associated with Fortran device 8 (IF the parameters in file IN9828 dictate that an old grid is to be read), and is provided. -An output file corresponding to Fortran device 9 should be defined, because you will usually want to save the values in the essential arrays for plotting, restarting, etc. These files can be very valuable, if they require runs of hours to produce! Under MVS (JCL) one must predefine the file type and size. I define the device-9 dataset with the statement "//FT09F001 DD ...". This output file will receive "card-format" text (no more than 80 characters per logical record), and it must be BIG (about 6000 records per "report" issued by LARAMY. Call it whatever you like. ----------------------------------------------------- IN9828 This is the parameter-input file that subroutine READIN of program LARAMY will read from Fortran device 1. It defines the physical problem, the tactics for solving it, and the initial conditions. It also has some graphics-control parameters following at the end which are not read by LARAMY, but which are read by programs in the other packages. The file is in "card-format" (80 character lines), and some lines are comments. Do not remove the comments, or add any lines, or subtract any lines; this will misalign all the values! You will notice that most lines have a value followed by a comment as to its meaning. The value is read by the Fortran free-format interpreter ("READ (1,*) "), and can be anywhere on the left side of the record. (It must be separated from the comment by one or more spaces, however.) During your first runs, you should check the output carefully to see that the program has interpreted the values as you intended. (One error I once made was to leave a space within a number in exponential notation. I entered "730. E+05" and the computer read this as "730." and ignored the rest. I should have typed "730.E+05".) ----------------------------------------------------- CORD11N This file is a card-format text file (no more than 80 characters per logical record) of packed numbers, which is read from Fortran device 8 by subroutine GOON of program LARAMY and used to fill the essential arrays in LARAMY with initial conditions. (The same type of file is produced by subroutine TAPE of program LARAMY on Fortran device 9, to record the results of the run and allow graphics post- processing and/or restarting.) This particular file represents a guesstimate of the structure of North America 85 million years ago. The interior is plane-layered, thermally equilibrated, and floats isostatically near sea level. The western margin is an active subduction zone, with a 3-4 km high cordillera inland from the trench. The mantle lithosphere grid is offset from the crustal grid, and begins under the eastern flanks of this cordillera. Heat flow is elevated in the altiplano, and depressed in the forearc. If you don't believe it, OK. This is a sample file that can be used for starting up LARAMY and/or for trying out the various graphics packages in advance of your own runs. The file is human-readable, but just barely. ----------------------------------------------------- P9828A00 This file contains the "printed" output that I obtained on Fortran device 6 (default output) when I ran OLD_LARAMY with input datasets IN9828 and CORD11N. You might wish to see if your system can (approximately) reproduce these results before running your own problems. This file will also show you the low quality of the printer- plots, and possibly convince you that you also want to use my graphics post-processing packages (e.g., Laramy2AI)! Since the file has line lengths of up to 133 characters, you should either print it, or inspect it using a word processor that "windows" rather than "wrapping" long lines. Otherwise, the tables and plots will be too distorted to interpret. ----------------------------------------------------- T9828A00 This file contains the "tape" output that I obtained on Fortran device 9 when I ran OLD_LARAMY with input datasets IN9828 and CORD11N. It contains all the values of the essential arrays reported every 5 Ma during the duration of the calculation (17 reports). It can be used with any of my graphics programs to plot various views of the output. One can also refer to any one of the 17 "reports" and use it as an initial condition for restarting the compution, by connecting this file as an input file on Fortran device 8 when running LARAMY. ----------------------------------------------------- T9828END This file is just a fragment from the last timestep contained in file T9828A00 (see description above). Several of my graphics programs require this file in order to determine how the state lines in BASEMAP have changed over time, so they can superpose them on the plots. It has no other function. ----------------------------------------------------- GOOD LUCK !!!!!!!!