Version History of MLOC

Program to locate clusters of earthquakes, using the Hypocentroidal Decomposition (HD) algorithm of Jordan and Sverdrup (BSSA, 71, 1105-1130, 1981). The code is based on Ken Creager's LOC code for single event location. It was originally written in January 1989 by Eric Bergman and has been exclusively developed and maintained by him since then. In particular it has been extensively developed to determine calibrated earthquake locations with accurate estimates of uncertainties for all hypocentral parameters.

Latest changes, in inverse date order:

2020/4/22: Getting a segmentation fault with a cluster that had a lot of events with stations right on top of them. I think it was the section in subroutine delaz (mloclib_geog.f90) that handled zero epicentral distance. I was setting delta=0. when it got to be less than about 50 m. Setting the distance to a very small number (but non-zero) fixed the problem.

2020/4/10: In the Iwaki cluster I ran into a bug in mlocout_ttsprd.f90 where the number of P readings exceeded the hard-wired limit of 30,000. This didn't break anything but produced a large number of warnings. I just set the limit of those arrays to be equal to the parameter that sets the limit for number of readings used (ntmax1).

2019/12/18: I changed the functioning of the RADF command, so that reading (and using) agency and deployment fields to resolve station code conflicts is restricted to specified station codes, not universal. In practice there are seldom more than a couple cases in any given cluster where agency and deployment codes are needed. Making sure that the entire dataset had correct specification of agency and deployment (to match whatever was in the station files) when you only need it for a couple readings from a single station was a nightmare. In fact, the code now uses the agency and deployment fields along with the station code for every comparison, but in most cases the fields are simply blank. RADF just specifies that those fields, both for station files and phase readings, will actually be read for certain station codes. 

2019/12/5 (v10.5.0): Major changes in plotting, related to the adoption of GMT6, which has a nice new option for handling plotting of topography using the "@earth_relief_rru" option in grdcut. I was also running into trouble with the old .grd files when using GMTv6. Some of the GLOBE tiles were ill-constituted and are now breaking stricter requirements. Rather than try to fix my old .grd files I decided to abandon GLOBE and GINA as options and use ETOPO1 (the "01m" option from the GMT server) as the sole option for dem1. The new system also supports high-rez DEMs from SRTM but I have not yet experimented with that. ETOPO1 has a little less resolution than GLOBE but it still works well for basic plotting, and it will be easier to support since I don't need to serve the data file. The argument to the command "dem1" is simply "on" or "off" now. It appears that the current plotting code will still work in GMT5 (with one exception) so that is still permitted. The exception is the call to grdcut, and for that there is a special loop for GMT5 there, referencing the old ETOPO1 file in /tables/gmt/dem/ETOPO. I also fixed a long-time annoyance with the .kml file. The symbols only displayed correctly if the directory is still in the working directory, so it could find the image files in the tables/kml directory. Now each cluster directory acquires a directory called "_kml" containing the necessary images and the .kml file refers to those, so it displays correctly as long as the .kml file is in the same folder with the _kml directory.

2019/11/18 (v10.4.8): Running under macOS 10.15 (Catalina). My old installation of GMT 5 violated the new security requirements, so I had to install GMT again. Version 6.0.0 was just released so I went with that. I installed the macOS bundle (gmt-6.0.0-darwin-x86_64.dmg) from GitHub in the Applications folder. After making the needed changes to my .bash_profile:

export GMTHOME=/Applications/GMT-6.0.0.app/Contents/Resources
export PATH="/Applications/GMT-6.0.0.app/Contents/Resources/bin:$PATH"
export PROJ_LIB=$GMTHOME/share/proj6
export MAGICK_CONFIGURE_PATH=$GMTHOME/lib/GraphicsMagick-1.3.33/config

so that I could call GMT from scripts most of the plots work as before but the ones that use topography (e.g., _base.bash) don't work because of minor syntax changes to grdcut and pscoast. This version of mloc requires GMT v6.

2019/09/13: Fixed a minor annoyance in mlocset.f90 where, in a forward model run (STEP 0), the bounce point file would not be closed and a warning would be issued.

2019/08/27: Fixed a bug in plotting the relative depth phase data for a single station, using a f10.3 format descriptor to write an integer. It crashes the program. It was missed before because it only came up when things were really ugly.

2019/07/14: Fixed a bug in "station_check" to prevent it's actions when a station code is coming out of a supplemental station file; it is intended to apply only to stations in the master station file. Fixed a bug when using command "radf on" where a reading with a deployment field entry would become "missing" if that station-deployment did not exist in the supplemental station file. In many cases that station should be matched in the master station file, even if the deployment fields don't match (because of all the aliasing in the ADSL protocols).

2019/6/10: Changed "run_gmt_script" (in mloclib_gmt.f90) to add S-P plots (type 9) to the set of plots that are not moved to the _comcat folder when ComCat/GCCEL output (command ccat) is requested.

2019/3/26: Somewhere along the line I dropped the depth component of the PERT command from being applied in subroutine starting_hypocenter" (mlocio_mnf.f90). The latitude, longitude and origin time components were being handled correctly. I have added the depth perturbation component back and rewritten that section to simplify the code a bit.

2019/2/6: I modified the code to allow a relocation to be done when latitude and longitude are fixed, e.g., with only origin time as free parameter, or depth and origin time. Formerly, the program would crash in this situation because of a call to smatinv2 (through elips) for the confidence elips.

2019/1/24: Relative depth phase readings will now be automatically flagged with 'p' if their epicentral distance is less than 26¡ or if they are on the list of stations suspected of reporting bogus depth phases. The format of the data file for bogus depth phases (bdps.dat) is changed to the generic (type 3) format for supplemental station files, and the operational epoch fields in that format are interpreted as the date range during which bogus depth phase reports are suspected.

2019/1/22: I rewrote dem1_gmt_5 (mlocout_gmt5.f90) to do a better job with situations where a cluster crosses the boundary between the different grids of the GLOBE DEM. The old logic failed with a cluster that crosses the prime meridian, and I think it would have failed with a cluster that crosses the equator as well. The same thing still needs to be done for the GINA DEM.

2019/1/15: Fixed a bug in rdp_single_5 (mlocout_gmt.f90) for the case when an event was assigned the cluster default depth and it also had teleseismic depth phase data. It was attempting to plot the uncertainty in assigned depth, which for the cluster default depth is ±99.9 km.

2018/12/18 (v10.4.7): I added a command "damp" to control whether the damping of cluster vector changes is made or not. This damping was introduced in v10.4.4 but recent experiences have led me to think it may not be optimal in some circumstances. The default is off.

2018/11/27: I fixed a bug in the new code that listed stations from supplemental station files that are NOT used in the dataset.

2018/11/14: I added a fairly objective estimate of uncertainty for the depth inferred from relative depth phases. The asymmetric uncertainty in depth (at ~90% confidence level) is calculated using the zM test (e.g., Langley, 'Practical Statistics, Simply Explained', Dover Publications, 1970, p 152.), to find the range of depths around the preferred depth in which the z statistic, relative to the value of z at the preferred depth, is less than the value at which the null hypothesis (no difference in mean residual) would be rejected at the 10% level of probability. The search is extended over the full range of test depths, to capture multiple minima which are fairly common at greater depths. The results are printed in the .depth_phases file and plotted in the individual plots of relative depth phases (command RDPP).

2018/11/12: There is a completely new analysis of relative depth phases, mostly implemented in the subroutine "rdp_depth_test" in mloclib_tt.f90. Over a range of depths (presently hardwired as 5-60 km in 1-km increments) each relative depth phase is compared to the theoretical time of all three depth phases (pP-P, sP-P and pwP-P) and smallest residual is selected to add to the sum of squared errors for each event at each depth. Next, the minimum squared error over the depth range for each event is found and the corresponding depth is taken as the best estimate of depth. Finally, for the preferred depth, the residual of the observed relative depth phase against each of the three theoretical phases is printed, with other information that makes it easy to see which depth phase readings should be flagged and which ones should be renamed. Details are printed in the .depth_phases output file. It is important to realize that this estimate of depth does not depend at all on the phase names assigned to the depth phases in the MNF files, but the estimate of focal depth will be influenced by readings that have large residuals with all possible depth phases. It does not depend on the focal depth used in the relocation either. The depth phase plots for individual events (command RDPP) are completely new. Now they show the curve of misfit vs. depth. The RDPP command now also accepts the argument "all" to plot all events for which any relative depth phase data were found.

2018/11/08: I rewrote much of the code dealing with relative depth phases and bounce point corrections. A major change was to convert 'normal' depth phases (i.e., pP and sP) into relative phases when they are read in (read_mnf_13), as long as there is a first-arriving P arrival for them. This also happens for pwP arrivals, which are being handled correctly for the first time I think. This simplifies things a lot through the rest of the code, compared to the old strategy of carrying the original phase and making conversions for relative phase data in several places. There is a new output file (.bptc) when bounce point corrections are made. Corrections are made to orphan depth phases as well as relative depth phases, but the bounce point corrections are only listed for relative depth phases. The default for bounce point correction is "on". Relative depth phases are not allowed to be renamed by the phase identification algorithm; they can only be changed by editing the MNF files.

2018/11/06: Change to the format of the generic supplemental station file (isstn=3) to carry the agency and deployment codes (which are optional). There is a test to catch an attempt to read a file in the old format. The .stn file now carries a list of all entries from supplemental station files that are found in the data set, written in this format so it is easy to cut and paste into a cluster-specific station file.

2018/11/02 (v10.4.6): I made major changes to the search for and reporting of station code and coordinate conflicts in the station files (supplementary and master). The search for duplicate station entries is now done only once, after all supplemental stations files and the master station file have been read. Formerly, the check was made after each entry was read and it was checked against prior entries, but this is opposite of the way that mloc actually decides which coordinates to use. Starting from the top, checking for duplicates is only done for entries in the supplemental station files, and then the check is made against all subsequent entries jn the full station list (supplemental + master). There is now a check (reported in the .stn file) for supplemental station file entries that do not appear in the arrival time dataset. The list of station codes (in the dataset) for which coordinates were found now carries a note if there was a duplicate of some kind found for this code later in the full station list. Other details of the .stn file have been rearranged to improve clarity.

2018/10/30: I added information about author, mloc version, run number and date to the .comcat output file as comment lines. I fixed a bug in the plot routine focal_depth_histogram that did not handle the plot boundaries correctly in some cases.

2018/10/27: I modified the handling of the "cal" command so that multiple entries can be processed without having to comment out the ones you don't want for a given run. The calibration parameters are set by the last instance of the command. I also removed testing for old versions of the command (e.g., "cali") and different parameterizations of the uncertainties.

2018/9/28: There were two places where the psconvert routine in GMT was being called in shell scripts without a prefix "gmt "  This works on my installation but when Bob Herrmann was putting the mloc ecosystem in VirtualBox it did not work without the prefix. Now the calls to psconvert are "gmt psconvertÉ". I also replaced the ETOPO2 DEM with ETOPO1.

2018/8/11 (v10.4.5): I consolidated three commands (lgou, tpou and spou) that were used to output empirical travel times for specific phases into a single command (ttou) that can be used for any phase. The ttou command handles S-P phase exactly the same as the command spou did. Output for other, 'standard' phases is different. It includes more information than the old lgou or tpou output files, and should have everything one might want in developing custom travel time models.

2018/7/11: Fixed a bug in the subroutine 'write_mnf_14' (mlocio_mnf.f90) in writing the hypocenter line. For indirect calibration it still wrote the direct or uncalibrated hypocentral coordinates. This only affected .comcat files. I also added an in-line comment feature for command processing. Before any command is processed it is sent to subroutine 'parse' to extract arguments. Now, the argument string is first inspected for the presence of an exclamation point "!". If found, the last exclamation point and any following text are stripped off before actual command processing.

2018/6/29: Fixed some bugs related to the processing of differential time data. These were introduced in v10.4.0 when I made the changes to use deployment/network codes to resolve station code conflicts, but I haven't worked with differential time data since then, until now.

2018/6/19: I changed the convergence test to apply only to the final hypocenters, not the cluster vectors and hypocentroid separately. The actual criteria are the same as what was used for the cluster vectors previously: 0.5 km for lat, lon and depth, 0.1 s for OT. I kept the original test as an option, controlled by a variable "convergence_test_index" that is set in mloc.f90. In both cases the mean of the cluster vector changes is still subtracted to improve convergence behavior.

2018/6/17 (v10.4.4): I added code in mlocinv.f90 to deal with an occasional situation when the cluster vectors have a significantly non-zero mean in one or more parameters that is then corrected by the hypocentroid but the solution oscillates and prevents convergence. After many years of chasing various theories about how to prevent the cluster vectors from having this character I decided to just remove the mean of the cluster vectors before moving on to the hypocentroid. This seems to solve the problem. Why it happens is still a mystery but it must ultimately reflect a numerical instability or weakness in the dataset. The inversion can still fail to converge if individual events (or the hypocentroid) are unstable.

At some point after releasing the May 21 version I made an additional change to the code, fiddling with plotting (probably the EPA plot) and made a bad error that sent the "dcal_distance" variable to map_boundaries in km, not degrees. This caused the .dcal plot to blow up. Probably the epap plot would have blown up as well, but I did not check it before fixing the bug. So this version is functionally identical to 10.4.3. In chasing down the bug I decided to rewrite some of the routines in mlocout_gmt5. The main change is that I broke the direct calibration raypath plot out of the routine that did the base plot. The two plot types were different enough to deserve separate routines ("base_map" and "dcal_map") and that simplified some of the coding that led to the plotting bug. The map_boundaries routine is also substantially changed (simplified). Now, map_boundaries only calculates the limits of the cluster; the code for calculating what the boundaries of the map should be for different plots is now moved into the respective plotting routines.

2018/5/21: I was wrong. there was still a bug in the windowing. Readings outside the window would still have a weight of 0.01, even with "weig off". That is now fixed.

2018/5/1: When I was fixing the divide-by-zero bug (April 15 version) I over-did it and made a change that created a new bug in the windowing algorithm, such that nothing would get filtered by the window algorithm when "WEIG off" was used. I think it's fixed now.

2018/4/17: I worked on a cluster which caused me to want to be able to turn off the ttsprd term in the calculation of uncertainty for the hypocentroid, i.e., assuming the TT model is perfect. That's how most codes work, so it may be useful for making comparisons with other relocation codes. There is a new command PTTT (perfect theoretical travel times) that does this. Default behavior is the traditional: adding ttsprd to sdread for calculating hypocentroid vnhat in mlocinv. It did not seem to make much difference in the case I tested, but further experimentation is needed.

2018/4/15: I caught an obscure bug that would occasionally cause a crash with a floating point exception, just after returning from the calculation of the hypocentroid. It was a divide-by-zero or possibly an overflow condition in the calculation of vnhat, which involved a division by the weight related to windowing (variable "weight"). It was caused by allowing the calculation of weight (in mlocset) to go too small where it rolls down smoothly from 1 to 0 with a sine curve between wind1 and wind2. Round-off error may have contributed too. Another possible cause was when a reading falls off the end of a phase branch at some iteration and becomes unassociated (phase ID becomes "UNKNOWN"). Now the calculated weight is not allowed to be smaller than 0.1, unless it is explicitly set equal to zero (and flagged) because the residual is outside wind2. There is also a check in mlocinv just before the division to catch any zero-weight readings that somehow slip through.

2018/4/5 (v10.4.3): For EPA plots that cover larger areas (> 5¡) the ETOPO2 DEM is now used, regardless of the choice made with command 'dem1'. If GLOBE or GINA are used, the PDF file gets very large and you don't need that much resolution anyway. I removed all the manual and automatic one-minute-error correction code, and support for the associated data flags ('k' and 'l'). There are still checks for possible one-minute errors in subroutines read_mnf_13 and mloc_set, but they only issue a warning; they do not change anything. In read_mnf_13 the check is based on the 'resisc' variable, the residual reported by the reporting agency. The check in mlocset is based on the residual calculated in mloc after initial processing (phase ID, starting locations, etc). The preferred method for handling one-minute errors (if you don't want to let them get flagged as outliers) is to make a copy of the phase record in the MNF file, flag the original one with 't', change the minute value of the arrival time as needed in the new phase record, and change the Author code to the mloc_author code (from mloc.conf). It is also helpful (to prevent false warnings) to delete the 'resisc' field in the edited phase record if that is what triggered the warning. This is a bit tedious but it preserves the original data and provides the necessary documentation on the change. At the very least, change the author field if you edit the arrival time. I removed the code related to honoring Engdahl's flag ('*') because it was never being used.

