PROGRAM Seismicity ! ! Reads a seismic catalog, such as: ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! (1)International Seismological Centre (ISC) mb body-wave magnitude ! (~3 to ~6.8) catalog, 1964-1991 (complete 1971-1991), on CD-ROM; ! ! (2)Harvard Centroid Moment Tensor (CMT) moment-magnitude Mw ! catalog for large (>(4.8~5.7)) events, Jan 1977-Dec 2004, downloaded ! from ftp://saf.havard.edu/CMT/allorder.dek (3 Mb; 1977-1994) ! plus monthly .dek files through dec04.dek; assembled as my file ! with generic name CMT_dek_file (see assignment statement below for ! actual name, or check in Seismicity.ini). ! See type (10) below for the new .ndk format from GCMT. ! ! (3)JISSIG historical catalog of significant earthquakes since ! start of history in continents, plus large instrumental events ! (on land or at sea) up to 1994. ! ! (4)Nuclear Explosion Data Base from Prototype International ! Data Center (PIDC NEDB), 1945 - 1998. ! ! (5)Any EarthQuake Catalog (.eqc) file previously produced by this program. ! ! (6)Pacheco and Sykes [1992] complete catalog of 698 M_s >= 7.0 ! shallow events in 1900-1989. ! ! (7)Engdahl and Villasensor [2002] "centennial" catalog, 1900-1999, ! with many relocations (1918-1942, 1956-1999) and a variety of ! magnitudes (up to 8 types). 11,956 events, complete down to: ! Ms = 7.0 through 1929; 6.5 through 1963; 5.5 through 1999. ! ! (8)Ekstrom and Nettles [1997] catalog, in CMT .dek format, and using ! CMT methods, for 108 m > 6.0 events in 1976, from input file ! CMT_1976.dek. ! ! (9)ANSS catalog, in file ANSS.catalog.dat as supplied by ! Danijel Schorlemmer of the RELM project, after assigning ! independent-mainshock probabilities to each event ! (but not cutting the aftershocks). ! !(10)Global Centroid Moment Tensor (GCMT) moment-magnitude Mw ! catalog for large (>(4.8~5.7)) events, Jan 1976-Dec 2022. ! See type (2) above for the pre-2005 .dek format from Harvard CMT. ! !(11)ISC_GEM Catalog of Storchak et al. [2012 GEM Tech. Rep.]: ! Thorough revision of ISC 1900-2009, from original sources, ! with relocations, and moment magnitudes. Likely completeness ! thresholds are 7.5 from 1900 and 6.8 from 1918. ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! and extracts the desired spatial/ ! temporal/magnitude rectangle, and writes as a .eqc file, ! in format required for smoothing and scoring by OrbScore, ! or for inclusion in a map by FiniteMap or NeoKineMap. ! Optionally, produces an Adobe Illustrator (.ai) file with ! a simple map of the events selected. ! ! by Peter Bird, Dept. of Earth & Space Sciences, UCLA ! Last revision 2023.10.14 (adding new selection modes for existing .EQC input files) USE DAdobe_Illustrator ! Fortran 90 MODULE DAdobe_Illustrator, in Peter Bird's file DAbobe_Illustrator.f90 USE DMap_Projections ! Fortran 90 MODULE DMap_Projections, in Peter Bird's file DMap_Projections.f90 USE DMap_Tools ! Fortran 90 MODULE DMap_Tools, in Peter Bird's file DMap_Tools.f90 USE DFLIB, ARCQQ => ARC ! provided with Digital Visual Fortran: ! Using GETFILEINFOQQ, which provides names of files ! matching spec.s like "v*.out". Helps user select correct file. ! If no substitute is available on your system when you compile, ! just omit SUBROUTINE File_List (and any CALLs to it). ! However, not using ARC, because I have my own Arc; so I am ! renaming their ARC to ARCQQ to avoid conflicts. INTEGER :: add_hours, add_minutes REAL*8 :: add_time LOGICAL :: all_depths, all_mags, all_years REAL*8 :: alt_eq_mag INTEGER :: alt_eq_mag_int CHARACTER*80 :: ANSS_RELM_file = ' ' LOGICAL :: any_FPS, any_longitude CHARACTER*80 :: answer = ' ', bottom_line = ' ', CMT_dek_file = ' ', GCMT_ndk_file = ' ', CMT1976_file = ' ', eqc_path = ' ', ISC_GEM_file = ' ' CHARACTER*200 :: appended_data = ' ', long_line = ' ' CHARACTER*1 :: asol CHARACTER*1 :: CD_drive, c1 CHARACTER*4 :: c4 CHARACTER*5 :: c5 CHARACTER*7 :: c7 CHARACTER*8 :: c8 INTEGER :: col INTEGER, DIMENSION(200) :: commas ! for ISC_GEM, we need to interpret .csv format. CHARACTER*24 :: date_time_c24 INTEGER :: day REAL*8 :: d0, d1, dxp, dyp REAL*8 :: dep, DT INTEGER :: desired_tectonic_style CHARACTER*80 :: dig_file, dig_pathfile REAL*8 :: E_ELon REAL*8 :: e1_lat, e1_lon, e2_lat, e2_lon, e3_lat, e3_lon INTEGER :: e1_plunge, e1_azimuth, e2_plunge, e2_azimuth, e3_plunge, e3_azimuth REAL*8, DIMENSION(3) :: e1_b_uvec, e1_f_uvec, e2_b_uvec, e2_f_uvec, e3_b_uvec, e3_f_uvec CHARACTER*80 :: Engdahl_file = ' ' REAL*8 :: epicenter_elevation REAL*8 :: epicenter_x_m, epicenter_x_points, epicenter_y_m, epicenter_y_points INTEGER :: eq_year, eq_year_memory CHARACTER*3 :: eq_depth_c3 CHARACTER*8 :: eq_depth_c8 REAL*8 :: eq_depth_real INTEGER :: eq_depth_int CHARACTER*2 :: eq_generic_type_c2 CHARACTER*10 :: eq_moment_over_1E20Nm_c10 REAL*8 :: eq_moment_over_1E20Nm CHARACTER*2 :: eq_month, eq_month_memory, eq_day, eq_day_memory, eq_hour, eq_hour_memory, eq_minute, eq_minute_memory, eq_second CHARACTER*1 :: eq_tenths, eq_lat_flag, eq_lon_flag INTEGER :: eq_months, eq_days, eq_hours, eq_minutes, eq_seconds_int, eq_tenths_int INTEGER :: eq_lat_int, eq_lon_int, eq_mag_int INTEGER :: eq_tectonic_style REAL*8 :: eq_Nlat, eq_Nlat_memory, eq_Elon, eq_Elon_memory, eq_mag, eq_seconds_real CHARACTER*8 :: eq_Nlat_c8, eq_Elon_c8 INTEGER :: events INTEGER :: file_number, filled, first_number, first_year REAL*8 :: glat, glon CHARACTER*80 :: graphics_path = ' ' REAL*8, DIMENSION(:,:), ALLOCATABLE :: grid1 REAL*8 :: grd1_lon_min, grd1_d_lon, grd1_lon_max REAL *8 :: grd1_lat_min, grd1_d_lat, grd1_lat_max INTEGER :: grd1_ncols, grd1_nrows INTEGER :: greg LOGICAL :: got_dig INTEGER :: hr INTEGER :: i CHARACTER*6 :: icat CHARACTER*80 :: in_pathfile = ' ' LOGICAL :: in_ok INTEGER :: int_temp INTEGER :: ios CHARACTER*5 :: isol INTEGER :: last_byte, last_day, last_number, last_year REAL*8 :: lon_range REAL*8 :: mag1 INTEGER :: max_depth, min_depth REAL*8 :: max_epicenter_elevation REAL*8 :: max_Nlat, min_Nlat REAL*8 :: max_mag, min_mag REAL*8 :: min_epicenter_elevation INTEGER :: minu INTEGER :: minutes INTEGER :: mon INTEGER :: M0_exponent REAL*8 :: M0_argument, M0, Mw INTEGER :: m1, m2 REAL*8 :: m7_diam_points, m9_diam_points, m10_diam_points REAL*8 :: North_argument_radians INTEGER :: ntel REAL*8 :: offset_x_m, offset_x_points, offset_y_m, offset_y_points CHARACTER*80 :: old_eqc_file = ' ', new_eqc_file = ' ', new_eqc_pathfile = ' ', Pacheco_file = ' ', PIDC_file = ' ' LOGICAL :: plot_FPS, polygons, problem REAL*8 :: radians_per_point, radius, radius_points INTEGER :: row REAL*8 :: sec LOGICAL :: selecting, select_tectonic_style, select_epicenter_elevation INTEGER :: source REAL*8 :: start_Elon REAL*8 :: t1, t2, t3, t4, t5, t6, t7, t8, t9, t10 REAL*8 :: test_Elon CHARACTER*1 :: time_byte CHARACTER*80 :: top_line = ' ' REAL*8, DIMENSION(3) :: turn_1_uvec, turn_2_uvec, turn_3_uvec, turn_4_uvec, tvec CHARACTER*2 :: two_digits LOGICAL :: valid_FPS REAL*8 :: W_Elon LOGICAL :: want_map, want_titles, whole_Earth REAL*8 :: x1_points, x2_points, x_used_points, y1_points, y2_points, y_used_points REAL*8 :: xbp, ybp, xp, xps, xpt, xp0, xp1, xp2, xp3 REAL*8 :: yp, ypt, yp0, yp1, yp2, yp3 INTEGER :: year CHARACTER*4 :: year_text INTEGER :: yr REAL*8, DIMENSION(3) :: uvec, omega_uvec, result_uvec WRITE (*,"(/' --==++*** SEISMICITY ***++==---')") WRITE (*,"(/' Extracts seismicity catalog from source file(s),')") WRITE (*,"( ' produces an .EQC file, as needed by scoring programs')") WRITE (*,"( ' Kagan_2009_GJI_I_scores, pseudoCSEP, and OrbScore2;')") WRITE (*,"( ' and/or for plotting of seismicity in FiniteMap or NeoKineMap;')") WRITE (*,"( ' and optionally creates an Adobe Illustrator .AI map file.'/)") WRITE (*,"(/' Press [Enter] to initialize...'\)") READ (*,"(A)") c1 !Extract memories from Seismicity.ini, if present in current directory: OPEN (UNIT = 11, FILE = 'Seismicity.ini', STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios == 0) THEN ! Seismicity.ini was found problem = .FALSE. ! may change below READ (11,"(A)",IOSTAT=ios) eqc_path problem = problem.OR.(ios /= 0) READ (11,"(A)",IOSTAT=ios) graphics_path problem = problem.OR.(ios /= 0) READ (11,*) source problem = problem.OR.(ios /= 0) READ (11,"(A1)") CD_drive problem = problem.OR.(ios /= 0) READ (11,"(A)") CMT_dek_file problem = problem.OR.(ios /= 0) READ (11,"(A)") CMT1976_file problem = problem.OR.(ios /= 0) READ (11,"(A)") Pacheco_file problem = problem.OR.(ios /= 0) READ (11,"(A)") Engdahl_file problem = problem.OR.(ios /= 0) READ (11,"(A)") PIDC_file problem = problem.OR.(ios /= 0) READ (11,"(A)") old_eqc_file problem = problem.OR.(ios /= 0) READ (11,"(A)") GCMT_ndk_file problem = problem.OR.(ios /= 0) READ (11,"(A)") ISC_GEM_file problem = problem.OR.(ios /= 0) READ (11,*) all_depths problem = problem.OR.(ios /= 0) READ (11,*) min_depth problem = problem.OR.(ios /= 0) READ (11,*) max_depth problem = problem.OR.(ios /= 0) READ (11,*) whole_Earth problem = problem.OR.(ios /= 0) READ (11,*) start_Elon problem = problem.OR.(ios /= 0) READ (11,*) lon_range problem = problem.OR.(ios /= 0) READ (11,*) min_Nlat problem = problem.OR.(ios /= 0) READ (11,*) max_Nlat problem = problem.OR.(ios /= 0) READ (11,*) all_mags problem = problem.OR.(ios /= 0) READ (11,*) min_mag problem = problem.OR.(ios /= 0) READ (11,*) max_mag problem = problem.OR.(ios /= 0) READ (11,*) all_years problem = problem.OR.(ios /= 0) READ (11,*) first_year problem = problem.OR.(ios /= 0) READ (11,*) last_year problem = problem.OR.(ios /= 0) READ (11,"(A)") new_eqc_file problem = problem.OR.(ios /= 0) READ (11,*) want_map problem = problem.OR.(ios /= 0) READ (11,*) got_dig problem = problem.OR.(ios /= 0) READ (11,"(A)") dig_file problem = problem.OR.(ios /= 0) READ (11,*) polygons problem = problem.OR.(ios /= 0) READ (11,*) plot_FPS problem = problem.OR.(ios /= 0) READ (11,*) min_mag problem = problem.OR.(ios /= 0) READ (11,*) m7_diam_points problem = problem.OR.(ios /= 0) READ (11,*) m9_diam_points problem = problem.OR.(ios /= 0) READ (11,*) m10_diam_points problem = problem.OR.(ios /= 0) READ (11,*) minutes problem = problem.OR.(ios /= 0) CLOSE(11) IF (problem) THEN WRITE (*,"(/' ERROR: Bad data, bad format, or missing lines in Seismicity.ini.')") WRITE (*,"( ' The easiest way to recover from this is to:')") WRITE (*,"( ' (1) Print out Seismicity.ini')") WRITE (*,"( ' (2) Delete Seismicity.ini')") WRITE (*,"( ' (3) Restart Seismicity, and enter your choices manually.')") WRITE (*,"( ' Press [Enter] when ready...'\)") READ (*,"(A)") c1 STOP ' ' END IF ELSE ! no .ini file; use defaults CLOSE (11, IOSTAT = ios) ! just to be sure eqc_path = ' ' graphics_path = ' ' source = 1 CD_drive = 'E' CMT_dek_file = 'CMT_1977_through_2004.dek' CMT1976_file = 'CMT1976.dek' Pacheco_file = 'Pacheco_and_Sykes_1992.dat' Engdahl_file = 'Centennial.dat' PIDC_file = 'PIDC_nuclear.dat' old_eqc_file = ' ' GCMT_ndk_file = 'GCMT_1976_through_2022.ndk' ISC_GEM_file = 'isc-gem-cat.csv' all_depths = .TRUE. min_depth = -20 max_depth = 1000 whole_Earth = .TRUE. start_Elon = -180.0D0 lon_range = 360.0D0 min_Nlat = -90.0D0 max_Nlat = +90.0D0 all_mags = .TRUE. min_mag = 0.0D0 max_mag = 10.0D0 all_years = .TRUE. first_year = -5000 last_year = 9999 new_eqc_file = 'my_choice.eqc' want_map = .TRUE. got_dig = .TRUE. dig_file = 'WorldMap.dig' polygons = .FALSE. plot_FPS = .TRUE. m7_diam_points = 12.0D0 m9_diam_points = 12.0D0 m10_diam_points = 12.0D0 minutes = 300 END IF ! .ini file, or defaults? !-------------------------(Define Paths)----------------------------- WRITE (*,"(//' ----------------------------------------------------------------------'& &/' Setting PATHS to Seismicity and Graphics Files'& &//' Seismicity stores its memory in Seismicity.ini, which is'& &/' placed in the current directory when Seismicity is run.'& &/' Normally, this should be the directory holding Seismicity.exe.'& &//' However, users may wish to keep seismic catalog files (new and old)'& &/' and seismicity-map graphics files separate in other directories.'& &//' When entering the paths below, you may include computer and drive'& &/' identifiers according to the conventions of your system.'& &/' In DOS or Windows, paths should end in ''\''.'& &/' In Unix, paths should end in ''/''.'& &//' PLEASE TYPE PATH NAMES CAREFULLY; there is no way to validate or'& &/' correct them using standard Fortran 90; errors will crash'& &/' Seismicity (at least)!'& &/' ----------------------------------------------------------------------')") CALL DPrompt_for_String('What is the path to the hard-disk directory where seismic catalog files (both existing and new) are stored?',eqc_path,eqc_path) eqc_path = ADJUSTL(eqc_path) CALL DPrompt_for_String('What is the path to the hard-disk directory where .ai graphics files (both existing and new) are stored?',graphics_path,graphics_path) graphics_path = ADJUSTL(graphics_path) WRITE (*,"(' IT WILL NOT BE NECESSARY TO TYPE THESE PATHS AGAIN!')") !-------------------------(end of defining paths)-------------------- 1 WRITE (*,*) WRITE (*,"(' SELECT THE DATA SOURCE:')") WRITE (*,"(' (1) International Seismological Center (ISC)')") WRITE (*,"(' catalog of body-wave magnitudes (mb) from ~3')") WRITE (*,"(' up to saturation at ~6.8; includes many events')") WRITE (*,"(' with magnitude zero (undetermined).')") WRITE (*,"(' Version on CD-ROM covers 1964-1991; complete from 1971.')") WRITE (*,"(' (2) Harvard Centroid Moment Tensor (CMT) .dek-format catalog')") WRITE (*,"(' of centroid moment tensors from events of moment-magnitude')") WRITE (*,"(' from 4.8~5.6 (incomplete detection)) and 5.7~9.5 (complete')") WRITE (*,"(' detection); can be accurately converted to moment magnitudes.')") WRITE (*,"(' Can also be plotted as fault-plane solutions.')") WRITE (*,"(' Files downloaded from Harvard cover 1977-2004.')") WRITE (*,"(' See new .ndk format, and later years, under (10) below.')") WRITE (*,"(' (3) JISSIG historical catalog of significant earthquakes,')") WRITE (*,"(' with coverage to 3000 B.C. in Asia, and to beginning')") WRITE (*,"(' of history on other continents. (Short history at sea; no')") WRITE (*,"(' events listed in Antarctica. Some early events lack date')") WRITE (*,"(' or even location! Almost half have Mw = 0.0!) By Utsu (1996).')") WRITE (*,"(' (4) Nuclear Explosion Data Base from Prototype International')") WRITE (*,"(' Data Center, with ~2000 events 1945-1998. Many lack origin')") WRITE (*,"(' times; 46% lack mb magnitudes.')") CALL Pause() WRITE (*,"(' (5) any EarthQuake Catalog (.eqc) file created by this program.')") WRITE (*,"(' (6) Pacheco and Sykes [1992] complete catalog of 698 Ms >= 7.0')") WRITE (*,"(' shallow earthquakes, 1900-1989, with generic mechanism types.')") WRITE (*,"(' All events have moment-magnitude. No relocations performed.')") WRITE (*,"(' (7) Engdahl and Villasensor [2002] Centennial catalog, 1900-1999,')") WRITE (*,"(' with many relocations (1918-1942, 1956-1999) and a variety of')") WRITE (*,"(' magnitudes (up to 8 kinds). 11,956 events, complete down to:')") WRITE (*,"(' M_s = 7.0 through 1929; 6.5 through 1963; 5.5 through 1999.')") WRITE (*,"(' (8) Ekstrom and Nettles [1997] extended CMT catalog for 1976:')") WRITE (*,"(' centroid moment tensors from 108 events of moment magnitude')") WRITE (*,"(' above 6.00; can be accurately converted to moment magnitudes.')") WRITE (*,"(' Can also be plotted as fault-plane solutions.')") WRITE (*,"(' File obtained 2003.01.21 from Yan Y. Kagan.')") WRITE (*,"(' (9) ANSS catalog, in file ANSS.catalog.dat as supplied by')") WRITE (*,"(' Danijel Schorlemmer of the RELM project, after assigning')") WRITE (*,"(' independent-mainshock probabilities to each event')") WRITE (*,"(' (but not cutting the aftershocks).')") CALL Pause() WRITE (*,"(' (10) Global Centroid Moment Tensor (GCMT) .ndk-format catalog')") WRITE (*,"(' of centroid moment tensors from events of moment-magnitude')") WRITE (*,"(' from 4.8~5.6 (incomplete detection)) and 5.7~9.5 (complete')") WRITE (*,"(' detection); can be accurately converted to moment magnitudes.')") WRITE (*,"(' Can also be plotted as fault-plane solutions.')") WRITE (*,"(' Files downloaded from Harvard cover 1976-2022.')") WRITE (*,"(' For data from the old .dek format before 2005, see (2) above.')") WRITE (*,"(' (11) ISC_GEM Catalog of Storchak et al. [2012 GEM Tech. Rep.]:')") WRITE (*,"(' Thorough revision of ISC 1900-2009, from original sources,')") WRITE (*,"(' with relocations, and moment magnitudes. Likely completeness')") WRITE (*,"(' thresholds are 7.5 from 1900 and 6.8 from 1918.')") CALL DPrompt_for_Integer('SELECT BY INTEGER CODE:',source,source) IF ((source > 11).OR.(source < 1)) THEN WRITE (*,"(' ERROR: Select integer in range of 1 to 11 !!!')") CALL Pause() GO TO 1 END IF WRITE (*,"(' ')") SELECT CASE (source) CASE (1) WRITE (*,"(' Please load the CD-ROM labelled:')") WRITE (*,"(' Seismicity Catalogs, Volume 2: Global and Regional, 2150 B.C.-1996 A.D.')") WRITE (*,"(' that was published by the National Geophysical Data Center,')") WRITE (*,"(' Boulder, Colorado, USA in March 1996.')") WRITE (*,"(' When you have finished, answer the following question:')") CALL DPrompt_for_String ('What letter identifies your CD-ROM drive?', CD_drive, answer) CD_drive = answer(1:1) CASE (2) CALL DPrompt_for_String('What file contains the CMT .dek catalog?',CMT_dek_file,CMT_dek_file) CASE (3) WRITE (*,"(' Please load the CD-ROM labelled:')") WRITE (*,"(' Seismicity Catalogs, Volume 2: Global and Regional, 2150 B.C.-1996 A.D.')") WRITE (*,"(' that was published by the National Geophysical Data Center,')") WRITE (*,"(' Boulder, Colorado, USA in March 1996.')") WRITE (*,"(' When you have finished, answer the following question:')") CALL DPrompt_for_String ('What letter identifies your CD-ROM drive?', CD_drive, answer) CD_drive = answer(1:1) CASE (4) CALL DPrompt_for_String('What file contains the PIDC NEDB catalog of nuclear explosions?',PIDC_file,PIDC_file) CASE (5) !CALL File_List( any = .FALSE., & ! & basemap = .FALSE., & ! & catalog = .TRUE. , & ! & fegrid = .FALSE., & ! & forces = .FALSE., & ! & gridded_data = .FALSE., & ! & parameters = .FALSE., & ! & velocity = .FALSE. , & ! & suggested_file = old_eqc_file ) CALL DPrompt_for_String('What .eqc file should be read?',old_eqc_file,old_eqc_file) CASE (6) CALL DPrompt_for_String('What file contains the Pacheco and Sykes [1992] catalog?',Pacheco_file,Pacheco_file) CASE (7) CALL DPrompt_for_String('What file contains the Engdahl & Villasensor [2002] catalog?',Engdahl_file,Engdahl_file) CASE (8) CALL DPrompt_for_String('What file contains the Ekstrom & Nettles [1997] catalog?',CMT1976_file,CMT1976_file) CASE (9) CALL DPrompt_for_String('What file contains the ANSS catalog in RELM format?',ANSS_RELM_file,ANSS_RELM_file) CASE (10) CALL DPrompt_for_String('What file contains the GCMT .ndk catalog?',GCMT_ndk_file,GCMT_ndk_file) CASE (11) CALL DPrompt_for_String('What file contains the ISC_GEM catalog in .csv?',ISC_GEM_file,ISC_GEM_file) END SELECT ! source 10 IF (source == 4) THEN ! depth selection is inappropriate! all_depths = .TRUE. ELSE ! ask user WRITE (*,*) CALL DPrompt_for_Logical('Do you want all depths?',all_depths,all_depths) END IF ! need to ask user about depths? IF (all_depths) THEN min_depth = -20 max_depth = 1000 ELSE ! select depths IF ((source == 10).OR.(source == 2).OR.(source == 8)) THEN ! GCMT warning WRITE (*, "(' According to Goran Ekstrom (pers. comm., 2002),')") WRITE (*, "(' [G]CMT depths are measured from sea level in the PREM')") WRITE (*, "(' model which has 3 km ocean and 21.4 km crust; thus')") WRITE (*, "(' 10 km means 7 km below the water/rock interface.')") WRITE (*, "(' He also reports that depths were constrained to be at')") WRITE (*, "(' least 10 km in years 1977-1992, then later required to be')") WRITE (*, "(' at least 15 km. Depths which do not converge are often')") WRITE (*, "(' fixed at the NEIC standard shallow depth of 33 km.')") WRITE (*, "(' Therefore, depths have poor resolution, and a request for')") WRITE (*, "(' a depth subset of shallow earthquakes is likely to lead to')") WRITE (*, "(' sparse and unsatisfactory results.')") ELSE IF (source == 6) THEN ! Pacheco and Sykes [1992] warning WRITE (*, "(' WARNING: 388 (46%) of the events in this catalog')") WRITE (*, "(' have no depth information, and depth will default to 0.')") WRITE (*, "(' Actual determined depths range from 2 to 66 km.')") WRITE (*, "(' Therefore, depths have poor resolution, and a request for')") WRITE (*, "(' a depth subset of shallow earthquakes is likely to lead to')") WRITE (*, "(' sparse and unsatisfactory results.')") ELSE IF (source == 7) THEN ! Engdahl & Villasensor [2002] notice WRITE (*, "(' NOTE: Reported depths range 0-720 km, with 0.1 km precision.')") WRITE (*, "(' 78% of events are listed as 0-70 km (shallow).')") WRITE (*, "(' Of these, 6% defaulted to 0 km, 7% to 15 km, and 9% to')") WRITE (*, "(' 35 km. Many of these are the events that were not relocated')") WRITE (*, "(' (1900-17, 1943-55). The rest of the histogram looks realistic.')") END IF CALL DPrompt_for_Integer('Enter minimum depth in kilometers',min_depth, min_depth) CALL DPrompt_for_Integer('Enter maximum depth in kilometers',max_depth, max_depth) END IF ! all_depths, or not IF (max_depth < min_depth) THEN WRITE (*,"(' ERROR: Maximum depth must be >= minimum.')") GO TO 10 END IF ! impossible window WRITE (*,*) CALL DPrompt_for_Logical ('Do you want the entire Earth?', whole_Earth, whole_Earth) IF (whole_Earth) THEN start_Elon = -180.0D0 lon_range = 360.0D0 any_longitude = .TRUE. min_Nlat = -90.0D0 max_Nlat = +90.0D0 ELSE ! part of Earth IF (start_Elon > 180.0) THEN W_Elon = start_Elon - 360.00D0 ELSE W_Elon = start_Elon END IF E_Elon = W_Elon + lon_range CALL DPrompt_for_Real('Enter West-side longitude (East = + )', W_Elon,W_Elon) start_Elon = DEasting(W_Elon) ! in range 0.0 ... 359.999 CALL DPrompt_for_Real('Enter East-side longitude (East = + )', E_Elon,E_Elon) IF (MOD(ABS(E_Elon - W_Elon), 360.0D0) == 0.0D0) THEN lon_Range = 360.0D0 any_longitude = .TRUE. ELSE lon_Range = DEasting(E_Elon - W_Elon) any_longitude = .FALSE. END IF 14 CALL DPrompt_for_Real('Enter Southernmost latitude (North = +)', min_Nlat,min_Nlat) CALL DPrompt_for_Real('Enter Northernmost latitude (North = +)', max_Nlat,max_Nlat) IF (max_Nlat < min_Nlat) THEN WRITE (*,"(' ERROR: Maximum latitude must be >= minimum.')") GO TO 14 END IF ! impossible window END IF ! whole_Earth, or not 16 WRITE (*,*) CALL DPrompt_for_Logical('Do you want all magnitudes?',all_mags,all_mags) IF (all_mags) THEN min_mag = 0.0D0 max_mag = +999.9D0 !Note: Allow all events to go into .eqc file; this does ! NOT mean that all small ones must be plotted, or ! that legend will have too many samples. ELSE ! select magnitudes SELECT CASE (source) CASE (1) ! ISC WRITE (*,"(' In answering the following, remember that ISC reports')") WRITE (*,"(' body-wave magnitudes mb, which rarely exceed 6.8!')") WRITE (*,"(' At the low end, this catalog has many zero (undetermined) ')") WRITE (*,"(' magnitudes; real values begin with ~3.')") CASE (2) ! CMT .dek format WRITE (*,"(' In answering the following, remember that Harvard reports')") WRITE (*,"(' moment-magnitudes Mw for events ~4.8 and up,')") WRITE (*,"(' but that event detection is only complete for 5.7 and up.')") CASE (3) ! JISSIG CASE (4) ! nuclear explosions WRITE (*,"(' In answering the following, remember that magnitudes')") WRITE (*,"(' were not reported by PIDC for 46% of events;')") WRITE (*,"(' in these cases magnitude defaults to zero.')") CASE (5) ! old_eqc_file CASE (6) ! Pacheco and Sykes [1992] WRITE (*,"(' In answering the following, remember that this catalog')") WRITE (*,"(' is nominally for Ms of 7.0 and larger, but that actual')") WRITE (*,"(' selection magnitude will be computed from seismic moment.')") CASE (7) ! Engdahl and Villasensor [2002] WRITE (*,"(' In answering the following, remember that this catalog')") WRITE (*,"(' has multiple magnitudes per event. Selection will be')") WRITE (*,"(' according to the first entry, which will be mw (if')") WRITE (*,"(' available), or Ms, or mB, or mb, etc.')") WRITE (*,"(' This catalog has entries down to 5.5, but is only')") WRITE (*,"(' complete down to Ms = 7.0.')") CASE (8) ! Ekstrom & Nettles [1997] extended CMT catalog for 1976 WRITE (*,"(' In answering the following, remember that this 1976 catalog')") WRITE (*,"(' only includes moment magnitudes of 6.00 and up.')") CASE (10) ! GCMT .ndk format WRITE (*,"(' In answering the following, remember that GCMT reports')") WRITE (*,"(' moment-magnitudes Mw for events ~4.8 and up,')") WRITE (*,"(' but that event detection is only complete for 5.7 and up.')") CASE (11) ! ISC_GEM WRITE (*,"(' In answering the following, consider Michael [2014, BSSA]:')") WRITE (*,"(' Completeness threshold likely 6.8 from 1918,')") WRITE (*,"(' and higher than 7.5 before that.')") END SELECT CALL DPrompt_for_Real('Enter minimum magnitude',min_mag,min_mag) CALL DPrompt_for_Real('Enter maximum magnitude',max_mag,max_mag) END IF ! all_mags, or not IF (max_mag < min_mag) THEN WRITE (*,"(' ERROR: Maximum magnitude must be >= minimum.')") GO TO 16 END IF ! impossible window 18 WRITE (*,*) CALL DPrompt_for_Logical('Do you want all available years?',all_years,all_years) IF (all_years) THEN SELECT CASE (source) CASE (1) ! ISC first_year = 1964 !Must be accurate, as used to construct filenames. last_year = 1991 CASE (2) ! CMT .dek format first_year = 1977 last_year = 2004 CASE (3) ! JISSIG first_year = -5000 !No reason for accuracy; allow all data. last_year = 9999 CASE (4) ! PIDC NEBD first_year = 1945 !No reason for accuracy; allow all data. last_year = 9999 CASE (5) ! old_eqc_file first_year = -5000 !No reason for accuracy; allow all data. last_year = 9999 CASE (6) ! Pacheco and Sykes [1992] first_year = 1900 !Just a comment, since all data are in one file. last_year = 1989 CASE (7) ! Engdahl and Villasensor [2002] first_year = 1900 !Just a comment, since all data are in one file. last_year = 1999 CASE (8) ! Ekstrom & Nettles [1997] extended CMT catalog for 1976 first_year = 1976 last_year = 1976 CASE (9) ! ANSS first_year = -5000 !No reason for accuracy; allow all data. last_year = 9999 CASE (10) ! GCMT .ndk format first_year = 1976 last_year = 2022 CASE (11) ! ISC_GEM catalog first_year = 1900 last_year = 2009 END SELECT ELSE ! select years SELECT CASE (source) CASE (1) first_year = MAX(first_year, 1964) CASE (2) first_year = MAX(first_year, 1977) CASE (4) first_year = MAX(first_year, 1945) CASE (6) first_year = MAX(first_year, 1900) CASE (7) first_year = MAX(first_year, 1900) CASE (8) first_year = 1976 last_year = 1976 CASE (10) first_year = MAX(first_year, 1976) CASE (11) first_year = MAX(first_year, 1900) last_year = MIN(last_year, 2009) END SELECT CALL DPrompt_for_Integer('Enter first year',first_year,first_year) CALL DPrompt_for_Integer('Enter last year',last_year,last_year) END IF ! all_years, or not IF (last_year < first_year) THEN WRITE (*,"(' ERROR: Last year must be >= first year.')") GO TO 18 END IF ! impossible window IF (source == 5) THEN ! input .EQC file WRITE (*, *) CALL DPrompt_for_Logical('Select by tectonic style (normal, strike-slip, thrust)?', .FALSE., select_tectonic_style) IF (select_tectonic_style) THEN WRITE (*, *) CALL DPrompt_for_Integer('Choose: 1 = normal; 2 = strike-slip; 3 = thrust:', 0, desired_tectonic_style) END IF WRITE (*, *) CALL DPrompt_for_Logical('Select by elevation at epicenter?', .FALSE., select_epicenter_elevation) IF (select_epicenter_elevation) THEN WRITE (*, *) WRITE (*, "(' Note that database ETOPO5.GRD *must* be available in current directory!')") WRITE (*, *) CALL Pause() ! ------------------------------------------------------------------------------------------------------- OPEN (UNIT = 21, FILE = "ETOPO5.GRD", STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') READ (21, *, IOSTAT = ios) grd1_lon_min, grd1_d_lon, grd1_lon_max READ (21, *, IOSTAT = ios) grd1_lat_min, grd1_d_lat, grd1_lat_max grd1_ncols = 1 + NINT((grd1_lon_max - grd1_lon_min) / grd1_d_lon) grd1_nrows = 1 + NINT((grd1_lat_max - grd1_lat_min) / grd1_d_lat) ALLOCATE ( grid1(grd1_nrows, grd1_ncols) ) READ (21, *, IOSTAT = ios) ((grid1(i,j), j = 1, grd1_ncols), i = 1, grd1_nrows) CLOSE (21) WRITE (*, "(' File ETOPO5.GRD has been read and memorized.')") CALL Pause() ! ------------------------------------------------------------------------------------------------------- WRITE (*, *) CALL DPrompt_for_REAL('Minimum elevation (in meters above sea-level):', 0.0D0, min_epicenter_elevation) CALL DPrompt_for_REAL('Maximum elevation (in meters above sea-level):', 0.0D0, max_epicenter_elevation) ELSE min_epicenter_elevation = -99999.9D0 max_epicenter_elevation = 99999.9D0 END IF ELSE select_tectonic_style = .FALSE. select_epicenter_elevation = .FALSE. END IF 20 WRITE (*,*) CALL DPrompt_for_String('Enter name for output (.eqc) file:',new_eqc_file,new_eqc_file) new_eqc_pathfile = TRIM(eqc_path) // new_eqc_file OPEN (UNIT = 2, FILE = new_eqc_pathfile, STATUS = 'NEW', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: File already exists (or disk full)!')") GO TO 20 END IF ! problem opening out_file events = 0 WRITE (*,*) ! LOOP on input files, if more than one SELECT CASE (source) CASE (1) ! ISC first_number = MAX(first_year, 1964) ! to avoid crash from trying to open nonexistent files! last_number = last_year ! don't limit this, in case new CD-ROM comes out later; we will just jump out of loop CASE (2:11) ! CMT, JISSIG, NEDB, .eqc, Pacheco & Sykes [1992], Engdahl & Villasensor [2002], Ekstrom & Nettles [1997], GCMT, ISC_GEM first_number = 1 last_number = 1 END SELECT ! source any_FPS = .FALSE. ! unless changed during the data-reading filing: DO file_number = first_number, last_number SELECT CASE (source) CASE (1) ! ISC WRITE (year_text, "(I4)") file_number two_digits = year_text(3:4) in_pathfile = CD_drive // ':\Data\Globecat\ISC\' // two_digits & & // '\' // two_digits // '.dat' CASE (2) ! CMT .dek format, through end 2004 in_pathfile = TRIM(eqc_path) // CMT_dek_file CASE (3) ! JISSIG in_pathfile = CD_drive // ':\Data\GlobeCat\JISSIG\JISSIG.dat' CASE (4) ! NEDB from PIDC in_pathfile = TRIM(eqc_path) // PIDC_file CASE (5) ! .eqc in_pathfile = TRIM(eqc_path) // old_eqc_file CASE (6) ! Pacheco and Sykes [1992] in_pathfile = TRIM(eqc_path) // Pacheco_file CASE (7) ! Engdahl and Villasensor [2002] in_pathfile = TRIM(eqc_path) // Engdahl_file CASE (8) ! Ekstrom & Nettles [1997] extended CMT catalog for 1976 in_pathfile = TRIM(eqc_path) // CMT1976_file CASE (9) ! ANSS catalog reprocessed by Danijel Schorlemmer of RELM to assign mainshock-independence probabilities in_pathfile = TRIM(eqc_path) // ANSS_RELM_file CASE (10) ! GCMT .ndk format, 2005+ in_pathfile = TRIM(eqc_path) // GCMT_ndk_file CASE (11) ! ISC_GEM [Storchak et al., 2012 GEM Tech. Rep.] in_pathfile = TRIM(eqc_path) // ISC_GEM_file END SELECT ! source OPEN (UNIT = 1, FILE = in_pathfile, STATUS = 'OLD', IOSTAT = ios, & & PAD = 'YES') ! Padding is needed since some old .eqc files have FPS and/or other appended data; some not! IF (ios == 0) THEN WRITE (*,"(' Reading ',A)") TRIM(in_pathfile) ELSE ! trouble opening in_file WRITE (*,"(' ',A,' could not be found.')") TRIM(in_pathfile) WRITE (*,"(' Terminating loop on multiple (or single) input file(s).')") EXIT filing END IF !Waste any header lines: IF (source == 11) THEN ! ISC_GEM DO i = 1, 57 READ( 1, *) END DO END IF !LOOP on records in input file next_line: DO SELECT CASE (source) CASE (1) ! ISC READ (1,100,IOSTAT=ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_lat_int, eq_lat_flag, eq_lon_int, eq_lon_flag, & & eq_depth_int, eq_mag_int, alt_eq_mag_int ! Note: Before 1971, magnitude was usually in ! columns 61-63, followed by "ISC" in 64-66. ! From 1971 on, it is usually in columns 36-38, ! followed by "MB" in 39-40. 100 FORMAT(4X,I4,A2,A2, & & A2,A2,A2,A1, & & I5,A1,I6,A1, & & I3,I3,22X,I3) IF (ios /= 0) EXIT next_line eq_Elon = eq_lon_int/1000.0D0 IF (eq_lon_flag == 'W') eq_Elon = - eq_Elon eq_Nlat = eq_lat_int / 1000.0D0 IF (eq_lat_flag == 'S') eq_Nlat = - eq_Nlat eq_mag = eq_mag_int / 100.0D0 alt_eq_mag = alt_eq_mag_int / 100.0D0 IF (eq_mag <= 0.0D0) eq_mag = alt_eq_mag CASE (2, 8) ! CMT .dek format, or Ekstrom & Nettles [1997] extended CMT for 1976 READ (1,201,IOSTAT=ios) & ! 1st line of each block of 4 lines. & eq_month, eq_day, eq_year, & & eq_hour, eq_minute, eq_second, eq_tenths !Note: The origin time in this line was determined by NEIC, ISC, or PDE using body waves; ! the origin time will be revised by lag DT in lines below. ! The latitude, longitude, and depth in this line are also from NEIC, ISC, or PDE, ! and will not be used; the CMT centroid position will be read in line 2. 201 FORMAT(9X,A2,1X,A2,1X,I2,1X, & & A2,':',A2,':',A2,'.',A1) IF (ios /= 0) EXIT next_line IF (eq_year > 60) THEN eq_year = eq_year + 1900 ELSE eq_year = eq_year + 2000 END IF READ (1,202,IOSTAT=ios) & ! 2nd line of each block of 4 lines & DT, eq_Nlat, eq_Elon, eq_depth_real 202 FORMAT(33X,F6.1,4X,F7.2,5X,F8.2,5X,F6.1) IF (ios /= 0) EXIT next_line !--- add DT seconds to origin time, possibly adjusting minutes & hours! --------- READ (eq_second,"(I2)") eq_seconds_int READ (eq_tenths,"(I1)") eq_tenths_int eq_seconds_real = eq_seconds_int + 0.1D0 * eq_tenths_int + DT more_minutes: DO IF (eq_seconds_real >= 59.95D0) THEN eq_seconds_real = MAX(0.0D0, eq_seconds_real - 60.0D0) READ (eq_minute,"(I2)") eq_minutes eq_minutes = eq_minutes + 1 IF (eq_minutes >= 60) THEN eq_minutes = eq_minutes - 60 READ (eq_hour,"(I2)") eq_hours eq_hours = eq_hours + 1 WRITE (eq_hour,"(I2)") eq_hours !Note: Not allowing date to change (too confusing); in very rare cases hour may reach 24! END IF ! eq_minutes reached 60 WRITE (eq_minute,"(I2)") eq_minutes ELSE EXIT more_minutes END IF ! seconds have exceeded 60, or not END DO more_minutes less_minutes: DO IF (eq_seconds_real < 0.0D0) THEN eq_seconds_real = MIN(59.9D0, eq_seconds_real + 60.0D0) READ (eq_minute,"(I2)") eq_minutes eq_minutes = eq_minutes - 1 IF (eq_minutes < 0) THEN eq_minutes = eq_minutes + 60 READ (eq_hour,"(I2)") eq_hours eq_hours = eq_hours - 1 WRITE (eq_hour,"(I2)") eq_hours !Note: Not allowing date to change (too confusing); in very rare cases hour may reach -1! END IF ! eq_minutes reached 60 WRITE (eq_minute,"(I2)") eq_minutes ELSE EXIT less_minutes END IF ! seconds are negative, or not END DO less_minutes WRITE (c4,"(F4.1)") eq_seconds_real eq_second = c4(1:2) eq_tenths = c4(4:4) !-------------------------------------------------------------------------------- IF (eq_hour(1:1) == ' ') eq_hour(1:1) = '0' IF (eq_hour(2:2) == ' ') eq_hour(2:2) = '0' IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' IF (eq_minute(2:2) == ' ') eq_minute(2:2) = '0' IF (eq_second(1:1) == ' ') eq_second(1:1) = '0' IF (eq_second(2:2) == ' ') eq_second(2:2) = '0' eq_depth_int = NINT(eq_depth_real) READ (1,203,IOSTAT=ios) M0_exponent ! 3rd line of each block of 4 lines 203 FORMAT(12X,I2) IF (ios /= 0) EXIT next_line READ (1,204,IOSTAT=ios) & ! 4th line & e3_plunge, e3_azimuth, e2_plunge, e2_azimuth, e1_plunge, e1_azimuth, & & M0_argument ! e3 is extension, e1 is shortening, e2 neutral 204 FORMAT(3(8X,I2,1X,I3),2X,F5.2) IF (ios /= 0) EXIT next_line any_FPS = .TRUE. M0 = M0_argument * 10.0D0**M0_exponent IF (M0 > 0.0D0) THEN Mw = (DLOG10(M0) - 9.05D0 - 7.0D0) / 1.5D0 !Hanks and Kanamori (1979; J.G.R., 2348-2350) !The -7. term is to convert from dyne-cm !(in the exponent of the Harvard CMT catalog) !to N-m (required in the Hanks and Kanamori formula). eq_mag = Mw ELSE eq_mag = 0.0D0 END IF CASE (3) ! JISSIG READ (1, 300, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, time_byte, & & eq_Nlat, eq_Elon, & & eq_depth_c3, eq_mag !Note: JISSIG.fmt implies eq_mag is Mw. 300 FORMAT ( 7X, & & I5,1X,A2,A2, 1X, & & A2,A2,1X,A1, 1X, & & F6.2, 1X, F7.2, 1X, & & A3, 1X, F3.1) IF (ios /= 0) THEN IF (EOF(1)) THEN EXIT next_line ELSE ! bad line within file (e.g., at 1000 A.D.!) CYCLE next_line END IF ! at EOF, or not END IF ! READ went badly !Discard worthless events with no location IF ((eq_Nlat == 0.0D0).AND.(eq_Elon == 0.0D0)) CYCLE next_line !If date within year is missing, don't bother to refine the time! IF ((eq_month == ' ').AND.(eq_day == ' 0')) THEN eq_month = '00' eq_day = '00' eq_hour = '00' eq_minute = '00' ELSE ! date is given IF (eq_month(1:1) == ' ') eq_month(1:1) = '0' IF (eq_day(1:1) == ' ') eq_day(1:1) = '0' IF (eq_hour == ' ') eq_hour = ' 0' IF (eq_hour(1:1) == ' ') eq_hour(1:1) = '0' IF (eq_minute == ' ') eq_minute = ' 0' IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' !If time was given, but in local time, then !correct local time to universal time (roughly!) IF ((eq_hour /= '00').OR.(eq_minute /= '00')) THEN IF (time_byte == 'L') THEN ! local time add_time = -eq_Elon / 15.0D0 ! "sun time", not "time zones" add_hours = add_time ! integer = real truncation add_minutes = NINT(60.0D0 * (add_time - add_hours)) !create integer equivalents of ascii values READ (eq_minute, *) eq_minutes READ (eq_hour, * ) eq_hours READ (eq_day, *) eq_days READ (eq_month, *) eq_months eq_minutes = eq_minutes + add_minutes IF (eq_minutes >= 60) THEN eq_minutes = eq_minutes - 60 eq_hours = eq_hours + 1 END IF IF (eq_minutes < 0) THEN eq_minutes = eq_minutes + 60 eq_hours = eq_hours - 1 END IF eq_hours = eq_hours + add_hours IF (eq_hours >= 24) THEN eq_hours = eq_hours - 24 eq_days = eq_days + 1 END IF IF (eq_hours < 0) THEN eq_hours = eq_hours + 24 eq_days = eq_days - 1 END IF IF (eq_days <= 0) THEN ! go to previous month IF (eq_months > 1) THEN eq_months = eq_months - 1 ELSE ! current month is January, go to previous year eq_months = 12 eq_year = eq_year - 1 END IF last_day = Days_Long (eq_months, eq_year) eq_days = eq_days + last_day END IF last_day = Days_Long (eq_months, eq_year) IF (eq_days > last_day) THEN ! go to next month IF (eq_months < 12) THEN eq_months = eq_months + 1 ELSE ! current month is December; go to next year eq_months = 1 eq_year = eq_year + 1 END IF eq_days = eq_days - last_day ! at least 1 END IF !convert back from integer to ascii WRITE (eq_month,"(I2)") eq_months IF (eq_month(1:1) == ' ') eq_month(1:1) = '0' WRITE (eq_day,"(I2)") eq_days IF (eq_day(1:1) == ' ') eq_day(1:1) = '0' WRITE (eq_hour,"(I2)") eq_hours IF (eq_hour(1:1) == ' ') eq_hour(1:1) = '0' WRITE (eq_minute,"(I2)") eq_minutes IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' END IF ! Local time was reported END IF ! time was reported within the day END IF ! no date given, or date IS given !seconds and tenths are never given in this catalog: eq_second = '00' eq_tenths = '0' !convert left-justified depth string to integer: IF (eq_depth_c3 == ' ') THEN eq_depth_int = 0 ELSE READ (eq_depth_c3, *, IOSTAT = ios) eq_depth_int IF (ios /= 0) eq_depth_int = 0 ! handles " *" END IF CASE (4) ! NEDB from PIDC READ (1, 400, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Nlat, eq_lat_flag, eq_Elon, eq_lon_flag, & & eq_mag !Note: PIDC_nuclear.dat header implies eq_mag is mb. 400 FORMAT (16X, & & I4,5X,A2,1X,A2, 4X, & & A2,1X,A2,1X,A2,1X,A1, 3X, & & F6.3,A1,2X,F7.3,A1, 3X, & & F3.1) IF (ios /= 0) THEN IF (EOF(1)) THEN EXIT next_line ELSE ! bad line within file (e.g., at 1000 A.D.!) CYCLE next_line END IF ! at EOF, or not END IF ! READ went badly IF (eq_lon_flag == 'W') eq_Elon = -eq_Elon IF (eq_lat_flag == 'S') eq_Nlat = -eq_Nlat !Discard worthless events with no location IF ((eq_Nlat == 0.0D0).AND.(eq_Elon == 0.0D0)) CYCLE next_line eq_depth_int = 0 CASE (5) ! .eqc appended_data = ' ' READ (1, 500, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data ! Note: PAD = "YES" permits READ with no appended data 500 FORMAT (9X, & & I5,'.',A2,'.',A2, 1X, & ! using I5 for negative years (B.C.) & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & & A) IF (ios /= 0) EXIT next_line IF (eq_month(1:1) == ' ') eq_month(1:1) = '0' IF (eq_day(1:1) == ' ') eq_day(1:1) = '0' IF (eq_hour(1:1) == ' ') eq_hour(1:1) = '0' IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' IF (eq_second == " ") THEN eq_second = "00" ELSE IF (eq_second(1:1) == ' ') THEN eq_second(1:1) = '0' END IF IF (eq_tenths == ' ') eq_tenths = '0' !try to read FPS from first part of appended data: e1_plunge = 0; e1_azimuth = 0; e2_plunge = 0; e2_azimuth = 0; e3_plunge = 0; e3_azimuth = 0 READ (appended_data, 501, IOSTAT = ios) e1_plunge, e1_azimuth, e2_plunge, e2_azimuth, e3_plunge, e3_azimuth 501 FORMAT (3(1X,I2,1X,I3)) valid_FPS = (e1_plunge /= 0).OR.(e1_azimuth /= 0).OR. & & (e2_plunge /= 0).OR.(e2_azimuth /= 0).OR. & & (e3_plunge /= 0).OR.(e3_azimuth /= 0) IF (valid_FPS.AND.select_tectonic_style) THEN IF ((e1_plunge > e2_plunge).AND.(e1_plunge > e3_plunge)) THEN eq_tectonic_style = 1 ! normal ELSE IF ((e2_plunge > e1_plunge).AND.(e2_plunge > e3_plunge)) THEN eq_tectonic_style = 2 ! strike-slip ELSE IF ((e3_plunge > e1_plunge).AND.(e3_plunge > e2_plunge)) THEN eq_tectonic_style = 3 ! thrust END IF ELSE eq_tectonic_style = 0 ! undefined END IF any_FPS = any_FPS .OR. valid_FPS IF (select_epicenter_elevation) THEN row = 1 + NINT((grd1_lat_max - eq_Nlat) / grd1_d_lat) test_Elon = eq_Elon IF (test_Elon < grd1_lon_min) test_Elon = test_Elon + 360.0D0 col = 1 + NINT((test_Elon - grd1_lon_min) / grd1_d_lon) epicenter_elevation = grid1(row, col) END IF CASE (6) ! Pacheco and Sykes [1992] READ (1, 600, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, & & eq_Nlat_c8, eq_Elon_c8, & ! this method allows for either " -10" or " 57.09" & eq_depth_c8, & ! this method allows for depth of either " 12" or " 7.50" & eq_generic_type_c2, & & eq_moment_over_1E20Nm_c10 600 FORMAT ( 1X, & & I4,A2,A2, 4X, & & A2,A2, & & A, A, & ! this method allows for either " -10" or " 57.09" & A, 25X, & ! this method allows for depth of either " 12" or " 7.50" & A, 1X, & ! eq_generic_type_c2 & A) ! eq_moment_over_1E20Nm_c10 IF (ios /= 0) EXIT next_line IF (eq_hour(1:1) == ' ') eq_hour(1:1) = '0' IF (eq_hour(2:2) == ' ') eq_hour(2:2) = '0' IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' IF (eq_minute(2:2) == ' ') eq_minute(2:2) = '0' eq_second = "00" eq_tenths = '0' READ (eq_Nlat_c8, *) eq_Nlat READ (eq_Elon_c8, *) eq_Elon !test for additional lines concerning the previous earthquake: IF ((eq_year == eq_year_memory ).AND. & &(eq_month == eq_month_memory ).AND. & &(eq_day == eq_day_memory ).AND. & &(eq_hour == eq_hour_memory ).AND. & &(eq_minute == eq_minute_memory).AND. & &(eq_Nlat == eq_Nlat_memory ).AND. & &(eq_Elon == eq_Elon_memory )) CYCLE next_line READ (eq_depth_c8, *, IOSTAT = ios) eq_depth_real IF (ios /= 0) eq_depth_real = 0.0D0 eq_depth_int = NINT(eq_depth_real) READ (eq_moment_over_1E20Nm_c10, *) eq_moment_over_1E20Nm eq_mag = (DLOG10(eq_moment_over_1E20Nm * 1.0D20) - 9.05D0) / 1.5D0 valid_FPS = .FALSE. ! some events have generic mechanism (n, s, t, r, c, ts, rs, ns) ! but this is not the same as a plottable fault plane solution! !Remember this earthquake, in case subsequent lines repeat it (to make room for alternate moments): eq_year_memory = eq_year eq_month_memory = eq_month eq_day_memory = eq_day eq_hour_memory = eq_hour eq_minute_memory = eq_minute eq_Nlat_memory = eq_Nlat eq_Elon_memory = eq_Elon CASE (7) ! Engdahl & Villasensor [2002]; see Appendix 1 for FORMAT READ (1,700,IOSTAT=ios) icat, asol, isol, & & yr, mon, day, & & hr, minu, sec, & & glat, glon, & & dep, & & greg, ntel, & & mag1 700 FORMAT(A6,A1,A5, & & I2,I3,I3, & & 1X,I3,I3,F6.2, & & 1X,F8.3,F8.3, & & F6.1, & & I4,I4, & & F4.1) IF (ios /= 0) EXIT next_line !translate their variables to mine: WRITE (eq_month, "(I2)") mon IF (eq_month(1:1) == ' ') eq_month(1:1) = '0' WRITE (eq_day, "(I2)") day IF (eq_day(1:1) == ' ') eq_day(1:1) = '0' eq_year = 1900 + yr WRITE (eq_hour, "(I2)") hr IF (eq_hour(1:1) == ' ') eq_hour(1:1) = '0' WRITE (eq_minute, "(I2)") minu IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' WRITE (c4, "(F4.1)") sec IF (c4(1:1) == ' ') c4(1:1) = '0' eq_second = c4(1:2) eq_tenths = c4(4:4) eq_Nlat = glat eq_Elon = glon eq_depth_int = NINT(dep) eq_mag = mag1 !CASE (8) is covered together with CASE (2) above. CASE (9) ! ANSS catalog processed by Danijel Schorlemmer of RELM to add independent-mainshock probabilities READ (1, 500, IOSTAT = ios) t1, t2, t3, t4, t5, t6, t7, t8, t9, t10 eq_year = INT(t3) ! I int_temp = NINT(t4) ! month, temporarily in integer form WRITE (eq_month, "(I2)") int_temp ! A2 IF (eq_month(1:1) == ' ') eq_month(1:1) = '0' int_temp = NINT(t5) ! day, temporarily in integer form WRITE (eq_day, "(I2)") int_temp ! A2 IF (eq_day(1:1) == ' ') eq_day(1:1) = '0' int_temp = NINT(t8) ! hour, temporarily in integer form WRITE (eq_hour, "(I2)") int_temp ! A2 int_temp = NINT(t9) ! minute, temporarily in integer form WRITE (eq_minute, "(I2)") int_temp ! A2 IF (eq_minute(1:1) == ' ') eq_minute(1:1) = '0' eq_second = "00" ! A2; not supplied in file eq_tenths = '0' ! A1; not supplied in file eq_Elon = t1 ! F eq_Nlat = t2 ! F eq_depth_int = NINT(t7) ! I eq_mag = t6 ! F WRITE (appended_data, "(F10.8)") t10 ! A (first 10 bytes of possible 200) IF (ios /= 0) EXIT next_line CASE (10) ! GCMT .ndk format, used in 2005+ READ (1,1001,IOSTAT=ios) & ! 1st line of each block of 5 lines. & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths !Note: The origin time in this line was determined by NEIC, ISC, or PDE using body waves; ! the origin time will be revised by lag DT in lines below. ! The latitude, longitude, and depth in this line are also from NEIC, ISC, or PDE, ! and will not be used; the GCMT centroid position will be read in line 3. 1001 FORMAT(5X,I4,1X,A2,1X,A2,1X, & & A2,':',A2,':',A2,'.',A1) IF (ios /= 0) EXIT next_line READ (1, *) ! 2nd line of each block of 5 lines has no important info (for me) READ (1,1002,IOSTAT=ios) & ! 3rd line of each block of 5 lines & DT, eq_Nlat, eq_Elon, eq_depth_real 1002 FORMAT(9X,F9.1,4X,F7.2,5X,F8.2,5X,F6.1) IF (ios /= 0) EXIT next_line !--- add DT seconds to origin time, possibly adjusting minutes & hours! --------- READ (eq_second,"(I2)") eq_seconds_int READ (eq_tenths,"(I1)") eq_tenths_int eq_seconds_real = eq_seconds_int + 0.1 * eq_tenths_int + DT more_minutes_10: DO IF (eq_seconds_real >= 59.95D0) THEN eq_seconds_real = MAX(0.0D0, eq_seconds_real - 60.0D0) READ (eq_minute,"(I2)") eq_minutes eq_minutes = eq_minutes + 1 IF (eq_minutes >= 60) THEN eq_minutes = eq_minutes - 60 READ (eq_hour,"(I2)") eq_hours eq_hours = eq_hours + 1 WRITE (eq_hour,"(I2)") eq_hours !Note: Not allowing date to change (too confusing); in very rare cases hour may reach 24! END IF ! eq_minutes reached 60 WRITE (eq_minute,"(I2)") eq_minutes ELSE EXIT more_minutes_10 END IF ! seconds have exceeded 60, or not END DO more_minutes_10 less_minutes_10: DO IF (eq_seconds_real < 0.0) THEN eq_seconds_real = MIN(59.9D0, eq_seconds_real + 60.0D0) READ (eq_minute,"(I2)") eq_minutes eq_minutes = eq_minutes - 1 IF (eq_minutes < 0) THEN eq_minutes = eq_minutes + 60 READ (eq_hour,"(I2)") eq_hours eq_hours = eq_hours - 1 WRITE (eq_hour,"(I2)") eq_hours !Note: Not allowing date to change (too confusing); in very rare cases hour may reach -1! END IF ! eq_minutes reached 60 WRITE (eq_minute,"(I2)") eq_minutes ELSE EXIT less_minutes_10 END IF ! seconds are negative, or not END DO less_minutes_10 WRITE (c4,"(F4.1)") eq_seconds_real eq_second = c4(1:2) eq_tenths = c4(4:4) !-------------------------------------------------------------------------------- eq_depth_int = NINT(eq_depth_real) READ (1,1003,IOSTAT=ios) M0_exponent ! 4th line in each block of 5 lines 1003 FORMAT(I2) IF (ios /= 0) EXIT next_line READ (1,1004,IOSTAT=ios) & ! 5th and final line of in each block of 5 lines & e3_plunge, e3_azimuth, e2_plunge, e2_azimuth, e1_plunge, e1_azimuth, & & M0_argument ! e3 is extension, e1 is shortening, e2 neutral 1004 FORMAT(3X,3(8X,1X,I2,1X,I3),1X,F7.3) IF (ios /= 0) EXIT next_line any_FPS = .TRUE. M0 = M0_argument * 10.0D0**M0_exponent IF (M0 > 0.0D0) THEN Mw = (DLOG10(M0) - 9.05D0 - 7.0D0) / 1.5D0 !Hanks and Kanamori (1979; J.G.R., 2348-2350) !The -7. term is to convert from dyne-cm !(in the exponent of the GCMT catalog) !to N-m (required in the Hanks and Kanamori formula). eq_mag = Mw ELSE eq_mag = 0.0D0 END IF CASE (11) ! ISC_GEM catalog, in .csv format: READ (1, 1101, IOSTAT = ios) long_line ! probably not more than ~168 bytes... 1101 FORMAT (A) IF (ios /= 0) EXIT next_line !Build a list of byte-indeces of all commas: commas = 0 ! whole INTEGER vector; clearing memory last_byte = LEN_TRIM(long_line) filled = 0 DO i = 1, last_byte IF (long_line(i:i) == ',') THEN filled = filled + 1 IF (filled > 200) THEN ! Weird! Expecting only ~23 at most. WRITE (*, "(' ERROR: More than 200 commas in line of .csv file. Raise limit and recompile!')") CALL Pause() STOP END IF commas(filled) = i END IF END DO !Parse the necessary fields: ! INTEGER::eq_year CHARACTER*2::eq_month,eq_day,eq_hour,eq_minute,eq_second CHARACTER*1::eq_tenths date_time_c24 = long_line(1:commas(1)) ! expecting wasted ' 's at the beginning and the end c5 = date_time_c24(1:5) READ (c5, *) eq_year eq_month = date_time_c24(7:8) ! where leading '0' is already in-place eq_day = date_time_c24(10:11) ! where leading '0' is already in-place eq_hour = date_time_c24(13:14) ! where leading '0' is already in-place eq_minute = date_time_c24(16:17) ! where leading '0' is already in-place eq_second = date_time_c24(19:20) ! where leading '0' is already in-place !Note that fraction-of-second is given to 2 digits; I will truncate to one: eq_tenths = date_time_c24(22:22) !REAL :: eq_Nlat, eq_Elon, eq_depth_real c7 = long_line((commas(1)+1):(commas(2)-1)) READ (c7, *) eq_Nlat ! either an integer, or real with up to 3 decimal places c8 = long_line((commas(2)+1):(commas(3)-1)) READ (c8, *) eq_Elon ! using range -180. to +180. c5 = long_line((commas(7)+1):(commas(8)-1)) ! depth is sometimes presented as integer, and sometimes as real: "437.8". IF (LEN_TRIM(c5) == 0) THEN eq_depth_real = 0.0D0 ELSE READ (c5, *) eq_depth_real END IF eq_depth_int = NINT(eq_depth_real) c5 = long_line((commas(10)+1):(commas(11)-1)) ! m = M_w is usually in "6.54" or "5.6" format, but could be "10.12" someday! READ (c5, *) eq_mag !Note that moment-tensor is given in form of 6 components. These are copied from GCMT in 1976-2009. !Prior to that, most moment-tensors are blank; information begins to appear for selected events in the 1960s. !This program is not currently reading the info, however. any_FPS = .FALSE. ! NOT quite correct; this could be improved someday. END SELECT !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF ((eq_depth_int >= min_depth).AND.(eq_depth_int <= max_depth)) THEN IF (any_longitude.OR.((DEasting(eq_Elon - start_Elon) <= lon_range)) .AND. & &(eq_Nlat >= min_Nlat).AND.(eq_Nlat <= max_Nlat)) THEN IF ((eq_mag >= min_mag).AND.(eq_mag <= max_mag)) THEN IF ((eq_year >= first_year).AND.(eq_year <= last_year)) THEN !FOUND ONE! events = events + 1 SELECT CASE (source) CASE (1) ! ISC: WRITE (2, 1100) file_number, & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag 1100 FORMAT ('ISC(',I4,')', 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) CASE (2, 8) ! CMT .dek format, or Ekstrom & Nettles [1997] extended CMT catalog for 1976: WRITE (2, 1200) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & ! and, for this catalog only: & e1_plunge, e1_azimuth, e2_plunge, e2_azimuth, e3_plunge, e3_azimuth 1200 FORMAT ("Harv. CMT", 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & ! this would be the end for most catalogs; & 3(1X,I2,1X,I3)) ! add 21 more bytes with FPS CASE (3) ! JISSIG WRITE (2, 1300) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag 1300 FORMAT ('JISSIG ', & ! no space after label, & I5,'.',A2,'.',A2, 1X, & ! because date is 5 bytes! & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) CASE (4) ! NEDB from PIDC WRITE (2, 1400) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag 1400 FORMAT ("NucExplDB", 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) CASE (5) ! input was from old_eqc_file selecting = (.NOT.select_tectonic_style).OR.(eq_tectonic_style == desired_tectonic_style) selecting = selecting.AND.((.NOT.select_epicenter_elevation).OR. & & ((epicenter_elevation >= min_epicenter_elevation).AND.(epicenter_elevation <= max_epicenter_elevation))) IF (selecting) THEN WRITE (2, 1500) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & TRIM(appended_data) ! if any (might include FPS, PB2001 segment, etc. 1500 FORMAT ('.eqc file', & ! no space after label, & I5,'.',A2,'.',A2, 1X, & ! because date is 5 bytes! & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & & A) END IF CASE (6) ! Pacheco and Sykes [1992]; append generic type code at end (n, s, t, r, c, ns, ts, rs, u) WRITE (2, 1600) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & eq_generic_type_c2 1600 FORMAT ('PachecoSy' 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, 1X, & & A) CASE (7) ! Engdahl and Villasensor [2002] Centennial catalog WRITE (2, 1700) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag 1700 FORMAT ('Centennia', 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) !CASE (8) is combined with CASE (2) above CASE (9) ! ANSS catalog processed by RELM to add mainshock-independence probabilities WRITE (2, 1900) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & TRIM(appended_data) ! mainshock-independence probability (F10.8) 1900 FORMAT ('ANSS/RELM', & ! no space after label, & I5,'.',A2,'.',A2, 1X, & ! because date is 5 bytes! & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & & 1X,A) CASE (10) ! GCMT .ndk format WRITE (2, 2000) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & ! and, for this catalog only: & e1_plunge, e1_azimuth, e2_plunge, e2_azimuth, e3_plunge, e3_azimuth 2000 FORMAT ("GlobalCMT", 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & ! this would be the end for most catalogs; & 3(1X,I2,1X,I3)) ! add 21 more bytes with FPS CASE (11) ! ISC-GEM: WRITE (2, 2100) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag 2100 FORMAT ("ISC-GEM ", 1X, & & I4,'.',A2,'.',A2, 1X, & & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2) END SELECT END IF ! year is OK END IF ! magnitude is OK END IF ! location is OK END IF ! depth is OK END DO next_line CLOSE (1, IOSTAT = ios) END DO filing ! different files (or only one) CLOSE (2) WRITE (*,"(/' ',I7,' events were found.')") events IF (events > 0) THEN WRITE (*,*) CALL DPrompt_for_Logical('Do you want to make an .ai map plot?',want_map,want_map) IF (want_map) THEN CALL DPrompter (xy_mode = .FALSE., lonlat_mode = .TRUE., path_out = graphics_path) !base map WRITE (*,*) CALL DPrompt_for_Logical('Do you have a basemap (.dig file) to display?',got_dig,got_dig) IF (got_dig) THEN 3000 CALL DPrompt_for_String('Enter name of .dig file:',dig_file,dig_file) dig_pathfile = TRIM(graphics_path) // dig_file OPEN (UNIT = 1, FILE = dig_pathfile, STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found (in this directory).')") GO TO 3000 END IF ! problem opening file CLOSE(1) CALL DSet_Stroke_Color ('foreground') CALL DSet_Line_Style (1.5D0, .FALSE.) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='gray______') CALL DPrompt_for_Logical('Is this file composed ENTIRELY of closed polygons?',polygons, polygons) CALL DPlot_Dig (7, dig_pathfile, polygons, & & 1, in_ok) END IF ! got_dig WRITE (*,*) SELECT CASE (source) CASE (1) ! ISC WRITE (*,"(' In answering the following, remember that ISC reports')") WRITE (*,"(' body-wave magnitude mb, which rarely exceeds 6.8,')") WRITE (*,"(' and that there is a gap in the catalog between ~3 and 0 (unknown).')") CALL DPrompt_for_Real('What is the smallest magnitude to plot (w/ diameter = 2 points)?',min_mag,min_mag) CALL DPrompt_for_Real('What diameter (in points) for magnitude 7.0?',m7_diam_points,m7_diam_points) d1 = MAX((m7_diam_points - 2.0D0), 0.0D0)/MAX((7.0D0 - min_mag), 0.001D0) CASE (2, 8, 10) ! CMT .dek format, or Ekstrom & Nettles [1997] CMT catalog for 1976, or GCMT .ndk format: CALL DPrompt_for_Logical('Do you want to plot fault-plane-solutions as stereographic projections of the lower focal hemisphere?',plot_FPS,plot_FPS) WRITE (*,*) WRITE (*,"(' In answering the following, remember that Harvard reports')") WRITE (*,"(' moment-magnitudes Mw beginning ~5.5 and going up to ~9.6.')") WRITE (*,"(' Diameter of symbol will be a linear function of magnitude.')") CALL DPrompt_for_Real('What is the smallest magnitude to plot (w/ diameter = 2 points)?',min_mag,min_mag) CALL DPrompt_for_Real('What diameter (in points) for magnitude 9.0?',m9_diam_points,m9_diam_points) d1 = MAX((m9_diam_points - 2.0D0), 0.0D0)/MAX((9.0D0 - min_mag), 0.001D0) CASE (3) ! JISSIG WRITE (*,"(' In answering the following, remember that JISSIG reports')") WRITE (*,"(' almost half of the magnitudes as zero. Also, smallest')") WRITE (*,"(' non-zero Mw is 2, but catalog is very incomplete below 7.')") WRITE (*,"(' Diameter of symbol will be a linear function of magnitude.')") CALL DPrompt_for_Real('What is the smallest magnitude to plot (w/ diameter = 2 points)?',min_mag,min_mag) CALL DPrompt_for_Real('What diameter (in points) for Mw magnitude 10?',m10_diam_points,m10_diam_points) d1 = MAX((m10_diam_points - 2.0D0), 0.0D0)/MAX((10.0D0 - min_mag), 0.001D0) CASE (4) ! NEDB WRITE (*,"(' In answering the following, remember that PIDC reports')") WRITE (*,"(' 46% of magnitudes as zero; others are mb 2.0 to 7.4.')") WRITE (*,"(' Diameter of symbol will be a linear function of magnitude.')") CALL DPrompt_for_Real('What is the smallest magnitude to plot (w/ diameter = 2 points)?',min_mag,min_mag) CALL DPrompt_for_Real('What diameter (in points) for magnitude 7?',m7_diam_points,m7_diam_points) d1 = MAX((m7_diam_points - 2.0D0), 0.0D0)/MAX((7.0D0 - min_mag), 0.001D0) CASE (5) ! .eqc file IF (any_FPS) THEN CALL DPrompt_for_Logical('Do you want to plot fault-plane-solutions as stereographic projections of the lower focal hemisphere?',plot_FPS,plot_FPS) END IF WRITE (*,"(' Diameter of symbol will be a linear function of magnitude.')") CALL DPrompt_for_Real('What is the smallest magnitude to plot (w/ diameter = 2 points)?',min_mag,min_mag) CALL DPrompt_for_Real('What diameter (in points) for magnitude 9.0?',m9_diam_points,m9_diam_points) d1 = MAX((m9_diam_points - 2.0D0), 0.0D0)/MAX((9.0D0 - min_mag), 0.001D0) CASE (6) ! Pacheco and Sykes [1992] WRITE (*,"(' In answering the following, remember that Pacheco and Sykes [1992]')") WRITE (*,"(' report only magnitudes from about 7.0 up to about 9.6.')") WRITE (*,"(' Diameter of symbol will be a linear function of magnitude.')") CALL DPrompt_for_Real('What is the smallest magnitude to plot (w/ diameter = 2 points)?',min_mag,min_mag) CALL DPrompt_for_Real('What diameter (in points) for Mw magnitude 10?',m10_diam_points,m10_diam_points) d1 = MAX((m10_diam_points - 2.0D0), 0.0D0)/MAX((10.0D0 - min_mag), 0.001D0) !CASE (8) is combined with CASE (2) above !CASE (10) is combined with CASE (2) above END SELECT d0 = 2.3D0 - d1 * min_mag !Note: extra 0.3 point of radius is to compensate for the ! overlap of the 0.6-point white outline into the interior. OPEN (UNIT = 2, FILE = new_eqc_pathfile, STATUS = 'OLD', & & PAD = 'YES') ! padding required because appended_data (perhaps beginning with a FPS) may not be there CALL DSet_Line_Style (0.6D0, .FALSE.) CALL DBegin_Group rereading: DO READ (2, 4000, IOSTAT = ios) & & eq_year, eq_month, eq_day, & & eq_hour, eq_minute, eq_second, eq_tenths, & & eq_Elon, eq_Nlat, & & eq_depth_int, eq_mag, & & appended_data 4000 FORMAT (9X, & & I5,'.',A2,'.',A2, 1X, & ! read with I5 in case of -3000 (B.C.) & A2,':',A2,':',A2,'.',A1, 1X, & & F8.3, 1X, F7.3, 1X, & & I3, F6.2, & & A ) ! Note: OPEN (..., PAD = "YES") permits READ even if no appended_data IF (ios /= 0) EXIT rereading READ (appended_data, *, IOSTAT = ios) e1_plunge, e1_azimuth, e2_plunge, e2_azimuth, e3_plunge, e3_azimuth IF (ios == 0) THEN valid_FPS = (e1_plunge /= 0).OR.(e1_azimuth /= 0).OR. & & (e2_plunge /= 0).OR.(e2_azimuth /= 0).OR. & & (e3_plunge /= 0).OR.(e3_azimuth /= 0) ELSE valid_FPS = .FALSE. END IF radius_points = 0.5D0 *(d0 + d1 * eq_mag) IF (radius_points >= 1.0D0) THEN ! large enough to plot IF (valid_FPS.AND.plot_FPS.AND.(radius_points >= 3.0D0)) THEN ! plot as FPS ! (1) Plot a small cross to mark position if FPS circle must be pulled aside: CALL DLonLat_2_Uvec (eq_Elon, eq_Nlat, uvec) radians_per_point = (3.527777D-4)*(mp_scale_denominator*DConformal_Deflation(uvec))/mp_radius_meters ! (page m / point)*(local scale)/(planet radius) radius = 3.0D0 * radians_per_point ! each arm of cross CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DNew_L45_Path (5, result_uvec) CALL DTurn_To (azimuth_radians = Pi, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DGreat_to_L45 (result_uvec) CALL DEnd_L45_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DTurn_To (azimuth_radians = Pi_over_2, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DNew_L45_Path (5, result_uvec) CALL DTurn_To (azimuth_radians = -Pi_over_2, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DGreat_to_L45 (result_uvec) CALL DEnd_L45_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL DBegin_Group ! for this one FPS symbol ! (2) Find Northward direction at epicenter, and express as ! an argument (counterclockwise from right, in radians): CALL DLonLat_2_Uvec (eq_Elon, eq_Nlat, uvec) radians_per_point = (3.527777D-4)*(mp_scale_denominator*DConformal_Deflation(uvec))/mp_radius_meters ! (page m / point)*(local scale)/(planet radius) radius = radius_points * radians_per_point CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DProject (uvec = uvec, x = epicenter_x_m, y = epicenter_y_m) CALL DProject (uvec = result_uvec, x = offset_x_m, y = offset_y_m) CALL DMeters_2_Points (epicenter_x_m,epicenter_y_m, epicenter_x_points,epicenter_y_points) CALL DMeters_2_Points (offset_x_m,offset_y_m, offset_x_points,offset_y_points) North_argument_radians = DATan2F((offset_y_points - epicenter_y_points), & &(offset_x_points - epicenter_x_points)) ! (3) Plot a white background circle (even for slide copy!): CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='white_____') CALL DLonLat_2_Uvec (eq_Elon, eq_Nlat, uvec) radians_per_point = (3.527777D-4)*(mp_scale_denominator*DConformal_Deflation(uvec))/mp_radius_meters ! (page m / point)*(local scale)/(planet radius) radius = radius_points * radians_per_point CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DNew_L45_Path (5, result_uvec) CALL DSmall_to_L45 (uvec, result_uvec) ! complete small circle CALL DEnd_L45_Path (close = .TRUE., stroke = .FALSE., fill = .TRUE.) ! (4) Save state of module Map_Projections: CALL DSave_mp_State () ! (5) Reset Map_Projections to show a tiny world at right location and size: ! NOTE: Since projection-plane (x,y) system is arbitrary, I will set it as ! equal to the page-points system (except that it is in meters instead of points): ! centered at lower left corner of page; +x to right; +y up; dimensions ! are those of the physical page space. CALL DSet_Zoom (scale_denominator = 1.0D0, & & x_center_meters = ai_window_xc_points / 2834.65D0, & & y_center_meters = ai_window_yc_points / 2834.65D0, & & xy_wrt_page_radians = 0.0D0) CALL DSet_Stereographic (radius_meters = 0.5D0 * radius_points / 2834.65D0, & ! factor 0.5 counters stereographic blowup of outer circle & projpoint_uvec = (/ -0.01745241D0, 0.0D0, 0.9998477D0 /), & & x_projpoint_meters = epicenter_x_points / 2834.65D0, & & y_projpoint_meters = epicenter_y_points / 2834.65D0, & & y_azimuth_radians = North_argument_radians - Pi_over_2) ! (6) Plot two black sectors on (front) side of little world. ! NOTE: Little world is seen from ~North pole ! (actually, from 89N, 180 E to prevent degeneracy), with ! its Greenwich meridian pointing to N on the big Earth, ! so that if 1.0*plunge is used as a North latitude, and ! -1.0*azimuth is used as a longitude, points plot correctly on ! the lower focal hemisphere. Points with negative ! plunge will not be seen, as they will be on the back side. CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='black_____') e1_lon = -1.0D0 * e1_azimuth e2_lon = -1.0D0 * e2_azimuth e3_lon = -1.0D0 * e3_azimuth e1_lat = 1.0D0 * e1_plunge e2_lat = 1.0D0 * e2_plunge e3_lat = 1.0D0 * e3_plunge CALL DLonLat_2_Uvec (lon = e1_lon, lat = e1_lat, uvec = e1_f_uvec) ! front or visible end CALL DLonLat_2_Uvec (lon = e2_lon, lat = e2_lat, uvec = e2_f_uvec) ! front or visible end CALL DLonLat_2_Uvec (lon = e3_lon, lat = e3_lat, uvec = e3_f_uvec) ! front or visible end !To prevent topological problems during drafting, adjust these three axes !to be exactly perpendicular to each other! Preserve e2_f_uvec exactly, !since this is the one that comes directly from data. CALL DCross (e1_f_uvec, e2_f_uvec, tvec) ! replacing e3, now perp. to e2 IF (tvec(3) < 0.0D0) tvec = -tvec CALL DMake_Uvec (tvec, e3_f_uvec) CALL DCross (e2_f_uvec, e3_f_uvec, tvec) ! replacing e1, now perp. to both IF (tvec(3) < 0.0D0) tvec = -tvec CALL DMake_Uvec (tvec, e1_f_uvec) e1_b_uvec = -e1_f_uvec ! back end of e1 axis; invisible e2_b_uvec = -e2_f_uvec ! back end of e2 axis; invisible e3_b_uvec = -e3_f_uvec ! back end of e3 axis; invisible tvec = e3_f_uvec + e1_b_uvec CALL DMake_uvec (tvec, turn_1_uvec) ! pole of 1st small circle arc tvec = e3_f_uvec + e1_f_uvec CALL DMake_uvec (tvec, turn_2_uvec) ! pole of 2nd small circle arc turn_3_uvec = -turn_1_uvec ! pole of 3rd small circle turn_4_uvec = -turn_2_uvec ! pole of 4th small circle CALL DNew_L45_Path (5, e2_f_uvec) CALL DSmall_To_L45 (pole_uvec = turn_1_uvec, to_uvec = e2_b_uvec) ! front to back CALL DSmall_To_L45 (pole_uvec = turn_2_uvec, to_uvec = e2_f_uvec) ! back to front CALL DEnd_L45_Path (close = .TRUE., stroke = .FALSE., fill = .TRUE., retro = .TRUE.) CALL DNew_L45_Path (5, e2_f_uvec) CALL DSmall_To_L45 (pole_uvec = turn_3_uvec, to_uvec = e2_b_uvec) ! front to back CALL DSmall_To_L45 (pole_uvec = turn_4_uvec, to_uvec = e2_f_uvec) ! back to front CALL DEnd_L45_Path (close = .TRUE., stroke = .FALSE., fill = .TRUE., retro = .TRUE.) ! (7) Reset (saved) state of module Map_Projections CALL DRestore_mp_State () ! (8) Plot the outer circle of lower focal hemisphere CALL DSet_Stroke_Color ('foreground') CALL DLonLat_2_Uvec (eq_Elon, eq_Nlat, uvec) radians_per_point = (3.527777D-4)*(mp_scale_denominator*DConformal_Deflation(uvec))/mp_radius_meters ! (page m / point)*(local scale)/(planet radius) radius = radius_points * radians_per_point CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DNew_L45_Path (5, result_uvec) CALL DSmall_to_L45 (uvec, result_uvec) ! complete small circle CALL DEnd_L45_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) CALL DEnd_Group ! for this one FPS symbol ELSE ! plot as solid dot ! EQs have black fill with white outline (to separate points) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CALL DSet_Stroke_Color ('background') CALL DLonLat_2_Uvec (eq_Elon, eq_Nlat, uvec) radians_per_point = (3.527777D-4)*(mp_scale_denominator*DConformal_Deflation(uvec))/mp_radius_meters ! (page m / point)*(local scale)/(planet radius) radius = radius_points * radians_per_point CALL DTurn_To (azimuth_radians = 0.0D0, base_uvec = uvec, & & far_radians = radius, & ! inputs & omega_uvec = omega_uvec, result_uvec = result_uvec) CALL DNew_L45_Path (5, result_uvec) CALL DSmall_to_L45 (uvec, result_uvec) ! complete small circle CALL DEnd_L45_Path (close = .TRUE., stroke = .TRUE., fill = .TRUE.) END IF ! FPS symbol or solid-dot symbol END IF ! large enough to plot END DO rereading CALL DEnd_Group CLOSE(2) ! parallels and meridians CALL DSet_Line_Style (0.6D0, .TRUE., 2.0D0, 5.0D0) CALL DSet_Stroke_Color ('foreground') WRITE (*,*) CALL DPrompt_for_Integer('How many minutes between parallels and meridians?',minutes,minutes) CALL DGraticule (minutes) !numbered margin: CALL DLonLat_Frame (minutes) !titles IF (ai_toptitles_reserved) THEN WRITE (*,*) CALL DPrompt_for_Logical('Do you want to add 2 title lines?',.TRUE.,want_titles) IF (want_titles) THEN CALL DPrompt_for_String('Enter top line (or a space for none):',' ',top_line) CALL DPrompt_for_String('Enter bottom line (or a space for none):',' ',bottom_line) CALL DTop_Titles (top_line, & & bottom_line) END IF ! want_titles END IF ! ai_top_titles_reserved !sample EQ magnitudes in the margin IF (ai_bottomlegend_reserved .OR. ai_rightlegend_reserved) THEN m1 = NINT(min_mag) SELECT CASE (source) CASE (1) ! ISC m2 = 7 !Note: ISC reports body-wave mb, which saturates at ~6.8. CASE (2, 8, 10) ! CMT, or Ekstrom & Nettles [1997] extended CMT for 1976, or GCMT .ndk format m2 = 9 CASE (3) ! JISSIG m2 = 9 CASE (4) ! NEDB m2 = 7 CASE (5) ! eqc m2 = 9 END SELECT CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') ! EQs have black fill with white outline (to separate points) CALL DSet_Stroke_Color ('background') CALL DSet_Line_Style (0.6D0, .FALSE.) CALL DBegin_Group IF (ai_bottomlegend_reserved) THEN CALL DReport_BottomLegend_Frame (x1_points, x2_points, y1_points, y2_points) x_used_points = 0.0D0 yp = (y1_points + y2_points) / 2.0D0 DO i = m1, m2 radius_points = 0.5D0 * (d0 + d1 * i) xp = x1_points + x_used_points + radius_points + 6.0D0 CALL DCircle_on_L12 (1, xp,yp, radius_points, .FALSE., .TRUE.) ypt = yp - radius_points - 12.0D0 WRITE (c1, "(I1)") i CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & c1) x_used_points = x_used_points + 2.0D0 * radius_points + 6.0D0 END DO IF (any_FPS) THEN ! sample thrust and normal in bottom legend CALL DBegin_Group xp = x1_points + x_used_points + radius_points + 6.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='white_____') CALL DCircle_on_L12 (1, xp,yp, radius_points, .FALSE., .TRUE.) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='black_____') CALL Cats_Eye (xp, yp, radius_points) CALL DSet_Stroke_Color ('foreground') CALL DCircle_on_L12 (1, xp,yp, radius_points, .TRUE., .FALSE.) ypt = yp - radius_points - 12.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & 'thrust') x_used_points = x_used_points + 2.0D0 * radius_points + 6.0D0 CALL DEnd_Group CALL DBegin_Group xp = x1_points + x_used_points + radius_points + 6.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='black_____') CALL DCircle_on_L12 (1, xp,yp, radius_points, .FALSE., .TRUE.) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='white_____') CALL Cats_Eye (xp, yp, radius_points) CALL DSet_Stroke_Color ('foreground') CALL DCircle_on_L12 (1, xp,yp, radius_points, .TRUE., .FALSE.) ypt = yp - radius_points - 12.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CALL DL12_Text (1, xp, ypt, 0.0D0, & & 12, 0.5D0, 0.0D0, & & 'normal') x_used_points = x_used_points + 2.0D0 * radius_points + 6.0D0 CALL DEnd_Group END IF ! sample FPS's needed in bottom legend ELSE ! ai_rightlegend_reserved CALL DReport_RightLegend_Frame (x1_points, x2_points, y1_points, y2_points) y_used_points = 0.0D0 radius_points = 0.5D0 * (d0 + d1 * m2) xp = x1_points + radius_points DO i = m1, m2 radius_points = 0.5D0 * (d0 + d1 * i) yp = y2_points - y_used_points - radius_points - 6.0D0 CALL DCircle_on_L12 (1, xp,yp, radius_points, .FALSE., .TRUE.) xpt = xp + radius_points + 6.0D0 WRITE (c1, "(I1)") i CALL DL12_Text (1, xpt, yp, 0.0D0, & & 12, 0.0D0, 0.4D0, & & c1) y_used_points = y_used_points + 2.0D0 * radius_points + 6.0D0 END DO IF (any_FPS) THEN ! sample thrust and normal in right legend CALL DBegin_Group yp = y2_points - y_used_points - radius_points - 6.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='white_____') CALL DCircle_on_L12 (1, xp,yp, radius_points, .FALSE., .TRUE.) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='black_____') CALL Cats_Eye (xp, yp, radius_points) CALL DSet_Stroke_Color ('foreground') CALL DCircle_on_L12 (1, xp,yp, radius_points, .TRUE., .FALSE.) xpt = xp + radius_points + 6.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CALL DL12_Text (1, xpt, yp, 0.0D0, & & 12, 0.0D0, 0.4D0, & & 'thrust') y_used_points = y_used_points + 2.0D0 * radius_points + 6.0D0 CALL DEnd_Group CALL DBegin_Group yp = y2_points - y_used_points - radius_points - 6.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='black_____') CALL DCircle_on_L12 (1, xp,yp, radius_points, .FALSE., .TRUE.) CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='white_____') CALL Cats_Eye (xp, yp, radius_points) CALL DSet_Stroke_Color ('foreground') CALL DCircle_on_L12 (1, xp,yp, radius_points, .TRUE., .FALSE.) xpt = xp + radius_points + 6.0D0 CALL DSet_Fill_or_Pattern (use_pattern=.FALSE., color_name='foreground') CALL DL12_Text (1, xpt, yp, 0.0D0, & & 12, 0.0D0, 0.4D0, & & 'normal') y_used_points = y_used_points + 2.0D0 * radius_points + 6.0D0 CALL DEnd_Group END IF ! sample FPS's needed in right legend END IF ! bottom, or right, legend reserved CALL DEnd_Group END IF ! either bottom or right legend reserved CALL DEnd_Page END IF ! want_map ELSE ! no events found WRITE (*,"(' No events in this catalog matched your selection criteria.')") WRITE (*,"(' Press [Enter] to end...'\)") READ (*,"(A)") c1 END IF ! events > 0, or zero !Save memories in Seismicity.ini, in current directory: OPEN (UNIT = 11, FILE = 'Seismicity.ini') WRITE (11,"(A)") TRIM(eqc_path) WRITE (11,"(A)") TRIM(graphics_path) WRITE (11,"(I12,' = source')") source WRITE (11,"(A1,11X,' = CD_drive')") CD_drive WRITE (11,"(A)") TRIM(CMT_dek_file) WRITE (11,"(A)") TRIM(CMT1976_file) WRITE (11,"(A)") TRIM(Pacheco_file) WRITE (11,"(A)") TRIM(Engdahl_file) WRITE (11,"(A)") TRIM(PIDC_file) WRITE (11,"(A)") TRIM(old_eqc_file) WRITE (11,"(A)") TRIM(GCMT_ndk_file) WRITE (11,"(A)") TRIM(ISC_GEM_file) WRITE (11,"(11X,L1,' = all_depths')") all_depths WRITE (11,"(I12,' = min_depth')") min_depth WRITE (11,"(I12,' = max_depth')") max_depth WRITE (11,"(11X,L1,' = whole_Earth')") whole_Earth WRITE (11,"(F12.4,' = start_Elon')") start_Elon WRITE (11,"(F12.4,' = lon_range')") lon_range WRITE (11,"(F12.4,' = min_Nlat')") min_Nlat WRITE (11,"(F12.4,' = max_Nlat')") max_Nlat WRITE (11,"(11X,L1,' = all_mags')") all_mags WRITE (11,"(F12.4,' = min_mag')") min_mag WRITE (11,"(F12.4,' = max_mag')") max_mag WRITE (11,"(11X,L1,' = all_years')") all_years WRITE (11,"(I12,' = first_year')") first_year WRITE (11,"(I12,' = last_year')") last_year WRITE (11,"(A)") TRIM(new_eqc_file) WRITE (11,"(11X,L1,' = want_map')") want_map WRITE (11,"(11X,L1,' = got_dig')") got_dig WRITE (11,"(A)") TRIM(dig_file) WRITE (11,"(11X,L1,' = polygons')") polygons WRITE (11,"(11X,L1,' = plot_FPS')") plot_FPS WRITE (11,"(F12.2,' = min_mag')") min_mag WRITE (11,"(F12.2,' = m7_diam_points')") m7_diam_points WRITE (11,"(F12.2,' = m9_diam_points')") m9_diam_points WRITE (11,"(F12.2,' = m10_diam_points')") m10_diam_points WRITE (11,"(I12,' = minutes')") minutes CLOSE (11) CONTAINS SUBROUTINE Cats_Eye (xp, yp, radius_points) !Creates a horizontal lens in the margin; used to avoid !repeating same code 4x in fault-plane-solution explanation IMPLICIT NONE REAL*8, INTENT(IN) :: xp, yp, radius_points REAL*8 :: xp0,xp1,xp2,xp3,yp0,yp1,yp2,yp3 xp0 = xp - radius_points xp1 = xp - 0.4D0 * radius_points ! adjust? xp2 = xp + 0.4D0 * radius_points ! adjust? xp3 = xp + radius_points yp0 = yp yp1 = yp + 0.6D0 * radius_points ! adjust? yp2 = yp1 yp3 = yp CALL DNew_L12_Path (1, xp0, yp0) CALL DCurve_to_L12 (xp1,yp1,xp2,yp2,xp3,yp3) xps = xp1 xp1 = xp2 xp2 = xps xp3 = xp0 yp1 = yp - (yp2 - yp) yp2 = yp1 CALL DCurve_to_L12 (xp1,yp1,xp2,yp2,xp3,yp3) CALL DEnd_L12_Path (close = .TRUE., stroke = .FALSE., fill = .TRUE.) END SUBROUTINE Cats_Eye INTEGER FUNCTION Days_Long (month, year) IMPLICIT NONE INTEGER, INTENT(IN) :: month, year SELECT CASE (month) CASE ( 1) ! January Days_Long = 31 CASE ( 2) ! February IF (MOD(year, 4) == 0) THEN ! leap year Days_Long = 29 ELSE Days_Long = 28 END IF CASE ( 3) ! March Days_Long = 31 CASE ( 4) ! April Days_Long = 30 CASE ( 5) ! May Days_Long = 31 CASE ( 6) ! June Days_Long = 30 CASE ( 7) ! July Days_Long = 31 CASE ( 8) ! August Days_Long = 31 CASE ( 9) ! September Days_Long = 30 CASE (10) ! October Days_Long = 31 CASE (11) ! November Days_Long = 30 CASE (12) ! December Days_Long = 31 END SELECT END FUNCTION Days_Long ! SUBROUTINE File_List( any, & ! & basemap, & ! & catalog, & ! & fegrid, & ! & forces, & ! & gridded_data, & ! & parameters, & ! & velocity, & ! & suggested_file ) ! ! Reports a list (on default device) of filenames of the ! ! (EXACTLY ONE!) type requested. ! ! ! ! Usage of CHARACTER*(*), INTENT(INOUT) :: suggested_file ! ! depends on how many files (of specified type) are ! ! found in the current path_in directory: ! ! * If none are found, suggested_file is unchanged (it may ! ! be a correct file name in some other directory). ! ! * If one file is found, suggested_file is changed to its name. ! ! * If multiple files are found: ! ! -if suggested_file is one of them, it is unchanged. ! ! -if suggested_file is not one, it is changed to ' '. ! ! ! ! Uses GETFILEINFOQQ of module DFLIB.F90 ! ! (DIGITAL Visual Fortran 5.0). ! IMPLICIT NONE ! LOGICAL, INTENT(IN) :: any, basemap, catalog, fegrid, forces, & ! & gridded_data, parameters, velocity ! CHARACTER*(*), INTENT(INOUT) :: suggested_file ! CHARACTER*1 :: first_letter ! CHARACTER*70 :: line = ' ', old_name ! CHARACTER*80 :: string0, string1, string2, & ! & use_path_in ! temporary version, may be changed ! CHARACTER*255 :: files ! INTEGER :: count, full_to, handle, old_result, result ! LOGICAL :: duplicate, matched !! TYPE file$info ! this type as defined in DFLIB.F90 !! INTEGER(4) creation !! INTEGER(4) lastwrite !! INTEGER(4) lastaccess !! INTEGER(4) length !! INTEGER(4) permit !! CHARACTER(255) name !! END TYPE file$info ! TYPE (FILE$INFO) info ! this type as defined in DFLIB.F90 ! use_path_in = TRIM(eqc_path) ! initialized as = global value !10 count = 0 ! matched = .FALSE. ! until we find a file == suggested_file ! IF (any) THEN ! WRITE (*,"(/' Here are all the files in the input directory:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.*' ! ! ELSE IF (basemap) THEN ! WRITE (*,"(/' The following appear to be basemap (.dig) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.DIG' ! ! ELSE IF (catalog) THEN ! WRITE (*,"(/' The following appear to be EarthQuake Catalog (.eqc) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.EQC' ! ! ELSE IF (fegrid) THEN ! WRITE (*,"(/' The following appear to be FE grid (.feg) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.FEG' ! ELSE IF (forces) THEN ! WRITE (*,"(/' The following appear to be nodal force (f___.out) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.OUT' ! must filter below to exclude velocities ! ELSE IF (gridded_data) THEN ! WRITE (*,"(/' The following appear to be gridded data (.grd) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.GRD' ! ELSE IF (parameters) THEN ! WRITE (*,"(/' The following appear to be parameter (.in) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.IN' ! ELSE IF (velocity) THEN ! WRITE (*,"(/' The following appear to be velocity (v___.out) files:')") ! files = TRIM(use_path_in) // & ! defined in FiniteMap above ! & '*.OUT' ! (must filter below to exclude forces) ! ELSE ! RETURN ! no files of any kind are wanted! ! END IF ! full_to = 0 ! keeps track of use of line ! handle = FILE$FIRST ! flag constant, defined in DFLIB as -1 ! old_result = -999 ! old_name = 'undefined' ! all_files: DO ! result = GETFILEINFOQQ (TRIM(files), info, handle) ! !check for duplicate return of last file (a bug in GETFILEINFOQQ): ! IF (result >= 1) THEN ! duplicate = (result == old_result) .AND. (info.name(1:result) == TRIM(old_name)) ! old_name = info.name(1:result) ! ELSE ! duplicate = .FALSE. ! old_name = ' ' ! END IF ! old_result = result ! !- - - - - - - - - - - - - - - - - - - ! IF (handle == FILE$ERROR) RETURN ! defined in DFLIB as -3 ! IF ((result == 0).OR.duplicate) THEN ! no (new) matching files found ! IF (full_to > 0) THEN ! WRITE (*,"(' ',A)") TRIM(line) ! GOTO 100 ! ELSE IF (count == 0) THEN ! WRITE (*,"(' No such files in the input directory.')") ! CALL DPrompt_for_String('What path shall we search (for this file only)?',use_path_in,use_path_in) ! GO TO 10 ! ELSE ! count > 0, but line empty ! GOTO 100 ! END IF ! END IF ! first_letter = info.name(1:1) ! IF (forces.AND.((first_letter == 'V').OR.(first_letter == 'v'))) THEN ! IF (handle == FILE$LAST) THEN ! GOTO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! IF (velocity.AND.((first_letter == 'F').OR.(first_letter == 'f'))) THEN ! IF (handle == FILE$LAST) THEN ! GOTO 100 ! ELSE ! CYCLE all_files ! END IF ! END IF ! !If we've gotten this far, we have a qualified file! ! count = count + 1 ! string0 = TRIM(suggested_file) ! CALL DUpper_Case(string0) ! string1 = info.name(1:result) ! string2 = string1 ! CALL DUpper_Case(string2) ! matched = matched .OR. (string0 == string2) ! IF ((full_to + 2 + result) > 70) THEN ! line would overflow ! WRITE (*,"(' ',A)") TRIM(line) ! full_to = 0 ! line = ' ' ! line = info.name(1:result) ! full_to = result ! ELSE ! line can accept this name ! IF (full_to == 0) THEN ! no leading spaces ! line = info.name(1:result) ! full_to = result ! ELSE ! use 2 leading spaces ! line = TRIM(line) // ' ' // info.name(1:result) ! full_to = full_to + 2 + result ! END IF ! END IF ! IF (handle == FILE$LAST) THEN ! IF (full_to > 0) WRITE (*,"(' ',A)") TRIM(line) ! GOTO 100 ! END IF ! END DO all_files ! 100 IF (count == 1) THEN ! collector point, replacing "RETURN" ! ! so that we can adjust suggested_file(?) ! suggested_file = TRIM(string1) ! ELSE IF (count > 1) THEN ! IF (.NOT.matched) THEN ! suggested_file = ' ' ! END IF ! END IF ! END SUBROUTINE File_List SUBROUTINE Pause () IMPLICIT NONE WRITE (*, "(/' Press [Enter] to continue...')") READ (*, *) RETURN END SUBROUTINE Pause END PROGRAM Seismicity