2018/4/3: At least one iteration is now required before convergence is permitted, in order to generate the eci values (cluster vector residuals) for cleaning.

2018/3/23 (v10.4.2): Fixed a bug in the station elevation correction (CORR), which caused crustal shear phases (Sg, Sb, Sn) to be corrected with the P velocity, not the shear velocity. I also added a new command (SECV) to let the user specify the velocities used for elevation corrections, and I added extra output in several places (e.g., the .summary file) to emphasize that the choice of making the station elevation correction (or not) also has consequences for the reference plane for focal depth when doing direct calibration. I made some adjustments to the plotting of the 5-km reference circle and hypocentroid confidence ellipse in the base plot. The positioning is done a little more carefully and the hypocentroid confidence ellipse is only drawn for calibrated clusters. I reduced the number of FYI messages that are routinely displayed, such as notices from "station_check" and those about annotations. These are still available if verbose screen display is selected (vscr on). FYIs about the number of stations that fail the date range or are missing from the stations file(s) are only displayed now if the number is non-zero. Added the date of the run to the information in the beginning of the .summary file, and I wrote the info about mloc version, run number, author and run date into comment lines in the beginning of .comcat files for forensic purposes.

2018/3/19: In the EPAP plots, the plotting of coastlines in blue has been removed because they interfere with the visibility of the blue circles.

2018/3/1: Ezgi Karasšzen modified mlocout_gmt5.f90 to add the hypocentroid confidence ellipse to the base plot, combined with the 5 km reference circle. She also rewrote the algorithm for placing the reference circle on the plot.

2018/2/13: After running, all GMT scripts are now saved in a subdirectory '_gmt_scripts'. None are deleted. This declutters the main cluster directory to some extent. I also decided to eliminate the dichotomy between "title" (first thing entered interactively) and the text string supplied to the OUTP command. At one time, early in mloc's development, I thought there might be value in being able to use different terms for them, but I never found that to be true in practice and the evolution of the code and my practice led to the expectation that they would be the same. Now and then I'd mis-type and have to clean up a mess of files that were named in two different ways. So I removed the OUTP command and renamed "title" to "basename". I added the '_rdp_summary' plot to the list of plots that are not added to the ComCat plots file, but I am planning to make major changes to the plots related to depth and relative depth phases. 

2018/2/12: Many more changes to the plots, adding better annotation and legends. Certain plots (epa, rdpp, tt5e, tt5s) are now prevented from appearing in the ComCat plot folder.

2018/2/9: There is now a subdirectory for ComCat output, named '_comcat'. It carries two files, the usual '.comcat' data file and a PDF that carries all plots for that run. The exact contents of the PDF file depend on what types of plots are requested by the user. It carries any of the plots called for using the PLTT command plus the basemap and the dcal raypath plot. It will also carry single-station tt5 plots, as well as any 'rdpp', 'epap' and single event tt5 plots. I have also started improving the labeling of all the plot styles so they are easier to understand for people who are not familiar with mloc output.

2018/1/31 (v10.4.1) I made major changes to the plotting for relative depth phases. Formerly, calling for a "tt8" plot with PLTT would cause an individual depth phase plot (stored in a subfolder) to be created for every event in the cluster for which any depth phase data was found. This was wasteful, since usually we are only interested in the data for one or a few events. Now this command creates a summary plot of depth phases for all events. I tried a new style of plot for this, in which all the relative depth phase times from all events are plotted as time vs distance, with reference lines for the theoretical times of pP-P and sP-P at a range of depths (10-50, increments of 10 km). I'm not sure I will keep the particular style of plot for the summary though. I think a histogram of the current event depths may be more useful in this context. The previous style of depth phase plot for individual events is still available, through the new command RDPP (Relative Depth Phase Plot). The command takes an event name and creates a plot for that event only, but still stored in a subdirectory, as before. Currently the command can be given up to 20 times; the number of instances is controlled by a parameter set in mloc.inc, as usual. I am thinking about changing the style of these plots to the style now used for the summary tt8 plot, but that'll wait for a later update. I also dropped support for GMT4 in this version, as changes to the current environment (GMT5) were not being propagated back to the GMT4 version and they were becoming too different. I also changed the file format for standard plots from postscript to PDF.

2018/1/19: I fixed a mistake in the the summary plot (tt1) where the color was wrong for the S phase.

2018/1/10: I added support for 5-character station codes for SEISAN-formatted files in the routine (redsta) that reads supplemental station files. Normally the station code starts in column 3, but if column 2 is non-blank it reads columns 2:6 for the station code. This was needed for a dataset of OBS data that I received from Frederik Tilmann. Nothing else changes.

2018/1/2: I modified the action of the SKIP command so that it could be used to skip (flag with "s") readings that are read in with a blank phase code. The phase code "blank" (without quotes) in the SKIP command is interpreted this way. The need for this arises because recent ISC datasets often have many secondary readings with no phase ID which usually get identified as late-arriving Sn or Sg phases. They need to be flagged but it's very tedious to do it by hand. Use of this feature will probably kill some readings that are associated to perfectly reasonable phases and may even be critical for a good relocation, so it will be highly advisable to carefully review the consequences of using this feature. Editing the phase line in an mnf file to put the correct (or at least a reasonable non-blank) phase code in columns 24:31 will bring back any desired readings.

2017/11/22: The format of .lres files has changed, adding agency and deployment fields, and therefore the utility program "lres" has changed as well. The previous version of lres will not work with the latest .lres files. The change was not needed for lres itself (because it only cares about the line number in .mnf file), but because the deployment field in the .lres file may be needed when deciding which station-phase groups to investigate with the utility program "rstat". I also added a new command RADF to control whether agency and deployment fields are read. The fields are still compared internally, but RADF can be used to cause those fields to be read in as blank from the .mnf files and the station files, so the behavior of MLOC will be just as it was before v10.4.0, i.e., using only station codes to assign station coordinates. The default is to NOT use agency and deployment fields.

2017/11/17 (v10.4.0). This version is a major update that supports the agency and deployment/network fields of IASPEI's ADSLC station nomenclature to provide an improved method to handle station code conflicts. It is no longer necessary (and it was never desirable) to append extra characters to station codes for this purpose. The agency and deployment fields were present in the MNF format from the beginning, but now they are read and made available throughout MLOC. The latest master station file format also supports these fields and they are now filled in. The majority of them are defined with agency 'ISC' and deployment 'IR', but there are exceptions, the most common being FDSN network codes. If FDSN network codes are used in the deployment field, the agency should be "FDSN". The agency field is used mainly for informational purposes in mloc, and shows up in certain output files, such as ".comcat" files and the ".stn" file. The deployment code plays a very important role, however, as it is concatenated with the station code to provide unambiguous identification of stations, where formerly station codes alone were used. It is important to note that when the master station file is read in, the agency and deployment are assigned for all entries as "ISC" and "IR", even though some entries have other agency/deployment codes, because it facilitates the de-conflicting process. "ISC" and "IR"" are legitimate (and the most general) agency and deployment codes for any station in the master station list. In the event data files, ISC data will normally come through with those default agency/deployment codes (if they have any codes at all). If a reading in an MNF event file has a blank entry for deployment, the processing (for the purpose of assigning station coordinates) is unchanged from previous versions. If the deployment field is non-blank, it is required to match both the station code and deployment in the station list (any supplemental files plus master file). Among the various station file formats supported by mloc, only the master station file format (type 0) and NEIC format (type 5) carry deployment codes. If the master station file format is used for a supplemental station file, different (non-"ISC"/"IR") agency and deployment codes would normally be used (otherwise the station should be in the master station file); they are not changed when read in MLOC, thus making it possible to distinguish these stations from the conflicting master file entries with the same station code. By using deployment codes, it is now possible to keep readings from multiple stations that use the same station code in a cluster. The format of several output files is changed as a result but in most cases it is of no consequence. One exception is ".rderr" files, which are incompatible in both directions between v10.4.0 and earlier versions of mloc. There is a new version of the utility program "rstat" also, due to the change in format of the .phase_data files.

2017/10/13: I changed the output for ".comcat" files, such that the coordinates are printed for all stations used, whether they are taken from the main station file or supplemental station files.

2017/10/11: This release is mainly related to a change in the MNF format for event files, which is now v1.3.3. There is a new record type, "ID record" or "evid record" to explicitly carry one or more event IDs. The most commonly-encountered are NEIC (ComCat) and ISC evids. See MNF documentation for details of the format. Along with this change is a more careful handling of MNF version information from the format record of an MNF file. The codes that create MNF files by conversion from other formats now insert a 3-part version code (e.g. "1.3.3") rather than the 2-part version ("1.3") used previously. Prior to v1.3.3 an event ID code (evid) could be carried at the end of the event record line. Only one evid could be carried for an event (except in comment lines) and there was no field for identifying the source of the evid. That functionality still exists for an MNF file with format "1.3", although this version of mloc will issue a warning in such cases. If the format is 1.3.3, however, evids (and their source code) can only be read from ID records. Although the MNF format v1.3.3 allows up to 40 characters to be used for an evid, only the first 10 characters are used by mloc. Their only use in mloc is to assist in associating events correctly with entries in an hdf file from a previous run, but they can be extremely useful when integrating an mloc cluster with a catalog database. The first ID record encountered will be taken as preferred, or the preference can be controlled explicitly with the standard usage code ("=") in column 3.

2017/09/20 (v10.3.4). I added a new command ÔnsmdÕ (standing for NEIC station metadata) to control a new tool for managing station codes and coordinates. If, after searching the master station file and any supplemental station files, there are still station codes in the dataset for which no coordinates have been found, a further search can be made through a file containing the entire list of station coordinates used by the NEIC metadata server. This file, named Ôneic_stn.datÕ, is stored in the tables/stn directory. The filename and path can be changed in the main program initialization section. The file is in the ÕNEICÕ format (isstn=5) and the search output is a list of matched entries in the .log file, from which it is easy to construct a new supplemental station file for a given cluster. It is important to understand that any station coordinates found through this search will not be used in the present run; you must create a supplemental station file in order to use this information in a later run. It is impractical to use neic_stn.dat directly as a supplemental station file because it carries over 200,000 entries. This command is useful for clusters that contain a significant amount of data that was acquired from the NEIC ComCat server, because many of the station codes in use by US regional networks have not been registered with the International Registry at the ISC and there are also many code conflicts, with the IR and also between different regional networks. The default is off.

2017/09/09: When the former master station format is used for a supplemental station file (isstn = 9) the sta_author variable was being set blank. Now it uses the supplemental station file index number, like other supplemental station files. I also changed the format of the list of station used in the .stn file so that the first 53 characters of each line are the Simplified Format (isstn = 3), to facilitate making a new supplemental station file for a few events out of a much larger file. The station author identifier is printed in column 55 and beyond, so it wonÕt be read by mloc if the whole line is used for a supplemental station file.

2017/08/08: Corrections to mlocout_gmt5.f90 to add Ògmt Ò in front of more of the GMT commands, especially the gmtset command and the commands like grdcut used for the plotting of topography. Versions of GMT since v5.2 require them, but older versions did not. It seems to be backward-compatible, however, so the new form of the commands does not break older versions. I also added a test in the plot for near-source TT to handle the case where teleseismic distance limits are used for the hypocentroid. Now it defaults to a maximum of 3.0 degrees in that case.

2017/07/07: Several minor bug fixes and sanity checks related to supplemental station files, to deal with a few issues discovered by mloc users. Added a check to make sure that a format index (isstn) has been supplied in the first line of supplemental station files. The variable is initialized to zero by the compiler so the code would assume type 0 format (new master file format) if no index was read (only a return character in line 1). This would work if the file actually was in that format, but it is a sloppy way to do things. Now there is an explicit initialization of the variable to 99 and a warning if no valid value for isstn is read from the supplemental station file. Also added a check that the SSTN command is not being used to call the master station file. I modified the format of the .stn output file for the section that lists station code conflicts and discrepancies about coordinates when there are supplemental station files. The order is reversed so that the first entry encountered (i.e., the one that will be used by mloc) is listed first. I also fixed a bug in which the ÒresultÓ field was not re-initialized after a non-null result so that if there were multiple entries of the same station in the supplemental file with slightly different coordinates (which often happens when Will Yeck pulls coordinates from ComCat) and the first one has a non-null result, the result field for the second instance would carry the result of the previous entry if it had a null result. I also modified redsta.f90 to check for those multiple stations code entries and issue a warning.

20170525: Changes to the calibration modules (cal_shift.f90 and dcal.f90) and mlocout_bloc.f90 by Ezgi Karasšzen to refine the use mloc calibrated locations as priors with BayesLoc. ÒI changed the default cluster depth error to 5 km (it was 10). Location error is now taken as the radius of the equivalent circle of the confidence ellipse, in km. There were some additional things that were outputted and not required by the Bayesloc prior input file, so I took those out (e.g. depth code, calibration level).Ó 

2017/02/24 (v10.3.3). Added optional fields for the date-on and date-off to the format for SEISAN station files (type 2). This format is often used for temporary networks with station codes that conflict with codes in the main station file, so use of the date-on and date-off feature helps resolve conflicts.

2017/01/18: I changed the near-source plot (tt5) so that it handles Pb, Pn, Sb, and Sn phases in addition to Pg and Sg. I usually donÕt use Pb or Sb but it could happen, and Pn and Sn are needed when an event is below the Moho, which often occurs in subduction zone clusters.

2016/11/02 (v10.3.2). I re-arranged the format of the .comcat file, adding some standard information (number of events, date range, calibration type, etc) as comment lines, so this info does not need to be repeated in the commentary file. I also moved all the code that writes the .comcat file into a separate subroutine (mlocout_comcat). Previously, it was spread out between the main program and mlocset but with the new content the code was getting cumbersome enough to deserve its own package.

2016/05/30: I updated mlocout_gmt5 so that all GMT scripts are written with the new (in GMT5) syntax where Ôgmt Ô should precede the module call (e.g., Ôgmt psxyÉÕ instead of ÔpsxyÉÕ. The old syntax still works in at least some installations, but apparently not all.

2016/05/10 (v10.3.1). I cleaned up the processing of differential time data in mnf v1.5 format. Was not keeping track of line numbers correctly when there were events in the data file that were not in the cluster, so the use of rstat to find out which line to flag was screwed up. Even if there were no missing events the value of ÒidiffÓ in output files was off by 1 because it did not count the format record, and it would be off also if there are any comment lines in the MNF file. I fixed a bug in the differential time code where it did not calculate TT for the Lg phase. Added the ÒskippÓ function in read_mnf_15 to handle some weirdo phases (e.g., Rg) that have shown up in differential time datasets. I removed the designation for ÒprimaryÓ readings. It was somewhat useful when I used the FFB format with its concept of reading groups, but it is often irrelevant or incorrect with the MNF format and the datasets I am getting from the ISC now. These frequently lose the order of first-arriving and secondary phases.

2016/05/05. I disabled the testing for duplicate readings in the code that reads differential times data. When you have a dataset that includes readings for the same station-phase-events triplet from two different authors the second one would get flagged as a duplicate. I think duplicate testing is not really needed in this case. I also set a default value of 0.1 s if no uncertainty is supplied in the differential time dataset, and I implemented a minimum reading error of 0.05 s for differential time data.

2016/03/29. I added a new command ÕshclÕ (set hypocentroid convergence limit) to allow the user to change the convergence criteria for the hypocentroid. This is useful in certain situations when the hypocentroid gets into an oscillatory state, typically between origin time and epicenter, and fails to converge. This is especially common when using S-P data to constrain the hypocentroid. Having a converged solution is helpful in exploring diagnostics to try to suppress the oscillation, or it may be concluded that a somewhat larger-than-normal convergence limit is sufficient for the cluster. I also did some clean-up to the mlocout_bloc.f90 module, especially with the algorithm for depth uncertainty. No change to the ÔblocÕ format though.

2016/03/20. Major rewrite of the input for differential time data, based on the new MNF v1.5 format. All differential time data must be in this format, and the subroutine for reading it (read_mnf_15) is in mlocio_mnf.f90.

2016/03/10 (v10.3.0). I added EzgiÕs code to implement the ÔblocÕ command to control creation of a special output file (.bloc) for import into BayesLoc. I also fixed a minor bug in mlocio_mnf.f90 (in subroutine write_mnf_14) where the hypocenter line for the Ô.comcatÕ output file still carried the usage code in column 3. That was dropped from the MNF v1.41 format specification because there is only a single hypocenter line in this format.

2016/01/19. I added a dashed line at the epicentral distance limit for the hypocentroid to the Õnear sourceÕ travel time plot (tt4).

2016/01/06. I fixed two minor bugs. I added code to catch cases where the old-style parameterization of the ÒcalÓ command (giving terms of the covariance matrix) might be encountered. This sometimes happens when I re-use an old command file when updating a cluster. Since those often were of the form ÔX 0 XÕ, i.e., circular, it now causes a crash in cal_shift.f90 when a confidence ellipse with zero semi-axis length is processed. Non-circular confidence ellipses given in the old parameterization will probably not crash but they will give bogus results. I may need to add more code to check for this. I also fixed a bug with the code for bounce point topography corrections for depth phases. It would crash in iupcor for an event with zero depth. Now a check is made and the correction is not made for zero focal depth, even though a depth phase would still exist for non-zero topography in that case.

2015/12/16 (v10.2.9). I added Bob EngdahlÕs code to calculate corrections to teleseismic depth phases for the topography/bathymetry at the bounce point, and also calculate a travel time for the pwP phase. This is only done when water depth is greater than 1.5 km. swP is not supported. Topography/bathymetry is taken from a version of the ETOPO5 DEM that is stored in the /tables/tau-p directory. The application of these corrections is controlled by a new command ÔbptcÕ (default ÔonÕ).

2015/12/11. Ezgi Karasšzen wrote code to implement a standard plot of empirical path anomalies, using the new command ÒepapÓ (which also has an alias ÒezgiÓ). Multiple plots can be made, for different phases, different distance ranges, and with or without raypaths delineated. The symbol size and color at the station indicates the sign and magnitude of the empirical path anomaly. I took the ÕshelÕ command out and made the selection of a shell script for GMT scripts an optional element of the configuration file.

2015/12/8 (v10.2.8). I redesigned the mloc configuration file to make it more flexible. Now it uses keywords at the beginning of each line to define what parameter is being set. Therefore the order of lines does not matter now. Also I set default values for all of the configuration parameters in mloc.f90. If you set up the ÒtablesÓ directory structure as recommended, many of these can be left as the default values. In that case, itÕs likely that only mloc_path and author would need to be set in mloc.conf. I generalized the gmt_path parameter to one for the color palette tables (cpt_path) and one for the DEMs (dem_path). Within the parent DEM folder there should be separate subfolders for the global DEMs that are supplied with mloc (GLOBE, GINA, ETOPO).

I split the ÔtopoÕ command into two new commands, Ôdem1Õ and Ôdem2Õ to permit a custom high-resolution DEM to be used for plotting topography in the plots (i.e., the baseplot and related plots) that cover a smaller area. The global DEMs that IÕve been using are fine for the _dcal plot that covers a larger area, but often are too coarse for the others. The new Ôdem1Õ command duplicates the functionality of the old ÔtopoÕ command. If a Ôdem2Õ command is issued the referenced DEM will be used for the baseplot and related plots. If Ôdem1Õ is given and Ôdem2Õ is not given, everything will be plotted as in the past, with the selected global DEM.

2015/12/5. I added code to mlocout_hdf.f90 to write a list of simplified hypocentral parameters to make it easier to import into a document. There is substantial amount of new code to implement a new command (EPAP) to make a special plot in GMT that shows empirical path anomalies for selected stations, but it is not finished yet. 

20150831. Added a check in mlocinv.f90 for the case when the number of connected readings for the cluster vector of an event is less than the number of free parameters. This usually only happens when first setting up a cluster, when the starting location for an event is particularly bad or the depth is changed significantly without adjusting the OT. Sometimes unreviewed events from the ISC Bulletin have this characteristic, when nearly all readings have big residuals (probably mis-associated). Whatever the cause, it leads to a floating point exception in the SVD.

20150731 (v10.2.7). I fixed a bug in cal_shift.f90, regarding the radius of doubt in indirect calibration. Two estimates of the radius of doubt are made, one based on a test of the null hypothesis that all residual calibration shift vectors have zero length, and the other on coverage statistics, i.e., given the number of calibration events, is the number of events that are ÔcoveredÕ by the 90% confidence ellipse statistically reasonable? In both tests, the radius of doubt is incremented until the test is ÔpassedÕ but in the coverage test, if there are fewer that 10 calibration events then ÔpassingÕ the coverage test means 100% coverage. So the radius of doubt was being over-estimated. Now, if there are less than 10 events the radius of doubt is taken from the test of zero-length residual shift vectors. If there are 10 or more calibration events it goes the same as before.

2015/7/20. I fixed a bug in subroutine delaz (in mloclib_geog.f90) which would throw an exception if the two test locations were identical. The chances of this are normally very small but it is guaranteed to happen in mlocout_summary.f90 when a single calibration event is used and the input location is the calibration location (in the calculation of ÔCalibrated-INPUTÕ). I ran into this with a cluster that was using a PNE as calibration event.

2015/7/5 (v10.2.6). I fixed a bug in the tt5 plots of individual events, where the theoretical TT curves would be drawn for the hypocentroid depth. Now they are drawn for the focal depth of that event. The routine Òwrite_mnf_14Ó is changed to write the v1.41 format, which carries the TT residual with 2 decimal places of precision. I noticed that both write_mnf_13 and write_mnf_14 did not consider the case of indirect calibration when writing the TT residual (dts) and that is now fixed.

2015/6/20 (v10.2.5). I made it possible for Lg phases to be re-identified as Sg and vice versa. There is now a minimum epicentral distance for the Lg phase, with a default of 2.5¡. It can be specified as a third argument to the ÔlgttÕ command.

2015/5/31 (v10.2.4). I tried to clean up some of the logic around the filtering of readings for various reasons. There were some anomalies in the output files for Òbad dataÓ for .dcal_phase_data files that I am trying to fix. The windowing algorithm is now modified so that for distances less than ÒwindloclimÓ the window is expanded by a factor of 2. This helps keep important readings in the direct calibration problem when the starting locations are poor. The WIND command is altered so that windloclim can optionally be specified as a third parameter. Default value is 1.2 degrees. I fixed the subroutine Òxsec_gmt_5Ó in mlocout_gmt5.f90 so it now works correctly with GMT5.

2015/5/25 (v10.2.3). I removed the choice to do some plots with reduced velocity. The choice to use data flags was replaced by the new command ÔflagÕ (default is to use data flags). Cleaned up some of the user interface, e.g., adding a prompt (Ò:Ó) for additional interactive commands.

2015/5/21. Minor change to the .kml format. I added the event number to the Òclick dataÓ information box for each event, in addition to using it for the Placemark name. I made some changes to the format of the .sp file (S-P data) to make it easier to use awk to extract values,including moving the fcode flag to the end of the line.

2015/5/5. I cleaned up the analysis and display of duplicate stations in mloclib_stations.f90. No change to the actual processing.

2015/4/17 (v10.2.2). I changed the phase reading usage flags for .comcat files after talking to Harley. Readings that were used for either the hypocentroid or cluster vectors (or both) are still Ò+Ó. Outliers are still flagged ÒxÓ. Readings that were not flagged as outliers but were not used in the relocation are flagged Ò-Ò. This can be the case for singlet station-phase readings, or readings outside the distance range defined for the cluster vectors, or readings that were skipped (SKIP command) for some reason. As before, duplicate readings are not recorded in the .comcat file, nor are readings with problematic phase names or readings associated with stations for which coordinates are missing or stations with suspected timing problems.

2015/4/16. I rewrote some of the code around the CCAT command. The argument to the command is now the filename (assumed to be in the data directory) of the commentary file. The former Òcomcat_idÓ is now simply the title of the run. I also fixed a bug in write_mnf_14.f90 that caused the event ID to be mis-printed for event numbers greater than 99.

2015/2/26 (v10.2.1). I added the command ÔoldrÕ to make a new output file (Ô.oldrÕ), similar in format to the .phase_data file, but carrying only readings with epicentral distance in a defined range (OLDR = Output Limited Distance Range). This is intended mainly to aid in analyzing phase IDs around the Pg-Pn cross-over distance, but it could be helpful in other contexts as well. I also fixed a bug that resulted from my recent modifications to the function ÔbetaiÕ. The revised code would call ÔoopsÕ and stop if a = 0 in the call to betai. I did not think that could happen with mloc, but in fact it could, from the subroutine rdbt_test2 in cal_shift.f90, when there is perfect coverage of calibration events. The fix is made in rdbt_test2, where that case is handled explicitly now.

2015/2/19 (v10.2). This is the first version of mloc that supports GMT5. Testing was done with v5.1.1 as installed by fink. GMT4 is still supported as well. GMT5 is now the default in mloc.f90, but the choice is really made in the userÕs mloc.conf file. No further development will be done in support of GMT4. I made major changes in the way events can be killed. The commands ÔmembÕ and ÔkillÕ can both kill events now. Use ÔmembÕ with the argument ÔkillÕ to kill a single event. Additional text can be included after the ÔkillÕ argument to explain the reason. Use the ÔkillÕ command to kill blocks of events, with the ÔonÕ and ÔoffÕ arguments (i.e., Ôkill onÕ and Ôkill offÕ). This is not backward-compatible; command files that use the kill command in the traditional manner must be rewritten, but it is a trivial exercise. Multiple blocks of events can be killed. The two methods can be combined in the same command file and the kill-on/kill-off block can include events killed by the memb command. I added the use of an external file that carries phase-specific default reading errors (maximum 20). The data file is named Òpsdre.datÓ and lives in a folder named ÒspreadÓ in the ÒtablesÕ directory. The format can be found in mloc.f90, where it is read. If it does not exist, processing proceeds as before. These values can still be over-ridden by reading an empirical reading error file (.rderr) from a previous run, and they are subject to the other forms of over-ride as well (see subroutine ÔreaderrÕ in mloclib.f90). In this version I fixed errors (introduced on 2015/1/3 in mloc v10.1.2) in my implementation of EngdahlÕs flagging of certain stations that are known to have timing problems during certain time periods (subroutine Ôstation_checkÕ in mloclib_stations.f90). I had the logic exactly reversed on the time periods for stations FCC, TUC, WIT and YKA. Oops!

2015/2/17 (v10.1.3). To make it easier for students to use mloc, I am trying to use gfortran to build it. I installed the latest stable version of gfortran (v4.9). I am hoping that it will work acceptably with my code now; I have heard that it has been improved significantly in recent versions. The v10.1.2 source code compiled successfully, but with many warnings that revealed poor coding practices, especially in the area of controlling type conversions and doing tests of equality between floating point numbers. I have fixed many of them in this version, but there are still a few remaining in packages (e.g., dsvd2 and libtau) in which it is not easy to evaluate the tolerances that are needed. With this release I am reorganizing the development of mloc a bit. I will try to keep a single codebase that can be compiled with several compilers. For the moment I will work with gfortran and Absoft, but I hope to bring back the Intel compiler too. Executables should be flagged with the corresponding compiler (mloc_g, mloc_a, mloc_i) with version numbers appended (e.g., mloc_g1013). There are no major changes in the code from v10.1.2 with this release.

2015/1/3. I added a subroutine Òstation_checkÓ in package mloclib_stations.f90 to correct various station-specific problems, mostly related to known timing problems during certain periods (which are flagged with fcode = ÕtÕ) and cases where a station was moved a significant distance without registering a new station code (which was fixed at a later date). In this case the incoming station code (stname) is changed accordingly. Most of the supporting information comes from Bob Engdahl. This subroutine is called just after reading in a data file in subroutine read_mnf_13.

2014/12/4. I extended the ÔrelsÕ command to take a third argument for the distance range in which to apply the specified reading errors. Usually you would only want to use these within a fairly short epicentral distance that is being used for direct calibration, but you may have Pg and Sg readings out to several hundred km.

2014/12/3 (v10.1.2). Added command Ôtt5eÕ which makes a local distance TT plot using only the readings from the specified event. This should help with analyzing problem events. The command can be issued multiple times; the current limit is 10. Plots are placed in a subdirectory. I also arranged things so that when you do a forward model (step 0) you will get the various flavors of local distance TT plot (tt5, tt5e, tt5s) if you have requested any. I find that forward modeling is usually done when some event(s) is giving trouble, and these plots help reveal the nature of the problem.

2014/11/26. Fixed a bug that would cause a crash in the function ÔellipÕ if the calibration shift in indirect calibration caused the focal depth of an event to go negative. Now the depth is reset to 1.0 km in that case, and a warning is issued.

2014/11/25 (v10.1.1). Added command ÔsubcÕ to select a subset of events that meet criteria related to suitability for direct calibration (number of readings within a certain distance and connectivity). The main body of a command file for the selected events is written in the .log file, to facilitate analysis of the subcluster.

2014/11/14. Fixed a bug in mlocio_mnf that did not cause a problem under the Absoft compiler, but which crashed hard under the Intel compiler. It only occurred when verbose screen output was requested (vscr on), and was caused by writing more than 132 characters to the ÔmsgÕ variable when reading a hypocenter line. There was a chance of a similar problem occurring when dealing with annotations, but I changed the code to avoid it.

2014/11/11 (v10.1.0). I added the command Òtt5sÓ which makes a local distance TT plot (tt5) for a single station. This is helpful for investigating phase ID issues for stations at distances near the Pg/Pn and Sg/Sn cross-overs. The command can be issued multiple times; the current limit is 10. Plots are placed in a subdirectory. I removed the command ÒrdepÓ and associated code, because the concept of revising starting depths according to relative depth phases did not prove to be useful.

2014/11/4: Slight tweak to the algorithm in mlocio_mnf.f90 that matches up events to an HDF file. In the old version, there was a problem if a new event was added to the cluster (thus, not in the HDF file) that has an OT less than 99 seconds different from one of the existing events, it would consider that a match and use the wrong hypocenter as a starting point. That resulted in all the readings being flagged as gross outliers and then the location would blow up with a floating point exception. Now the window is ten seconds for a match. It is still possible for this bug to occur, but less likely.

2014/10/21 (v10.0.9). I discovered and fixed a bug in cal_shift.f90 that led to an incorrect confidence ellipse for indirect calibration. The augmented GT covariance matrix (agcv) was not being fully defined. The epicenter cross terms (agcv(1,2) = agcv(2,1)) where being left initialized at 0. so the final confidence ellipse always had zero rotation (azimuth). Also fixed a bug in the subroutine Òwrite_mnf_13Ó that caused .datf_mnf files to have incorrect hypocentral parameters. In the case of indirect calibration it was still writing the origin time, latitude, longitude and focal depth from the standard hypocentroid estimation (either uncalibrated or direct calibration), i.e., the unshifted hypocenter. The uncertainties in hypocentral parameters were correct however.

2014/9/23 (v10.0.8). I changed the syntax of the cal_ command, so that it takes a confidence ellipse and uncertainties for depth and OT rather than elements of a covariance matrix. I also added output to the .log file that formats the corresponding cal_ commands for all events in a calibrated relocation, either from direct or indirect calibration. This facilitates using these events as calibration events in a future relocation.

2014/9/21. I fixed a bug that prevented the format record from being written to .dat0 and .datf files.

2014/8/26 (10.0.7). I added a new command CCAT to control the creation of a new output file (Ô.comcatÕ) that will be used to import the results into the USGS/NEIC COMCAT catalog server. At this point the output file is identical to a .datf file. The modifications needed for COMCAT will be added in future versions.

2014/7/1. I cleaned up the handling of arrival time precision in mlocio_mnf.f90, so that the input data are checked for consistency (and iptim is set appropriately if there is inconsistency) and output data (write_mnf) are written with the correct number of decimal places. Before everything was printed with three decimal places.

2014/6/26 (10.0.6). I changed the MARE command so that the order of arguments is local distance then other distances, and added a third parameter for reading error of teleseismic depth phases. I also changed the output of the depth_phases file so that readings out to 100 km are printed, because these readings (if direct crustal arrivals) can help constrain focal depth through changes in OT that are coupled with different assumptions about depth.

2014/6/20 (10.0.5). I rewrote much of the code related to setting starting locations. Functionality (in terms of which sources are preferred) has changed. The major change is that I took out the logic that would try to set up an indirect calibration if events classified as GTXX were read. The code was getting too complex, and this situation does not arise often enough to warrant automatic handling. GTXX events are noted with an FYI when read but calibration events can only be declared in a command file with CAL_ commands. There were some other changes as well, especially concerned with setting focal depth. Many variables were renamed, to maintain a standard terminology. I also took out support for MNF v1.1 and v1.2 data formats. I didnÕt want to take the time to modify the subroutine read_mnf_11_12 to be consistent with the changes made to subroutine read_mnf_13 and itÕs better to convert such files to MNF v1.3.

2014/6/18 (10.0.4). New version of the SKIP command allows specification of station code, phase code, and author code to be skipped. Commands AUSK and PHSK are gone. This allows more control of which readings are skipped. Wildcard still makes it easy to duplicate the effect of the former individual commands.

2014/6/14. I fixed a minor bug in subroutine hdf_match (in mlocio_mnf.f90) that would cause a floating point exception sometimes when using locations from a previous run (command rhdf). The problem was that a test variable (adt) was not being initialized. If an event was matched in the HDF file it would be set, but if the event was not matched, there was still a formatted write statement that referenced it. Sometimes it would work, but sometimes it would blow up. Now it is initialized but I also changed the write statement so that it is not referenced. It makes more sense to reference the actual origin time in any case.

2014/5/28 (v10.0.3). I removed the option tt_corr = 2 from command ÒcorrÓ, because I no longer use EngdahlÕs master station file that carried the patch correction data. I never used the patch corrections in any case, as they are not relevant to the way things are done in mloc. I also changed the output of the phase_data file. The column that used to carry ÒTT SPRDÓ (ttsprd(iev,i) now carries ÒSTA CORRÓ (stacs(iev,i). The TT spread values can easily be looked up in the .ttsprd file, but there wasnÕt any place to inspect the station corrections and I thought the station corrections were more important, especially for data being used for direct calibration. The format of the line is exactly the same. I fixed a minor bug in the logic to note possibly incorrectly flagged readings (Ò?Ó) in the Òbad dataÓ section of a .phase_data file, by including the baseline offset of the TT curve for the appropriate branch.

2014/5/13 (v10.0.2). Discovered a poor decision about the way I wrote the code to use the values of TT spread and offset from a previous run. Until now, if weighting by reading error was turned off (Òweig offÓ), the value of spread would default to 1.0 and offset would default to zero for all phases and distances, even if a .ttsprd file had been read. This can cause problems when doing direct calibration with weighting turned off, as is typical during early stages of analysis. If there is a large TT offset for teleseismic P, for example, it is easy to flag readings that should have been kept, leading to a biased teleseismic dataset. I think this may be why certain clusters have had problems with convergence. Now, if a .ttsprd file has been read those values of offset and phase-specific spread will be applied even if data weighting is turned off. If a .ttsprd file is not read, the code uses the default values of phase-specific spread specified in subroutine ttsig (all TT offsets = zero in this case).

2014/5/6 (v10.0.1). Cleaned up the handling of station files (subroutine redsta in mloclib_stations.f90). New master station file uses isstn = 0. This format can also be used for a supplemental station file. The former master station file (and supplemental station files using that format) is supported as isstn = 9. Therefore, older supplemental station files using isstn = 0 will need to be edited to change the value of isstn if they are to be used again. All station files now have the concept of comment lines (Ô#Õ in column 1).

2014/4/28 (v10.0.0) Added support for the new master station list format (isstn = 9). There is still support for the previous master station list, but the file needs to be edited to add a Ò0Ó in the first column of the first line, to tell mloc which format to use. As before the choice of master station files is made in the configuration file. This version re-introduces support for differential time data (command ÒdiffÓ), but the format is not settled yet. In this version it is reading (mlocio_diff.f90) a format I made for a specific data set from Mike Begnaud. There is still work to be done in deciding how to represent differential time data in the MNF format. Even if the data is not read from an MNF-formatted file, there is the question about how it is represented in a .datf file.

2014/4/23. Fixed a bug in mlocinv that caused the array bound of ntq to be exceeded when there are a large enough number of instances of a particular station-phase. The array limit had been set to the maximum number of events but there are now data sets that have multiple instances of the same station phase in many events and that original limit can be exceeded. I created a new limit parameter in mloc.inv (ntqmax) that can be set to a multiple of the maximum number of events. Presently set to nevmax*2. I also added a test to check that the new limit is not being exceeded.

2014/3/24. Added commands ÒtpouÓ and ÒtpttÓ to handle T-phases for location, in the same way that Lg is handled. 

2014/3/24 (v9.9.9). Making the transition from GMT4 to GMT5, will support both versions for a while. I added a line to the configuration file to specify which version is to be used. Then I pulled a few utilities out of mlocout_gmt.f90 that are independent of the GMT version and put them in a new library package called Òmloclib_gmt.f90Ó and created a near-copy of the modified Òmlocout_gmt.f90Ó, called Òmlocout_gmt5.f90Ó. The only difference (for the moment) is that all routines in mlocout_gmt5 have a Ò5Ó appended to their name. Obviously, the actual codes in mlocout_gmt5 will need to be modified to conform to GMT5 standards, but until I finish that task, I can still run with GMT4. I also made some improvements to the plot for relative depth phases. It now automatically adjusts the dimensions of the plot to suit the depth of the event, and I plot the theoretical TT curves for the minimum and maximum epicentral distance, in addition to the average.

2014/3/18. I fixed the code dealing with phase names to handle P* and S* properly - change to Pg and Sg. Previously, it was possible (and common) for P* readings to be re-identified as Sg, because the code that decides what is a P wave didnÕt recognize P* at all.

2014/3/12. I moved the earthquake icons for kml files into a subdirectory of the tables directory (tables/kml).

2014/3/3. I added the code value ÒtÓ to flag readings for which the timing is suspect. This is mainly used for the readings that are converted to S-P.

2014/1/22. I made many small edits to clean up the logic and clarity of the code related to depth constraint. Fixed a bug that printed 0.1 for depth uncertainties (instead of blank) for events that had no entry for depth uncertainties in the input file or the command file.

2014/1/17. If no magnitude record is declared as the preferred by Ò=Ò usage flag, the first encountered magnitude record will be used.

2013/12/9 (v9.9.8). I added support to take calibration events from the input files if they are in MNF 1.3 format and have a calibration code of GT00 or GT01. As per the GTCNU standard, I do not support the concept of GTX with X > 1. Such events are now read as calibration events and written back to the .datf file with the same GTXX hypocentral parameters. Obviously, the fields for OT uncertainty, depth uncertainty and epicenter confidence ellipse must be supplied in order for the location parameters to be used for location calibration. When such events are written back to the .datf file, the mloc hypocenter is added as a new hypocenter record, but the GTXX hypocenter is kept as the preferred one. If the "calX" command is used interactively or in a command file to declare the same event as a calibration event (possibly with different location parameters or covariance), the parameters in the input file are over-ridden for relocation calibration, but the original calibration parameters from the input file are written back to the .datf file. There is no way for mloc to declare a GT-class event; they can only be brought in from the input file. Even if an event is calibrated in mloc with uncertainties that would qualify as GT00 or GT01, the software will declare it as a CH00 or CH01 event in all output files, including the .datf file. Conversion to GT-class must be done outside mloc.

2013/12/4 (v9.9.7). Improved support for MNF v1.3, especially in the hypocenter record, where I've added the fields for uncertainty in origin time and the epicenter confidence ellipse parameters in the .datf file. I altered the MNF v1.3 format as concerns the order and formatting of the fields for the confidence ellipse. I also added a new command LONR which allows the user to specify the desired longitude range (-180 to 180 or 0 to 360). The default is -180 to 180, but the other option is needed for a cluster that spans the antimeridian (180¡ from the prime meridian). I fixed a couple of bugs related to depth codes, setting the starting depth, and writing MNF v1.3 output to a .datf file. If RHDF is used to get starting locations from an hdf file, depths that have the code 'c' (cluster default depth) are over-ridden by the DEPC value in the current command file, but depths with other depth codes are honored, unless they are explicitly over-ridden by a DEPX command. When mloc is run with a free depth solution, the depth code written to the .datf file (hypocenter record) is 'm', not the depth code used to constrain starting depth. I fixed a minor formatting bug in the magnitude record written to a .datf file, and a bug relating to the evid and event_comment fields of the event record when the input data file is in MNF v1.1 or v1.2 format.

2013/12/2. I added the near-source travel-time plot (type 4).

2013/11/26. The .hdf file was carrying the number of readings lost due to large residual, i.e., those that show up in the 'bad' section of the .phase_data file with 'PRES'. What I want here is the number of readings flagged as outliers (code = 'x') so I changed the code to do that. After all the cleaning, most of the 'PRES' events should be flagged with 'x'.

2013/11/26. After dropping the 'free parameter' flags in v9.9.6, I realized I needed a way to indicate in the .hdf files whether or not a solution was a free depth solution, so I added a "free depth" flag immediately after the depth constraint code. It is 'f' for free depth, blank otherwise. I also changed the output of magnitude to leave the field blank if there is no magnitude, rather than writing '0.0'. I rewrote the code related to the travel time plots, changing the association between index and plot type to something more rational. I also added code for a new plot for near-source data, but I have not yet put in the plot code. So there are now 9 plot types.

2013/11/23 (v9.9.6). I made a change to the HDF file format, replacing the old fields concerning free parameters, calibration status and location level with the recently-developed "calibration code" (GTCNU). I changed the output of the default values for depth uncertainties in the .hdf files and .datf file: now the fields are left blank, instead of printing '99.9' or whatever the default value happens to be. I also modified the subroutine that writes MNF format to the .datf file so that all the old hypocenter, depth, and magnitude lines are carried forward. The new hypocenter line (with '=' usage flag, of course) is written first. The usage flag for any subsequent hypocenter and depth records is reset to blank. I also rewrote mlocio_mnf.f90 to handle the issue of input data files having different MNF formats more gracefully. In subroutine stafind.f90 I added the flag 'm' (i.e., for fcode) for readings from stations that could not be found in the station file. They cannot have their phase re-identified (how could they?) and will no longer be printed in the .phase_data file and other places (except for the list of missing stations in the 'bad' section and in the .stn file), where they only caused confusion because of the arbitrary epicentral distance and azimuth that had to be assigned to them. They are written out to the .datf file, with the epicentral distance and azimuth taken from the input file, so that info is not lost.

2013/11/20 (v9.9.5). Implementing more support for reading and writing MNF v1.3, including the new depth records. The code can read data files in MNF version 1.1, 1.2, or 1.3, but the output file ".datf" is only in v1.3. The file ".dat0" is in the original input file format. I did a substantial rewrite on code throughout the program that deals with setting the starting depth. This aspect was getting very messy, and benefited greatly from reorganization. I may apply the same principles to the latitude, longitude and origin time variables. I added the "mloc_author", the person doing the analysis (EBergman in my case), to the mloc.conf file so other users could easily apply their names to their results (I certainly don't want them to be attributed to me).

2013/11/5 (v9.9.4). I cleaned up and expanded the usage of the 'mloc.conf' configuration file. The previous reference to the 'table_path' has been replaced by specific entries containing the pathname for the tau-p directory, the ellipticity directory and the master station file. The master station file is referenced explicitly to allow easy usage of a user's own master station file, rather than my 'eab_stn.dat' file, because each user is likely to edit their master station file to some extent to solve problems with on-off dates, as well adding newly deployed stations. In that case the file should be renamed to avoid confusion with my original file (which I will continue to modify myself). These changes provide more flexibility in setting up the program's file structure in a new user's directory structure.

2013/10/31 (v9.9.3). I added a new command "bdps" to help deal with the problem that has recently come to light regarding the reporting of teleseismic depth phases. The ISC discovered that certain Chinese stations have been reporting bogus depth phase readings that are generated from the theoretical arrival times according to some standard earth model and the preliminary location (and depth) in the PDE. The full dimensions of the problem, as to how many stations (it's probably not limited to Chinese stations) have done this and over what time period, are not yet well known, but I wanted a mechanism to flag (at least) depth phases from stations that are suspect in this regard. I introduced an optional station file, which is specified by the "bdps" command, that carries the list of stations suspected of this activity. For convenience I copy the station entries from the master station list for this purpose, and I keep this file in the tables/stn directory. For the moment I only check station code, but it's useful to have the other station information handy as well. In the future I may implement a check of the date-on and date-off fields (which would be altered to reflect the date range in which bogus depth phases have been reported) to facilitate the processing. For now the usage of this information shows up in the ".depth_phases" file, as an asterisk next to suspect stations, and in the type 6 travel-time plots (_tt6), in which the symbols for suspect stations are plotted at a smaller size.

2013/10/26. I added the confidence ellipse of calibration locations to the plots where calibration locations were plotted (base plot and selected event plots). The 'x' and the ellipse of the calibration locations are in red.

2013/10/21. I added another TT plot, (type 8) for S-P data plotted as a function of epicentral distance. Curves for depths of 5, 10, 15 and 20 km are plotted. I also rewrote much of the code for the type 6 plot, for teleseismic relative depth phases. The plot is done over focal depth, with all times reduced to the theoretical P arrival. Theoretical times are plotted for 60¡ epicentral distance, because they don't change much with distance. Individual plots are made for each event, because they have different focal depths, and all plots are stored in a subdirectory named "_tt6".

2013/10/16. I added the command "mdou" to control creation of an output file (.map_dat) containing epicenter data in a format that can be easily imported into a GMT script for further plotting. I added a 7th TT plot type to the "pltt" command, which displays the entire teleseismic P branch, plus Pdiff and PcP, everything reduced to theoretical P time.

2013/10/15 (v9.9.2). I added the command 'spou' to create a new kind of output file containing S-P data, for my project with Vera Schulte-Pelkum. I also made S-P readings immune to the automated determination of duplicate readings in mlocio_mnf.f90 because S-P data rarely arrive in such a way that they are duplicates. Readings that appear to be duplicates, according to the residual, are nearly always separate samples. It is still possible to explicitly flag an S-P reading as a duplicate.

2013/10/14 (v9.9.1). I modified the format of the hypocenter line in the puke file, replacing the old "GTX" nomenclature for calibration level with the new GTCNU nomenclature. Also replaced the single estimate of depth uncertainty with the asymmetric "plus/minus" values.

2013/09/30 (v9.9.0). I moved the specification of the environmental variables mloc_path, table_path, gmt_path and dirsym into a separate configuration file named "mloc.conf" which must reside in the same directory with the mloc executable. These values were formerly hard-coded into mloc.f90, mostly in the block data section, but having them in a separate file makes it easier to hand out an executable to someone with them having to edit the source code and recompile.

2013/09/29 (v9.8.8). Further modification to support MNF v1.2 output, especially the GTCNU system of nomenclature for calibration and location accuracy. I also rewrote the "cali" command so that there are several flavors, controlled by the 4th character of the command, similar to the "dep" command. The different flavors relate to which hypocentral parameters are considered to be calibrated. This was needed in order to compose the correct calibration code (GTCNU).

2013/09/24 (v9.8.7). Began adding support for reading and writing the latest version of the MNF format (v1.2). I renamed the code module "mlocio_mnf" and renamed the read and write routines to remove references to MNF versions. Version checking is just handled in the routines. As far as read_mnf is concerned there are no changes between v1.1 and v1.2, so it only requires accepting a format line with v1.2 specified. The subroutine write_mnf is  still putting out MNF v1.1.

2013/06/20. Fixed a bug in redsta.f90 with the type 4 station file (CSB) was dividing seconds by 36000 instead of 3600 because I'd copied the code from the FFB station format (type 1) which gives secs*10 in i3 format.

2013/06/17. I rewrote the section in mlocio_mnf_11 that pulls starting locations from an HDF file, to improve robustness against consecutive events with OTs very close to each other. I took out the code related to archiving warning messages and flushing them at certain points. I rewrote stafind.f90 to eliminate the possibility of over-writing the bounds of the arrays dealing with missing stations.

2013/05/22 (v9.8.6). I wanted to compile the code with gfortran for use in the workshop but it failed with a segmentation fault that I traced to the hyposat_loc package. In the process of fixing the problem I converted the package to f90 and I turned all double precision variables back to real because it looked like there were some conversion problems between the two types that were causing the segmentation fault. This did fix the problem, so I now get identical results from the Absoft and gfortran versions.

2013/05/13 (v9.8.5). Fixed some bugs relating to relative depth phases.

2013/04/25 (v9.8.3). Support for MNF v1.1.

2013/04/22. I modified the hdf file format to carry an evid for each event and rewrote the algorithm that matches an event with a previous hdf file, so that it tries to make the match using evids first and only falls back on the date/time match if no evid has been supplied (evid = 0). This provides a much more reliable method of association.

2013/04/16 (v9.8.1). I changed the format of the PUKE file significantly to add additional information in a form that will serve people who want to get the entire data set for various purposes, but who do not need to deal with the MNF data files. It's not as specialized as the tomography file format but it could be processed for tomography pretty easily.

2013/04/11. I removed the alternate station code functionality (command 'astn') because the problem of unregistered, conflicting station codes can be handled more effectively now within the MNF format. The unregistered code can be carried in the ADSLC portion of the phase line, and the corrected station code is read from the 'functional' part of the phase line. The 6th character of the station code can be used to de-conflict stations codes if necessary. In the case where stations have been registered at ISC but the agency still uses the old, conflicting code for internal processing and their own bulletins, the switch in station code should be made when converting to MNF format. Therefore there is also no longer a need for a flag to prevent the use of alternate station codes in mloc. I also removed the 'rsrc' command and related code because the 8-character 'author' field in MNF provides adequate room to document the source of a reading. I removed the mlocio_xcorr package and the 'diff' command because I intend to support differential time data in the MNF format, not as a separate file. The variables needed for differential time processing (e.g., idiff) are still in place. I need to design a special phase line format for differential time data. The line will need to carry info on both events.

2013/04/09. Fixed a long-standing bug in mlocset, in which the travel-time plots would only be made if the command 'datf' had been issued.

2013/04/04 (v9.8.0). With this version, support for the FFB format is discontinued. I have designed a new format specifically for use in mloc, MLOC Native Format, or MNF format. Data files have the filename suffix 'mnf'. MNF is based on IMS 1.0 format. It has many advantages over FFB, in particular being far easier to read, write, and edit. It also has the capability to carry multiple hypocenters for a given event, while designating one as the preferred hypocenter, and it has the ability to carry magnitude estimates far more conveniently, also with flexibility on designating one as preferred. There is only a single format for all standard phase readings and differential phase readings like 'S-P' and 'pP-P'. A new line type will be defined to handle differential time data, i.e., time difference of the same station-phase from two events. MNF supports the full 'IP style' addressing (agency.deployment.station.location.channel). It also carries evid, orid and arrid values, as well as an 8-character field for "author" of each pick. There are two fields for station code, one for the original code of the reporting agency and one for a code to be read into mloc, thus allowing for de-confliction of station codes without losing the original value. Similarly, phase name is carried in two fields, one of which is kept unchanged for documentation purposes. There are flags defined that govern the usage of the reading (e.g., flagged as an outlier) and also on the station code and phase name fields that can be changed for processing. I have written conversion programs to convert to MNF from FFB and the flavor of IMS 1.0 used by NEIC. FFB files can no longer be read in directly. I also removed the 'sing' command (sequential single event processing) and all associated code structure. This style of processing has not proven to be useful and it complicates the program structure significantly.

2013/03/22 (v9.7.1). Until this version there was always a test made in "model_init" regarding focal depth as a free parameter. The number of depth phases and near phases (defined as epicentral distance less than 1.5 focal depths) that were counted during the reading of the input data file were compared with limits to see if there is enough data to allow depth to be a free parameter. I removed that test in this version because I am starting to work more often with free depth solutions and I want to control this decision more explicitly, with commands 'frec' and 'freh'. I cleaned up some of the code that deals with documenting how depth is controlled.

2013/02/24. I added two commands for enhanced plotting. "star" is used to designate that the current event's location should be highlighted by plotting a solid star in red. "stat" is used to plot solid blue triangles to indicate the locations of seismic stations in the various map plots.

2013/02/22 (v9.7.0). I changed the format of the HDF files to carry more precision for latitude and longitude because some of the recent clusters with dense local networks (Prague, Oklahoma for example) needed it. With only three decimal places (~100 m resolution) a high-resolution plot of the locations was pixilated. So I'm carrying five decimal places (~ 1 m resolution) now. I also increased the precision of depth and the semi-axis lengths of confidence ellipses from one to two decimal places (10 m resolution), and changed the precision of the confidence ellipse area from unitary to one decimal place (0.1 km**2) to avoid displaying confidence ellipse areas of "0" for really well-located events. This required a change in the subroutine "hdf_read" to be able to read the various flavors of HDF files correctly, and a change in mloc.inc to make the "hdf_line" character variable longer.

2013/02/21 (v9.6.5). I rewrote mlocout_gmt to provide additional control over plotting of the maps. The naming of the plot files is now abstracted to mlocset, where a filename suffix is part of the call to mlocout_gmt. In addition there are logical variables in mlocout_gmt that control whether event numbers, confidence ellipses, and the relocation vectors are plotted. There is finer control over the relocation vectors through the "vect" command. If confidence ellipses are not plotted, the locations are indicated by open circles with diameter 1 km. This is useful for a cluster like Trinidad, Colorado where there is a very dense cluster with many events but many of the confidence ellipses are fairly large, obscuring the pattern of the epicenters. I added commands "eplt" and "splt" to control creation of confidence ellipse plots (no numbers or relocation vectors) and seismicity map plots (same symbol size for all events, no numbers) with standard elements. Further customization is easy through hand-editing the .bash files, or adding a new call to mlocout_gmt in the mlocset.



******************************************************************************************
Below this line, the version history is given in increasing date order:
******************************************************************************************
Written January 26, 1989 by Eric Bergman. Modified to use tau-p software, 9/20/94, by eab. Huge number of undocumented modifications during 1998-2004. Modified to use geocentric coordinates internally, 9/29/04 by eab.

Version 4.1 (2/3/05)

has the following major changes: Built with Absoft Pro Fortran v9.0. No longer using IMSL libraries. Turned off phase association, which was not really functional anyway. New format for covariance matrix output. Significant change to calculation of confidence ellipses for cluster vectors see mlocinv for documentation.

Version 4.11 (2/11/05)

added subroutine 'croux' for robust estimate of reading errors.

Version 4.2 (2/21/05)

Major reorganization to rd96b, some functions broken out as subroutines or functions attached to the rd96b module. I also broke out readerr as a separate module, broke out stafind into a subroutine in mloclib.

I added a "verb" (verbose) command to control extra output, especially to the application window. This is needed in single-event mode to avoid huge amounts of non-critical output.

The biggest change is that I added phase association, which now runs before each iteration except the third. The tools for handling phase association with custom determinations for median offset, spread, and relative observability for each phase are included but not yet being used. It cant be used unless you have a consistent description for all phases of interest. For now I use default values (median offset = 0., spread = 3.0 s, relative observability = 1.0). I am using the probabilistic approach to deciding which association to make when two or more phases are close together, with a Buland-style combo PDF (gaussians plus cauchy).

Basic rules for phase association: 1) With a couple of exceptions, initial phase records can only be associated as first-arriving phase. 2) Secondary phases can only be associated as secondary phases. 3) Depths phases can only be associated as depth phases. 4) Depth phase type cannot change (e.g., pP cannot become sP) 5) Non-depth phases cannot be associated as depth phases. Further work is needed on association of depth phases, a la Engadhl.

Version 4.3

incorporates the new shifting algorithm when there are reference events. This is done in subroutine refshift. There are two modes, the original average shift with no weighting and no consideration of GT uncertainties (mode 0), and one that includes GT covariance in the stimate of shift vector uncertainty and also uses weighted means to estimate the best shift vector (mode 1). The choice of mode is made the first time a reference command is invoked.

Version 4.4

Includes a depth-determination strategy in single event mode. A range (9 max) of trial depths is defined in the routine "model_init" based on the depth given in the input file. A full free-depth solution is found at each trial depth. A starting origin time for each trial depth is estimated, and phase-association is done before each new start depth. The median of the converged depths is used for a final run. I changed the rules about phase association for depth phases. p-type depth phases can become s-types and vice versa. It is also possible for depth phases to be changed to non-depth phases and vice-versa.

Version 4.5

Includes the new output file ".saic" that is used by Istvan Bondar for importing results into the database.

It also has the gmtout subroutine that creates a GMT script for plotting.

Version 4.6

Removes most of the v4.4 depth-determination strategy, which didn't work. I now set starting depths in a separate program (clean).

Version 4.7 (6/2/2005)

changes the nature of the .dat file output, to reflect the actual data and flags used in the run. Some readings (e.g., "LR") are not carried through, and the flags in input files may or may not be used. Now the program always creates the original-style .dat file of event data as it is read in, but this is now the ".dat0" file. This can be useful as an audit trail for the progression of flagging during a series of runs and for the phase identification changes that may take place during a run. Whether or not a .dat file is created is still an option, using the command "datf".

In this version I also added some extra info at the end of each reading line in the .phase_data file (mlocout2), giving the original phase ID in the input file and the index and number of readings in the station block. These are used by the program rstat3 when analyzing cluster residuals and trying to decide what to edit in the input files. I also reversed the order of iphas and nphase output to the LRES and XDAT files, so that it reads more intuitively ("ith of j readings"), and changed the programs lres and xdat accordingly.

I changed the meaning of the "dept p" command, which used to mean "use the pP-P depth info in the ISC Bulletin". Now it is used for single event mode, meaning to check the variable psedep in the hypocenter continuation line to see if the solution had a fixed depth solution (psedep = 99). This determines whether depth will be a free parameter for this event.

I also added a check for open azimuth, which in single event mode is used to determine if the epicenter should be fixed (openaz > 180 degrees) or not.

I added a hook to be able to turn off phase association for individual phases. The flag is a "!" in the last (8th) character place of the phase name.

Version 4.8 (6/8/05)

Changes the phase association algorithm to do reading groups all at once, instead of phase-by-phase. This allows for some useful rules, such as prohibiting association as a depth phase unless the parent phase has been associated.

Version 4.9 (7/25/05)

Changes the use of the "PRES" command and the two parameters (prmax1 and prmax2). They are still multipliers for ttsprd for each phase reading, but now they describe the transition from full weight to zero weight, as opposed to sharp cutoff (that changed with iterations) before. This should help stabilize the inversion by making it less sensitive to data falling in and out of the problem on subsequent iterations, thereby make flagging less critical. There's a new variable "weight" to carry this value, which is based on the path-corrected absolute residual "dts". If dts is less than prmax1*ttsprd, weight = 1.0. If dts is greater than prmax2*ttsprd, weight = 0. In between is a 1.-cosine reduction. Weight is recalculated in  mlocset before the inversion at each iteration. This weight is applied to vnhat in mlocinv for both cluster and hypocentroid. Obviously, prmax2 must be equal to or larger than prmax1. Same for all iterations. Default values are 3.0 and 4.0. I also changed mlocout2 (phase_data) to display weight in place of the station correction.

In this version I also made a small fix in subroutine croux for the case n=3. In the original version, it behaved badly in the case of two (of three) values being very close together, leading to unrealistically small estimate of reading error. It would ignore the "outlier" completely. After discussion with Christophe Croux I modified the algorithm and the finite-sample correction factor for n=3 to give a more sensible result. See the subroutine for details. It should now be safe to use Sn as the standard estimate of scale for reading errors, no need for the "largest of three" option.

I wrote subroutine "ttsprdout" to give some output on the observed spread of readings for different phases, but it's not so obvious what to do with it, since it combines reading error and systematic variations between ak135 and the real Earth. I need to do a more careful analysis, including effect of changing epicentral distance.

Version 5.0 (8/15/05)

Most of the improved statistics for GT shifts. See refshift2 for extended discussion. At this point the main remaining task is to implement a better method to estimate the "radius of doubt". Here I use the RMS of the reduced shift vectors (average subtracted), but that overestimates, too sensitive to outliers.

Version 5.1 (11/15/05

Implemented a statistical test to get radius of doubt. Based on the test of all cluster vectors having zero length in J&S. Null hypothesis here is that all reduced (or residual) GT shift vectors (subtract average) have zero length. Code adapted from the old mloctst.f program. Since all covariance matrices have already been scaled to reflect varying nf, mf, squared error, etc, the critical value is just 1.0 and the confidence level is 90%. If the null hypothesis cannot be rejected, it means that the individual estimates of the shift vector are mutually consistent at this level of confidence and no extra term is needed for the estimate of uncertainty of the optimal GT shift vector. However, if the null hypothesis CAN be rejected, the length of each reduced GT shift vector is reduced by an increment (0.2 km at the moment) and it's tested again. Vector lengths cannot go negative. Once the null hypothesis cannot be rejected, the accumulated increments that have been shaved off the reduced shift vectors give the radius of doubt (equivalent to asking what is the radius of a circle within which all the reduced cluster vectors could be found, taking into account their uncertainty of course). This is converted to CV components, and added to the CV of the optimal shift vector (agcv). The area of the associated ellipse, converted to a radius, gives the GT level (90%) of the cluster, the basic level of uncertainty of the calibration of the cluster.

Version 5.2 (12/3/05)

Fixed bug in refshift2 for the case of only one reference event. In that case, set rdbt=0. and skip the call to subroutine rdbt_test.

Along the same lines, bypass the call to subroutine croux in the case of one reference event, to get the spread of depth and OT shifts, and set those values to 0.

I fixed a bug at the end of refshift2, where it did the wrong thing to get the final uncertainty of depth and OT for reference events. It was adding depgt**2 and otgt**2 instead of depsn**2 and otsn**2. This had the effect of adding the uncertainty of the reference location in twice.

Version 5.21 (12/12/05)

Changed phase_assoc_grp so that Engdahl's pwP phase identifications would not be changed.

Version 5.22 (12/13/05)

Changed phase_assoc_grp so that Engdahl's PKPpre phase identifications would not be changed.

I also reduced (from 10 s to 4 s) the amount that a primary phase needs to be late in order for it to be identified as a secondary phase. This came up with PKPbc phases being stronger and more easily observed than PKPdf in the range where PKPdf is the first-arriving phase. Need to do some more work on this logic.

Version 5.3 (1/10/06)

Adds the capability to specify a custom crustal model for travel times of local phases. This is invoked through the TABL command, which has been rewritten. The algorithm for this is taken directly from Johannes Schweitzer's hyposat program.

version 5.4 (3/3/06)

For local distance readings, ttsig now has a minimum value, which I've set to be equal to the ttsig for that phase at 1 degree. This is needed so that readings are not dropped out on the first iteration because the starting location is a little too far off.

Don't reassociate depth phases if the event depth is fixed.

Changed PRES command to WIND (window) to reflect the use of these parameters to establish the window of acceptable residuals for all phases.

I converted the old PKPL and RAYP commands into parallel commands (HLIM and CLIM) to set epicentral distance ranges for cluster vectors and hypocentroid. In each case, up to three distance ranges for use can be set.

It's only an internal change but I got rid of the old ifltrh and ifltrc variables, and replaced them with a set of logical variables to keep track of which readings would be used for cluster vectors and/or hypocentroid, and why. This should be cleaner and make it easier to handle the logic when the datasets for cluster vector and hypocentroid are different.

I changed HYPL to PHYP, and all it does is to set a logical variable to use P phases only, for the hypocentroid.

New output file (.hypo_phase_data) to hold local-distance readings when they are used for the hypocentroid.

version 5.41 (3/27/06)

If a phase falls off the end of the branch in the TT tables after we're done with phase association, it is renamed as "PUNK" or "SUNK", depending on the type of the prior phase name, and regardless of whether phase re-association is permitted in this run. We don't want to lose the information on type of phase. In the .saic output file, the residual for such readings is set to -999.0, as per Istvan's request. In the .phase_data file, however, I still print the full travel time of such readings in the column for residuals.

version 5.5 (3/28/06)

The empirical travel time spread info is now output into a separate file ".ttsprd". These spreads are now calculated only from readings that pass the cluster or hypocentroid filter. There is an option to use the default values as before (ttsig) or to read values of ttsprd from one of these output files (ttsig2). In addition to the spread of each phase, the ttsprd file contains an estimate of the offset of each phase set from the reference TT tables. In applying the window (WIND command) to set weights, each reading will have the offset for that phase subtracted first. This helps with the situation where some phases have baseline offsets relative to the phase(s) used to estimate the hypocentroid, and it is statistically valid because the window width is based on an estimate of the spread of residuals (Croux's Sn) that is independent of central tendency. The offsets are not used to build the residual data for inversion of cluster vectors or hypocentroid, only in the application of the weight windows. Default estimates of spread have zero offset (for now).

version 5.6 (4/5/06)

No significant changes to function but major reorganization of program modules. Combined all output subroutines in a single module. Eliminated libsun.f, some routines were not needed and the rest were combined into libtau.f. Started cleaning up the code in libtau.f. Declared most logical unit numbers in the block_init portion of mloc.f, rather than allocating them at various places in the program using lunit.

version 5.7 (4/14/06)

Major changes to some commands: dept, lat, long, and time. In each case the command always requires at least one parameter, an integer 0 (fixed) or 1 (free). If the command is given before any events are declared, it applies to all events unless over-ridden by the same command issued during the event declaration phase of the command file. Also, the meaning of the second parameter (if given) depends on whether it is a "pre-event" call, or specific to an event. Before any event is declared, it is treated as a perturbation to the starting value, either read from the input file or specified directly in the event declaration phase. In the event declaration phase, the second parameter is taken as the starting value, which may be modified further if a universal perturbation was declared earlier. For the time command there are three parameters (hour-minute-sec) that would follow the mandatory first parameter. There are limits to the possible perturbation of each parameter (see the code).

I added a new command sstn to make use of a supplemental station file easier. The command is followed by the pathname of the file, but the pathname must start with "t" under the assumption that the supplementary station file is stored in the "tables" directory. I created a "stn" subdirectory in tables to hold such files.

Similary, the tabl command now has the option to give a pathname to a crustal model for local-distance travel times in a command file. The file must be in the tables directory, and I created a subdirectory "crust" to hold such files. It's still possible to select crustal models interactively with the "l" parameter.

All these changes are motivated by situations that arise when using local-distance data to estimate the hypocentroid directly. It is very useful to be able to perturb the starting parameters of all events by an amount similar to the mislocation vector determined in a previous run. Otherwise, there can be problems with the starting model (based on teleseismic data) being too far away and losing most or all of the local distance data to high residuals. Also, if ttsprd data are to be used from an earlier run to correct regional and teleseismic data, those offsets will be relative to the final (unbiased) location. So the ability to start the inversion close to the final solution helps a lot with convergence.

I changed ttsig and ttsig2 so that there is less increase in ttsig at iteration 0: 1.2 instead of 2.

It makes no sense to calculate reference shift if the relocation does not converge, so I added appropriate logic to mlocset.

I set gmtout so that it now plots the augmented confidence ellipses (including GT uncertainty) if reference data has been used.

I modified rdttsig, ttsig, and ttsig2 so that a minimum TT spread of 0.7 s would be used. If crustal readings (Pg, Sg, etc) at fairly short distances are being used for cluster vectors, but the hypocentroid is estimated as usual from teleseismic P, there can be rather large station-specific offsets for the crustal phases from the ak135 model that causes problems keeping the data in the problem. You cannot assume a small TT spead for stations at near distances unless these data are actually being used for the hypocentroid.

version 5.8 (5/4/06)

I redid some of the commands related to event parameters. The DEPT, LAT, LONG, and TIME commands were too complex. I added two new commands to do some of the tasks. FREE is used to declare free and fixed parameters for all events at once (if issued before any events are defined) or event-by-event. PERT is used to define a perturbation that is applied to one or more parameters for all events. Now, the only function of DEPT, LAT, LONG, and TIME is to set a starting value for the parameter. If issued before any events are defined, applies to all events. Otherwise, applies only to current event.

version 5.81 (5/9/06)

Added one more iteration, and made it easier to change in the future by using a parameter itmax (now = 4) ste in mloc.inc, instead of hard-wiring a specific value throughout the code. Ever since I started doing phase re-identification, and especially since the new WIND command with variable weighting on the shoulders, convergence has been a problem.

Made some minor fomatting changes to output files, ~.phase_data and ~.summary. Changed the logic that controls what is written into ~.xdat, so that only readings with large residuals that are currently unflagged and delta greater than 2.5 degrees are written.

Changed weighting scheme so that readings at local distance range (delta < 2.5 degrees) are always given full weight. It is impossible to control the weighting of these readings with the usual (WIND + TTSPRD) logic, even when an offset has been read from a ttsprd file. It has to be done manually, because the spread of theoretical travel times at local distances is outweighed by other factors, such as the bias from locating the hypocentroid with teleseismic P readings.

Added current (re-identified) phase name to end of line in ~.lres to make it easier to search.

version 5.82 (6/14/06)

Minor changes that I forget.

Version 5.83 (6/23/06)

Fixed a few minor bugs in the way phase names are cleaned up and the setting of the 'p' flag for skipped phases.

Added a new command 'plot', which is used to determine which events are plotted in the GMT script. In any case, a GMT script plotting all events is always created. If the 'plot' command is issued after the MEMB command for an event, then it is selected for a second GMT script of selected events.

Version 5.9 (6/29/06)

Many changes is rd96b related to processing input files. I'm making more consistent use of the "p" flag for records discounted because of the phase, and the "d" flag for duplicate readings. To deal with duplicate reading groups, I am keeping track of a reading group index so that LRES can find the correct readings to flag. I introduced the "d" flag for duplicate readings within an event, which have the same station, phase, and arrival time, and are almost always in different reading groups. Only one of such readings is kept in the problem. If station and phase are identical, but the times are different, I keep both readings, as independent samples.

I'm now catching PFAKE readings in type 5 records with the numerical phase ID 111. These are kept in the problem because its easier to keep track of reading group indices if I don't skip any, but they all get the "p" flag and are not used.

For secondary phases, I'm catching many of the amplitude readings by numerical phase ID = 117. Other kinds of junk records are caught with the "zero-day" test, where they have a blank day field. These are skipped.

Version 5.92 (7/6/06)
I changed the .rderr file, to carry only Sn as a measure of spread, and also the average, as measure of central tendency (for empirical path anomalies). 

Version 5.93 (7/8/2006)

Moved output of the details of calculation of GT shift in refshift2 from the window display into a file ".cal". The ".summary_ref" file is mostly superfluous now.

Cleaned up an unfortunate choice of terminology from way back, from "phase association" to "phase identification". Phase association is the process of deciding what event a reading is "associated" with. I renamed the subroutine phase_assoc_grp to phreid, renamed all the "assoc"-based variables, and edited comments to be clearer. There is no change to functioning of the program EXCEPT that the "asoc" command is now "phid". I modified the command parser to accept the old form of the command, with a warning.

Version 5.94 (7/10/2006)

Rewrote the command processing structure, cleaning up a few loose ends at the same time. New 'mlochelp' subroutine for the help command. No substantive changes in how things work, but it will be easier to add or subtract or change commands in the future. All the gotos and line numbers are gone. Now it is just a long set of if-then-else options.

Version 5.95 (8/4/2006)

Added a new command (FREH) to control free parameters for the hypocentroid. The old "FREE" command is renamed "FREC", and it only controls free parameters for cluster vectors. This was done because there are many cases in which there is insufficient connectivity for depth phases to have depth as a free parameter for cluster vectors, but there are still enough depth phases to allow a free-depth soltuion for the hypocentroid.

When a phase is renamed as PUNK or SUNK because we can't find a theoretical travel time for it, we now set fcode = 'p'. If we're using all phases (not PHYP) for the hypocentroid (or if there are several of these for the same station for cluster vector), these will blow up the inversion if they are not flagged, because their derivatives are 0.

The PHYP command now will restrict the phases used for the hypocentroid to P, Pg, Pb, and Pn. It formerly accepted only P, but I needed the equivalent functionality for the case when local or regional network data are used, and the epicentral distance range (hlim) can be used to control which P-type phases are actually used.

Reorganized the treatment of readings for which no phase name is supplied, or for which the supplied phase name is nonsensical. Now using "UNKNOWN ", and "UNKNOWNP" and "UNKNOWNS" in cases (e.g., PPP) where the phase is known but is not part of the ak135 phase set. Istvan did not like blank phase names in his data sets for RCA.

Version 5.96 (8/31/2006)

I changed the output of the .hdf file, such that the column that used to carry the sample variance (error/ndata) now holds the shatsqci value, normalized sample variance, using the full Bayesian formula (see mlocinv). So these values represent the spread of individual events around the overall shatsqc that is printed in the .summary file. If I've chosen K correctly, the spread will be close to the standard deviation predicted from Fig 3 of Jordan and Sverdrup (1981). For K = 3, the SD should be around 0.35. I am estimating the spread (with Sn) and printing it out just below the final value of shatsqc in the .summary file. Large deviations on the high side may indicate a need for additional cleaning.

I added a new command (pltt) to control creation of three kinds of travel-time plots (GMT scripts): a summary chart and special ones for the PKP caustic and for local-regional data. The local-regional one can be done with reduced travel times if desired.

I changed the behavior of the code that makes the GMT script for the standard plot. When there are calibration events, their positions are plotted with black crosses (no ellipses yet). The event ellipses are just the relative location ellipses, but the events are plotted where the calibrated cluster would put them. Thus, an event with calibration data is not plotted exactly at the calibration point, nor is it given the calibration location confidence ellipse, as it was before.

Version 6.0 (9/25/2006)

The main change with v6.0 is the addition of capability to use differential time data, specifically from Ben Kohl's waveform cross-correlation studies. I introduce two dummy arrival times for each differential arrival time datum. The first one is based on the theoretical arrival time for the first station-phase-event. Each station-phase pair is given an index based on its record number in the special data file that is read by the new "diff" command. Regular bulletin data all have an index = 0, so association goes on as usual. This index is used internally to keep the differential time pairs unique. The datum is number of seconds between the two waveforms, assuming they occurred on the same day. This is added to the first dummy arrival time to generate the second dummy arrival time. I also read Ben's estimate of uncertainty, which ranges from about 0.03 to 0.09 s. These data are so accurate that I had to introduce a second error term for the cluster vectors to account for errors in the assumed velocity model for the differential theoretical times. Value of 0.05 to 0.10 s seem to be about right. A new command "cvtt" (crustal velocity travel time) controls this. Minor changes to the phase_data file to add the differential time index to all readings. Differential time data in the phase_data file can be identified two ways, by the non-zero index number in columns 142-145, or by the presence of an asterisk in column 1, in front of the station code. To make parsing the phase_data file easier I changed the write format so that all lines are padded with zeros out to column 145.

The subroutine that reads in differential time data is called "rd_diff.f". This is specific to the file format from Ben Kohl. 

I fixed a bug concerning the supplemental station code file, so that the program chooses the station from the supplemental file if there is a conflict with the main station data file.

Version 6.1 (10/4/2006)

I added the "rsrc" command to keep track of reading sources.

Version 6.2 (10/12/2006)

Identical to v6.1, only it is compiled on an Intel Mac using Absoft patch with Pro Fortran v9.0. Executable runs under Rosetta.

I set windloclim = 0. in mloc.f because experience showed that forcing full weight for short-range readings caused more trouble than it saved. If there were any large outliers at short range it would usually blow up the inversion. It was particularly troublesome when using a custom crustal model, as phases could disappear suddenly and leave a previously good reading looking like a huge outlier.

I made a number of small changes to improve the handling of PKP phases and pwP phases. pwP phases now get the "p" flag, until I implement the code to calculate travel times for them. The logic of phreid and related code is changed so that if a phase is identified as PKPdf and it's residual is less than 3 seconds, it will be treated like a first-arriving phase, even if Pdiff is theoretically the first-arriving phase at that distance. Too many PKPdf readings were being changed to PKiKP.

I added the idiff index for differential time data to the .saic output format.

Version 7.0 (11/15/2006)

Major code conversion to compile with Fortran 90. I had trouble with libtau.f90 - it would compile, but I had a problem reading the ak135.hed file. For now I am still using the libtau.f version, compiling with f77.

Version 7.01 (11/17/2006)

Replaced "mloc.inc" with an equivalent module "mloc_declarations", but we don't use named common blocks any longer. The tau-p codes are still using include files for now. 

Version 7.03 (11/20/2006) Minor changes to the main program and command processing. Some comands can now be toggled on and off (e.g., verb, phid, phyp).

Version 7.1 (11/30/2006) New command "rhdf" to read starting locations from an HDF file.

Version 7.11 (12/5/2006) Added commands "lgou" and "lgtt" to better control handling of Lg phases. I moved the output section from mlocout3 to mlocout2, where it's part of the .phase_data loops, instead of the PUKE file output.

Version 7.12 (12/28/06) Trying to fix the problem with local velocity structure that was introduced when I started using declaration statements. I went back to using ttloc.inc and the F77 version of hyposat_loc.f.

Version 7.13 (1/8/07) Fixed a bug introduced in v7.1 in the way an HDF file was used to set the starting location. The problem was that it was being done after reading input files, so phase identification could be wrong if the HDF file changed the starting depth significantly.

Version 7.2 (1/19/2007) I finally had to undo much of the coding upgrade done since v7.0 to get hyposat_loc to work properly. I'm not sure if it was the f90 coding or the use of the v9.0 Absoft compiler on the Intel processor that was the problem.

Version 7.21 (1/23/2007) Added the ability (in rd_loc_mod) to have comment lines in the files that define local crustal models. Comments are needed to keep track of the source of the different models.

Version 7.22 (1/25/2007) Added command "fmap" to plot fault maps (or other digitized lines) in GMT. Eliminated the .summary_cal file by adding the unique parts of it into the standard .summary file.

Version 7.25 (2/1/2007) Added the ability to define multiple selected-event plots by specifiying an argument for the "plot" command. Maximum 9 selected-event plots. Events can belong to more than one set of selected events.

Made the "sstn" command more flexible by allowing for an index in the first line that sets the format of the file. It supports the original format (same as stn.dat) and the ISC FFB format for station data, in geographic coordinates.

Version 7.3 (2/8/07) Deleted the "data" command. Data file formats are now taken from the filename of the data files. The suffix ".ffb" (ISC FFB format) is the only one fully implemented; the code is backward compatible to the old-style naming convention without a suffix. The suffix ".isf" (ISC ISF format) will be implemented in a later version. This will allow clusters to be composed of events in different data formats. By adding support for formats used by NEIC and other rapid-response agencies, it will be easier to quickly assemble clusters for important recent events.

Version 7.31 (2/26/07) I removed the bit of code in rd96b in which depth phases would not be re-identified if depth is fixed. I modified "model_init" so that starting parameters set by the "lat", "long", "dept", or "time" commands take precedence over an HDF file.

Version 7.4 (4/4/07) I fixed a bug in mloc.f90 which caused the last command in a command file to be processed twice. It doesn't matter for most commands, but when "cali" was the last command it would increment the number of calibration events again. I also improved the coding in cal_shift.f90 to handle larger numbers of calibration events. the limit was 50 in several places, but the Pahute cluster had 52. Now the limit on number of calibration events is the same as the number of events in the cluster.

Starting with the 7.4 release, all version 7 releases are being built on the Mac Mini because v9.0 of the Absoft Pro Fortran compiler requires an older version of Xcode (2.2 or less) than is required for the v10 compiler on the MacBook Pro. There are still problems getting mloc to compile under v10. When I get things working with v10 of the compiler on the MacBook Pro, all version 8 releases will be built there.

Version 8.0 (4/10/07) First version built with Absoft Pro Fortran v10. There were problems getting the tau-p software to work. I finally used the libtau.f95 code, which required some additional modifications. I moved all the named common blocks into '.inc' files and added several 'block data' subprograms to initialize variable in those common blocks.

I was not able to get the program running reliably when built under Absoft Pro Fortran v10 on my MacBook Pro. There were crashes in the main program and also the MRWE environment. Until further notice, I moved development of mloc onto the Mac Mini (G4) running Absoft Pro Fortran v9.0 under and older version of XTools. The versions built on the Mac Mini run fine on my MacBook Pro.

Version 8.1 (4/14/07) Added the new command "cvff" (cluster vector fudge factor) to apply an additional uncertainty of cluster vectors that I thought was revealed by tests with Pahute and Rainier Mesa nukes. When GT0 locations are used I found that the radius of doubt for the calibration is about 1.0 km, even though these arrival time data sets are extremely clean and high quality. I thought the best explanation is that the arrival time data, even though they have been cleaned to mimic a gaussian distribution, still have non-gaussian aspects that cause additional bias in relative locations. Unmodeled crustal velocity effects may contribute as well. Whatever the causes, I concluded that  1.0 km is probably a minimum value for this unmodeled uncertainty in any cluster. The cvff command specifies a radius (km) of additional uncertainty that is added to every cluster vector's covariance matrix to help account for this non-gaussian bias. Most of the work is done in mlocinv. However further testing and debugging revealed that I was calculating the radius of doubt incorrectly. the problem showed up with the Pahute cluster when I had 52 calibraion events, and exceptionally large number. I was assuming kcrit = 1 but it needed to be dependent on the number of calibration events. With the corrected algorithm, the radius of doubt for Pahute is zero. So I don't think there is a need to add additional uncertainty routinely but I left the cvff command in place for special situations.

In this version I also fixed a bug in "shorten" that caused the radius of doubt loop to blow up when a residual cluster vector that had already been zeroed out is trimmed again.

Version 8.2 (5/6/07) Fixed a bug in mlocout_ttsprd, in which differential time data would be used to estimate ttsprd. This is incorrect because differential time data is put into mloc with dummy arrival times based on the theoretical Earth model. I noticed this with the niijima1 cluster because it had so much Pn differential time data that it was screwing up the estimate of ttsprd for Pn, making it smaller than it should be because all the xcorr data have arrival times close to what is predicted by ak135. 

Version 8.3 (8/22/07) I finally got the program working with the new Absoft compiler (v10). The big change is the addition of the "dcal" command for direct calibration, and associated files. The output file of local-distance data is changed from "hypo_phase_data" to "dcal_phase_data". The main effect of the dcal command is to change the .hdf file output. If direct calibration is specified by the dcal command, the uncertainty of the hypocentroid is added to that of the cluster vectors to get the confidence ellipse printed in the .hdf file - in this case the name is changed to .hdf_dcal. You can still use local data (or whatever) to locate the hypocentroid and leave "dcal" off, in which case you get the traditional .hdf file with confidence ellipses based only on the cluster vector covariance matrices. There is another new file (.dcal) which summarizes the cumulative uncertainties. However, it doesn't really show much that is not in the .hdf_dcal file. I may add more of the details of the covariance matrix manipulation.

(8/28/07) I added station coordinates and epicentral distance from the hypocentroid to the .rderr file to aid in plotting empirical path anomalies. I also changed the code (mlocout_phase_data) so that the .rderr file is based on the calibrated (shifted) residuals in the case where calibration events are used. For direct calibration or no calibration, it's the same as before.

(8/29/07) Changed the format of .rderr files: the first line now carries the hypocentroid used for plotting empirical path anomalies. If the program tries to read an old-format file, it'll bomb out and use default reading errors. This can be fixed by adding a line that contains "Hypocentroid" in the first 12 characters. Actual values for the hypocentroid are not read by MLOC. I also changed the algorithm to assign 'p' flag to any depth phase with surface focus distance less than 28 degrees, following Bob's strategy. This is done only in mlocio_ffb, reading in the data and then re-identifying phases relative to the starting locations. It's still possible for a phase to be re-identified as a depth phase at regional distance after the first iteration, but this is a pretty rare occurence and can be dealt with manually. I fixed a bug with the .hdf_dcal file in which the uncertainties of depth and OT were not being increased by the uncertainty of the hypocentroid parameters. I also commented out printing of the .dcal file, as it displays nothing that cannot be found in the .hdf_dcal file.

(9/19/07) I added logic to put the 'p' flag on local and regional multiples (PgPg, PbPb, PnPn, and equivelent S phases). These are mostly junk. This is done only when the data files are read in, after phase identification relative to the starting location. It's possible for some of these phases to pop up after iteration 1 and subsequent phase reidentification.

Version 8.4 (9/20/07) Fixed a bug in the code that generated the .rderr file. In mlocset I wasn't making the automatic 1-minute error correction for the residuals used for this output file. It is still possible for a reading to get the automatic correction in the main program but not get it here, because the calibration shift could shift it outside the +/- 5 second window. I saw this in the first test I made. In that case, inspection of the .rderr file reveals the problem and manual application of the 'k' or 'l' flag fixes it.

(10/2/07) In mlocout_hdf I added code to calculate the GT level (for direct or indirect calibration) and put it in column 6 of .hdf_cal and .hdf_dcal files, where Bob expects to find it. GT level is calculated as nearest whole number to the larger semi-axis length. Only values up to GT9 are permitted, so calibrated events with larger uncertainty have a blank entry.

(10/8/07) Modifications to mlocout_puke to bring the calculation of GT level into agreement with recent changes to the hdf files. I also extended the format to allow inclusion of parameters that are relevant to both the hypocentroid and cluster vector for hypocenter lines (e.g., open azimuth) and phase reading lines (e.g., defining phase).

(10/15/07) I added code to mlocio_ffb to prevent phase re-identification for PKPbc, as a primary or secondary readings. This over-rides the logic in phreid that allows re-identification of PKPbc when it's the primary reading but not the first theoretical arrival. PKPbc is a such a strong phase that it's usually incorrect to re-identify it as one of the nearby (but weaker) PKP branches.

(10/26/07) I added command "ctyp" to control how uncertainty is calculated for calibration events in indirect calibration mode. Up until now, the calibration uncertainty was used, but in some cases it may be preferable to do it like the non-calibration events: add the calibration uncertainty to the cluster vector uncertainty. A third option takes the best of the other two, as measured by the length of the semi-major axis.

Version 8.5 (11/30/07) This version is functionally identical to v8.4, except that it runs in the Terminal, not with MRWE. I was experiencing sudden crashes under MRWE since upgrading the OS to Leopard (10.5) and Jeremy at Absoft suggested running in the Terminal. The crashes caused the MRWE window to disappear without warning so it seems likely that the problem was there. This only requires changes to two subroutines in mloclib, the ones that read the .ttsprd and .rderr files. It's a little more cumbersome to do the file input this way. As a result of giving up the nice MRWE interface for selecting files, I added two commands (TFIL and RFIL) to control selection of .ttsprd and .rderr files.

In early December I added the ability to use differential phase data (e.g., S-P, the only one implemented so far). Each differential phase type is handled separately (for theoretical time, derivatives, etc), and is flagged by the phase name, ("S-P").

Version 8.6 (12/21/07) is identical to 8.5 except it is being built under the Absoft 10.1
compiler.

Version 8.7 (2/1/2008) has modifications to cal_shift.f90 to calculate the coverage statistic for each event and the overall percentage of coverage when indirect calibration is used. I fixed a serious bug in the "shorten" subroutine that is used in the radius of doubt calculation. I also determined that the previous approach to calculating kcrit90 for radius of doubt was incorrect, leading to an over-estimate, and therefore underestimating rdbt. I'm still not confident that I understand the right way to determine the critical value for this test, but I think it has to be less than or equal to the number of calibration events. At this time I am setting kcrit90 to the number of calibration events.

(2/7/08) I completely redid the radius of doubt test. The new philosophy is based on a statistical test of whether the observed number of "uncovered" residual shift vectors is consistent with the assumed confidence level of the ellipses. The cumulative binomial probability distribution provides the probability of the observed number (or more) of "uncovered" residual shift vectors. The test is done on a covariance matrix that is the sum of the HDC relative locations, the calibration event uncertainties, and the uncertainty of the inferred calibration shift (sgcv = scv + gcv). If the test fails (probability of the observed number less than 90%) a small increment is added to these covariance matrices to inflate the ellipses, and the test is repeated. When the probability of the observed number of uncovered rsidual shift vectors exceeds 90%, the added component to the covariance matrices is converted to a radius of doubt, which is added to the the final uncertainty of all events.

(5/13/2008) I have switched development of mloc from the Absoft v10.1 compiler to gfortran v4.2.3. The first production version, mloc8.7g, has some slight differences in the code base from the last Absoft-based version (mloc8.7), but there should be no difference in functionality and it gives identical output on the test cluster (tangshan9.6). Future gfortran-based versions will drop the "g" suffix.

Version 8.8 (5/29/08). Current development is with the Intel compiler. In this version, I added the "topo" command to specify plotting of topography in the GMT .command file. It uses the GLOBE 1-km DEM, but I only have the G tile available at this time (0-90E,0-50N).

(6/11/2008) I added "astn" command to handle alternate station codes, for the case where we get data with "internal" codes that have been registered at the ISC with a different code. Like ITSN stations in Iran. The codes are fixed upon reading in the data files, but the data files themselves are unchanged.

(7/8/2008) I modified the "sstn" command to allow for multiple supplemental station files, up to a limit set by the variable "nsuppstnmax". I also found and fixed a bug when using relative phase data (e.g., S-P) in indirect calibration mode. It was miscalculating residuals relative to the calibrated location. The problem only showed up in the .rderr file.

(8/20/2008) I modified the .rderr format, adding the azimuth from hypocentroid to station for each station-phase. I needed this in order to do a 3-D plot of anomaly vs. location in cylindrical coordinates.

(10/3/2008) I modified mlocout_gmt to handle topo maps with tiles of the GLOBE dataset other than the G tile. So far, only the H tile is added, but the structure is there to add others very easily.

(1/10/2009) I split up the "verb" command into "vscr" and "vlog" to control extra output to the screen and log file, respectively. I added command "rels" to allow explicit values (generally smaller than the previous default values) to be assigned for reading errors of direct phases (Pg, Pb, Sg, Sb). This is mainly needed when doing a single event location of a calibration event and we can't get empirical reading errors from cluster analysis. It's important to try to get the reading errors right (on average) in order to have the statistics be consistent (as reflected by the SHAT value). At the same time I made some modifications to how the ttsprd and windowing is done for local distance phases (epicentral distance less than the parameter "windloclim", set to 1.8 degrees. I fiddled with the exact algorithm to calculate ttsprd but in the end did not change the values much. TT spreads in this distance range are a little smaller than they used to be, but it shouldn't have much effect. Within this distance range, the window is fixed at (4.0, 6.0 s) because the normal procedure of scaling the window from ttsprd fails with the very small ttsprd values that are appropriate for short distances.

(2/4/2009) I added the depth of the calibrated hypocentroid to the .rderr file's first line. This is needed when using the .rderr file as input for a program that will compute travel times (or residuals) through a model, and it needs the source depth.

I added the "step" command to better control the number of iterations the program will run. This is helpful when a run is blowing up after a few iterations. You can stop it early and see what's going on. "step 0" is more or less equivalent to the old "fwd" command to just run a forward model and get the residuals. I rearranged things a little bit in mlocset to avoid calling mlocinv when we just want the forward model.

Version 8.9 (2/28/2009) I changed the version number because I started using v11.0 of the Intel Fortran compiler with this release. I don't think this will effect the way the program runs.

Version 8.93 (5/18/2009) Added the "tomo" command to output files suited for input to tomography codes.

Version 9.0 (5/19/2009) Still unsatisfied with the radius of doubt calculation, I brought back the algorithm from v8.1 that implements a test of the null hypothesis that all residual cluster vectors have zero length. I renamed it "rdbt_test1" and for now, it runs alongside the test based on coverage statistics ("rdbt_test2), which is still the source of the estimate of rdbt. The main problem with rdbt_test2 is that in most cases we have only a couple calibration events so the 90% criterion on coverage is actually a 100% criteria, which inflates the estimate of rdbt. I want to do some further testing before I decide which one I want to use as standard. I need to think some more about the "shorten" algorithm and whether this is really correct, once a residual shift vector has been reduced to zero. I also rearranged the logging so that "verbose" logging is no longer the default and the things you get in standard and verbose mode are more balanced. Standard runs (non-verbose) provide a reasonable amount of info, but you can get a lot more by turning on verbose logging with VLOG.

(5/23/2009) Removed the SLIM command - it turned out to be not useful. Added the PPRI command to allow input of phase names that cannot be changed. This was mainly motivated by work in Iran, where Sn is rarely observed and I wanted the input Sg phase identifications to be kept, even beyond the theoretical distance where Sn should come in first. Fixed things so the STEP command works as desired: it sets the number of inversions (mlocinv) that are run, not the number of iterations.

6/20/2009) Added code to write the local velocity model (if one is used) to the log file after I ran into a problem where I couldn't reproduce an older run because I had modified the local velocity model and kept the same filename.

Version 9.1 (10/7/2009) This started out to be exactly the same as version 9.0 but it is compiled under OS X 10.6.1 (Snow Leopard) using the version 11.1 of the Intel fortran compiler. I had two loader errors:
ld: warning: for symbol _brkc_ tentative definition of size 152976 from ./build/mloc.o is is smaller than the real definition of size 152968 from ./build/libtau.o
ld: warning: for symbol _stat_ tentative definition of size 46368000 from ./build/phreid.o is being replaced by a real definition of size 368 from /usr/bin/ifort-11.1-base/lib/libifport.a(stat.o)
I fixed the second one by simply renaming the named common block "stat" to "station" in mloc.inc. I fixed the second one (without really understanding what I did) by cleaning up the code related to tau-p include files. libtau.inc was only called in the main program, and was no longer needed. It duplicated all of ttlim.inc and carried other unnecessary stuff. I got rid of the remaining reference to iottim in the main program, which was the only reason I ever had that include file there. I also ran the compiler with "-warn all" which caught a few minor problems and led me to go through and systematically use "implicit none" in all codes, and get rid of unused variables.

Version 9.1.1 (10/12/2009) I worked more on cal_shift and found a bug in my code for the radius of doubt based on coverage statistics. I was formulating the hypothesis test incorrectly and getting values of rdbt that were too large. It needs to be done on the lower tail of the cumulative distribution, finding the probability of the observed number of uncovered events (or more) and rejecting the null hypothesis only if the probability is less than 10%. Before it was being done at the 90% cumulative probability, not 10%. Testing with the Pahute cluster (pahute7), the two methods (coverage statistics and test for zero-length residual GT shift vectors) give similar results and I am more confident that I'm doing it right. I'm very confident about the test based on coverage statistics, still a little unsure about the critical value for the test based on zero-length residual GT shift vectors. But it can't be very far off.

Version 9.1.2 (11/3/2009) Added a new supplemental station file format ("MLOC" format), for decimal geographic coordinates and changed the logging so that all supplemental stations are written to the log file.

Version 9.1.3 (12/9/2009) There was an update to XCode and that required uninstalling the Intel compiler and reinstalling it. I got rid of the -warn all flag in the Makefile too - it creates a lot of extra files that aren't needed except for debugging. So there should be no functional difference between 9.1.2 and 9.1.3, but you never know.

Version 9.1.4 (3/18/2010) I modified the subroutine POSPHAS so that when the theoretical first-arriving S phase is Sn, it also considers the possibility of Sb and Sg phases. In the Iran project it is very common for Sg or Sb to be read instead of Sn as the first S phase. After I tested this it seems to create as many problems as it solves so I commented the changes out. I may take another look at this later. It may be better to just make a blanket rule that first-arriving S phases in a certain distance range are Sg, not Sn. So at this point, v9.1.4 is functionally identical to v9.1.3.

Version 9.1.5 (4/9/2010). I added a new graphic output file (actually several files if the "plot i" option is used to create '.sel' files) that is created when doing direct calibration. It shows all known stations in the region in direct calibration distance as "x"'s, and also shows the raypaths that are actually used in direct calibration. Stations for which only S-P data are available are shown as open circles. regular seismic stations as filled triangles. Circles are drawn at 1.0¡ and 2.0¡ from the hypocentroid.

Version 9.1.6 (4/13/2010) I added the "skip" command to allow flagging ("s") of all readings from a specified station. This is mainly intended for cases where it is suspected that a particular station has systematic timing errors and you want to quickly take it out of the problem without having to hand-edit all the input data files. I also found a bug in the code, related to mlocout_rderr. I had dimensioned the "qres" array to nevmax in several subroutines under the assumption that there would not be more than one reading for a specific station-phase per event, and I had even hard-coded a limit of 200 for "ndatamax" in one place. But now we are getting data sets in which there are many repeated measurements (not duplicates) of the same station-phase for one event. The zarand9 cluster broke that array bound. So I rewrote the code to put that array bound (n_qres_max) in mloc.inc and set it to 400. I also rewrote the "phyp" command so that it does not turn off S-P readings. I added the "phsk" command to allow skipping certain phases.

Version 9.1.7 (4/29/2010). I added another type of travel time plot (tt5), for the crustal S phases, Sn, and Lg. This one should normally be done as a reduced TT plot, in which case the reduction velocity is the typical Lg group velocity, 3.5 km/s (31.7 sec/deg). I revisited the issue with identification of crustal S phases that was first addressed in 9.1.4. I think I got it right this time, but there are still issues of identification as Sg vs Lg.

Version 9.1.8 (10/9/2010) I finally decided to try a different tack as far as phase ID. The old method was getting too complicated with all sorts of conflicting rules, and it was not doing a good job for crustal phases, especially around cross-over distances. I'm going to a more traditional approach: the list of candidate phase IDs for each arrival is simply the list of phases in the theoretical TT model for that distance and depth. The probability of each candidate phase is calculated using the same method as before and the highest probability wins. There are hooks to ensure that phases types (P or S) stay the same, and depth phases are allowed only to be identified as depth phases (so mainly, pP and sP can switch places). There is no restriction on multiple IDs of the same phase in the same reading group. This methodology is implemented in subroutine "phreid2". I also modified the function "prob" to use the values for spread and offset for a given phase from a .ttsprd file if one has been read in. There is still no distance-dependence and the value of "relative observability" is still always 1.0 for all phases. The results indicate that I should really do distance-dependent spreads for this. So I removed the use of spread values from a ttsprd file, but I kept the use of the offset. I changed the "phid" command to take an integer argument to select between the two types of phase identification.

10/15/2010 I modified the algorithm for calculating travel times when there is a local crustal model. Instead of a complete break at some distance, the local model is now used to calculate all crustal phases to whatever distance is specified in the .cr file. No crustal phases will be returned with theoretical travel times and derivatives beyond this distance. The global model is used to calculate non-crustal phases at all distances, which brings back some teleseismic phases that are observed at short epicentral distances. This does a pretty good job now, but there is a problem extending the hyposat software as far as I'd like to get the back-branch Pn out to around 20 degrees. At some distance in the 15-20 range the program bombs out because of the way it's coded. I have some tips from Johannes about where the problem is, so I may be able to fix it. For now I can go out to 17¡ safely, but 19¡ is too far.

10/23/2010 I removed the use of offset from a .ttsprd file in subroutine "prob". It turns out that even this causes problems (e.g., in direct calibration) if it is not treated as distance dependent.

2/8/2011 Version 9.2 changes the way the tomo3 data sets are calculated. Instead of referencing all empirical path anomalies to the cluster hypocentroid, I now keep track of which events contribute to a particular station-phase empirical path anomaly and calculate a hypocentroid for those events. Obviously, this is of most importance for geographically large clusters. 

3/23/2011 Version9.2.1: I made a single change in phreid2 so that a primary reading that matches the theoretical first-arriving phase will not have its name changed. The problem came up when doing direct calibration where the teleseismic P arrivals are consistently several seconds later than the theoretical arrival time, and they were being re-identified as PcP in the distance range where PcP comes in soon behind P.

3/28/2011 Version 9.2.2: Bob gave me a new stn.dat file which has a completely changed format, so the code has been altered to deal with that. Supplemental station data format types are unchanged. In addition to format changes, the coordinates are now given in geographic coordinates, which is more convenient, but it required some changes to the code that processes the main station file. The new format also includes dates (in Julian date format: YYYYDDD) for start and end of operation, so I added code to use this information (if it is present) in assigning station coordinates. There are stations which have moved substantially without changing station code.

3/31/2011 Version 9.3: This update was triggered when I noticed anomalies in the _all and _all_dcal plots, soon after updating from GMT 4.5.1 to 4.5.6. Ellipses and circles (which are plotted as a special case of an ellipse) were plotting at half the size they should. I learned that they changed the function of the -SE option for plotting ellipses at v4.5.2, so that it takes the full axis lengths instead of the semi-axes. So I modified all the related code in mlocout_gmt.f90. At the same time I fixed a long-standing issue where the S-P symbols were not plotted correctly in the tt3 plots (I just needed to implement the same code I'd written to do this in tt4 plots). I also finally implemented the "vect" command which controls which relocation vectors are plotted.

4/11/2011 Version 9.3.1: I modified the "map" command to allow multiple digital fault map files to be read. Current limit is 4.

6/20/2011 Added code to redsta (mloclib) to check for duplicate station codes in supplemental station files.

6/24/2011: v9.3.2. I made some additional output to the log file, to explicitly list each station actually found in the dataset, with coordinates. See subroutine "stations_used" in mloclib.f90. I added a lot more logging for various station-related issues.

7/21/2011: v9.3.3. Minor changes. I raised the limit on number of supplemental station files to 8. I also took out the portion of code in the "zero epicenter distance" loop that resets variables. It is just a warning now. It seemed to be causing a bug in certain obscure situations.

8/1/2011: v9.3.4. I created a new log file for station data.

9/17/2011: v9.3.5. I rewrote the mlocout_gmt subroutines so that the GMT scripts can be created for different shells. I implemented bash in addition to the original csh. I did a lot of reorganization of the code to make it more compact and easier to edit.

9/27/2011 v9.4.0. I started improving the handling of depth phase data. Adjusting focal depths to fit observed depth phases is dangerous, because in a calibrated location the teleseismic P data often has a baseline offset. I added a new plot (type 6) for displaying relative depth phase data (pP-P and sP-P), and changed the .depth_phase output file to show relative depth phase residuals when they can be calculated. It still shows the raw depth phase residuals otherwise.

12/20/2011 v9.4.1. I changed the format of .hdf, .hdf_dcal, and .hdf_cal files to provide additional information, remove useless legacy info, and to make sure data columns are separated by spaces, to aid in parsing. The format is still similar to Bob's original HDF format but has many differences which are documented in the code. In particular I now have a facility to annotate the source for setting starting depths - see the code for the DEPT command. I also added the 'anno' command to allow annotations for events. The text is limited to 20 characters, which is appended to the line for each event in the various HDF files. As of 1/8/2012 I added the depth-source information after the depth in the ".depth_phases" and ".dcal_phase_data" files. It helps in reviewing depth-related information in those files to know what the source of assumed depth was. 

1/10/2012 v9.4.2. I added a new version of the standard "_all" plot, in which event numbers and location shift vectors are not plotted. That is, only confidence ellipses are shown. It has the filename suffix "_all_ell". This is useful for seismotectonic interpretation when there are many events close together. It is hardwired in for now, but it would be easy to add a command to control whether this plot is made.

1/11/2012 I added some of the output data that had been left blank in the HDF files: Number of stations for hypocentroid, cluster vectors, and flagged, nearest and farthest stations for cluster vectors, and largest open azimuth for cluster vectors.

1/23/2012 I added near-source readings to the depth_phases file so that all the information used to constrain depth is in one place.

1/25/2012 v9.4.3 I fixed a small but important plotting bug. Under indirect calibration the TT plots would be made with the corrected travel time (based on the shifted OT and epicenter) but they were plotted at the unshifted epicentral distance. Now the symbols are plotted at the shifted epicentral distance. It can make a noticeable difference in some plots, especially TT2 and TT4.

2/8/2012 v9.4.4 I made a change to the algorithm for automatically detecting duplicate readings. It was being defeated too often by recent ISC data files because it required a match to station, phase name and exact time. There are many duplicates coming through the ISC with different phase names (e.g., same times but one called Pg and one called Pb), and also with arrival times that vary slightly because they've been rounded off differently (two decimal places vs one, or sometimes even to the nearest second) at some point of the processing. Multiple readings of the same station-phase are not a problem if they are truly independent readings, but duplicate readings screw up the statistics for empirical reading errors. I was having to spend too much time hunting these things down and fixing them manually. Now the algorithm is more draconian. Each reading is checked against all previous readings as before, but now the only criteria for a match are the station name and the arrival time reading, with a range of ±0.1 s. Only readings with a blank value for fcode are tested, so if there is any other flag set it will stay the same.

2/17/2012 v9.4.5 I added the subroutine "mlocout_kml.f90" to create a .kml file of the epicenters for display in Google Earth. It codes depth ranges by color and icon sizes by magnitude.

2/20/2012 v9.4.6 I added the FPLT command to better control the plotting of fault maps. Now you can specify the options to GMT's psxy -Sf command to indicate the type of faulting. These are very complex options so I just made it possible for the user to put in the entire -Sf command string for each fault file to meet their needs. I also made a few minor changes to some output files, at Reza's suggestion, to make it easier to read them and work with them. Nothing substantive. I added a header line that defines the data columns in the .depth_phases file and I added comments to the TT4 file to indicate what the different plotting segments were about.

3/8/2012 v9.4.7 I added the XSEC command to plot cross-sections (depth profiles). Multiple cross-sections can be defined. Cross-sections are defined by the lat-lon of the end-points, a depth limit (minimum is always 0), and a width (km) within which events will be gathered for the cross-section.

5/10/2012 v9.4.8 I enhanced the "dep " command to optionally take one or two more values for an assigned uncertainty in depth. If one value is given, it's taken as a symmetric depth uncertainty. If two are given the first is taken as the "plus" (deeper) uncertainty, the second is taken as the "minus" (shallower) uncertainty. These values are not used in the relocation, but are carried in the HDF files for future reference. As a result, the format of the HDF files is changed to carry two values for depth uncertainty.

6/7/2012. I modified "mlocout_summary.f90" to calculate and display (for each event) the information on how much the solution differed from the hypocenter provided from the input file (which is not necessarily the same as the hypocenter used for the start of the inversion). I also added a line at the beginning of the .depth_phases file, giving the median of constrained depths (i.e., not from input file or cluster default). This should be helpful for setting cluster default depth.

6/10/2012 v9.5.0. I added the command "mech" which allows information on focal mechanism or moment tensor to be associated with an event. The argument is just a text string, so the user has full control over what parameters are supplied and the format. This information is appended to the relocated hypocenter data in an output file with '.focal_mech" filename suffix. The intended use is to create a file that will be easy to read as input for plotting focal mechanisms in GMT, but it could be used for other purposes.

7/11/2012 v9.5.1. I added the "kill" command to simplify excluding events from the cluster. Instead of commenting out 3+ lines, you just need to insert the kill command before the 'membÕ command. This also makes it easier to keep track of events that have been excluded, but I have not yet written code to do that. The command only kills the next event.

7/31/2012 v9.5.2. I changed two default parameters to better handle the Louisa, Virginia cluster, which has an unusual amount of high quality data at very short epicentral distances and well-resolved events less than 10 s apart. I set the minimum reading error (rderrmin) to 0.09 s from 0.15 s and I set the time window for distinguishing events (adt_limit) in an hdf file to 7.0 s from 15 s.

8/6/2012 v9.5.3. I added a new command "mare" to control minimum allowed readings errors. It takes separate arguments for local distance phases and all others. I set the default values to 0.15 s for non-local phases, and 0.10 s for local phases (Pg, Pb, Sg, Sb).

8/27/2012. Because my new MacBook Pro runs Mountain Lion (10.8) I would have to buy a new license for the Intel compiler to continue development. I may still do that but I decided to try the Absoft compiler again. It was surprisingly easy to get the code to compile and run in the Absoft environment. I only had to fix a few lines that were too long before the continuation symbol- the Intel compiler is evidently more lenient about line lengths. I ran both versions on tabas8.37. There are many very small differences, which I think come out of the TT calculations, on the order of 1/100 s. Slight differences in the numerical routines, I suspect. The relocation results are very close, and I'm not sure what I could do to make them any closer. The only significant difference is in the output of the working directory for various files (e.g., TT model) in the .summary file. The Absoft compiler gives the relative pathname, whereas the Intel compiler gives the full pathname. I don't think this has any major repercussions. So it appears that I can continue development with Absoft. I think it would be good to get the Intel license and carry on development in parallel. Each compiler flags different things so running the code through both of them makes for a more robust code base.

8/27/2012 v9.5.4. Fixed bugs in mlocout_gmt.f90 and mloc.f90 that prevented the plotting of faults.

9/3/2012 v9.5.5. Major rewrite of mloc.f90 and the help and messaging systems. I put the processing of all commands into separate subroutines and moved them into a new source code file called 'mloc_commands.f90'. I also moved the help system for each command into the same subroutine so they are self-contained and easier to document. I also added subroutines 'warnings' and 'fyi' to handle the output of messages throughout the program in a consistent manner. No changes to the relocation algorithm.

9/6/2012. Added printout of events with large open azimuth in mlocout_hdf.

9/12/2012 v9.5.6. Added the command "clec" to control the convergence criteria for cluster vector epicenters.

9/27/2012 v9.6. Changed the way eigenvalues are truncated to try to solve problems with Maule clusters. Until now, the smallest M eigenvalues were zeroed out, where M is the number of free parameters for a single event's cluster vector (usually M = 3). But only the smallest eigenvalue is really small, usually 10-5 or so down from all other eigenvalues. After a few trials I also changed the way convergence is calculated. Instead of applying individual criteria to the changes in hypocentroid and each cluster vector I apply them to the sum of cluster vector + hypocentroid. This leads to faster convergence when the cluster vectors oscillating with the hypocentroid, but I think it is just masking a more fundamental problem.

9/30/2012 v9.6.1. Those experiments in v9.6 with zeroing eigenvalues and convergence criteria do not seem useful, so I took them out and went back to the old algorithms in this version. I also took out the "clef" command, as this is not a helpful way to deal with convergence problems.

10/4/2012 v9.6.2. I added the command "tikh" to implement Tikhonov Regularization for the estimation of perturbations to the cluster vectors.

10/7/2012 I added some extra logging for the kill command to carry the comments about why an event was killed.

10/31/2012 I fixed a bug in the code for plotting cross-sections (it1 not defined) that was not fatal under the Intel compiler but caused a crash under the Absoft compiler. The cross-section plots under Intel were not correct but in most cases the error was too small to notice.

11/12/2012 Added a new station file format, for Kevin Mackey's station list at MSU.

11/16/2012 v9.6.3. I added plotting for residual calibration shift vectors (in blue), controlled by the VECT command.

11/21/2012 v9.6.4. I added two new subdirectories to the tables directory, to organize things a bit better. One ("ellipticity") is for the file (tau.table) used for ellipticity corrections, and one ("tau-p") for files related to tau-p calculations. I also moved Bob's master station file into the "stn" subdirectory and renamed it "ere_stn.dat". The only change to the code is fixing a few hard-wired pathnames to reflect the new organization. Obviously, older versions of bloc will not run properly with the new organization of the tables directory unless copies of the affected files are left in the old locations with the correct names. This would be easy to do.

11/23/2012. I added some array-boundary checking loops in mlocout_ttsprd. I was hitting the limit with the wenchuan4 dataset and getting a segmentation fault.