PROGRAM Slippery ! Processes tabular data about geologic features offset by faulting !(in a tab-delimited ASCII flat-file, one line per datum), and: ! *estimates the probability-density and cumulative-probability functions for offset at the fault trace; ! *convolves with plausible changes in elastic offset to add uncertainty; ! *estimates the probability-density and cumulative-probability functions for age of offset feature; ! *convolves these to estimate probability-density and cumulative-probability functions for the long-term offset rate; ! *accumulates a regional pdf for offset-rate of each type (D, L, N, R, P, T); ! *computes a combined pdf and cpf for long-term offset rate on each fault section; ! *outputs traditional scalar measures of long-term offset rate (L) in an ASCII table: ! ! l:P(L 0. DOUBLE PRECISION, DIMENSION(nbins) :: pd ! probability density (per unit of x) in each bin DOUBLE PRECISION, DIMENSION(nbins) :: xlist ! independent variable value at the center of each bin DOUBLE PRECISION, DIMENSION(0:nbins) :: rxlist ! independent variable value at the right of each bin LOGICAL :: normalized ! only T if integral of this pdf is 1.0000D0 END TYPE pdf TYPE cpf ! cumulative probability function CHARACTER*80 :: name DOUBLE PRECISION :: x_min, x_max ! limits on range of independent variable LOGICAL :: logarithmic ! should never be T unless x_min > 0. DOUBLE PRECISION, DIMENSION(0:nbins) :: P ! probability that value is <= this x. P(0) = 0.0D0 for convenience. DOUBLE PRECISION, DIMENSION(0:nbins) :: xlist ! independent variable value at the right of each bin LOGICAL :: normalized ! only T if P(nbins) = 1.0000D0 END TYPE cpf TYPE (pdf) :: age_pdf, coseismic_C_pdf TYPE (pdf) :: inheritance_time_pdf TYPE (pdf) :: L_D_prior_pdf, L_L_prior_pdf, L_N_prior_pdf, L_P_prior_pdf, L_R_prior_pdf, L_T_prior_pdf, L_X_prior_pdf TYPE (pdf) :: L_D_posterior_pdf, L_L_posterior_pdf, L_N_posterior_pdf, L_P_posterior_pdf, L_R_posterior_pdf, L_T_posterior_pdf TYPE (pdf) :: offset_D_pdf, offset_E_pdf, offset_F_pdf, offset_F_prior_pdf, rate_pdf TYPE (cpf) :: age_cpf, coseismic_C_cpf TYPE (cpf) :: inheritance_time_cpf TYPE (cpf) :: L_D_prior_cpf, L_L_prior_cpf, L_N_prior_cpf, L_P_prior_cpf, L_R_prior_cpf, L_T_prior_cpf, L_X_prior_cpf TYPE (cpf) :: L_D_posterior_cpf, L_L_posterior_cpf, L_N_posterior_cpf, L_P_posterior_cpf, L_R_posterior_cpf, L_T_posterior_cpf TYPE (cpf) :: offset_D_cpf, offset_E_cpf, offset_F_cpf, offset_F_prior_cpf, rate_cpf TYPE (pdf), DIMENSION(storehouse_size) :: rate_pdf_store TYPE (cpf), DIMENSION(storehouse_size) :: rate_cpf_store !================================================================================================== three_halves_plate_rate = 1.5D0 * plate_rate D_dip_radians = D_dip_degrees / 57.29577951D0 L_dip_radians = L_dip_degrees / 57.29577951D0 N_dip_radians = N_dip_degrees / 57.29577951D0 P_dip_radians = P_dip_degrees / 57.29577951D0 R_dip_radians = R_dip_degrees / 57.29577951D0 T_dip_radians = T_dip_degrees / 57.29577951D0 WRITE (*, "( ' PROGRAM Slippery')") WRITE (*, "(/' Processes tabular data about geologic features offset by faulting')") WRITE (*, "( ' (in a tab-delimited ASCII flat-file, one line per datum), and:')") WRITE (*, "(/' *estimates the probability-density and cumulative-probability functions for offset at the fault trace;')") WRITE (*, "( ' *convolves with plausible changes in elastic offset to add uncertainty;')") WRITE (*, "( ' *estimates the probability-density and cumulative-probability functions for age of offset feature; ')") WRITE (*, "( ' *convolves these to estimate probability-density and cumulative-probability functions for the long-term offset rate;')") WRITE (*, "( ' *accumulates a regional pdf for offset-rate of each type (D, L, N, R, P, T);')") WRITE (*, "( ' *computes a combined pdf and cpf for long-term offset rate on each fault section;')") WRITE (*, "( ' *outputs traditional scalar measures of long-term offset rate (L) in an ASCII table:')") WRITE (*, "(/' l:P(L "Inherited" at start, L_R_posterior_pdf%name = "New posterior distribution,Long-term rate,type-R offsets (mm/a)" ! and still not exceeding L_T_posterior_pdf%name = "New posterior distribution,Long-term rate,type-T offsets (mm/a)" ! 80 bytes L_D_posterior_pdf%x_min = zero L_L_posterior_pdf%x_min = zero L_N_posterior_pdf%x_min = zero L_P_posterior_pdf%x_min = zero L_R_posterior_pdf%x_min = zero L_T_posterior_pdf%x_min = zero L_D_posterior_pdf%x_max = plate_rate L_L_posterior_pdf%x_max = plate_rate * 0.5D0 ! kludge; last factor specific to California! L_N_posterior_pdf%x_max = plate_rate * DTAN(N_dip_radians) * 0.1D0 ! kludge; last factor specific to California! L_P_posterior_pdf%x_max = plate_rate * 0.5D0 ! kludge; last factor specific to California! L_R_posterior_pdf%x_max = plate_rate L_T_posterior_pdf%x_max = plate_rate * DTAN(T_dip_radians) * 0.2D0 ! kludge; last factor specific to California! ! N.B. ALSO adjust final factors for L_X_prior_pdf's, above! CALL Fill_pdf_xlist (L_D_posterior_pdf) CALL Fill_pdf_xlist (L_L_posterior_pdf) CALL Fill_pdf_xlist (L_N_posterior_pdf) CALL Fill_pdf_xlist (L_P_posterior_pdf) CALL Fill_pdf_xlist (L_R_posterior_pdf) CALL Fill_pdf_xlist (L_T_posterior_pdf) DO i = 1, nbins L_D_posterior_pdf%pd(i) = zero L_L_posterior_pdf%pd(i) = zero L_N_posterior_pdf%pd(i) = zero L_P_posterior_pdf%pd(i) = zero L_R_posterior_pdf%pd(i) = zero L_T_posterior_pdf%pd(i) = zero END DO !display all 6 prior long-term offset-rate distributions (based on plate_rate): CALL Normalize_pdf (L_D_prior_pdf) CALL Normalize_pdf (L_L_prior_pdf) CALL Normalize_pdf (L_N_prior_pdf) CALL Normalize_pdf (L_P_prior_pdf) CALL Normalize_pdf (L_R_prior_pdf) CALL Normalize_pdf (L_T_prior_pdf) CALL Integrate_pdf_to_cpf (L_D_prior_pdf, L_D_prior_cpf) CALL Integrate_pdf_to_cpf (L_L_prior_pdf, L_L_prior_cpf) CALL Integrate_pdf_to_cpf (L_N_prior_pdf, L_N_prior_cpf) CALL Integrate_pdf_to_cpf (L_P_prior_pdf, L_P_prior_cpf) CALL Integrate_pdf_to_cpf (L_R_prior_pdf, L_R_prior_cpf) CALL Integrate_pdf_to_cpf (L_T_prior_pdf, L_T_prior_cpf) CALL CLEARSCREEN ($GCLEARSCREEN) CALL Plot_3_Distributions(L_D_prior_pdf, L_D_prior_cpf, L_L_prior_pdf, L_L_prior_cpf, L_N_prior_pdf, L_N_prior_cpf, .TRUE.) CALL CLEARSCREEN ($GCLEARSCREEN) CALL Plot_3_Distributions(L_P_prior_pdf, L_P_prior_cpf, L_R_prior_pdf, L_R_prior_cpf, L_T_prior_pdf, L_T_prior_cpf, .TRUE.) 10 CALL CLEARSCREEN ($GCLEARSCREEN) WRITE (*, "(' Select which kind of input file you wish to process:')") WRITE (*, "(' 1 = Peter Bird''s format (e.g., f_Gorda-Cal-Nev.txt)')") WRITE (*, "(' 2 = Sue Perry''s PaleoSites/CombinedEvs format (e.g., f_WGCEP_2007_02.txt)')") WRITE (*, "(' Enter integer code (1 or 2) for type of input file:')") READ (*, *) input_file_type IF ((input_file_type < 1).OR.(input_file_type > 2)) THEN WRITE (*, "(' ERROR: Your integer response was out of the legal range.')") CALL Pause() GO TO 10 END IF IF (input_file_type == 2) THEN WRITE (*, *) WRITE (*, "(' *** IMPORTANT NOTE: ***')") WRITE (*, "(' To prepare the spreadsheet for analysis, you must:')") WRITE (*, "(' (1) Save a copy under a temporary name;')") WRITE (*, "(' (2) Sort rows by FAULT_NAME as primary index; and')") WRITE (*, "(' by MEASURED_COMPONENT_OF_SLIP as secondary index;')") WRITE (*, "(' (3) Delete right-hand columns AR~AY;')") WRITE (*, "(' (4) Save As...Text (tab-delimited).')") CALL Pause() END IF 100 CALL CLEARSCREEN ($GCLEARSCREEN) WRITE (*, "(' Enter [Drive:][\path\]filename for tab-delimited ASCII flat-file data base:')") READ (*, "(A)") input_pathfile OPEN (UNIT = 1, FILE = TRIM(input_pathfile), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not find: ',A)") TRIM(input_pathfile) CALL Pause() GO TO 100 END IF IF (input_file_type == 2) READ (1, *) ! to get past the header line WRITE (*, "(/)") CALL Prompt_for_Logical("Do you want program to pause, for [Enter], so you can inspect offset PDFs?:", & & .TRUE.,pause_to_see_offset) WRITE (*, "(/)") CALL Prompt_for_Logical("Do you want program to pause, for [Enter], so you can inspect datum rate PDFs?:", & & .TRUE.,pause_to_see_rate) WRITE (*, "(/)") CALL Prompt_for_Logical("Do you want program to pause, for [Enter], so you can inspect COMBINED rate PDFs?:", & & .TRUE.,pause_to_see_combined_rate) 120 WRITE (*, "(/////)") WRITE (*, "(' Enter [Drive:][\path\]filename for OUTPUT of revised data base:')") READ (*, "(A)") output_pathfile OPEN (UNIT = 2, FILE = TRIM(output_pathfile), STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not create: ',A)") TRIM(output_pathfile) CALL Pause() GO TO 120 END IF 130 WRITE (*, "(/)") WRITE (*, "(' Enter [Drive:][\path\]filename for OUTPUT of disagreements between rates on same fault:')") READ (*, "(A)") del_pathfile OPEN (UNIT = 3, FILE = TRIM(del_pathfile), STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not create: ',A)") TRIM(del_pathfile) CALL Pause() GO TO 130 END IF !Header line for OUTPUT tab-delimited flat-file ASCII data base; !Note that 1st header differs for the 2 file types: IF (input_file_type == 1) THEN ! Peter Bird's format: line = "F0000" // tab // "Sense" // tab // "Fault_name" // tab // & & "Grade" // tab // "Reference" // tab // & & "Kilometers" // tab // "sigma_(km)" // tab // & & "Since_(Ma)" // tab // "Before_(Ma)" // tab // & & "l:P(L 0) show_string_c40(40:40) = 'M' show_offset_km_c19 = editors_offset_km_c9 show_offset_sigma_km_c9 = editors_offset_sigma_km_c9 ! (currently not displayed) show_began_Ma_c18 = editors_began_Ma_c9 show_ended_Ma_c17 = editors_ended_Ma_c9 show_neorate_mmpa_c16 = editors_neorate_mmpa_c9 show_neorate_sigma_mmpa_c9 = editors_neorate_sigma_mmpa_c9 ! (currently not displayed) show_offset_km_c19 = ADJUSTR(show_offset_km_c19) show_offset_sigma_km_c9 = ADJUSTR(show_offset_sigma_km_c9) show_began_Ma_c18 = ADJUSTR(show_began_Ma_c18) show_ended_Ma_c17 = ADJUSTR(show_ended_Ma_c17) show_neorate_mmpa_c16 = ADJUSTR(show_neorate_mmpa_c16) show_neorate_sigma_mmpa_c9 = ADJUSTR(show_neorate_sigma_mmpa_c9) WRITE (*, "(1X,'SenseV',1X,'Fault',35X,1X,8X,'Offset_(km)',1X,8X,'Began_(Ma)',1X,7X,'Ended_(Ma)',1X,5X,'Rate_(mm/a)' & &/1X, A5, A1,1X, A40, 1X,A19, 1X,A18, 1X,A17, 1X,A16)") & ! 121 bytes after 1X & Findex, offset_type_c1, show_string_c40, & & show_offset_km_c19, & & show_began_Ma_c18, show_ended_Ma_c17, & & show_neorate_mmpa_c16 ! show datum line in raw ASCII form: show_string_c40 = datum_name(1:40) IF (LEN_TRIM(datum_name(41:50)) > 0) show_string_c40(40:40) = 'M' datum_offset_km_c20 = ADJUSTL(datum_offset_km_c20) show_offset_km_c19 = datum_offset_km_c20(1:19) IF (LEN_TRIM(datum_offset_km_c20(20:20)) > 0) show_offset_km_c19(19:19) = 'M' show_offset_km_c19 = ADJUSTR(show_offset_km_c19) datum_offset_sigma_km_c20 = ADJUSTL(datum_offset_sigma_km_c20) show_offset_sigma_km_c9 = datum_offset_sigma_km_c20(1:9) IF (LEN_TRIM(datum_offset_sigma_km_c20(10:20)) > 0) show_offset_sigma_km_c9(9:9) = 'M' show_offset_sigma_km_c9 = ADJUSTR(show_offset_sigma_km_c9) ! (currently not displayed) datum_began_Ma_c20 = ADJUSTL(datum_began_Ma_c20) show_began_Ma_c18 = datum_began_Ma_c20(1:18) IF (LEN_TRIM(datum_began_Ma_c20(19:20)) > 0) show_began_Ma_c18(18:18) = 'M' show_began_Ma_c18 = ADJUSTR(show_began_Ma_c18) datum_ended_Ma_c20 = ADJUSTL(datum_ended_Ma_c20) show_ended_Ma_c17 = datum_ended_Ma_c20(1:17) IF (LEN_TRIM(datum_ended_Ma_c20(18:20)) > 0) show_ended_Ma_c17(17:17) = 'M' show_ended_Ma_c17 = ADJUSTR(show_ended_Ma_c17) datum_neorate_mmpa_c20 = ADJUSTL(datum_neorate_mmpa_c20) show_neorate_mmpa_c16 = datum_neorate_mmpa_c20(1:16) IF (LEN_TRIM(datum_neorate_mmpa_c20(17:20)) > 0) show_neorate_mmpa_c16(16:16) = 'M' show_neorate_mmpa_c16 = ADJUSTR(show_neorate_mmpa_c16) datum_neorate_sigma_mmpa_c20 = ADJUSTL(datum_neorate_sigma_mmpa_c20) show_neorate_sigma_mmpa_c9 = datum_neorate_sigma_mmpa_c20(1:9) IF (LEN_TRIM(datum_neorate_sigma_mmpa_c20(10:20)) > 0) show_neorate_sigma_mmpa_c9(9:9) = 'M' show_neorate_sigma_mmpa_c9 = ADJUSTR(show_neorate_sigma_mmpa_c9) ! (currently not displayed) WRITE (*, "(1X,3X,A1,2X,1X, A40,1X,A19,1X,A18,1X,A17,1X,A16)") & ! 121 bytes after 1X & grade, show_string_c40, & & show_offset_km_c19, & & show_began_Ma_c18, show_ended_Ma_c17, & & show_neorate_mmpa_c16 IF (warning_count /= 0) THEN ! routine computation is impossible DO i = 1, warning_count WRITE (*, "(' *ERROR**> ',A)") TRIM(warning_messages(i)) END DO IF (pause_to_see_rate) THEN WRITE (*, "(/' NO'/' GRAPHS'/' WILL'/' BE'/' SHOWN'/)") CALL Pause() END IF IF (warning_messages(1)(1:3) == "EOF") THEN ! N.B. If this message occurs, it is always in position (1). CALL Finish_up_fault() ! compute disagreements and merge pdfs into best-estimate pdf, etc. EXIT data ! jump out of main loop entirely ELSE CYCLE data ! do nothing with this line, proceed on to the next END IF END IF !------------------------------------------------------------------------------- useful_for_bootstrap = .TRUE. ! (until set .FALSE. by any of 3 criteria below) IF (age_is_upper_limit) useful_for_bootstrap = .FALSE. IF (F_is_lower_limit) useful_for_bootstrap = .FALSE. IF (input_file_type == 2) THEN ! don't allow double-counting of offset features copied to "alt 2" or "alt 3" traces! DO j = 1, (LEN_TRIM(fault_name) - 4) IF ((fault_name(j:(j+4)) == "alt 2").OR.(fault_name(j:(j+4)) == "alt 3")) useful_for_bootstrap = .FALSE. END DO END IF ! age of offset feature, Delta_t: age_pdf%name = "age = Delta t (Ma)" IF (age_is_Gaussian) THEN !Gaussian (truncated and normalized): IF (sigma_Delta_t_stated) THEN !use values of Delta_t_nominal, sigma_Delta_t ELSE IF (age_range_stated) THEN Delta_t_nominal = (age_nominal_high + age_nominal_low) / 2.0D0 sigma_Delta_t = (age_nominal_high - age_nominal_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: sigma_Delta_t = Implied_Sigma(Delta_t_string) END IF age_pdf%x_min = MAX(minimum_age_Ma, Delta_t_nominal - 3.72D0 * sigma_Delta_t) !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. age_pdf%x_max = Delta_t_nominal + 3.72D0 * sigma_Delta_t CALL Fill_pdf_xlist(age_pdf) !----begin equation (21) ----------------------- prefix = 1.0D0 / (sigma_Delta_t * root_2_Pi) DO i = 1, nbins t_prime = age_pdf%xlist(i) age_pdf%pd(i) = prefix * DEXP(-(t_prime - Delta_t_nominal)**2 / (2.0D0 * sigma_Delta_t**2)) END DO !----end equation (21) ------------------------- CALL Normalize_pdf (age_pdf) ! to adjust for the clipping of the tails IF (age_from_raw_14C .OR. age_from_raw_nuclides) THEN IF (age_from_raw_14C) THEN median_inheritance_Ma = 51.0D-6 ! 51 years ELSE IF (age_from_raw_nuclides) THEN median_inheritance_Ma = 0.020 ! 20,000 years END IF omega_inheritance = -DLOG(0.5D0) / median_inheritance_Ma !create (non-negative) distribution for inheritance time, in Ma (equation 23): inheritance_time_pdf%name = "inheritance (Ma)" inheritance_time_pdf%x_min = minimum_age_Ma inheritance_time_pdf%x_max = 9.21D0 / omega_inheritance ! excluding only 0.0001 of the tail on the right CALL Fill_pdf_xlist(inheritance_time_pdf) DO i = 1, nbins inheritance_time_pdf%pd(i) = omega_inheritance * DEXP(-omega_inheritance * inheritance_time_pdf%xlist(i)) END DO inheritance_time_pdf%normalized = .TRUE. CALL Integrate_pdf_to_cpf (inheritance_time_pdf, inheritance_time_cpf) !re-create age_pdf with lower limit, to allow for inheritance: age_pdf%x_min = MAX(minimum_age_Ma, Delta_t_nominal - 3.72D0 * sigma_Delta_t - 9.21D0 / omega_inheritance) age_pdf%x_max = Delta_t_nominal + 3.72D0 * sigma_Delta_t CALL Fill_pdf_xlist(age_pdf) !----begin equation (21) ----------------------- prefix = 1.0D0 / (sigma_Delta_t * root_2_Pi) DO i = 1, nbins t_prime = age_pdf%xlist(i) age_pdf%pd(i) = prefix * DEXP(-(t_prime - Delta_t_nominal)**2 / (2.0D0 * sigma_Delta_t**2)) END DO !----end equation (21) ------------------------- CALL Normalize_pdf (age_pdf) ! to adjust for the clipping of the tails !---- begin equation (23) ---- DO i = 1, nbins ! values of t_prime t_prime = age_pdf%xlist(i) sum = zero ! initialing integral DO j = 1, nbins ! integration steps on t axis: t = inheritance_time_pdf%xlist(j) t_double_prime = t + t_prime lab_pd = Get_One_pd(age_pdf, t_double_prime) d_t = age_pdf%rxlist(j) - age_pdf%rxlist(j-1) sum = sum + inheritance_time_pdf%pd(j) * lab_pd * d_t END DO ! integration steps age_pdf%pd(i) = sum END DO ! values of t_prime CALL Normalize_pdf (age_pdf) !---- end equation (23) ------ END IF ! age_from_raw_14C .OR. age_from_raw_nuclides ELSE IF (age_has_two_bounds) THEN !---- consider lower limit on age: ------------------------------------ IF (sigma_Delta_t_min_stated) THEN !use values of Delta_t_min, sigma_Delta_t_min ELSE IF (Delta_t_min_range_stated) THEN Delta_t_min = (Delta_t_min_high + Delta_t_min_low) / 2.0D0 sigma_Delta_t_min = (Delta_t_min_high - Delta_t_min_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: sigma_Delta_t_min = Implied_Sigma(Delta_t_min_string) END IF !---- consider upper limit on age: ------------------------------------ IF (sigma_Delta_t_max_stated) THEN !use values of Delta_t_max, sigma_Delta_t_max ELSE IF (Delta_t_max_range_stated) THEN Delta_t_max = (Delta_t_max_high + Delta_t_max_low) / 2.0D0 sigma_Delta_t_max = (Delta_t_max_high - Delta_t_max_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: sigma_Delta_t_max = Implied_Sigma(Delta_t_max_string) END IF age_pdf%x_min = MAX(minimum_age_Ma, Delta_t_min - 3.72D0 * sigma_Delta_t_min) !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. age_pdf%x_max = Delta_t_max + 3.72D0 * sigma_Delta_t_max CALL Fill_pdf_xlist(age_pdf) !----- begin equation (25) ------------------------------ DO i = 1, nbins t_prime = age_pdf%xlist(i) age_pdf%pd(i) = (one / (Delta_t_max - Delta_t_min)) * & & 0.5D0 * (one + ERF (REAL((t_prime - Delta_t_min) / (1.414213562D0 * sigma_Delta_t_min)))) * & & 0.5D0 * ERFC(REAL((t_prime - Delta_t_max) / (1.414213562D0 * sigma_Delta_t_max))) END DO !------- end equation (25) ------------------------------ CALL Normalize_pdf(age_pdf) ELSE IF (age_is_lower_limit) THEN !---- consider lower limit on age: ------------------------------------ IF (sigma_Delta_t_min_stated) THEN !use values of Delta_t_min, sigma_Delta_t_min ELSE IF (Delta_t_min_range_stated) THEN Delta_t_min = (Delta_t_min_high + Delta_t_min_low) / 2.0D0 sigma_Delta_t_min = (Delta_t_min_high - Delta_t_min_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: READ (Delta_t_min_string, *) Delta_t_min sigma_Delta_t_min = Implied_Sigma(Delta_t_min_string) END IF !---- fixed upper limit on age: ------------------------------------ Delta_t_max = basement_age_Ma sigma_Delta_t_max = sigma_basement_age_Ma age_pdf%x_min = MAX(minimum_age_Ma, Delta_t_min - 3.72D0 * sigma_Delta_t_min) !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. lambda = -DLOG(0.5D0) / Delta_t_min ! lambda is in units of (Ma)**(-1). age_pdf%x_max = MIN((Delta_t_max + 3.72D0 * sigma_Delta_t_max), & & (Delta_t_min + 9.21D0 / lambda)) CALL Fill_pdf_xlist(age_pdf) !----- begin equation (25) ------------------------------ DO i = 1, nbins t_prime = age_pdf%xlist(i) age_pdf%pd(i) = 0.5D0 * (one + ERF (REAL((t_prime - Delta_t_min) / (1.414213562D0 * sigma_Delta_t_min)))) * & & lambda * MIN(one, (DEXP(-lambda * (t_prime - Delta_t_min)))) * & & 0.5D0 * ERFC(REAL((t_prime - Delta_t_max) / (1.414213562D0 * sigma_Delta_t_max))) END DO !------- end equation (25) ------------------------------ CALL Normalize_pdf(age_pdf) ELSE IF (age_is_upper_limit) THEN !---- consider upper limit on age: ------------------------------------ IF (sigma_Delta_t_max_stated) THEN !use values in Delta_t_max, sigma_Delta_t_max ELSE IF (Delta_t_max_range_stated) THEN Delta_t_max = (Delta_t_max_high + Delta_t_max_low) / 2.0D0 sigma_Delta_t_max = (Delta_t_max_high - Delta_t_max_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: READ (Delta_t_max_string, *) Delta_t_max sigma_Delta_t_max = Implied_Sigma(Delta_t_max_string) END IF age_pdf%x_min = minimum_age_Ma age_pdf%x_max = Delta_t_max + 3.72D0 * sigma_Delta_t_max !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. CALL Fill_pdf_xlist(age_pdf) !----- begin equation (27) ------------------------------ DO i = 1, nbins t_prime = age_pdf%xlist(i) age_pdf%pd(i) = (one / Delta_t_max) * & & 0.5D0 * ERFC(REAL((t_prime - Delta_t_max) / (1.414213562D0 * sigma_Delta_t_max))) END DO !------- end equation (27) ------------------------------ CALL Normalize_pdf(age_pdf) END IF ! age_is_Gaussian, age_has_two_bounds, age_is_lower_limit, OR age_is_upper_limit !transfer age (Delta_t) info to cpf format: CALL Integrate_pdf_to_cpf (age_pdf, age_cpf) !-------------------------------------------------------- ! decide nature of offset (D, L, N, P, R, T): !-------------------------------------------------------- ! Fault-trace offset, F: offset_F_pdf%name = "Fault-trace offset (km)" !different methods for best-estimate, upper-limit, lower-limit, or zero offsets: IF (F_is_Gaussian) THEN !Gaussian (truncated and normalized): Delta_F_nominal, sigma_Delta_F are required. IF (sigma_F_stated) THEN !use values of Delta_F_nominal, sigma_Delta_F ELSE IF (F_range_stated) THEN Delta_F_nominal = (Delta_F_nominal_high + Delta_F_nominal_low) / 2.0D0 sigma_Delta_F = (Delta_F_nominal_high - Delta_F_nominal_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: sigma_Delta_F = Implied_Sigma(Delta_F_string) END IF offset_F_pdf%x_min = MAX(zero, Delta_F_nominal - 3.72D0 * sigma_Delta_F) !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. offset_F_pdf%x_max = Delta_F_nominal + 3.72D0 * sigma_Delta_F CALL Fill_pdf_xlist(offset_F_pdf) offset_F_pdf%normalized = .FALSE. !---- begin equation (8) --------------------- prefix = 1.0D0 / (sigma_Delta_F * root_2_Pi) DO i = 1, nbins f_prime = offset_F_pdf%xlist(i) offset_F_pdf%pd(i) = prefix * DEXP(-(f_prime - Delta_F_nominal)**2 / (2.0D0 * sigma_Delta_F**2)) END DO CALL Normalize_pdf (offset_F_pdf) ! to adjust for the clipping of the tails !---- end equation (8) ---------------------------- ELSE IF (F_is_upper_limit) THEN !erfc, complementary error function (truncated and normalized): Delta_F_max, sigma_Delta_F_max are required: IF (sigma_F_stated) THEN !use values of Delta_F_max, sigma_Delta_F_max ELSE IF (F_range_stated) THEN Delta_F_max = (Delta_F_max_high + Delta_F_max_low) / 2.0D0 sigma_Delta_F_max = (Delta_F_max_high - Delta_F_max_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: sigma_Delta_F_max = Implied_Sigma(Delta_F_string) END IF offset_F_pdf%x_min = zero !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. offset_F_pdf%x_max = Delta_F_max + 3.72D0 * sigma_Delta_F_max CALL Fill_pdf_xlist(offset_F_pdf) offset_F_pdf%normalized = .FALSE. !---- begin equation (9) --------------------- prefix = 1.0D0 / Delta_F_max DO i = 1, nbins f_prime = offset_F_pdf%xlist(i) offset_F_pdf%pd(i) = (1.0D0 / Delta_F_max) * (0.5D0 * ERFC(REAL((f_prime - Delta_F_max) / (1.414213562D0 * sigma_Delta_F_max)))) END DO CALL Normalize_pdf (offset_F_pdf) ! to adjust for the clipping of the tails !---- end equation (9) ---------------------------- ELSE IF (F_is_lower_limit) THEN !erf * prior; Delta_F_min, sigma_Delta_F_min are required: IF (sigma_F_stated) THEN !use values of Delta_F_min, sigma_Delta_F_min ELSE IF (F_range_stated) THEN Delta_F_min = (Delta_F_min_high + Delta_F_min_low) / 2.0D0 sigma_Delta_F_min = (Delta_F_min_high - Delta_F_min_low) / 3.0D0 ELSE ! use implied standard deviation based on decimal position of right-most non-zero digit: sigma_Delta_F_min = Implied_Sigma(Delta_F_string) END IF !First, create offset_F_prior_pdf: !find appropriate prior, and call it "X": IF (offset_type_c1 == 'D') THEN L_X_prior_pdf = L_D_prior_pdf ELSE IF (offset_type_c1 == 'L') THEN L_X_prior_pdf = L_L_prior_pdf ELSE IF (offset_type_c1 == 'N') THEN L_X_prior_pdf = L_N_prior_pdf ELSE IF (offset_type_c1 == 'P') THEN L_X_prior_pdf = L_P_prior_pdf ELSE IF (offset_type_c1 == 'R') THEN L_X_prior_pdf = L_R_prior_pdf ELSE IF (offset_type_c1 == 'T') THEN L_X_prior_pdf = L_T_prior_pdf END IF offset_F_prior_pdf%x_min = MAX(zero, Delta_F_min - 3.72D0 * sigma_Delta_F_min) !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. offset_F_prior_pdf%x_max = L_X_prior_pdf%x_max * age_pdf%x_max CALL Fill_pdf_xlist(offset_F_prior_pdf) offset_F_prior_pdf%normalized = .FALSE. !---- begin equation (11) ---- DO i = 1, nbins ! values of l_prime f_prime = offset_F_prior_pdf%xlist(i) sum = zero ! initialing integral DO j = 1, nbins ! integration steps on l_prime axis: l_prime = L_X_prior_pdf%xlist(j) IF (l_prime > 0.0D0) THEN t_prime = f_prime / l_prime ! (most of the time) t_pd = Get_One_pd(age_pdf, t_prime) d_l_prime = L_X_prior_pdf%rxlist(j) - L_X_prior_pdf%rxlist(j-1) sum = sum + (L_X_prior_pdf%pd(j) * t_pd / l_prime) * d_l_prime END IF END DO ! integration steps offset_F_prior_pdf%pd(i) = sum END DO ! values of f_prime CALL Normalize_pdf (offset_F_prior_pdf) !---- end equation (11) ------ !Second, use this offset_F_prior_pdf to compute offset_F_pdf: offset_F_pdf%x_min = MAX (zero, Delta_F_min - 3.72D0 * sigma_Delta_F_min) !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. offset_F_pdf%x_max = offset_F_prior_pdf%x_max CALL Fill_pdf_xlist(offset_F_pdf) offset_F_pdf%normalized = .FALSE. !---- begin equation (12) --------------------- DO i = 1, nbins f_prime = offset_F_pdf%xlist(i) offset_F_pdf%pd(i) = Get_One_pd (offset_F_prior_pdf, f_prime) * & & 0.5D0 * (1.0D0 + ERF(REAL((f_prime - Delta_F_min) / (1.414213562D0 * sigma_Delta_F_min)))) END DO CALL Normalize_pdf (offset_F_pdf) ! to adjust for the clipping of the tails !---- end equation (12) ---------------------------- ELSE IF (F_is_zero) THEN IF (sigma_F_stated) THEN ! use truncated, normalized Gaussian distribution: offset_F_pdf%x_min = zero offset_F_pdf%x_max = 3.72D0 * sigma_Delta_F !N.B. Range of +- 3.72 sigmas covers range from P = 0.0001 to P = 0.9999 in a Gaussian distribution. CALL Fill_pdf_xlist(offset_F_pdf) offset_F_pdf%normalized = .FALSE. !---- begin equation (8, modified by Delta_F_nominal = zero) --------------------- prefix = 1.0D0 / (sigma_Delta_F * root_2_Pi) DO i = 1, nbins f_prime = offset_F_pdf%xlist(i) offset_F_pdf%pd(i) = prefix * DEXP(-f_prime**2 / (2.0D0 * sigma_Delta_F**2)) END DO CALL Normalize_pdf (offset_F_pdf) ! to adjust for the clipping of the tails !---- end equation (8, modified by Delta_F_nominal = zero) ---------------------------- ELSE IF (F_range_stated) THEN ! use upper limit in a uniform ("boxcar" distribution): offset_F_pdf%x_min = zero offset_F_pdf%x_max = Delta_F_nominal_high CALL Fill_pdf_xlist(offset_F_pdf) DO i = 1, nbins offset_F_pdf%pd(i) = 1.0D0 / Delta_F_nominal_high END DO offset_F_pdf%normalized = .TRUE. ELSE ! pure, unadorned "0" in table; infer a minimal level of geologic imprecision (epsilon): offset_F_pdf%x_min = zero offset_F_pdf%x_max = epsilon ! on the order of 1.D-5 km = 0.01 m = 1 cm, as in 100 cracks of 100 micron offset. CALL Fill_pdf_xlist(offset_F_pdf) DO i = 1, nbins offset_F_pdf%pd(i) = 1.0D0 / epsilon END DO offset_F_pdf%normalized = .TRUE. END IF ! sigma_F_stated, or F_range_stated, or neither END IF ! F_is_Gaussian, OR F_is_upper_limit, OR F_is_lower_limit, or F_is_zero !transfer F info to cpf format: offset_F_cpf%name = offset_F_pdf%name CALL Integrate_pdf_to_cpf (offset_F_pdf, offset_F_cpf) !-------------------------------------------------------- !cryptic elastic offset, offset_E_pdf & offset_E_cpf: offset_E_cpf%name = "Elastic offset (km)" !determine medians of maximum slip & maximum offset in historic ground-breaking continental earthquakes, !following Figure 1 of my paper, which is based on data of Wells & Coppersmith [1994]: IF (offset_type_c1 == 'D') THEN median_coseismic_slip = 1.15D0 / 1000.D0 ! (converting m to km) median_coseismic_offset = median_coseismic_slip * DCOS(D_dip_radians) ! (converting slip to extension/opening) ELSE IF (offset_type_c1 == 'L') THEN median_coseismic_slip = 2.0D0 / 1000.D0 ! (converting m to km) median_coseismic_offset = median_coseismic_slip ! (sinistral strike-slip = slip) ELSE IF (offset_type_c1 == 'N') THEN median_coseismic_slip = 1.15D0 / 1000.D0 ! (converting m to km) median_coseismic_offset = median_coseismic_slip * DSIN(N_dip_radians) ! (converting slip to throw) ELSE IF (offset_type_c1 == 'P') THEN median_coseismic_slip = 1.3D0 / 1000.D0 ! (converting m to km) median_coseismic_offset = median_coseismic_slip * DCOS(P_dip_radians) ! (converting slip to convergence) ELSE IF (offset_type_c1 == 'R') THEN median_coseismic_slip = 2.0D0 / 1000.D0 ! (converting m to km) median_coseismic_offset = median_coseismic_slip ! (dextral strike-slip = slip) ELSE IF (offset_type_c1 == 'T') THEN median_coseismic_slip = 1.3D0 / 1000.D0 ! (converting m to km) median_coseismic_offset = median_coseismic_slip * DSIN(T_dip_radians) ! (converting slip to throw) END IF lambda = -DLOG(0.5D0) / median_coseismic_offset ! lambda is now in units of (km)**(-1). !create (non-negative) distributions for coseismic offset C (equation 15ab): coseismic_C_pdf%name = "Coseismic offset (km)" coseismic_C_cpf%name = coseismic_C_pdf%name coseismic_C_pdf%x_min = zero coseismic_C_pdf%x_max = 9.21D0 / lambda ! excluding only 0.0001 of the tail on the right CALL Fill_pdf_xlist(coseismic_C_pdf) coseismic_C_cpf%x_min = zero coseismic_C_cpf%x_max = 9.21D0 / lambda ! excluding only 0.0001 of the tail on the right CALL Fill_cpf_xlist(coseismic_C_cpf) coseismic_C_cpf%P(0) = zero DO i = 1, nbins coseismic_C_pdf%pd(i) = lambda * DEXP(-lambda * coseismic_C_pdf%xlist(i)) coseismic_C_cpf%P(i) = 1.0D0 - DEXP(-lambda * coseismic_C_cpf%xlist(i)) END DO coseismic_C_pdf%normalized = .TRUE. coseismic_C_cpf%normalized = .TRUE. !finished creating coseismic distributions (equation 15ab). !Next, create offset_F_prior_pdf & _cpf based on L_X_prior_pdf and on age_pdf (equation 11): IF (offset_type_c1 == 'D') THEN L_X_prior_pdf = L_D_prior_pdf ELSE IF (offset_type_c1 == 'L') THEN L_X_prior_pdf = L_L_prior_pdf ELSE IF (offset_type_c1 == 'N') THEN L_X_prior_pdf = L_N_prior_pdf ELSE IF (offset_type_c1 == 'P') THEN L_X_prior_pdf = L_P_prior_pdf ELSE IF (offset_type_c1 == 'R') THEN L_X_prior_pdf = L_R_prior_pdf ELSE IF (offset_type_c1 == 'T') THEN L_X_prior_pdf = L_T_prior_pdf END IF offset_F_prior_pdf%x_min = zero offset_F_prior_pdf%x_max = L_X_prior_pdf%x_max * age_pdf%x_max CALL Fill_pdf_xlist(offset_F_prior_pdf) offset_F_prior_pdf%normalized = .FALSE. !---- begin equation (11) ---- DO i = 1, nbins ! values of l_prime f_prime = offset_F_prior_pdf%xlist(i) sum = zero ! initialing integral DO j = 1, nbins ! integration steps on l_prime axis: l_prime = L_X_prior_pdf%xlist(j) IF (l_prime > 0.0D0) THEN t_prime = f_prime / l_prime ! (most of the time) t_pd = Get_One_pd(age_pdf, t_prime) d_l_prime = L_X_prior_pdf%rxlist(j) - L_X_prior_pdf%rxlist(j-1) sum = sum + (L_X_prior_pdf%pd(j) * t_pd / l_prime) * d_l_prime END IF END DO ! integration steps offset_F_prior_pdf%pd(i) = sum END DO ! values of f_prime CALL Normalize_pdf (offset_F_prior_pdf) !---- end equation (11) ------ CALL Integrate_pdf_to_cpf (offset_F_prior_pdf, offset_F_prior_cpf) !finished creating offset_F_prior_pdf & _cpf based on L_X_prior_pdf and on age_pdf. !Now, product of complements gives complement of desired cpf (equation 17) for all positive Delta E: offset_E_cpf%x_min = -coseismic_C_cpf%x_max offset_E_cpf%x_max = coseismic_C_cpf%x_max CALL Fill_cpf_xlist(offset_E_cpf) !-----begin equation (17) ------------------ DO i = 0, nbins e = offset_E_cpf%xlist(i) IF (e <= 0.0D0) THEN offset_E_cpf%P(i) = zero ! for now.... ELSE ! e > zero offset_E_cpf%P(i) = -((one - Get_One_cp(offset_F_prior_cpf, e)) * (one - Get_One_cp(coseismic_C_cpf, e)) - one) !N.B. These are only temporary "z_cpf" values which will be differentiated. Afterwards, they will be overwritten. END IF END DO !----- end equation (17) --------------------------------- !----- begin equation (18), determining z_pdf: CALL Differentiate_cpf_to_pdf (offset_E_cpf, offset_E_pdf) !----- end equation (18) --------------------------------- !finished filling in POSITIVE sides of offset_E_pfd & _cpf; !Now fill in negative side of offset_E_pdf (without worrying about balancing, normalizing, or _cpf): left_side: DO i = 1, nbins e = offset_E_pdf%xlist(i) IF (e <= 0.0D0) THEN offset_E_pdf%pd(i) = Get_One_pd (coseismic_C_pdf, -e) ELSE ! e > zero; end of loop EXIT left_side END IF END DO left_side !determine and apply factor "k" (equation 19) so that mean or expection of offset_E_pdf is zero: CALL No_Expectations_for_pdf(offset_E_pdf) ! (equation 19) CALL Normalize_pdf (offset_E_pdf) ! (equation 19) CALL Integrate_pdf_to_cpf (offset_E_pdf, offset_E_cpf) ! overwrites previous "z_cpf" values on positive side !-------------------------------------------------------- !Distant (far-field) total offset, Delta D offset_D_pdf%name = "Distant (far-field) offset (km)" offset_D_pdf%x_min = MAX(zero, offset_F_pdf%x_min + offset_E_pdf%x_min) offset_D_pdf%x_max = offset_F_pdf%x_max + offset_E_pdf%x_max CALL Fill_pdf_xlist(offset_D_pdf) offset_D_pdf%normalized = .FALSE. !---- begin equation (20) ---- DO i = 1, nbins ! values of d_prime d_prime = offset_D_pdf%xlist(i) sum = zero ! initialing integral !Use integration method based on (6) or on (7): whichever is least stiff and most accurate: IF (offset_F_pdf%x_max > offset_E_pdf%x_max) THEN ! normal case; large offset relative to one EQ DO j = 1, nbins ! integration steps e_prime = offset_E_pdf%xlist(j) f_prime = d_prime - e_prime F_pd = Get_One_pd(offset_F_pdf, f_prime) d_e_prime = offset_E_pdf%rxlist(j) - offset_E_pdf%rxlist(j-1) sum = sum + offset_E_pdf%pd(j) * F_pd * d_e_prime END DO ! integration steps ELSE ! special case (e.g. F_is_zero = .TRUE.) of very small Delta F relative to Delta E: DO j = 1, nbins ! integration steps f_prime = offset_F_pdf%xlist(j) e_prime = d_prime - f_prime E_pd = Get_One_pd(offset_E_pdf, e_prime) d_f_prime = offset_F_pdf%rxlist(j) - offset_F_pdf%rxlist(j-1) sum = sum + offset_F_pdf%pd(j) * E_pd * d_f_prime END DO ! integration steps END IF offset_D_pdf%pd(i) = sum END DO ! values of d_prime CALL Normalize_pdf (offset_D_pdf) !---- end equation (20) ------ CALL Integrate_pdf_to_cpf (offset_D_pdf, offset_D_cpf) !-------------------------------------------------------- CALL Plot_3_Distributions(offset_F_pdf, offset_F_cpf, offset_E_pdf, offset_E_cpf, offset_D_pdf, offset_D_cpf, pause_to_see_offset) !-------------------------------------------------------- CALL CLEARSCREEN ($GCLEARSCREEN) show_string_c40 = fault_name(1:40) IF (LEN_TRIM(fault_name(41:50)) > 0) show_string_c40(40:40) = 'M' show_offset_km_c19 = editors_offset_km_c9 show_offset_sigma_km_c9 = editors_offset_sigma_km_c9 ! (currently not displayed) show_began_Ma_c18 = editors_began_Ma_c9 show_ended_Ma_c17 = editors_ended_Ma_c9 show_neorate_mmpa_c16 = editors_neorate_mmpa_c9 show_neorate_sigma_mmpa_c9 = editors_neorate_sigma_mmpa_c9 ! (currently not displayed) show_offset_km_c19 = ADJUSTR(show_offset_km_c19) show_offset_sigma_km_c9 = ADJUSTR(show_offset_sigma_km_c9) show_began_Ma_c18 = ADJUSTR(show_began_Ma_c18) show_ended_Ma_c17 = ADJUSTR(show_ended_Ma_c17) show_neorate_mmpa_c16 = ADJUSTR(show_neorate_mmpa_c16) show_neorate_sigma_mmpa_c9 = ADJUSTR(show_neorate_sigma_mmpa_c9) WRITE (*, "(1X,'SenseV',1X,'Fault',35X,1X,8X,'Offset_(km)',1X,8X,'Began_(Ma)',1X,7X,'Ended_(Ma)',1X,5X,'Rate_(mm/a)' & &/1X, A5, A1,1X, A40, 1X,A19, 1X,A18, 1X,A17, 1X,A16)") & ! 121 bytes after 1X & Findex, offset_type_c1, show_string_c40, & & show_offset_km_c19, & & show_began_Ma_c18, show_ended_Ma_c17, & & show_neorate_mmpa_c16 ! show datum line in raw ASCII form: show_string_c40 = datum_name(1:40) IF (LEN_TRIM(datum_name(41:50)) > 0) show_string_c40(40:40) = 'M' datum_offset_km_c20 = ADJUSTL(datum_offset_km_c20) show_offset_km_c19 = datum_offset_km_c20(1:19) IF (LEN_TRIM(datum_offset_km_c20(20:20)) > 0) show_offset_km_c19(19:19) = 'M' show_offset_km_c19 = ADJUSTR(show_offset_km_c19) datum_offset_sigma_km_c20 = ADJUSTL(datum_offset_sigma_km_c20) show_offset_sigma_km_c9 = datum_offset_sigma_km_c20(1:9) IF (LEN_TRIM(datum_offset_sigma_km_c20(10:20)) > 0) show_offset_sigma_km_c9(9:9) = 'M' show_offset_sigma_km_c9 = ADJUSTR(show_offset_sigma_km_c9) ! (currently not displayed) datum_began_Ma_c20 = ADJUSTL(datum_began_Ma_c20) show_began_Ma_c18 = datum_began_Ma_c20(1:18) IF (LEN_TRIM(datum_began_Ma_c20(19:20)) > 0) show_began_Ma_c18(18:18) = 'M' show_began_Ma_c18 = ADJUSTR(show_began_Ma_c18) datum_ended_Ma_c20 = ADJUSTL(datum_ended_Ma_c20) show_ended_Ma_c17 = datum_ended_Ma_c20(1:17) IF (LEN_TRIM(datum_ended_Ma_c20(18:20)) > 0) show_ended_Ma_c17(17:17) = 'M' show_ended_Ma_c17 = ADJUSTR(show_ended_Ma_c17) datum_neorate_mmpa_c20 = ADJUSTL(datum_neorate_mmpa_c20) show_neorate_mmpa_c16 = datum_neorate_mmpa_c20(1:16) IF (LEN_TRIM(datum_neorate_mmpa_c20(17:20)) > 0) show_neorate_mmpa_c16(16:16) = 'M' show_neorate_mmpa_c16 = ADJUSTR(show_neorate_mmpa_c16) datum_neorate_sigma_mmpa_c20 = ADJUSTL(datum_neorate_sigma_mmpa_c20) show_neorate_sigma_mmpa_c9 = datum_neorate_sigma_mmpa_c20(1:9) IF (LEN_TRIM(datum_neorate_sigma_mmpa_c20(10:20)) > 0) show_neorate_sigma_mmpa_c9(9:9) = 'M' show_neorate_sigma_mmpa_c9 = ADJUSTR(show_neorate_sigma_mmpa_c9) ! (currently not displayed) WRITE (*, "(1X,3X,A1,2X,1X, A40,1X,A19,1X,A18,1X,A17,1X,A16)") & ! 121 bytes after 1X & grade, show_string_c40, & & show_offset_km_c19, & & show_began_Ma_c18, show_ended_Ma_c17, & & show_neorate_mmpa_c16 IF (warning_count /= 0) THEN ! routine computation is impossible DO i = 1, warning_count WRITE (*, "(' *ERROR**> ',A)") TRIM(warning_messages(i)) END DO CALL Silent_Pause() IF (warning_messages(1)(1:3) == "EOF") THEN ! N.B. When this message occurs, it is always in position (1). EXIT data ! jump out of main loop entirely ELSE CYCLE data ! do nothing with this line, proceed on to the next END IF END IF !-------------------------------------------------------- ! Long-term offset rate, L: rate_pdf%name = "Long-term offset rate (km/Ma, or mm/a)" rate_pdf%x_min = offset_D_pdf%x_min / age_pdf%x_max IF (age_pdf%x_min > zero) THEN rate_pdf%x_max = MIN((offset_D_pdf%x_max / age_pdf%x_min), three_halves_plate_rate) ELSE rate_pdf%x_max = three_halves_plate_rate END IF IF (rate_pdf%x_max > (0.5D0 * three_halves_plate_rate)) useful_for_bootstrap = .FALSE. CALL Fill_pdf_xlist(rate_pdf) !---- begin equation (29) ---- DO i = 1, nbins ! values of l_prime l_prime = rate_pdf%xlist(i) sum = zero ! initialing integral DO j = 1, nbins ! integration steps t_prime = age_pdf%xlist(j) d_prime = l_prime * t_prime D_pd = Get_One_pd(offset_D_pdf, d_prime) d_t_prime = age_pdf%rxlist(j) - age_pdf%rxlist(j-1) sum = sum + age_pdf%pd(j) * D_pd * t_prime * d_t_prime END DO ! integration steps rate_pdf%pd(i) = sum END DO ! values of l_prime CALL Normalize_pdf (rate_pdf) !---- end equation (29) ------ !transfer L info to cpf format: CALL Integrate_pdf_to_cpf (rate_pdf, rate_cpf) !-------------------------------------------------------- low_limit_c9 = ADJUSTL(ASCII9(REAL(Get_x_from_P(rate_cpf, 0.025D0)))) mode_c9 = ADJUSTL(ASCII9(REAL(Get_Mode(rate_pdf)))) ! mode median_c9 = ADJUSTL(ASCII9(REAL(Get_Median(rate_cpf)))) ! median mean_c9 = ADJUSTL(ASCII9(REAL(Get_Mean(rate_pdf)))) ! mean high_limit_c9 = ADJUSTL(ASCII9(REAL(Get_x_from_P(rate_cpf, 0.975D0)))) sigma_c9 = ADJUSTL(ASCII9(REAL(Get_Sigma(rate_pdf)))) ! standard deviation !Do not report "Unbounded" or "Unreasonable" numbers with digits; use code 'U' as described in paper: IF (Get_x_from_P(rate_cpf, 0.025D0) > (1.5D0 * plate_rate)) low_limit_c9 = 'U' IF (Get_Mode(rate_pdf) > (1.5D0 * plate_rate)) mode_c9 = 'U' IF (Get_Median(rate_cpf) > (1.5D0 * plate_rate)) median_c9 = 'U' IF (Get_Mean(rate_pdf) > (1.5D0 * plate_rate)) mean_c9 = 'U' IF (Get_x_from_P(rate_cpf, 0.975D0) > (1.5D0 * plate_rate)) high_limit_c9 = 'U' line = Findex // tab // offset_type_c1 // tab // TRIM(fault_name) // tab // & & TRIM(grade) // tab // TRIM(datum_name)// tab // & & TRIM(datum_offset_km_c20) // tab // TRIM(datum_offset_sigma_km_c20) // tab // & & TRIM(datum_began_Ma_c20) // tab // TRIM(datum_ended_Ma_c20) // tab // & & TRIM(low_limit_c9) // tab // TRIM(mode_c9) // tab // TRIM(median_c9) // tab // TRIM(mean_c9) // tab // TRIM(high_limit_c9) // tab // & & TRIM(sigma_c9) WRITE (2, "(A)") TRIM(line) !-------------------------------------------------------- !remember the PDF and CPF for this L datum, until all data for this fault segment have been completed: number_in_store = number_in_store + 1 IF (number_in_store > storehouse_size) THEN WRITE (*,"(' WARNING: Storage limit reached. Increase INTEGER PARAMETER storehouse_size.')") CALL Pause() number_in_store = storehouse_size ELSE ! normal case rate_pdf_store(number_in_store) = rate_pdf rate_cpf_store(number_in_store) = rate_cpf IF (age_is_Gaussian) THEN Delta_t_store(number_in_store) = Delta_t_nominal ELSE IF (age_has_two_bounds) THEN Delta_t_store(number_in_store) = (Delta_t_min + Delta_t_max) / 2.0D0 ELSE IF (age_is_lower_limit) THEN Delta_t_store(number_in_store) = Delta_t_min ELSE IF (age_is_upper_limit) THEN Delta_t_store(number_in_store) = Delta_t_max END IF !But, Delta_t used in computation of alpha (for combined rate computation) by be less: IF (F_is_zero) THEN Delta_t_for_alpha_store(number_in_store) = minimum_age_Ma ELSE Delta_t_for_alpha_store(number_in_store) = Delta_t_store(number_in_store) END IF !because an unbroken overlap formation is a valid datum since 5 Ma (or less) !just as much as it is a valid datum since (e.g.) 35 Ma. grade_store(number_in_store) = grade datum_name_store(number_in_store) = datum_name !set logical values which will contribute to decision of asterisk_on_data_count? positive_offset_query(number_in_store) = .NOT.F_is_zero END IF !-------------------------------------------------------- CALL Plot_3_Distributions(offset_D_pdf, offset_D_cpf, age_pdf, age_cpf, rate_pdf, rate_cpf, pause_to_see_rate) !-------------------------------------------------------- IF (useful_for_bootstrap) THEN IF (offset_type_c1 == 'D') THEN CALL Add_pdf_to_accumulator (rate_pdf, one, L_D_posterior_pdf) D_count = D_count + 1 ELSE IF (offset_type_c1 == 'L') THEN CALL Add_pdf_to_accumulator (rate_pdf, one, L_L_posterior_pdf) L_count = L_count + 1 ELSE IF (offset_type_c1 == 'N') THEN CALL Add_pdf_to_accumulator (rate_pdf, one, L_N_posterior_pdf) N_count = N_count + 1 ELSE IF (offset_type_c1 == 'P') THEN CALL Add_pdf_to_accumulator (rate_pdf, one, L_P_posterior_pdf) P_count = P_count + 1 ELSE IF (offset_type_c1 == 'R') THEN CALL Add_pdf_to_accumulator (rate_pdf, one, L_R_posterior_pdf) R_count = R_count + 1 ELSE IF (offset_type_c1 == 'T') THEN CALL Add_pdf_to_accumulator (rate_pdf, one, L_T_posterior_pdf) T_count = T_count + 1 END IF END IF ! useful_for_bootstrap END DO data !=============== END LOOP ON LINES OF INPUT FLAT-FILE... =================================================== CLOSE (UNIT = 1, DISP = "KEEP") ! input_pathfile = tab-delimited ASCII data base CLOSE (UNIT = 2, DISP = "KEEP") ! output_pathfile = tab-delimited ASCII data base, with new MEDIAN offset-rates and their SIGMAs. CLOSE (UNIT = 3, DISP = "KEEP") ! del_pathfile = list of disagreements between long-term offset rates on the same fault segment. CALL CLEARSCREEN ($GCLEARSCREEN) WRITE (*, "(' Job completed.')") WRITE (*, "(' Press [Enter] to see the 6 posterior rate pdfs/cpfs, and to save them to a file....')") CALL Silent_Pause() !display all 6 POSTERIOR long-term offset-rate distributions (and indicate "new" or "inherited"): IF (D_count > 0) THEN WRITE (c6, "(I6)") D_count; c6 = ADJUSTL(c6) L_D_posterior_pdf%name = TRIM(L_D_posterior_pdf%name)//" {N="//TRIM(c6)// "}" CALL Normalize_pdf (L_D_posterior_pdf) ELSE L_D_posterior_pdf = L_D_prior_pdf END IF IF (L_count > 0) THEN WRITE (c6, "(I6)") L_count; c6 = ADJUSTL(c6) L_L_posterior_pdf%name = TRIM(L_L_posterior_pdf%name)//" {N="//TRIM(c6)// "}" CALL Normalize_pdf (L_L_posterior_pdf) ELSE L_L_posterior_pdf = L_L_prior_pdf END IF IF (N_count > 0) THEN WRITE (c6, "(I6)") N_count; c6 = ADJUSTL(c6) L_N_posterior_pdf%name = TRIM(L_N_posterior_pdf%name)//" {N="//TRIM(c6)// "}" CALL Normalize_pdf (L_N_posterior_pdf) ELSE L_N_posterior_pdf = L_N_prior_pdf END IF IF (P_count > 0) THEN WRITE (c6, "(I6)") P_count; c6 = ADJUSTL(c6) L_P_posterior_pdf%name = TRIM(L_P_posterior_pdf%name)//" {N="//TRIM(c6)// "}" CALL Normalize_pdf (L_P_posterior_pdf) ELSE L_P_posterior_pdf = L_P_prior_pdf END IF IF (R_count > 0) THEN WRITE (c6, "(I6)") R_count; c6 = ADJUSTL(c6) L_R_posterior_pdf%name = TRIM(L_R_posterior_pdf%name)//" {N="//TRIM(c6)// "}" CALL Normalize_pdf (L_R_posterior_pdf) ELSE L_R_posterior_pdf = L_R_prior_pdf END IF IF (T_count > 0) THEN WRITE (c6, "(I6)") T_count; c6 = ADJUSTL(c6) L_T_posterior_pdf%name = TRIM(L_T_posterior_pdf%name)//" {N="//TRIM(c6)// "}" CALL Normalize_pdf (L_T_posterior_pdf) ELSE L_T_posterior_pdf = L_T_prior_pdf END IF CALL Integrate_pdf_to_cpf (L_D_posterior_pdf, L_D_posterior_cpf) ! (sets x-scale of cpf, but won't try to normalize if pdf == 0.0) CALL Integrate_pdf_to_cpf (L_L_posterior_pdf, L_L_posterior_cpf) CALL Integrate_pdf_to_cpf (L_N_posterior_pdf, L_N_posterior_cpf) CALL Integrate_pdf_to_cpf (L_P_posterior_pdf, L_P_posterior_cpf) CALL Integrate_pdf_to_cpf (L_R_posterior_pdf, L_R_posterior_cpf) CALL Integrate_pdf_to_cpf (L_T_posterior_pdf, L_T_posterior_cpf) CALL Plot_3_Distributions(L_D_posterior_pdf, L_D_posterior_cpf, L_L_posterior_pdf, L_L_posterior_cpf, L_N_posterior_pdf, L_N_posterior_cpf, .TRUE.) CALL CLEARSCREEN ($GCLEARSCREEN) CALL Plot_3_Distributions(L_P_posterior_pdf, L_P_posterior_cpf, L_R_posterior_pdf, L_R_posterior_cpf, L_T_posterior_pdf, L_T_posterior_cpf, .TRUE.) CALL CLEARSCREEN ($GCLEARSCREEN) CALL Write_6_pdfs (L_D_posterior_pdf, L_L_posterior_pdf, L_N_posterior_pdf, L_P_posterior_pdf, L_R_posterior_pdf, L_T_posterior_pdf) CONTAINS SUBROUTINE Add_pdf_to_accumulator (pdf1, weight, accumulator_pdf) IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 DOUBLE PRECISION, INTENT(IN) :: weight ! dimensionless; of order 1.D0 TYPE (pdf), INTENT(INOUT) :: accumulator_pdf DOUBLE PRECISION :: x INTEGER :: i IF (weight /= 0.0D0) THEN DO i = 1, nbins x = accumulator_pdf%xlist(i) accumulator_pdf%pd(i) = accumulator_pdf%pd(i) + weight * Get_One_pd (pdf1, x) END DO END IF END SUBROUTINE Add_pdf_to_accumulator CHARACTER*9 FUNCTION ASCII9(x) ! Returns a right-justified 9-byte (or shorter) version of a ! floating-point number, in "human" format, with no more ! than 3 significant digits. IMPLICIT NONE REAL, INTENT(IN) :: x CHARACTER*9 :: temp9 CHARACTER*19 :: temp19 INTEGER :: j, k1, k9, zeros LOGICAL :: punt REAL :: x_log DOUBLE PRECISION :: y IF (isNan(x)) THEN ASCII9="NaN " RETURN ELSE IF (x == 0.0) THEN ASCII9=" 0" RETURN ELSE IF (x > 0.0) THEN punt = (x >= 999999999.5).OR.(x < 0.0000100) ELSE ! x < 0.0 punt = (x <= -99999999.5).OR.(x > -0.000100) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp9,'(1P,E9.2)') x !consider possible improvements, from left to right: IF (temp9(3:5) == '.00') THEN ! right-shift over it temp19(6:9) = temp9(6:9) temp19(4:5) = temp9(1:2) temp19(1:3) = ' ' temp9 = temp19(1:9) ELSE IF (temp9(5:5) == '0') THEN ! right-shift over it temp19(6:9) = temp9(6:9) temp19(2:5) = temp9(1:4) temp19(1:1) = ' ' temp9 = temp19(1:9) END IF IF (temp9(7:7) == '+') THEN ! right-shift over it temp19(8:9) = temp9(8:9) temp19(2:7) = temp9(1:6) temp19(1:1) = ' ' temp9 = temp19(1:9) END IF IF (temp9(8:8) == '0') THEN ! right-shift over it temp19(9:9) = temp9(9:9) temp19(2:8) = temp9(1:7) temp19(1:1) = ' ' temp9 = temp19(1:9) END IF ASCII9 = temp9 ELSE ! can represent without exponential notation x_log = LOG10(ABS(x)) zeros = Int_Below(x_log) - 2 y = (10.D0**zeros) * NINT(ABS(x) / (10.D0**zeros)) IF (x < 0.0) y = -y WRITE (temp19,"(F19.8)") y ! byte 11 is the '.' !Avoid results like "0.7400001" due to rounding error! IF (temp19(18:19) == '01') temp19(18:19) = '00' !Find first important byte from right; change 0 -> ' ' k9 = 10 ! (if no non-0 found to right of .) right_to_left: DO j = 19, 12, -1 IF (temp19(j:j) == '0') THEN temp19(j:j) = ' ' ELSE k9 = j EXIT right_to_left END IF END DO right_to_left !put leading (-)0 before . of fractions, if it fits IF (x > 0.0) THEN IF (temp19(10:11) == ' .') temp19(10:11) = '0.' ELSE ! x < 0.0 IF (k9 <= 17) THEN IF (temp19(9:11) == ' -.') temp19(9:11) = '-0.' END IF END IF k1 = k9 - 8 ASCII9 = temp19(k1:k9) END IF ! punt, or not END FUNCTION ASCII9 DOUBLE PRECISION FUNCTION Del_of_2_cpfs(cpf1, cpf2) !Computes difference measure: 0.5D0 <= del <= 1.0D0 IMPLICIT NONE TYPE (cpf), INTENT(IN) :: cpf1, cpf2 DOUBLE PRECISION P_X_lt_y, P_Z_lt_y, sup, test, y INTEGER :: k sup = 0.0D0 DO k = 0, nbins !test at points which are natural/easy in cpf1: y = cpf1%xlist(k) P_X_lt_y = cpf1%P(k) P_Z_lt_y = Get_One_cp(cpf2, y) test = P_X_lt_y - 2.0D0 * P_X_lt_y * P_Z_lt_y + P_Z_lt_y sup = MAX(sup, test) !also test at points which are natural/easy in cpf2: y = cpf2%xlist(k) P_Z_lt_y = cpf2%P(k) P_X_lt_y = Get_One_cp(cpf1, y) test = P_X_lt_y - 2.0D0 * P_X_lt_y * P_Z_lt_y + P_Z_lt_y sup = MAX(sup, test) END DO Del_of_2_cpfs = sup Del_of_2_cpfs = MAX(MIN(Del_of_2_cpfs, 1.0D0), 0.5D0) ! Any values outside this range are due to imperfect sup{} operator. END FUNCTION Del_of_2_cpfs SUBROUTINE Differentiate_cpf_to_pdf (cpf1, pdf1) IMPLICIT NONE TYPE (cpf), INTENT(IN) :: cpf1 TYPE (pdf), INTENT(OUT) :: pdf1 INTEGER :: i !DON'T ASSUME that pdf was initialized; do it HERE! pdf1%name = cpf1%name pdf1%x_min = cpf1%x_min pdf1%x_max = cpf1%x_max CALL Fill_pdf_xlist(pdf1) DO i = 1, nbins pdf1%pd(i) = (cpf1%P(i) - cpf1%P(i-1)) / (cpf1%xlist(i) - cpf1%xlist(i-1)) END DO CALL Normalize_pdf (pdf1) ! just for safety END SUBROUTINE Differentiate_cpf_to_pdf LOGICAL FUNCTION Digit(c1) !Returns .TRUE. if c1 == one of {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, .} IMPLICIT NONE CHARACTER(1), INTENT(IN) :: c1 Digit = .FALSE. IF (c1 == '0') THEN Digit = .TRUE. RETURN END IF IF (c1 == '1') THEN Digit = .TRUE. RETURN END IF IF (c1 == '2') THEN Digit = .TRUE. RETURN END IF IF (c1 == '3') THEN Digit = .TRUE. RETURN END IF IF (c1 == '4') THEN Digit = .TRUE. RETURN END IF IF (c1 == '5') THEN Digit = .TRUE. RETURN END IF IF (c1 == '6') THEN Digit = .TRUE. RETURN END IF IF (c1 == '7') THEN Digit = .TRUE. RETURN END IF IF (c1 == '8') THEN Digit = .TRUE. RETURN END IF IF (c1 == '9') THEN Digit = .TRUE. RETURN END IF IF (c1 == '.') THEN Digit = .TRUE. RETURN END IF END FUNCTION Digit SUBROUTINE Dump_3_PDFs_for_Excel (pdf1, pdf2, pdf3) IMPLICIT NONE TYPE (pdf) :: pdf1, pdf2, pdf3 CHARACTER*2 :: c2 CHARACTER*12 :: file_name = "Excel_01.dat" INTEGER :: i INTEGER :: Excel_file_index = 1 SAVE WRITE (c2, "(I2)") Excel_file_index IF (c2(1:1) == ' ') c2(1:1) = '0' file_name(7:8) = c2 OPEN (UNIT = 99, FILE = file_name) WRITE (99, "(A)") TRIM(pdf1%name) DO i = 1, nbins WRITE (99, "(2ES12.4)") pdf1%xlist(i), pdf1%pd(i) END DO ! N.B. This format is wasteful of CR/LF bytes, but it permits the PDF to be ! captured in Excel for easy plotting of the pdf as an x/y piecewise-linear graph or bar graph. WRITE (99, "(A)") TRIM(pdf2%name) DO i = 1, nbins WRITE (99, "(2ES12.4)") pdf2%xlist(i), pdf2%pd(i) END DO ! N.B. This format is wasteful of CR/LF bytes, but it permits the PDF to be ! captured in Excel for easy plotting of the pdf as an x/y piecewise-linear graph or bar graph. WRITE (99, "(A)") TRIM(pdf3%name) DO i = 1, nbins WRITE (99, "(2ES12.4)") pdf3%xlist(i), pdf3%pd(i) END DO ! N.B. This format is wasteful of CR/LF bytes, but it permits the PDF to be ! captured in Excel for easy plotting of the pdf as an x/y piecewise-linear graph or bar graph. CLOSE (99) Excel_file_index = Excel_file_index + 1 END SUBROUTINE Dump_3_PDFs_for_Excel SUBROUTINE Fill_cpf_xlist (cpf1) !Uses cpf attributes %x_min and %x_max to decide whether %logarithmic, !and then fills in the nbins (global parameter) x-values at the right edges of bins. !Also, the %x_min value is entered in %xlist(0). IMPLICIT NONE TYPE (cpf), INTENT(INOUT) :: cpf1 DOUBLE PRECISION :: fraction, high_log, log_x, low_log, x INTEGER :: i IF (cpf1%x_max <= cpf1%x_min) THEN WRITE (*, "(' ERROR in Fill_cpf_xlist: %x_min ',ES9.2,' >= %x_max ',ES9.2)") cpf1%x_min, cpf1%x_max WRITE (*, "(' ',A)") TRIM(cpf1%name) CALL Pause() STOP END IF cpf1%logarithmic = .FALSE. ! initializing; may change below IF (cpf1%x_min > 0.0D0) THEN IF (cpf1%x_max / cpf1%x_min > 10.0D0) THEN ! logarithmic cpf1%logarithmic = .TRUE. low_log = DLOG(cpf1%x_min) high_log = DLOG(cpf1%x_max) cpf1%xlist(0) = cpf1%x_min DO i = 1, nbins fraction = (1.0D0 * i) / (1.0D0 * nbins) log_x = low_log + fraction * (high_log - low_log) x = DEXP(log_x) cpf1%xlist(i) = x END DO END IF END IF IF (.NOT.cpf1%logarithmic) THEN ! fill in linear scale cpf1%xlist(0) = cpf1%x_min DO i = 1, nbins fraction = (1.0D0 * i) / (1.0D0 * nbins) x = cpf1%x_min + fraction * (cpf1%x_max - cpf1%x_min) cpf1%xlist(i) = x END DO END IF END SUBROUTINE Fill_cpf_xlist SUBROUTINE Fill_pdf_xlist (pdf1) !Uses pdf attributes %x_min and %x_max to decide whether %logarithmic, !and then fills in the nbins (global parameter) x-values at the centers of bins. IMPLICIT NONE TYPE (pdf), INTENT(INOUT) :: pdf1 DOUBLE PRECISION :: fraction, high_log, log_x, low_log, x INTEGER :: i IF (pdf1%x_max <= pdf1%x_min) THEN WRITE (*, "(' ERROR in Fill_pdf_xlist: %x_min ',ES9.2,' >= %x_max ',ES9.2)") pdf1%x_min, pdf1%x_max WRITE (*, "(' ',A)") TRIM(pdf1%name) CALL Pause() STOP END IF pdf1%logarithmic = .FALSE. ! initializing; may change below IF (pdf1%x_min > 0.0D0) THEN IF (pdf1%x_max / pdf1%x_min > 10.0D0) THEN ! logarithmic pdf1%logarithmic = .TRUE. low_log = DLOG(pdf1%x_min) high_log = DLOG(pdf1%x_max) pdf1%rxlist(0) = pdf1%x_min DO i = 1, nbins fraction = (i - 0.5D0) / (1.0D0 * nbins) log_x = low_log + fraction * (high_log - low_log) x = DEXP(log_x) pdf1%xlist(i) = x fraction = (1.0D0 * i) / (1.0D0 * nbins) log_x = low_log + fraction * (high_log - low_log) x = DEXP(log_x) pdf1%rxlist(i) = x END DO pdf1%rxlist(nbins) = pdf1%x_max END IF END IF IF (.NOT.pdf1%logarithmic) THEN ! fill in linear scale pdf1%rxlist(0) = pdf1%x_min DO i = 1, nbins fraction = (i - 0.5D0) / (1.0D0 * nbins) x = pdf1%x_min + fraction * (pdf1%x_max - pdf1%x_min) pdf1%xlist(i) = x fraction = (1.0D0 * i) / (1.0D0 * nbins) x = pdf1%x_min + fraction * (pdf1%x_max - pdf1%x_min) pdf1%rxlist(i) = x END DO pdf1%rxlist(nbins) = pdf1%x_max END IF END SUBROUTINE Fill_pdf_xlist SUBROUTINE Finish_up_fault() ! Compute disagreements and merge pdfs into combined best-estimate pdf, ! write both to output, empty the storehouse, etc. ! Note that this is only a SUB because it can be called from two different places !(where new fault title line is read, and where EOF is read), ! not because any name-substitution is used. ! In fact, there are no arguments, as all variables are in the global main program! IMPLICIT NONE CHARACTER(1), PARAMETER :: tab = CHAR(9) ! special "HT" tab character in ASCII sequence CHARACTER(1) :: c1, grade CHARACTER(2) :: c2 CHARACTER(3) :: c3 CHARACTER(7) :: del_c7 CHARACTER(9) :: low_limit_c9, mode_c9, median_c9, mean_c9, high_limit_c9, sigma_c9, & & show_offset_km_c9, show_offset_sigma_km_c9, show_began_Ma_c9, & & show_ended_Ma_c9, show_neorate_mmpa_c9, show_neorate_sigma_mmpa_c9 CHARACTER(10) :: Delta_Delta_t_c10 CHARACTER(21) :: file_name = "Excel_01_combined.dat" CHARACTER(40) :: combined_data_name CHARACTER(500) :: line DOUBLE PRECISION, PARAMETER :: one = 1.0D0 DOUBLE PRECISION :: age_Ma, alpha_1, alpha_2, alpha_3, alpha_i, del, Delta_Delta_t, L_max, L_plate, l_prime, multiplier INTEGER :: i, j INTEGER :: Excel_file_index = 1 TYPE (pdf) :: L_combined_pdf, datum_1_pdf, datum_2_pdf, datum_3_pdf, X_pdf TYPE (pdf) :: normalized_pdf1_pdf2, normalized_pdf2_pdf3, normalized_pdf1_pdf3, normalized_pdf1_pdf2_pdf3 TYPE (cpf) :: L_combined_cpf, X_cpf, Z_cpf SAVE IF (old_offset_type_c1 == 'D') THEN L_X_prior_pdf = L_D_prior_pdf ELSE IF (old_offset_type_c1 == 'L') THEN L_X_prior_pdf = L_L_prior_pdf ELSE IF (old_offset_type_c1 == 'N') THEN L_X_prior_pdf = L_N_prior_pdf ELSE IF (old_offset_type_c1 == 'P') THEN L_X_prior_pdf = L_P_prior_pdf ELSE IF (old_offset_type_c1 == 'R') THEN L_X_prior_pdf = L_R_prior_pdf ELSE IF (old_offset_type_c1 == 'T') THEN L_X_prior_pdf = L_T_prior_pdf ELSE ! might by '?', for example, if input_file_type == 2, and warning_count > 0. RETURN END IF !--------------------------------------------------------------------------------------------------------- IF (number_in_store >= 2) THEN ! compute and output del's: !Compute disagreements for all possible pairs: DO i = 2, number_in_store DO j = 1, (i-1) X_cpf = rate_cpf_store(i) Z_cpf = rate_cpf_store(j) del = Del_of_2_cpfs(X_cpf, Z_cpf) WRITE (del_c7, "(F7.4)") del ! (N.B. Printing more digits would just expose values < 0.5 due to imperfect sup{} operator.) Delta_Delta_t = ABS(Delta_t_for_alpha_store(i) - Delta_t_for_alpha_store(j)) ! in Ma, non-negative !N.B. Delta_t_for_alpha_store is different from Delta_t_store because it is set to minimum_age_Ma ! in cases of pinning overlap assemblages (F_is_zero = .TRUE.). ! This is because an unfaulted overlap assemblage was NEVER faulted at ANY time; ! it is a stronger constraint than merely the nominal constraint of zero AVERAGE offset rate; ! it is a constraint of zero offset rate at all times since it was laid down. Delta_Delta_t = MAX(Delta_Delta_t, minimum_age_Ma) ! Avoid "0" which cannot be plotted on a log scale. WRITE (Delta_Delta_t_c10, "(ES10.2)") Delta_Delta_t line = del_c7 // tab // old_Findex // tab // old_offset_type_c1 // tab // TRIM(old_fault_name) // tab // & & TRIM(ADJUSTL(Delta_Delta_t_c10)) // tab // & & grade_store(i) // tab // TRIM(datum_name_store(i)) // tab // & & grade_store(j) // tab // TRIM(datum_name_store(j)) WRITE (3, "(A)") TRIM(line) ! to disagreements data base on UNIT 3. END DO END DO END IF ! number_in_store >= 2; can compute del's !--------------------------------------------------------------------------------------------------------- !------------------------------------------------------------------------------ !Survey list of data in store, checking for the special case of !at least one invalid pin !(because there is at least one younger postive offset). !If so, set asterisk_on_data_count = .TRUE., !and then reduce the data count by eliminating all invalid (older) pin data! at_least_1_positive_offset = .FALSE. ! unless... IF (number_in_store >= 1) THEN DO i = 1, number_in_store IF (positive_offset_query(i)) at_least_1_positive_offset = .TRUE. END DO END IF IF (at_least_1_positive_offset) THEN at_least_1_invalid_pin = .FALSE. ! unless... IF (number_in_store >= 2) THEN DO i = 1, number_in_store ! evaluate invalid_pin_query(i)... invalid_pin_query(i) = .FALSE. ! unless... IF (.NOT.positive_offset_query(i)) THEN ! datum i is a pin DO j = 1, number_in_store IF (j /= i) THEN invalid_pin_query(i) = invalid_pin_query(i) .OR. & & (positive_offset_query(j) .AND. & & (Delta_t_store(j) < Delta_t_store(i))) END IF ! j /= i (duh!) END DO ! checking all data j for invalidating conflict END IF ! datum i is a pin IF (invalid_pin_query(i)) at_least_1_invalid_pin = .TRUE. END DO ! considering all data i, to determine invalid_pin_query(i) END IF ! at least 2 data in store asterisk_on_data_count = at_least_1_invalid_pin ELSE asterisk_on_data_count = .FALSE. END IF ! at_least_one_positive_offset, or not IF (asterisk_on_data_count) THEN !Go through storehouse and delete any data for which invalid_pin_query() = .TRUE.: total_deletions = 0 ! just initializing... DO i = 1, number_in_store IF (invalid_pin_query(i)) total_deletions = total_deletions + 1 END DO ! total deletions count is now complete DO j = 1, total_deletions !find index of datum to delete checking_validity: DO i = 1, number_in_store IF (invalid_pin_query(i)) EXIT checking_validity ! this should always occur! END DO checking_validity !i now points to the datum which is to be deleted: IF (i < number_in_store) THEN ! rows must be moved DO k = i, (number_in_store - 1) ! k points to datum that will be replaced with next-higher row datum_name_store(k) = datum_name_store(k + 1) grade_store(k) = grade_store(k + 1) Delta_t_store(k) = Delta_t_store(k + 1) Delta_t_for_alpha_store(k) = Delta_t_for_alpha_store(k + 1) rate_pdf_store(k) = rate_pdf_store(k + 1) rate_cpf_store(k) = rate_cpf_store(k + 1) positive_offset_query(k) = positive_offset_query(k + 1) END DO ! k = datum to be overwritten with next-higher row END IF ! rows must be moved number_in_store = number_in_store - 1 ! because last datum is now EITHER invalid or duplicate END DO ! j = 1, total_deletions END IF ! asterisk_on_data_count !------------------------------------------------------------------------------ !Compute combined pdf, selecting appropriate equation based on number of data (remaining): IF (number_in_store == 0) THEN !Best-estimate pdf is just the prior for this offset class. L_combined_pdf = L_X_prior_pdf ELSE IF (number_in_store == 1) THEN !Best-estimate pdf is mostly the one datum, with a small fraction of the prior: grade = grade_store(1) age_Ma = Delta_t_for_alpha_store(1) alpha_1 = Model_alpha(grade, age_Ma) datum_1_pdf = rate_pdf_store(1) L_combined_pdf%x_max = MAX(datum_1_pdf%x_max, L_X_prior_pdf%x_max) L_combined_pdf%x_min = L_combined_pdf%x_max * 1.D-08 CALL Fill_pdf_xlist(L_combined_pdf) DO i = 1, nbins l_prime = L_combined_pdf%xlist(i) L_combined_pdf%pd(i) = alpha_1 * Get_one_pd(L_X_prior_pdf, l_prime) + & & (one - alpha_1) * Get_one_pd(datum_1_pdf, l_prime) END DO ELSE IF (number_in_store == 2) THEN !datum_1: grade = grade_store(1) age_Ma = Delta_t_for_alpha_store(1) alpha_1 = Model_alpha(grade, age_Ma) datum_1_pdf = rate_pdf_store(1) !datum_2: grade = grade_store(2) age_Ma = Delta_t_for_alpha_store(2) alpha_2 = Model_alpha(grade, age_Ma) datum_2_pdf = rate_pdf_store(2) !normalized_pdf1_pdf2: normalized_pdf1_pdf2%x_max = MAX(datum_1_pdf%x_max, datum_2_pdf%x_max) normalized_pdf1_pdf2%x_min = normalized_pdf1_pdf2%x_max * 1.0D-08 CALL Fill_pdf_xlist(normalized_pdf1_pdf2) DO i = 1, nbins l_prime = normalized_pdf1_pdf2%xlist(i) normalized_pdf1_pdf2%pd(i) = Get_one_pd(datum_1_pdf, l_prime) * & & Get_one_pd(datum_2_pdf, l_prime) END DO CALL Normalize_pdf(normalized_pdf1_pdf2) !L_combined_pdf: L_combined_pdf%x_max = MAX(datum_1_pdf%x_max, datum_2_pdf%x_max, L_X_prior_pdf%x_max) L_combined_pdf%x_min = L_combined_pdf%x_max * 1.0D-08 CALL Fill_pdf_xlist(L_combined_pdf) DO i = 1, nbins l_prime = L_combined_pdf%xlist(i) L_combined_pdf%pd(i) = alpha_1 * alpha_2 * Get_one_pd(L_X_prior_pdf, l_prime) + & & (one - alpha_1) * alpha_2 * Get_one_pd(datum_1_pdf, l_prime) + & & alpha_1 * (one - alpha_2) * Get_one_pd(datum_2_pdf, l_prime) + & & (one - alpha_1) * (one - alpha_2) * Get_one_pd(normalized_pdf1_pdf2, l_prime) END DO ELSE IF (number_in_store == 3) THEN !datum_1: grade = grade_store(1) age_Ma = Delta_t_for_alpha_store(1) alpha_1 = Model_alpha(grade, age_Ma) datum_1_pdf = rate_pdf_store(1) !datum_2: grade = grade_store(2) age_Ma = Delta_t_for_alpha_store(2) alpha_2 = Model_alpha(grade, age_Ma) datum_2_pdf = rate_pdf_store(2) !datum_3: grade = grade_store(3) age_Ma = Delta_t_for_alpha_store(3) alpha_3 = Model_alpha(grade, age_Ma) datum_3_pdf = rate_pdf_store(3) !normalized_pdf1_pdf2: normalized_pdf1_pdf2%x_max = MAX(datum_1_pdf%x_max, datum_2_pdf%x_max) normalized_pdf1_pdf2%x_min = normalized_pdf1_pdf2%x_max * 1.0D-08 CALL Fill_pdf_xlist(normalized_pdf1_pdf2) DO i = 1, nbins l_prime = normalized_pdf1_pdf2%xlist(i) normalized_pdf1_pdf2%pd(i) = Get_one_pd(datum_1_pdf, l_prime) * & & Get_one_pd(datum_2_pdf, l_prime) END DO CALL Normalize_pdf(normalized_pdf1_pdf2) !normalized_pdf2_pdf3: normalized_pdf2_pdf3%x_max = MAX(datum_2_pdf%x_max, datum_3_pdf%x_max) normalized_pdf2_pdf3%x_min = normalized_pdf2_pdf3%x_max * 1.0D-08 CALL Fill_pdf_xlist(normalized_pdf2_pdf3) DO i = 1, nbins l_prime = normalized_pdf2_pdf3%xlist(i) normalized_pdf2_pdf3%pd(i) = Get_one_pd(datum_2_pdf, l_prime) * & & Get_one_pd(datum_3_pdf, l_prime) END DO CALL Normalize_pdf(normalized_pdf2_pdf3) !normalized_pdf1_pdf3: normalized_pdf1_pdf3%x_max = MAX(datum_1_pdf%x_max, datum_3_pdf%x_max) normalized_pdf1_pdf3%x_min = normalized_pdf1_pdf3%x_max * 1.0D-08 CALL Fill_pdf_xlist(normalized_pdf1_pdf3) DO i = 1, nbins l_prime = normalized_pdf1_pdf3%xlist(i) normalized_pdf1_pdf3%pd(i) = Get_one_pd(datum_1_pdf, l_prime) * & & Get_one_pd(datum_3_pdf, l_prime) END DO CALL Normalize_pdf(normalized_pdf1_pdf3) !normalized_pdf1_pdf2_pdf3: normalized_pdf1_pdf2_pdf3%x_max = MAX(datum_1_pdf%x_max, datum_2_pdf%x_max, datum_3_pdf%x_max) normalized_pdf1_pdf2_pdf3%x_min = normalized_pdf1_pdf2_pdf3%x_max * 1.0D-08 CALL Fill_pdf_xlist(normalized_pdf1_pdf2_pdf3) DO i = 1, nbins l_prime = normalized_pdf1_pdf2_pdf3%xlist(i) normalized_pdf1_pdf2_pdf3%pd(i) = Get_one_pd(datum_1_pdf, l_prime) * & & Get_one_pd(datum_2_pdf, l_prime) * & & Get_one_pd(datum_3_pdf, l_prime) END DO CALL Normalize_pdf(normalized_pdf1_pdf2_pdf3) !L_combined_pdf: L_combined_pdf%x_max = MAX(datum_1_pdf%x_max, datum_2_pdf%x_max, datum_3_pdf%x_max, L_X_prior_pdf%x_max) L_combined_pdf%x_min = L_combined_pdf%x_max * 1.0D-08 CALL Fill_pdf_xlist(L_combined_pdf) DO i = 1, nbins l_prime = L_combined_pdf%xlist(i) L_combined_pdf%pd(i) = alpha_1 * alpha_2 * alpha_3 * Get_one_pd(L_X_prior_pdf, l_prime) + & & (one - alpha_1) * alpha_2 * alpha_3 * Get_one_pd(datum_1_pdf, l_prime) + & & alpha_1 * (one - alpha_2) * alpha_3 * Get_one_pd(datum_2_pdf, l_prime) + & & alpha_1 * alpha_2 * (one - alpha_3) * Get_one_pd(datum_3_pdf, l_prime) + & & alpha_1 * (one - alpha_2) * (one - alpha_3) * Get_one_pd(normalized_pdf2_pdf3, l_prime) + & & (one - alpha_1) * alpha_2 * (one - alpha_3) * Get_one_pd(normalized_pdf1_pdf3, l_prime) + & & (one - alpha_1) * (one - alpha_2) * alpha_3 * Get_one_pd(normalized_pdf1_pdf2, l_prime) + & & (one - alpha_1) * (one - alpha_2) * (one - alpha_3) * Get_one_pd(normalized_pdf1_pdf2_pdf3, l_prime) END DO ELSE ! number_in_store > 3; use general asymptotic formula: !Find L_max to be a bound on l_prime in L_combined_pdf: L_max = L_X_prior_pdf%x_max DO i = 1, number_in_store X_cpf = rate_cpf_store(i) L_max = MAX(L_max, X_cpf%x_max) END DO L_plate = plate_rate ! unless.... IF (old_offset_type_c1 == 'N') L_plate = plate_rate * DTAN(N_dip_radians) IF (old_offset_type_c1 == 'T') L_plate = plate_rate * DTAN(T_dip_radians) !L_max = MIN(L_max, L_plate) !Initialize L_combined_pdf: L_combined_pdf%x_max = L_max L_combined_pdf%x_min = L_max * 1.0D-08 ! forcing logarithmic sampling, which is CRITICAL when some of the data show ~1 cm/10 Ma rates! CALL Fill_pdf_xlist(L_combined_pdf) !Set values to those of the first term in the extended product: grade = grade_store(1) age_Ma = Delta_t_for_alpha_store(1) alpha_1 = Model_alpha(grade, age_Ma) X_pdf = rate_pdf_store(1) DO j = 1, nbins l_prime = L_combined_pdf%xlist(j) L_combined_pdf%pd(j) = (alpha_1 / L_max) + & & (one - alpha_1) * Get_one_pd(X_pdf, l_prime) END DO CALL Normalize_pdf(L_combined_pdf) ! (probably not necessary) !Create and multiply other terms (without storing them), one by one... DO i = 2, number_in_store grade = grade_store(i) age_Ma = Delta_t_for_alpha_store(i) alpha_i = Model_alpha(grade, age_Ma) X_pdf = rate_pdf_store(i) DO j = 1, nbins l_prime = L_combined_pdf%xlist(j) L_combined_pdf%pd(j) = L_combined_pdf%pd(j) * & & ((alpha_i / L_max) + & & (one - alpha_i) * Get_one_pd(X_pdf, l_prime)) END DO ! j = 1, nbins CALL Normalize_pdf(L_combined_pdf) END DO ! i = 2, number_in_store END IF ! number_in_store == 0; 1; 2; 3; or more CALL Normalize_pdf(L_combined_pdf) ! just in case... CALL Integrate_pdf_to_cpf (L_combined_pdf, L_combined_cpf) WRITE (c3, "(I3)") number_in_store IF (asterisk_on_data_count) THEN IF (input_file_type == 1) THEN combined_data_name = "combined offset rate (" // TRIM(ADJUSTL(c3)) // "* GCN data)" ELSE IF (input_file_type == 2) THEN combined_data_name = "combined offset rate (" // TRIM(ADJUSTL(c3)) // "* WGCEP data)" ELSE combined_data_name = "combined offset rate (" // TRIM(ADJUSTL(c3)) // "* ??? data)" END IF ELSE IF (input_file_type == 1) THEN combined_data_name = "combined offset rate (" // TRIM(ADJUSTL(c3)) // " GCN data)" ELSE IF (input_file_type == 2) THEN combined_data_name = "combined offset rate (" // TRIM(ADJUSTL(c3)) // " WGCEP data)" ELSE combined_data_name = "combined offset rate (" // TRIM(ADJUSTL(c3)) // " ??? data)" END IF END IF L_combined_pdf%name = TRIM(combined_data_name) !----- output combined rate to file: ------------------------------- low_limit_c9 = ADJUSTL(ASCII9(REAL(Get_x_from_P(L_combined_cpf, 0.025D0)))) mode_c9 = ADJUSTL(ASCII9(REAL(Get_Mode(L_combined_pdf)))) ! mode median_c9 = ADJUSTL(ASCII9(REAL(Get_Median(L_combined_cpf)))) ! median mean_c9 = ADJUSTL(ASCII9(REAL(Get_Mean(L_combined_pdf)))) ! mean high_limit_c9 = ADJUSTL(ASCII9(REAL(Get_x_from_P(L_combined_cpf, 0.975D0)))) sigma_c9 = ADJUSTL(ASCII9(REAL(Get_Sigma(L_combined_pdf)))) ! standard deviation !Do not report "Unbounded" or "Unreasonable" numbers with digits; use code 'U' as described in paper: IF (Get_x_from_P(L_combined_cpf, 0.025D0) > (1.5D0 * plate_rate)) low_limit_c9 = 'U' IF (Get_Mode(L_combined_pdf) > (1.5D0 * plate_rate)) mode_c9 = 'U' IF (Get_Median(L_combined_cpf) > (1.5D0 * plate_rate)) median_c9 = 'U' IF (Get_Mean(L_combined_pdf) > (1.5D0 * plate_rate)) mean_c9 = 'U' IF (Get_x_from_P(L_combined_cpf, 0.975D0) > (1.5D0 * plate_rate)) high_limit_c9 = 'U' line = old_Findex // tab // old_offset_type_c1 // tab // TRIM(old_fault_name) // tab // & & 'S' // tab // TRIM(combined_data_name)// tab // & ! N.B. "Grade" entry is set to 'S' for "Slippery". & tab // tab // & ! N. B. Offset and offset sigma are blank. & tab // tab // & ! N. B. began_Ma and ended_Ma are blank. & TRIM(low_limit_c9) // tab // TRIM(mode_c9) // tab // TRIM(median_c9) // tab // TRIM(mean_c9) // tab // TRIM(high_limit_c9) // tab // & & TRIM(sigma_c9) WRITE (2, "(A)") TRIM(line) !----- output combined rate to screen: ------------------------------- CALL CLEARSCREEN ($GCLEARSCREEN) show_string_c40 = old_fault_name(1:40) IF (LEN_TRIM(old_fault_name(41:50)) > 0) show_string_c40(40:40) = 'M' show_offset_km_c19 = old_editors_offset_km_c9 show_offset_sigma_km_c9 = old_editors_offset_sigma_km_c9 ! (currently not displayed) show_began_Ma_c18 = old_editors_began_Ma_c9 show_ended_Ma_c17 = old_editors_ended_Ma_c9 show_neorate_mmpa_c16 = old_editors_neorate_mmpa_c9 show_neorate_sigma_mmpa_c9 = old_editors_neorate_sigma_mmpa_c9 ! (currently not displayed) show_offset_km_c19 = ADJUSTR(show_offset_km_c19) show_offset_sigma_km_c9 = ADJUSTR(show_offset_sigma_km_c9) show_began_Ma_c18 = ADJUSTR(show_began_Ma_c18) show_ended_Ma_c17 = ADJUSTR(show_ended_Ma_c17) show_neorate_mmpa_c16 = ADJUSTR(show_neorate_mmpa_c16) show_neorate_sigma_mmpa_c9 = ADJUSTR(show_neorate_sigma_mmpa_c9) WRITE (*, "(1X,'SenseV',1X,'Fault',35X,1X,8X,'Offset_(km)',1X,8X,'Began_(Ma)',1X,7X,'Ended_(Ma)',1X,5X,'Rate_(mm/a)' & &/1X, A5, A1,1X, A40, 1X,A19, 1X,A18, 1X,A17, 1X,A16)") & ! 121 bytes after 1X & old_Findex, old_offset_type_c1, show_string_c40, & & show_offset_km_c19, & & show_began_Ma_c18, show_ended_Ma_c17, & & show_neorate_mmpa_c16 !combined offset-rate PDF/CPF rectangle occupies 2.05 < y < 2.75, 0.5 < x < 3.5 (out of 3 x 4 screen) CALL Plot_1_Distribution(L_combined_pdf, L_combined_cpf, 0.5D0, 3.5D0, 2.05D0, 2.75D0) IF (pause_to_see_combined_rate) THEN READ (*, "(A)") c1 IF ((c1 == 'x').OR.(c1 == 'X')) THEN WRITE (c2, "(I2)") Excel_file_index IF (c2(1:1) == ' ') c2(1:1) = '0' file_name(7:8) = c2 OPEN (UNIT = 99, FILE = file_name) WRITE (99, "(A)") TRIM(L_combined_pdf%name) DO i = 1, nbins WRITE (99, "(2ES12.4)") L_combined_pdf%xlist(i), L_combined_pdf%pd(i) END DO ! N.B. This format is wasteful of CR/LF bytes, but it permits the PDF ! captured in Excel for easy plotting of the pdf as an x/y piecewise-linear graph or bar graph. END IF END IF !----- output is now complete for this fault ------------------------- !reset the store to EMPTY: number_in_store = 0 END SUBROUTINE Finish_up_fault DOUBLE PRECISION FUNCTION Get_Mean (pdf1) IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 DOUBLE PRECISION :: dx, sum INTEGER :: i sum = 0.0D0 ! initializing integral DO i = 1, nbins dx = pdf1%rxlist(i) - pdf1%rxlist(i-1) sum = sum + pdf1%pd(i) * pdf1%xlist(i) * dx END DO Get_Mean = sum END FUNCTION Get_Mean DOUBLE PRECISION FUNCTION Get_x_from_P (cpf1, P) IMPLICIT NONE TYPE (cpf), INTENT(IN) :: cpf1 DOUBLE PRECISION, INTENT(IN) :: P DOUBLE PRECISION :: fraction, P_high, P_low, x_left, x_right INTEGER :: i, i_left, i_right LOGICAL :: success IF (P <= 0.0D0) THEN Get_x_from_P = cpf1%x_min ELSE IF (P >= 1.0D0) THEN Get_x_from_P = cpf1%x_max ELSE i_left = -1; i_right = -1 ! should be changed during search below search : DO i = 1, nbins IF ((cpf1%P(i-1) <= P).AND.(cpf1%P(i) >= P)) THEN i_right = i i_left = i - 1 EXIT search END IF END DO search success = (i_right > 0) IF (success) THEN x_left = cpf1%xlist(i_left) x_right = cpf1%xlist(i_right) P_low = cpf1%P(i_left) P_high = cpf1%P(i_right) IF (P_high > P_low) THEN fraction = (P - P_low) / (P_high - P_low) ELSE fraction = 0.0D0 END IF Get_x_from_P = x_left + fraction * (x_right - x_left) ELSE Get_x_from_P = 0.0D0 END IF END IF END FUNCTION Get_x_from_P DOUBLE PRECISION FUNCTION Get_Median (cpf1) IMPLICIT NONE TYPE (cpf), INTENT(IN) :: cpf1 Get_Median = Get_x_from_P (cpf1, 0.5D0) END FUNCTION Get_Median DOUBLE PRECISION FUNCTION Get_Mode (pdf1) !Finds value of independent variable (x) which is associated with peak in pdf1. !Caution: IF pdf1 is flat-topped, model is not well-defined; this routine returns first high value. IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 DOUBLE PRECISION :: largest INTEGER :: i, i_best largest = 0.0D0 ! initializing search i_best = 1 DO i = 1, nbins IF (pdf1%pd(i) > largest) THEN i_best = i largest = pdf1%pd(i) END IF END DO Get_Mode = pdf1%xlist(i_best) END FUNCTION Get_Mode DOUBLE PRECISION FUNCTION Get_One_cp (cpf1, x) !extracts one value of cumulative probability (P) from a cumulative probability function, !for a value of the independent variable (x) which is not necessarily a grid point. IMPLICIT NONE TYPE (cpf), INTENT(IN) :: cpf1 DOUBLE PRECISION, INTENT(IN) :: x DOUBLE PRECISION :: fraction INTEGER :: i IF (x < cpf1%x_min) THEN Get_One_cp = 0.0D0 ELSE IF (x > cpf1%x_max) THEN Get_One_cp = 1.0D0 ELSE IF (cpf1%logarithmic) THEN i = NINT(0.5D0 + nbins * (DLOG(x) - DLOG(cpf1%x_min)) / (DLOG(cpf1%x_max) - DLOG(cpf1%x_min))) ELSE ! linear xlist i = NINT(0.5D0 + nbins * (x - cpf1%x_min) / (cpf1%x_max - cpf1%x_min)) END IF i = MAX(1, MIN(i, nbins)) fraction = (x - cpf1%xlist(i-1)) / (cpf1%xlist(i) - cpf1%xlist(i-1)) fraction = MAX(0.0D0, MIN(fraction, 1.0D0)) Get_One_cp = fraction * cpf1%P(i) + (1.0D0- fraction) * cpf1%P(i-1) END IF END FUNCTION Get_One_cp DOUBLE PRECISION FUNCTION Get_One_pd (pdf1, x) IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 DOUBLE PRECISION, INTENT(IN) :: x INTEGER :: i IF (x < pdf1%x_min) THEN Get_One_pd = 0.0D0 ELSE IF (x > pdf1%x_max) THEN Get_One_pd = 0.0D0 ELSE IF (pdf1%logarithmic) THEN i = NINT(0.5D0 + nbins * (DLOG(x) - DLOG(pdf1%x_min)) / (DLOG(pdf1%x_max) - DLOG(pdf1%x_min))) ELSE ! linear xlist i = NINT(0.5D0 + nbins * (x - pdf1%x_min) / (pdf1%x_max - pdf1%x_min)) END IF i = MAX(1, MIN(i, nbins)) Get_One_pd = pdf1%pd(i) END IF END FUNCTION Get_One_pd DOUBLE PRECISION FUNCTION Get_Sigma(pdf1) IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 DOUBLE PRECISION :: dx, mean, sum, variance, x INTEGER :: i !get the mean sum = 0.0D0 ! just initializing DO i = 1, nbins x = pdf1%xlist(i) dx = pdf1%rxlist(i) - pdf1%rxlist(i-1) sum = sum + pdf1%pd(i) * x * dx END DO mean = sum !get the variance sum = 0.0D0 ! just initializing DO i = 1, nbins x = pdf1%xlist(i) dx = pdf1%rxlist(i) - pdf1%rxlist(i-1) sum = sum + pdf1%pd(i) * (x - mean)**2 * dx END DO variance = sum !standard deviation is the square root of the variance Get_Sigma = DSQRT(variance) END FUNCTION Get_Sigma DOUBLE PRECISION FUNCTION Implied_Sigma (string) ! Takes a real number (in form of a text string) and finds its ! implied standard deviation, equal to one-half of the change ! in magnitude caused by a change of one in the right-most non-zero digit. ! Note that string must be readable as a * format number, so any ! special prefix and/or suffix characters such as <, >, ~, ? must be removed first. ! Any of these formats is acceptable: F or E or D or I. ! Position of number within the string is not important. ! Meaningful part of string must not be longer than 80 bytes, ! although padding with spaces on right does no harm. IMPLICIT NONE CHARACTER(*), INTENT(IN) :: string CHARACTER(1) :: c1, c1p CHARACTER(80) :: string_copy DOUBLE PRECISION :: nominal_value, perturbed_value INTEGER :: D_or_E_position, i, ios, length, rightmost_nonzero_position, start_position length = LEN_TRIM(string) IF (length > 80) THEN WRITE (*, "(' ERROR: String longer than 80 bytes sent to Implied_Sigma:',A)") TRIM(string) CALL Pause() STOP END IF READ (string, *, IOSTAT = ios) nominal_value IF (ios /= 0) THEN WRITE (*, "(' ERROR: Unreadable string sent to Implied_Sigma:',A)") TRIM(string) CALL Pause() STOP END IF string_copy = ADJUSTR(string) D_or_E_position = MAX(SCAN(string_copy, 'D'), SCAN(string_copy, 'd'), SCAN(string_copy, 'E'), SCAN(string_copy, 'e')) IF (D_or_E_position > 0) THEN start_position = D_or_E_position - 1 ELSE ! no D or E is present (F or I format): start_position = 80 END IF search: DO i = start_position, 1, -1 c1 = string_copy(i:i) IF ((c1 == '1').OR.(c1 == '2').OR.(c1 == '3').OR.(c1 == '4').OR.(c1 == '5').OR.(c1 == '6').OR.(c1 == '7').OR.(c1 == '8').OR.(c1 == '9')) THEN rightmost_nonzero_position = i EXIT search ! (also keeping value in c1) END IF END DO search IF (c1 == '1') c1p = '2' IF (c1 == '2') c1p = '3' IF (c1 == '3') c1p = '4' IF (c1 == '4') c1p = '5' IF (c1 == '5') c1p = '6' IF (c1 == '6') c1p = '7' IF (c1 == '7') c1p = '8' IF (c1 == '8') c1p = '9' IF (c1 == '9') c1p = '8' string_copy(rightmost_nonzero_position:rightmost_nonzero_position) = c1p READ (string_copy, *, IOSTAT = ios) perturbed_value IF (ios /= 0) THEN WRITE (*, "(' ERROR: Unreadable string created in Implied_Sigma:',A)") TRIM(string_copy) CALL Pause() STOP END IF Implied_Sigma = DABS(perturbed_value - nominal_value) * 0.5D0 END FUNCTION Implied_Sigma INTEGER FUNCTION Int_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL, INTENT(IN) :: x INTEGER :: i REAL :: y i = INT(x) IF (x >= 0.) THEN Int_Below = i ELSE ! x < 0. y = 1.*i IF (y <= x) THEN Int_Below = i ELSE ! most commonly Int_Below = i - 1 END IF END IF END FUNCTION Int_Below SUBROUTINE Integrate_pdf_to_cpf (pdf1, cpf1) IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 TYPE (cpf), INTENT(INOUT) :: cpf1 DOUBLE PRECISION :: dx, sum INTEGER :: i !DON'T assume that the cpf is properly initialized! cpf1%name = TRIM(pdf1%name) cpf1%x_min = pdf1%x_min cpf1%x_max = pdf1%x_max CALL Fill_cpf_xlist(cpf1) !actual integration sum = 0.0D0 cpf1%P(0) = 0.0D0 DO i = 1, nbins dx = pdf1%rxlist(i) - pdf1%rxlist(i-1) sum = sum + dx * pdf1%pd(i) cpf1%P(i) = sum END DO IF (sum > 0.0D0) CALL Normalize_cpf(cpf1) END SUBROUTINE Integrate_pdf_to_cpf DOUBLE PRECISION FUNCTION Model_alpha(grade, age_Ma) !Model based on analysis of f_Gorda-Cal-Nev_edited.txt, 2006.03.27 IMPLICIT NONE CHARACTER(1), INTENT(IN) :: grade ! A, B, or C DOUBLE PRECISION, INTENT(IN) :: age_Ma ! age of offset feature !Component of Model_alpha based on grade reflects fraction of incorrect or unrepresentative* data: ! (*because of variations in slip rate along the trace). IF (grade == 'A') THEN Model_alpha = 0.043D0 ! primary literature source ELSE IF (grade == 'B') THEN Model_alpha = 0.094D0 ! secondary literature source ELSE IF (grade == 'C') THEN Model_alpha = 0.145D0 ! tertiary literature source ELSE WRITE (*, "(' ERROR: Illegal grade (',A,') sent to Model_alpha.')") grade CALL Pause() STOP END IF !Other component of alpha, based on age_Ma, reflects data inapplicable to neotectonics: IF (age_Ma > 3.0D0) THEN Model_alpha = Model_alpha + 0.5 * MIN(ALOG10(REAL(age_Ma)/3.0), 0.87)/((1.0-0.05)*(1.0-0.043)) !Peaks at 0.521 (for grade A) at 22.23 Ma. END IF END FUNCTION Model_alpha SUBROUTINE No_Expectations_for_pdf (pdf1) !Weights negative side of pdf to achieve an overall mean or expectation of zero. !Note that this leaves a non-normalized pdf, so you should call Normalize_pdf(pdf1) afterward. IMPLICIT NONE TYPE (pdf), INTENT(INOUT) :: pdf1 DOUBLE PRECISION :: dx, factor, left_sum, right_sum INTEGER :: i IF ((pdf1%x_min >= 0.0D0).OR.(pdf1%x_max <= 0.0D0)) THEN WRITE (*, "(' ERROR: pdf send to No_Expectations_for_pdf must include both negative and positive sides.')") CALL Pause() STOP END IF left_sum = 0.0D0 ! initializing integral on negative side right_sum = 0.0D0 ! initializing integral on positive side DO i = 1, nbins dx = pdf1%rxlist(i) - pdf1%rxlist(i-1) IF (pdf1%xlist(i) < 0.0D0) THEN left_sum = left_sum + pdf1%pd(i) * pdf1%xlist(i) * dx ELSE right_sum = right_sum + pdf1%pd(i) * pdf1%xlist(i) * dx END IF END DO IF (left_sum >= 0.0D0) THEN WRITE (*, "(' ERROR: pdf send to No_Expectations_for_pdf has all zeros on negative side.')") CALL Pause() STOP END IF factor = right_sum / (-left_sum) adjusting: DO i = 1, nbins IF (pdf1%xlist(i) < 0.0D0) THEN pdf1%pd(i) =pdf1%pd(i) * factor ELSE EXIT adjusting END IF END DO adjusting END SUBROUTINE No_Expectations_for_pdf SUBROUTINE Normalize_cpf (cpf1) IMPLICIT NONE TYPE (cpf), INTENT(INOUT) :: cpf1 DOUBLE PRECISION :: inverse, sum INTEGER :: i sum = cpf1%P(nbins) IF (sum == 1.0D0) RETURN IF (sum == 0.0D0) THEN WRITE (*, "(' ERROR: All-zero cpf sent to Normalize_cpf.')") CALL Pause() STOP END IF inverse = 1.0D0 / sum DO i = 1, nbins cpf1%P(i) = cpf1%P(i) * inverse END DO cpf1%normalized = .TRUE. END SUBROUTINE Normalize_cpf SUBROUTINE Normalize_pdf (pdf1) !If pdf contains nothing but zeros, normalization is (silently) skipped. IMPLICIT NONE TYPE (pdf), INTENT(INOUT) :: pdf1 DOUBLE PRECISION :: dx, inverse, sum INTEGER :: i sum = 0.0D0 DO i = 1, nbins dx = pdf1%rxlist(i) - pdf1%rxlist(i-1) sum = sum + dx * pdf1%pd(i) END DO IF (sum == 1.0D0) RETURN IF (sum == 0.0D0) RETURN ! no error message, and no change to pdf inverse = 1.0D0 / sum DO i = 1, nbins pdf1%pd(i) = pdf1%pd(i) * inverse END DO pdf1%normalized = .TRUE. END SUBROUTINE Normalize_pdf SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause SUBROUTINE Plot_1_Distribution (pdf1, cpf1, x1, x2, y1, y2) !Plots a single pdf/cpf in box x1 < x < x2, y1 < y < y2): IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 TYPE (cpf), INTENT(IN) :: cpf1 DOUBLE PRECISION, INTENT(IN) :: x1, x2, y1, y2 DOUBLE PRECISION :: fraction, largest, wx, wy, x_left, x_median, x_right, zero_fraction CHARACTER(9) :: c9a, c9b, c9c, c9d, c9e, c9f CHARACTER(320) :: text_string INTEGER(2) :: axis_style = #FFFF, result INTEGER(4) :: i, ideg = 0, old_color REAL :: zoom TYPE (wxycoord) wt !determine whether to show full-width of pdf and cdf (zoom = 1.0), !or to only display the portion up to x_of_DF_high = 2 * x:cdf(x)=0.975 (zoom > 1.0) !----------------------------------------------------------------------------------- x_right = Get_x_from_P (cpf1, 0.975D0) zoom = MAX(1.0, ((pdf1%x_max - pdf1%x_min)/(2*x_right - pdf1%x_min))) !--------------------------------------------------------------------- !use simple 1-pixel vertical lines for pdf values: largest = 0.0D0 DO i = 1, nbins largest = MAX(largest, pdf1%pd(i)) END DO IF (largest > 0.0D0) THEN CALL SETLINESTYLE(axis_style) ! solid old_color = SETCOLORRGB(#FF0000) ! blue DO i = 1, nbins !horizontal position fraction = (pdf1%xlist(i) - pdf1%x_min) / (pdf1%x_max - pdf1%x_min) wx = x1 + zoom * fraction * (x2 - x1) IF (wx <= x2) THEN ! only plot if on-page !vertical position fraction = pdf1%pd(i) / largest wy = y1 + fraction * (y2 - y1) CALL MOVETO_W (wx, wy, wt) wy = y1 result = LINETO_W (wx, wy) END IF END DO END IF ! use white line for the cumulative probability: CALL SETLINESTYLE(axis_style) ! solid old_color = SETCOLORRGB(#FFFFFF) ! white wx = x1; wy = y1 CALL MOVETO_W (wx, wy, wt) DO i = 0, nbins !horizontal position fraction = (cpf1%xlist(i) - cpf1%x_min) / (cpf1%x_max - cpf1%x_min) wx = x1 + zoom * fraction * (x2 - x1) IF (wx <= x2) THEN ! only plot if on-page !vertical position fraction = cpf1%P(i) wy = y1 + fraction * (y2 - y1) result = LINETO_W (wx, wy) END IF END DO ! add white cross-bar at median: CALL SETLINESTYLE(axis_style) ! solid old_color = SETCOLORRGB(#FFFFFF) ! white x_median = Get_Median(cpf1) IF (cpf1%x_max > cpf1%x_min) THEN fraction = (x_median - cpf1%x_min) / (cpf1%x_max - cpf1%x_min) ELSE fraction = -0.1D0 END IF IF ((fraction > 0.0D0).AND.(fraction < 1.0D0)) THEN wx = x1 + zoom * fraction * (x2 - x1) wy = y1 + 0.6D0 * (y2 - y1) CALL MOVETO_W (wx, wy, wt) wy = y1 + 0.4D0 * (y2 - y1) result = LINETO_W (wx, wy) END IF ! label the median with a printed value CALL SETGTEXTROTATION(ideg) old_color = SETCOLORRGB(#FFFFFF) ! white result = SETFONT('t''Arial''h12w7') CALL OUTGTEXT(TRIM(ADJUSTL(ASCII9(REAL(Get_Median(cpf1)))))) ! add white cross-bar at x:P(X < x) = 0.025 CALL SETLINESTYLE(axis_style) ! solid old_color = SETCOLORRGB(#FFFFFF) ! white x_right = Get_x_from_P (cpf1, 0.025D0) IF (cpf1%x_max > cpf1%x_min) THEN fraction = (x_right - cpf1%x_min) / (cpf1%x_max - cpf1%x_min) ELSE fraction = -0.1D0 END IF IF ((fraction > 0.0D0).AND.(fraction < 1.0D0)) THEN wx = x1 + zoom * fraction * (x2 - x1) wy = y1 CALL MOVETO_W (wx, wy, wt) wy = y1 + 0.1D0 * (y2 - y1) result = LINETO_W (wx, wy) END IF ! label the median with a printed value CALL SETGTEXTROTATION(ideg) old_color = SETCOLORRGB(#FFFFFF) ! white result = SETFONT('t''Arial''h12w7') wy = y1 + 0.2D0 * (y2 - y1) CALL MOVETO_W (wx, wy, wt) CALL OUTGTEXT(TRIM(ADJUSTL(ASCII9(REAL(Get_x_from_P(cpf1, 0.025D0)))))) ! add white cross-bar at x:P(X < x) = 0.975 CALL SETLINESTYLE(axis_style) ! solid old_color = SETCOLORRGB(#FFFFFF) ! white x_right = Get_x_from_P (cpf1, 0.975D0) IF (cpf1%x_max > cpf1%x_min) THEN fraction = (x_right - cpf1%x_min) / (cpf1%x_max - cpf1%x_min) ELSE fraction = -0.1D0 END IF IF ((fraction > 0.0D0).AND.(fraction < 1.0D0)) THEN wx = x1 + zoom * fraction * (x2 - x1) wy = y2 CALL MOVETO_W (wx, wy, wt) wy = y1 + 0.9D0 * (y2 - y1) result = LINETO_W (wx, wy) END IF ! label the median with a printed value CALL SETGTEXTROTATION(ideg) old_color = SETCOLORRGB(#FFFFFF) ! white result = SETFONT('t''Arial''h12w7') CALL OUTGTEXT(TRIM(ADJUSTL(ASCII9(REAL(Get_x_from_P(cpf1, 0.975D0)))))) ! draw the green box: CALL SETLINESTYLE(axis_style) ! solid old_color = SETCOLORRGB(#00FF00) ! green wx = x1; wy = y2 CALL MOVETO_W (wx, wy, wt) wy = y1 result = LINETO_W (wx, wy) wx = x2 result = LINETO_W (wx, wy) wy = y2 result = LINETO_W (wx, wy) wx = x1 result = LINETO_W (wx, wy) !if linear scale, with negative values, then add vertical green bar at origin: IF (.NOT.pdf1%logarithmic) THEN IF (pdf1%x_min < 0.0D0) THEN IF (pdf1%x_max > pdf1%x_min) THEN zero_fraction = -pdf1%x_min / (pdf1%x_max - pdf1%x_min) wx = x1 + zero_fraction * (x2 - x1) CALL SETLINESTYLE(axis_style) old_color = SETCOLORRGB(#00FF00) ! green wy = y2 CALL MOVETO_W (wx, wy, wt) wy = y1 result = LINETO_W (wx, wy) END IF END IF END IF !label with name on top: CALL SETGTEXTROTATION(ideg) old_color = SETCOLORRGB(#FFFFFF) ! white result = SETFONT('t''Arial''h18w10') wx = x1; wy = y2 + 0.08D0 CALL MOVETO_W(wx, wy, wt) CALL OUTGTEXT(TRIM(pdf1%name)) !label independent variable values at bottom: c9a = ADJUSTL(ASCII9(REAL(pdf1%x_min))) ! min c9b = ADJUSTL(ASCII9(REAL(Get_Mode(pdf1)))) ! mode c9c = ADJUSTL(ASCII9(REAL(Get_Median(cpf1)))) ! median c9d = ADJUSTL(ASCII9(REAL(Get_Mean(pdf1)))) ! mean c9e = ADJUSTL(ASCII9(REAL(pdf1%x_max))) ! max c9f = ADJUSTL(ASCII9(REAL(Get_Sigma(pdf1)))) ! standard deviation IF (pdf1%logarithmic) THEN text_string = "min = "//TRIM(c9a)//", mode = "//TRIM(c9b)//", MEDIAN = "//TRIM(c9c)//", mean = "//TRIM(c9d)//", max = "//TRIM(c9e)//" (log sampling); S.D. = "//TRIM(c9f) ELSE text_string = "min = "//TRIM(c9a)//", mode = "//TRIM(c9b)//", MEDIAN = "//TRIM(c9c)//", mean = "//TRIM(c9d)//", max = "//TRIM(c9e)//" (linear sampling); S.D. = "//TRIM(c9f) END IF CALL SETGTEXTROTATION(ideg) old_color = SETCOLORRGB(#FFFFFF) ! white result = SETFONT('t''Arial''h16w8') wx = x1; wy = y1 CALL MOVETO_W(wx, wy, wt) CALL OUTGTEXT(TRIM(text_string)) END SUBROUTINE Plot_1_Distribution SUBROUTINE Plot_3_Distributions (pdf1, cpf1, pdf2, cpf2, pdf3,cpf3, pause) ! Fills screen (pre-blanked?) with plots of pdf/cdf of offset, pdf/cdf of age, and pdf/cdf of offset rate. IMPLICIT NONE TYPE (pdf) :: pdf1, pdf2, pdf3 TYPE (cpf) :: cpf1, cpf2, cpf3 LOGICAL, INTENT(IN) :: pause CHARACTER*1 :: c1 !offset graph rectangle occupies 2.05 < y < 2.75, 0.5 < x < 3.5 CALL Plot_1_Distribution(pdf1, cpf1, 0.5D0, 3.5D0, 2.05D0, 2.75D0) !age graph rectangle occupies 1.15 < y < 1.85, 0.5 < x < 3.5 CALL Plot_1_Distribution(pdf2, cpf2, 0.5D0, 3.5D0, 1.15D0, 1.85D0) !offset-rate graph rectangle occupies 0.25 < y < 0.95, 0.5 < x < 3.5 CALL Plot_1_Distribution(pdf3, cpf3, 0.5D0, 3.5D0, 0.25D0, 0.95D0) IF (pause) THEN READ (*, "(A)") c1 IF ((c1 == 'x').OR.(c1 == 'X')) THEN CALL Dump_3_PDFs_for_Excel (pdf1, pdf2, pdf3) END IF END IF END SUBROUTINE Plot_3_Distributions SUBROUTINE Prompt_for_Logical (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a logical value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! The same happens IF (mt_flashby), without waiting for the user. ! Note that prompt_text should usually end with '?'. ! It can be more than 70 bytes long if necessary, but cannot ! contain line breaks, tabs, or other formatting. ! Trailing blanks in prompt_text are ignored. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: prompt_text LOGICAL, INTENT(IN) :: default LOGICAL, INTENT(OUT) :: answer CHARACTER*1 :: inbyte CHARACTER*3 :: yesno INTEGER :: blank_at, bytes, ios, written LOGICAL :: finished bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) IF (default) THEN yesno = 'Yes' ELSE yesno = 'No' END IF written = 0 DO WHILE ((bytes - written) > 70) blank_at = written + INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 70 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(yesno) finished = .TRUE. ! unless changed below READ (*,"(A)") inbyte IF (LEN_TRIM(inbyte) == 0) THEN answer = default ELSE SELECT CASE (inbyte) CASE ('Y') answer = .TRUE. CASE ('y') answer = .TRUE. CASE ('T') answer = .TRUE. CASE ('t') answer = .TRUE. CASE ('R') answer = .TRUE. CASE ('r') answer = .TRUE. CASE ('O') answer = .TRUE. CASE ('o') answer = .TRUE. CASE ('N') answer = .FALSE. CASE ('n') answer = .FALSE. CASE ('F') answer = .FALSE. CASE ('f') answer = .FALSE. CASE ('W') answer = .FALSE. CASE ('w') answer = .FALSE. CASE DEFAULT WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") inbyte WRITE (*,"(' (Only the first letter of your answer is used.)')") WRITE (*,"(' To agree, enter Y, y, T, t, O, o, R, or r.')") WRITE (*,"(' To disagree, enter N, n, F, f, W, or w.')") WRITE (*,"(' Please try again:')") finished = .FALSE. END SELECT END IF ! a byte was entered END DO ! until finished END SUBROUTINE Prompt_for_Logical SUBROUTINE Read_Line_from_Birds_Table (title_line, Findex, offset_type_c1, fault_name, & & editors_offset_km_c9, editors_offset_sigma_km_c9, & & editors_began_Ma_c9, editors_ended_Ma_c9, & & editors_neorate_mmpa_c9, editors_neorate_sigma_mmpa_c9, & & grade, datum_name, & & datum_offset_km_c20, datum_offset_sigma_km_c20, & & datum_began_Ma_c20, datum_ended_Ma_c20, & & datum_neorate_mmpa_c20, datum_neorate_sigma_mmpa_c20, & & F_is_Gaussian, F_is_upper_limit, F_is_lower_limit, F_is_zero, & & sigma_F_stated, F_range_stated, & & Delta_F_nominal, sigma_Delta_F, Delta_F_nominal_high, Delta_F_nominal_low, Delta_F_string, & & Delta_F_max, sigma_Delta_F_max, Delta_F_max_high, Delta_F_max_low, & & Delta_F_min, sigma_Delta_F_min, Delta_F_min_high, Delta_F_min_low, & & age_is_Gaussian, age_has_two_bounds, age_is_lower_limit, age_is_upper_limit, & & sigma_Delta_t_stated, age_range_stated, age_from_raw_14C, age_from_raw_nuclides, & & Delta_t_nominal, sigma_Delta_t, age_nominal_high, age_nominal_low, Delta_t_string, & & sigma_Delta_t_min_stated, Delta_t_min_range_stated, Delta_t_min, sigma_Delta_t_min, & & Delta_t_min_high, Delta_t_min_low, Delta_t_min_string, & & sigma_Delta_t_max_stated, Delta_t_max_range_stated, Delta_t_max, sigma_Delta_t_max, & & Delta_t_max_high, Delta_t_max_low, Delta_t_max_string, & & warning_count, warning_messages) IMPLICIT NONE LOGICAL, INTENT(INOUT) :: title_line CHARACTER(5), INTENT(INOUT) :: Findex CHARACTER(1), INTENT(INOUT) :: offset_type_c1 CHARACTER(50), INTENT(INOUT) :: fault_name CHARACTER(9), INTENT(INOUT) :: editors_offset_km_c9, editors_offset_sigma_km_c9, & & editors_began_Ma_c9, editors_ended_Ma_c9,& & editors_neorate_mmpa_c9, editors_neorate_sigma_mmpa_c9 ! -------------------------------------------------------- CHARACTER(1), INTENT(OUT) :: grade CHARACTER(50), INTENT(OUT) :: datum_name CHARACTER(20), INTENT(OUT) :: datum_offset_km_c20, datum_offset_sigma_km_c20, & & datum_began_Ma_c20, datum_ended_Ma_c20,& & datum_neorate_mmpa_c20, datum_neorate_sigma_mmpa_c20 ! -------------------------------------------------------- LOGICAL, INTENT(OUT) :: F_is_Gaussian, F_is_upper_limit, F_is_lower_limit, F_is_zero LOGICAL, INTENT(OUT) :: sigma_F_stated, F_range_stated DOUBLE PRECISION, INTENT(OUT) :: Delta_F_nominal, sigma_Delta_F, Delta_F_nominal_high, Delta_F_nominal_low CHARACTER(20), INTENT(OUT) :: Delta_F_string LOGICAL, INTENT(OUT) :: sigma_Delta_t_max_stated, Delta_t_max_range_stated DOUBLE PRECISION, INTENT(OUT) :: Delta_F_max, sigma_Delta_F_max, Delta_F_max_high, Delta_F_max_low DOUBLE PRECISION, INTENT(OUT) :: Delta_F_min, sigma_Delta_F_min, Delta_F_min_high, Delta_F_min_low ! -------------------------------------------------------- LOGICAL, INTENT(OUT) :: age_is_Gaussian, age_has_two_bounds, age_is_lower_limit, age_is_upper_limit LOGICAL, INTENT(OUT) :: sigma_Delta_t_stated, age_range_stated, age_from_raw_14C, age_from_raw_nuclides DOUBLE PRECISION, INTENT(OUT) :: Delta_t_nominal, sigma_Delta_t, age_nominal_high, age_nominal_low CHARACTER(20), INTENT(OUT) :: Delta_t_string LOGICAL, INTENT(OUT) :: sigma_Delta_t_min_stated, Delta_t_min_range_stated DOUBLE PRECISION, INTENT(OUT) :: Delta_t_min, sigma_Delta_t_min DOUBLE PRECISION, INTENT(OUT) :: Delta_t_min_high, Delta_t_min_low CHARACTER(20), INTENT(OUT) :: Delta_t_min_string DOUBLE PRECISION, INTENT(OUT) :: Delta_t_max, sigma_Delta_t_max DOUBLE PRECISION, INTENT(OUT) :: Delta_t_max_high, Delta_t_max_low CHARACTER(20), INTENT(OUT) :: Delta_t_max_string ! -------------------------------------------------------- INTEGER, INTENT(OUT) :: warning_count ! >0 indicates non-routine return; EOF or missing data; no computation. CHARACTER(120), DIMENSION(2), INTENT(OUT) :: warning_messages ! (1) should begin with "EOF" if end-of-file found on device 1. CHARACTER(1), PARAMETER :: tab = CHAR(9) ! special "HT" tab character in ASCII sequence CHARACTER(1), PARAMETER :: quote = CHAR(34) ! special double-quotation-mark character in ASCII sequence CHARACTER(1) :: c1, left_c1, right_c1 CHARACTER(6) :: ios_c6 CHARACTER(50) :: left_c50, right_c50, temp_c50 CHARACTER(50), DIMENSION(8) :: column CHARACTER(500) :: line DOUBLE PRECISION, PARAMETER :: zero = 0.0D0 DOUBLE PRECISION, PARAMETER :: one = 1.0D0 DOUBLE PRECISION, PARAMETER :: D_dip_degrees = 20.0D0, & & L_dip_degrees = 73.0D0, & & N_dip_degrees = 55.0D0, & & P_dip_degrees = 20.0D0, & & R_dip_degrees = 73.0D0, & & T_dip_degrees = 20.0D0 DOUBLE PRECISION :: datum_ended_Ma, datum_ended_Ma_2, left_value, low_value, high_value, offset_magnifier, right_value, value DOUBLE PRECISION :: D_dip_radians, L_dip_radians, N_dip_radians, P_dip_radians, R_dip_radians, T_dip_radians INTEGER :: i, ios, j, n, n_right, n_tab, n1, n2 LOGICAL :: braces, paleotectonic D_dip_radians = D_dip_degrees / 57.29577951D0 L_dip_radians = L_dip_degrees / 57.29577951D0 N_dip_radians = N_dip_degrees / 57.29577951D0 P_dip_radians = P_dip_degrees / 57.29577951D0 R_dip_radians = R_dip_degrees / 57.29577951D0 T_dip_radians = T_dip_degrees / 57.29577951D0 ! N.B. Fault-title-line info is NOT initialized, so it can remain unchanged during a second CALL! ! All values associated with an offset datum are initialized to blank / zero / .FALSE.. grade = '' datum_name = '' datum_offset_km_c20 = '' datum_offset_sigma_km_c20 = '' datum_began_Ma_c20 = '' datum_ended_Ma_c20 = '' datum_neorate_mmpa_c20 = '' datum_neorate_sigma_mmpa_c20 = '' F_is_Gaussian = .FALSE. F_is_upper_limit = .FALSE. F_is_lower_limit = .FALSE. F_is_zero = .FALSE. sigma_F_stated = .FALSE. F_range_stated = .FALSE. Delta_F_nominal = zero sigma_Delta_F = zero Delta_F_nominal_high = zero Delta_F_nominal_low = zero Delta_F_string = '' sigma_Delta_t_max_stated = .FALSE. Delta_t_max_range_stated = .FALSE. Delta_F_max = zero sigma_Delta_F_max = zero Delta_F_max_high = zero Delta_F_max_low = zero Delta_F_min = zero sigma_Delta_F_min = zero Delta_F_min_high = zero Delta_F_min_low = zero age_is_Gaussian = .FALSE. age_has_two_bounds = .FALSE. age_is_lower_limit = .FALSE. age_is_upper_limit = .FALSE. sigma_Delta_t_stated = .FALSE. age_range_stated = .FALSE. age_from_raw_14C = .FALSE. age_from_raw_nuclides = .FALSE. Delta_t_nominal = zero sigma_Delta_t = zero age_nominal_high = zero age_nominal_low = zero Delta_t_string = '' sigma_Delta_t_min_stated = .FALSE. Delta_t_min_range_stated = .FALSE. Delta_t_min = zero sigma_Delta_t_min = zero Delta_t_min_high = zero Delta_t_min_low = zero Delta_t_min_string = '' Delta_t_max = zero sigma_Delta_t_max = zero Delta_t_max_high = zero Delta_t_max_low = zero Delta_t_max_string = '' warning_count = 0 ! unless changed below! warning_messages = '' ! each of 20 strings will be set to the null string READ (1, "(A)", IOSTAT = ios) line ! Unit 1 has already been OPENed for tab-delimited ASCII flat-file ! with one fault-name or one offset feature per line. (No other header lines.) IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) IF (ios == -1) THEN warning_messages(warning_count) = "EOF" ELSE WRITE (ios_c6, "(I6)") ios warning_messages(warning_count) = "ios = " // TRIM(ADJUSTL(ios_c6)) // " for line: " // TRIM(line) END IF RETURN ! to main program END IF get_columns: DO i = 1, 8 IF (LEN_TRIM(line) > 0) THEN n_tab = SCAN(line, tab) IF (n_tab > 0) THEN IF (n_tab > 1) THEN n_right = MIN(50, (n_tab - 1)) column(i) = line(1:n_right) !strip off quotation marks, if any DO j = 1, LEN_TRIM(column(i)) IF (column(i)(j:j) == quote) THEN column(i)(j:50) = column(i)((j+1):50) // ' ' END IF END DO line = ADJUSTL(line((n_tab + 1):500)) ELSE ! n_tab = 1 column(i) = '' line = ADJUSTL(line(2:500)) END IF ELSE ! n_tab = 0; no tab left in line column(i) = TRIM(line) line = '' END IF ELSE column(i) = '' END IF END DO get_columns c1 = column(1)(1:1) title_line = ((c1 == 'F').OR.(c1 == 'f')) IF (title_line) THEN Findex = column(1)(1:5) offset_type_c1 = column(1)(6:6) fault_name = TRIM(column(2)) editors_offset_km_c9 = TRIM(column(3)) editors_offset_sigma_km_c9 = TRIM(column(4)) editors_began_Ma_c9 = TRIM(column(5)) editors_ended_Ma_c9 = TRIM(column(6)) editors_neorate_mmpa_c9 = TRIM(column(7)) editors_neorate_sigma_mmpa_c9 = TRIM(column(8)) !N.B. All are returned left-justified. ELSE ! a datum line grade = TRIM(column(1)) datum_name = TRIM(column(2)) datum_offset_km_c20 = TRIM(column(3)) datum_offset_sigma_km_c20 = TRIM(column(4)) datum_began_Ma_c20 = TRIM(column(5)) datum_ended_Ma_c20 = TRIM(column(6)) datum_neorate_mmpa_c20 = TRIM(column(7)) datum_neorate_sigma_mmpa_c20 = TRIM(column(8)) !N.B. All are returned left-justified. !clean columns #3~8 of unwanted characters: '?'; '('; ')'; '-' --> '~'; any '~' NOT between digits !BUT still leave any '>' and '<' and '{' and '}' !!! DO i = 3, 8 !eliminate '?' n = SCAN(column(i), '?') IF (n > 0) THEN IF (n > 1) THEN ! usual case column(i) = column(i)(1:(n-1)) // column(i)((n+1):50) // ' ' ELSE ! leading '?'; not expected column(i) = column(i)(2:50) // ' ' END IF END IF !eliminate '?' again, just to be sure n = SCAN(column(i), '?') IF (n > 0) THEN IF (n > 1) THEN ! usual case column(i) = column(i)(1:(n-1)) // column(i)((n+1):50) // ' ' ELSE ! leading '?'; not expected column(i) = column(i)(2:50) // ' ' END IF END IF !eliminate '(' , as in ">(0.3-2) in many entries attributed to Stewart, 1998" n = SCAN(column(i), '(') IF (n > 0) THEN IF (n > 1) THEN column(i) = column(i)(1:(n-1)) // column(i)((n+1):50) // ' ' ELSE ! column(i) = column(i)(2:50) // ' ' END IF END IF !eliminate ')' , as in ">(0.3-2) in many entries attributed to Stewart, 1998" n = SCAN(column(i), ')') IF (n > 0) THEN IF (n > 1) THEN column(i) = column(i)(1:(n-1)) // column(i)((n+1):50) // ' ' ELSE ! column(i) = column(i)(2:50) // ' ' END IF END IF !change any "-" (which could be confused with negative sign) to '~' DO j = 1, LEN_TRIM(column(i)) IF (column(i)(j:j) == '-') column(i)(j:j) = '~' END DO !eliminate any '~' in leading position, or anywhere NOT surrounded by digits (e.g, ">~4"): n = SCAN(column(i), '~') IF (n > 0) THEN IF (n == 1) THEN column(i) = column(i)(2:50) // ' ' ELSE ! n > 1; embedded '~' left_c1 = column(i)((n-1):(n-1)) right_c1 = column(i)((n+1):(n+1)) IF ((.NOT.Digit(left_c1)).OR.(.NOT.Digit(right_c1))) THEN column(i) = column(i)(1:(n-1)) // column(i)((n+1):50) // ' ' END IF END IF END IF !repeat elimination of any '~' in leading position, or anywhere NOT surrounded by digits (e.g, ">~4"): n = SCAN(column(i), '~') IF (n > 0) THEN IF (n == 1) THEN column(i) = column(i)(2:50) // ' ' ELSE ! n > 1; embedded '~' left_c1 = column(i)((n-1):(n-1)) right_c1 = column(i)((n+1):(n+1)) IF ((.NOT.Digit(left_c1)).OR.(.NOT.Digit(right_c1))) THEN column(i) = column(i)(1:(n-1)) // column(i)((n+1):50) // ' ' END IF END IF END IF END DO ! cleaning columns #3~8 of unwanted characters (but leaving >, <, {, }, ~ in cases like 5~6). !Test for paleotectonic datum == positive ending time in column(6) !and re-cast as datum with zero displacement since age of the overlap assemblage! !Also test for blank ending time, combined with editor's positive ending time, and reject the datum. IF (LEN_TRIM(column(6)) == 0) THEN ! if this field empty, defer to editor's opinion about ending time. IF (LEN_TRIM(editors_ended_Ma_c9) == 0) THEN ! editor has no opinion, either. paleotectonic = .FALSE. ! assume fault could still be active. ELSE ! editor did express an opinion: READ (editors_ended_Ma_c9, *, IOSTAT = ios) value IF (ios /= 0) THEN WRITE (*, "(' ERROR: Illegible editors_ended_Ma_c9 = ',A)") TRIM(editors_ended_Ma_c9) CALL Pause() STOP END IF paleotectonic = (value > 0.0D0) IF (paleotectonic) THEN ! skip this datum, which does not contribute any information warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "No ending time, and editor considers fault to be inactive." RETURN END IF END IF ELSE IF (SCAN(column(6), '<') > 0) THEN paleotectonic = .FALSE. ! only an upper limit is available on the ending time; could be zero. ELSE IF (SCAN(column(6), '>') > 0) THEN n = SCAN(column(6), '>') IF (n == 1) THEN temp_c50 = column(6)(2:50) ELSE temp_c50 = column(6)(1:(n-1)) // column(6)((n+1):50) // ' ' END IF READ (temp_c50, *, IOSTAT = ios) value IF (ios /= 0) THEN WRITE (*, "(' ERROR: Illegible datum_ended_Ma = ',A)") TRIM(column(6)) CALL Pause() STOP ELSE IF (value > 0.0D0) THEN ! meaningful constraint; not a case of " rate < 0 / (>0)" paleotectonic = .TRUE. ! ending time is known to be > 0. age_is_lower_limit = .TRUE. ! Sic! Not a typo! Will use end time for displacement as start time for no-displacement. Delta_t_min_string = TRIM(temp_c50) Delta_t_min = value END IF END IF ELSE ! non-blank end time which is neither upper or lower limit; must READ the value(s) to decide on fault activity: n = SCAN(column(6), '~') IF (n == 0) THEN ! only one #; try to read: READ (column(6), *, IOSTAT = ios) datum_ended_Ma IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable number in ending-time column: " // TRIM(column(6)) RETURN ELSE ! READ succeeded paleotectonic = (datum_ended_Ma > zero) IF (paleotectonic) THEN ! Will use the end of displacement as the start of non-displacement: age_is_Gaussian = .TRUE. Delta_t_nominal = datum_ended_Ma Delta_t_string = TRIM(column(6)) ! can be used to infer sigma END IF END IF ! READ failed or succeeded ELSE ! range of ages for ending time; try to read two values !2nd value is read first, because its is needed to decide whether paleotectonic: right_c50 = column(6)((n+1):50) READ (right_c50, *, IOSTAT = ios) datum_ended_Ma_2 IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable number in ending-time column, after ~: " // TRIM(temp_c50) RETURN ELSE ! 1st READ, of 2nd number in ending-time, succeeded paleotectonic = (datum_ended_Ma_2 > zero) IF (paleotectonic) THEN ! Will use end of displacement as start of non-displacement: !Get first number in pair: left_c50 = column(6)(1:(n-1)) READ (left_c50, *, IOSTAT = ios) datum_ended_Ma IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable number in ending-time column, before ~: " // TRIM(temp_c50) RETURN ELSE ! 2nd READ, of first number in pair of ending-times, also succeeded age_has_two_bounds = .TRUE. Delta_t_max = MAX(datum_ended_Ma, datum_ended_Ma_2) Delta_t_min = MIN(datum_ended_Ma, datum_ended_Ma_2) IF (Delta_t_max == datum_ended_Ma) THEN ! e.g., "5.2~3" Delta_t_max_string = TRIM(left_c50) Delta_t_min_string = TRIM(right_c50) ELSE ! strings need to be swapped; e.g., ending time was "3~5.2" Delta_t_max_string = TRIM(right_c50) Delta_t_min_string = TRIM(left_c50) END IF IF (Delta_t_min <= 0.0D0) THEN ! not necessarily paleotectonic at all; correct mistakes: paleotectonic = .FALSE. age_has_two_bounds = .FALSE. END IF END IF ! 2nd READ, of 1st number in ending-time, failed or succeeded END IF ! 2nd number in ending-time implied paleotectonic END IF ! 1st READ, of 2nd number in ending-time, failed or succeeded END IF ! single value, or range of values for end time END IF IF (paleotectonic) THEN F_is_zero = .TRUE. !Note that overlap-assemblage age and type-of-age were set in logic above, !along each of the branches which contained " paleotectonic = .TRUE. " RETURN END IF !From here on, assume that the datum is neotectonic, so ending_Ma = zero: !Search for '>' or '<' in offset, set logical flag, then strip it out: IF (SCAN(column(3), '>') > 0) THEN F_is_lower_limit = .TRUE. n = SCAN(column(3), '>') ! > 0 IF (n == 1) THEN column(3) = column(3)(2:50) // ' ' ELSE ! n>1 column(3) = column(3)(1:(n-1)) // column(3)((n+1):50) // ' ' END IF ELSE IF (SCAN(column(3), '<') > 0) THEN F_is_upper_limit = .TRUE. n = SCAN(column(3), '<') ! > 0 IF (n == 1) THEN column(3) = column(3)(2:50) // ' ' ELSE ! n>1 column(3) = column(3)(1:(n-1)) // column(3)((n+1):50) // ' ' END IF END IF !Offset may look like any of these: "5", "5~6", "{5}", "{5~6}", "onstrained", " ". !Look for braces, {}, note the fact, set offset_magnifier, and then strip them out: n1 = SCAN(column(3), '{') n2 = SCAN(column(3), '}') braces = ((n1 > 0).AND.(n2>0)) IF (braces) THEN IF (offset_type_c1 == 'D') THEN !{throw} --> opening offset_magnifier = one / DTAN(D_dip_radians) ELSE IF (offset_type_c1 == 'L') THEN !{throw} --> sinistral strike-slip offset_magnifier = one F_is_lower_limit = .TRUE. ELSE IF (offset_type_c1 == 'N') THEN !{slip} --> throw offset_magnifier = DSIN(N_dip_radians) ELSE IF (offset_type_c1 == 'P') THEN !{throw} --> convergence offset_magnifier = one / DTAN(P_dip_radians) ELSE IF (offset_type_c1 == 'R') THEN !{throw} --> dextral strike-slip offset_magnifier = one F_is_lower_limit = .TRUE. ELSE IF (offset_type_c1 == 'T') THEN !{convergence} --> throw offset_magnifier = DTAN(T_dip_radians) END IF !strip out '{': IF (n1 == 1) THEN column(3) = column(3)(2:50) // ' ' ELSE ! n1>1 column(3) = column(3)(1:(n1-1)) // column(3)((n1+1):50) // ' ' END IF !strip out '}', recognizing that it has probably moved (so re-evaluate n2): n2 = SCAN(column(3), '}') IF (n2 == 1) THEN column(3) = column(3)(2:50) // ' ' ELSE ! n2>1 column(3) = column(3)(1:(n2-1)) // column(3)((n2+1):50) // ' ' END IF ELSE ! no braces offset_magnifier = one END IF ! braces, or not !Can now assume that braces are gone (but remember to apply "* offset_magnifier"). !Offset may look like any of these: "5", "5~6", "onstrained", " ". !Check for blank offset: IF (LEN_TRIM(column(3)) == 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Offset is blank." RETURN END IF !Offset may look like any of these: "5", "5~6", "onstrained". !If '~' is present in offset, try to read 2 offset values surrounding it: IF (SCAN(column(3), '~') > 0) THEN n = SCAN(column(3), '~') !Read left value: temp_c50 = column(3)(1:(n-1)) READ (temp_c50, *, IOSTAT = ios) left_value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable offset, before ~: " // TRIM(column(3)) RETURN END IF !Read right value: temp_c50 = column(3)((n+1):50) READ (temp_c50, *, IOSTAT = ios) right_value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable offset, after ~: " // TRIM(column(3)) RETURN END IF low_value = MIN(left_value, right_value) high_value = MAX(left_value, right_value) !********** GOT RANGE FOR OFFSET ************ F_range_stated = .TRUE. IF (F_is_lower_limit) THEN Delta_F_min_low = low_value * offset_magnifier Delta_F_min_high = high_value * offset_magnifier ELSE IF (F_is_upper_limit) THEN Delta_F_max_low = low_value * offset_magnifier Delta_F_max_high = high_value * offset_magnifier ELSE F_is_Gaussian = .TRUE. Delta_F_nominal_low = low_value * offset_magnifier Delta_F_nominal_high = high_value * offset_magnifier END IF !******************************************** ELSE ! no '~'; single value, or garbage: READ (column(3), *, IOSTAT = ios) value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable offset: " // TRIM(column(3)) RETURN END IF !********** GOT SINGLE-VALUED OFFSET ******** F_is_zero = (value == zero) IF (F_is_zero) THEN F_is_lower_limit = .FALSE. F_is_upper_limit = .FALSE. F_is_Gaussian = .FALSE. ELSE ! F is not zero: IF (F_is_lower_limit) THEN Delta_F_min = value * offset_magnifier ELSE IF (F_is_upper_limit) THEN Delta_F_max = value * offset_magnifier ELSE F_is_Gaussian = .TRUE. ! by elimination Delta_F_nominal = value * offset_magnifier END IF END IF Delta_F_string = TRIM(column(3)) !******************************************** END IF ! 2 offset values, or just one (or garbage) !Check whether sigma was provided for offset: IF (LEN_TRIM(column(4)) > 0) THEN READ (column(4), *, IOSTAT = ios) value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable sigma of offset: " // TRIM(column(4)) RETURN END IF sigma_F_stated = .TRUE. IF (F_is_Gaussian) THEN sigma_Delta_F = value * offset_magnifier ELSE IF (F_is_upper_limit) THEN sigma_Delta_F_max = value * offset_magnifier ELSE IF (F_is_lower_limit) THEN sigma_Delta_F_min = value * offset_magnifier ELSE IF (F_is_zero) THEN sigma_Delta_F = value * offset_magnifier END IF END IF ! something found in offset_sigma_km column ! = = = = = = = end of offset and its sigma; start of begin_Ma !!! ! begin_Ma may look like any of these: " ", "late Q.", "5", "6~5", "5~6", "<6", ">5". !Check for blank begin_Ma: IF (LEN_TRIM(column(5)) == 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Begin_Ma is blank." RETURN END IF ! begin_Ma may look like any of these: "late Q.", "5", "6~5", "5~6", "<6", ">5". !Search for '>' or '<' in Begin_Ma, set logical flag, then strip it out: IF (SCAN(column(5), '>') > 0) THEN age_is_lower_limit = .TRUE. n = SCAN(column(5), '>') ! > 0 IF (n == 1) THEN column(5) = column(5)(2:50) // ' ' ELSE ! n>1 column(5) = column(5)(1:(n-1)) // column(5)((n+1):50) // ' ' END IF ELSE IF (SCAN(column(5), '<') > 0) THEN age_is_upper_limit = .TRUE. n = SCAN(column(5), '<') ! > 0 IF (n == 1) THEN column(5) = column(5)(2:50) // ' ' ELSE ! n>1 column(5) = column(5)(1:(n-1)) // column(5)((n+1):50) // ' ' END IF END IF ! begin_Ma may look like any of these: "late Q.", "5", "6~5", "5~6". !If '~' is present in Begin_Ma, try to read 2 Begin_Ma values surrounding it: IF (SCAN(column(5), '~') > 0) THEN n = SCAN(column(5), '~') !Read left value (usually higher): left_c50 = column(5)(1:(n-1)) READ (left_c50, *, IOSTAT = ios) left_value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable Begin_Ma, before ~: " // TRIM(column(5)) RETURN END IF !Read right value (usually lower): right_c50 = column(5)((n+1):50) READ (right_c50, *, IOSTAT = ios) right_value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable Begin_Ma, after ~: " // TRIM(column(5)) RETURN END IF low_value = MIN(left_value, right_value) high_value = MAX(left_value, right_value) !high value should be on the left, and low value on right; if not, then swap text strings: IF (.NOT.(high_value == left_value)) THEN temp_c50 = left_c50 left_c50 = right_c50 right_c50 = temp_c50 END IF !********** GOT RANGE FOR Begin_Ma ************ IF (age_is_lower_limit) THEN Delta_t_min_range_stated = .TRUE. Delta_t_min_low = low_value Delta_t_min_high = high_value ELSE IF (age_is_upper_limit) THEN Delta_t_max_range_stated = .TRUE. Delta_t_max_low = low_value Delta_t_max_high = high_value ELSE age_has_two_bounds = .TRUE. Delta_t_max = high_value Delta_t_max_string = TRIM(left_c50) Delta_t_min = low_value Delta_t_min_string = TRIM(right_c50) END IF !N.B. Calling program, Slippery.f90, currently makes ! no use of age_from_raw_14C or age_from_raw_nuclides ! in this case, so no values are set. (Default = .FALSE.) !******************************************** ELSE ! no '~'; single value, or garbage: READ (column(5), *, IOSTAT = ios) value IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Unreadable Begin_Ma: " // TRIM(column(5)) RETURN END IF !******** GOT SINGLE-VALUED Begin_Ma ******** IF (age_is_lower_limit) THEN Delta_t_min = value Delta_t_min_string = TRIM(column(5)) ELSE IF (age_is_upper_limit) THEN Delta_t_max = value Delta_t_max_string = TRIM(column(5)) ELSE age_is_Gaussian = .TRUE. Delta_t_nominal = value Delta_t_string = TRIM(column(5)) END IF IF (value <= 0.03) THEN ! <= 30,000 years; see p. 126 of Yeats, Sieh, & Allen [1997] The Geology of Earthquakes. age_from_raw_14C = .TRUE. ELSE IF (value <= 0.6) THEN ! <= 600,000 years {rather arbitrary, but younger than Bishop tuff = 0.768 Ma} age_from_raw_nuclides = .TRUE. END IF !******************************************** END IF ! 2 Begin_Ma values, or just one (or garbage) END IF ! title line, or datum line? END SUBROUTINE Read_Line_from_Birds_Table SUBROUTINE Read_Line_from_PaleoSites_Table (old_Findex, old_offset_type_c1, & & new_group, Findex, offset_type_c1, fault_name, & & editors_offset_km_c9, editors_offset_sigma_km_c9, & & editors_began_Ma_c9, editors_ended_Ma_c9, & & editors_neorate_mmpa_c9, editors_neorate_sigma_mmpa_c9, & & source_letter, grade, datum_name, & & datum_offset_km_c20, datum_offset_sigma_km_c20, & & datum_began_Ma_c20, datum_ended_Ma_c20, & & datum_neorate_mmpa_c20, datum_neorate_sigma_mmpa_c20, & & F_is_Gaussian, F_is_upper_limit, F_is_lower_limit, F_is_zero, & & sigma_F_stated, F_range_stated, & & Delta_F_nominal, sigma_Delta_F, Delta_F_nominal_high, Delta_F_nominal_low, Delta_F_string, & & Delta_F_max, sigma_Delta_F_max, Delta_F_max_high, Delta_F_max_low, & & Delta_F_min, sigma_Delta_F_min, Delta_F_min_high, Delta_F_min_low, & & age_is_Gaussian, age_has_two_bounds, age_is_lower_limit, age_is_upper_limit, & & sigma_Delta_t_stated, age_range_stated, age_from_raw_14C, age_from_raw_nuclides, & & Delta_t_nominal, sigma_Delta_t, age_nominal_high, age_nominal_low, Delta_t_string, & & sigma_Delta_t_min_stated, Delta_t_min_range_stated, Delta_t_min, sigma_Delta_t_min, & & Delta_t_min_high, Delta_t_min_low, Delta_t_min_string, & & sigma_Delta_t_max_stated, Delta_t_max_range_stated, Delta_t_max, sigma_Delta_t_max, & & Delta_t_max_high, Delta_t_max_low, Delta_t_max_string, & & warning_count, warning_messages) IMPLICIT NONE LOGICAL, INTENT(INOUT) :: new_group CHARACTER(5), INTENT(INOUT) :: Findex, old_Findex CHARACTER(1), INTENT(INOUT) :: offset_type_c1, old_offset_type_c1 CHARACTER(50), INTENT(INOUT) :: fault_name CHARACTER(9), INTENT(INOUT) :: editors_offset_km_c9, editors_offset_sigma_km_c9, & & editors_began_Ma_c9, editors_ended_Ma_c9,& & editors_neorate_mmpa_c9, editors_neorate_sigma_mmpa_c9 ! -------------------------------------------------------- CHARACTER(1), INTENT(OUT) :: source_letter, grade CHARACTER(50), INTENT(OUT) :: datum_name CHARACTER(20), INTENT(OUT) :: datum_offset_km_c20, datum_offset_sigma_km_c20, & & datum_began_Ma_c20, datum_ended_Ma_c20,& & datum_neorate_mmpa_c20, datum_neorate_sigma_mmpa_c20 ! -------------------------------------------------------- LOGICAL, INTENT(OUT) :: F_is_Gaussian, F_is_upper_limit, F_is_lower_limit, F_is_zero LOGICAL, INTENT(OUT) :: sigma_F_stated, F_range_stated DOUBLE PRECISION, INTENT(OUT) :: Delta_F_nominal, sigma_Delta_F, Delta_F_nominal_high, Delta_F_nominal_low CHARACTER(20), INTENT(OUT) :: Delta_F_string LOGICAL, INTENT(OUT) :: sigma_Delta_t_max_stated, Delta_t_max_range_stated DOUBLE PRECISION, INTENT(OUT) :: Delta_F_max, sigma_Delta_F_max, Delta_F_max_high, Delta_F_max_low DOUBLE PRECISION, INTENT(OUT) :: Delta_F_min, sigma_Delta_F_min, Delta_F_min_high, Delta_F_min_low ! -------------------------------------------------------- LOGICAL, INTENT(OUT) :: age_is_Gaussian, age_has_two_bounds, age_is_lower_limit, age_is_upper_limit LOGICAL, INTENT(OUT) :: sigma_Delta_t_stated, age_range_stated, age_from_raw_14C, age_from_raw_nuclides DOUBLE PRECISION, INTENT(OUT) :: Delta_t_nominal, sigma_Delta_t, age_nominal_high, age_nominal_low CHARACTER(20), INTENT(OUT) :: Delta_t_string LOGICAL, INTENT(OUT) :: sigma_Delta_t_min_stated, Delta_t_min_range_stated DOUBLE PRECISION, INTENT(OUT) :: Delta_t_min, sigma_Delta_t_min DOUBLE PRECISION, INTENT(OUT) :: Delta_t_min_high, Delta_t_min_low CHARACTER(20), INTENT(OUT) :: Delta_t_min_string DOUBLE PRECISION, INTENT(OUT) :: Delta_t_max, sigma_Delta_t_max DOUBLE PRECISION, INTENT(OUT) :: Delta_t_max_high, Delta_t_max_low CHARACTER(20), INTENT(OUT) :: Delta_t_max_string ! -------------------------------------------------------- INTEGER, INTENT(OUT) :: warning_count ! >0 indicates non-routine return; EOF or missing data; no computation. CHARACTER(120), DIMENSION(20), INTENT(OUT) :: warning_messages ! (1) should begin with "EOF" if end-of-file found on device 1. CHARACTER(1), PARAMETER :: tab = CHAR(9) ! special "HT" tab character in ASCII sequence CHARACTER(1), PARAMETER :: quote = CHAR(34) ! special double-quotation-mark character in ASCII sequence CHARACTER(1) :: c1, left_c1, right_c1 CHARACTER(2) :: dominant_sense_c2, secondary_sense_c2 CHARACTER(6) :: ios_c6 CHARACTER(9) :: max_age_Ma_c9, max_offset_km_c9, min_age_Ma_c9, min_offset_km_c9, nominal_age_Ma_c9, nominal_offset_km_c9, & & max_end_time_Ma_c9, min_end_time_Ma_c9, nominal_end_time_Ma_c9 CHARACTER(50) :: left_c50, right_c50, temp_c50 CHARACTER(50), DIMENSION(51) :: column ! A ~ AY = 51 columns (not all used here) CHARACTER(500) :: line ! As of 2007.01.31, average line was 376 bytes. Here I round up. PAD = "YES" should be specified in the OPEN. DOUBLE PRECISION, PARAMETER :: zero = 0.0D0 DOUBLE PRECISION, PARAMETER :: one = 1.0D0 DOUBLE PRECISION, PARAMETER :: D_dip_degrees = 20.0D0, & & L_dip_degrees = 73.0D0, & & N_dip_degrees = 55.0D0, & & P_dip_degrees = 20.0D0, & & R_dip_degrees = 73.0D0, & & T_dip_degrees = 20.0D0 DOUBLE PRECISION :: datum_ended_Ma, datum_ended_Ma_2, left_value, low_value, high_value, & & offset_magnifier, rake_in_radians, right_value, value DOUBLE PRECISION :: D_dip_radians, L_dip_radians, N_dip_radians, P_dip_radians, R_dip_radians, T_dip_radians INTEGER :: i, ios, j, n, n_hyphen, n_right, n_tab, n1, n2, rake_in_degrees LOGICAL :: braces, neotectonic, paleotectonic REAL :: max_end_time, max_end_time_Ma, min_end_time, min_end_time_Ma, & & max_offset_km, max_offset_m, min_offset_km, min_offset_m, & & max_age, max_age_Ma, min_age, min_age_Ma, & & nominal_age, nominal_age_Ma, nominal_end_time, nominal_end_time_Ma, & & nominal_offset_m, nominal_offset_km, & & pivot_time_Ma D_dip_radians = D_dip_degrees / 57.29577951D0 L_dip_radians = L_dip_degrees / 57.29577951D0 N_dip_radians = N_dip_degrees / 57.29577951D0 P_dip_radians = P_dip_degrees / 57.29577951D0 R_dip_radians = R_dip_degrees / 57.29577951D0 T_dip_radians = T_dip_degrees / 57.29577951D0 ! All values associated with the new offset datum are initialized to blank / zero / .FALSE., to aid in debugging, ! EXCEPT for offset_type_c1, which "rides" at its own value until a new valid offset type is read. ! This prevents artificially breaking the grouping of offsets with the same type on the same fault section, ! even though a defective line may intervene in the PaleoSites database. ! Instead, offset_type_c1 is only initialized to '' when a new fault section is detected. new_group = .FALSE. Findex = '' fault_name = '' editors_offset_km_c9 = '' editors_offset_sigma_km_c9 = '' editors_began_Ma_c9 = '' editors_ended_Ma_c9 = '' editors_neorate_mmpa_c9 = '' editors_neorate_sigma_mmpa_c9 = '' source_letter = '' datum_name = '' datum_offset_km_c20 = '' datum_offset_sigma_km_c20 = '' datum_began_Ma_c20 = '' datum_ended_Ma_c20 = '' datum_neorate_mmpa_c20 = '' datum_neorate_sigma_mmpa_c20 = '' F_is_Gaussian = .FALSE. F_is_upper_limit = .FALSE. F_is_lower_limit = .FALSE. F_is_zero = .FALSE. sigma_F_stated = .FALSE. F_range_stated = .FALSE. Delta_F_nominal = zero sigma_Delta_F = zero Delta_F_nominal_high = zero Delta_F_nominal_low = zero Delta_F_string = '' sigma_Delta_t_max_stated = .FALSE. Delta_t_max_range_stated = .FALSE. Delta_F_max = zero sigma_Delta_F_max = zero Delta_F_max_high = zero Delta_F_max_low = zero Delta_F_min = zero sigma_Delta_F_min = zero Delta_F_min_high = zero Delta_F_min_low = zero age_is_Gaussian = .FALSE. age_has_two_bounds = .FALSE. age_is_lower_limit = .FALSE. age_is_upper_limit = .FALSE. sigma_Delta_t_stated = .FALSE. age_range_stated = .FALSE. age_from_raw_14C = .FALSE. age_from_raw_nuclides = .FALSE. Delta_t_nominal = zero sigma_Delta_t = zero age_nominal_high = zero age_nominal_low = zero Delta_t_string = '' sigma_Delta_t_min_stated = .FALSE. Delta_t_min_range_stated = .FALSE. Delta_t_min = zero sigma_Delta_t_min = zero Delta_t_min_high = zero Delta_t_min_low = zero Delta_t_min_string = '' Delta_t_max = zero sigma_Delta_t_max = zero Delta_t_max_high = zero Delta_t_max_low = zero Delta_t_max_string = '' warning_count = 0 ! unless changed below! warning_messages = '' ! entire list of 20 strings will be set to the null string READ (1, "(A)", IOSTAT = ios) line ! Unit 1 has already been OPENed for tab-delimited ASCII flat-file ! with one offset feature per line. (Header line has already been read, in Main.) IF (ios /= 0) THEN warning_count = MIN(20, warning_count + 1) IF (ios == -1) THEN warning_messages(warning_count) = "EOF" ELSE WRITE (ios_c6, "(I6)") ios warning_messages(warning_count) = "ios = " // TRIM(ADJUSTL(ios_c6)) // " for line: " // TRIM(line) END IF RETURN ! to main program END IF get_columns: DO i = 1, 51 IF (LEN_TRIM(line) > 0) THEN n_tab = SCAN(line, tab) IF (n_tab > 0) THEN IF (n_tab > 1) THEN n_right = MIN(50, (n_tab - 1)) column(i) = line(1:n_right) !strip off quotation marks, if any DO j = 1, LEN_TRIM(column(i)) IF (column(i)(j:j) == quote) THEN column(i)(j:50) = column(i)((j+1):50) // ' ' END IF END DO line = ADJUSTL(line((n_tab + 1):500)) ELSE ! n_tab = 1 column(i) = '' line = ADJUSTL(line(2:500)) END IF ELSE ! n_tab = 0; no tab left in line column(i) = TRIM(line) line = '' END IF ELSE column(i) = '' END IF END DO get_columns !Extract Findex: Findex = "WG000" ! template; but some 0's will be overwritten. IF (LEN_TRIM(column(3)) == 1) THEN ! 1-digit number Findex(5:5) = column(3)(1:1) ELSE IF (LEN_TRIM(column(3)) == 2) THEN ! 2-digit number Findex(4:5) = column(3)(1:2) ELSE IF (LEN_TRIM(column(3)) == 3) THEN ! 3-digit number Findex(3:5) = column(3)(1:3) ELSE ! too many digits warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Too many digits (>3) in WG_FAULT_ID, in column 3." Findex(1:5) = column(3)(1:5) END IF !--------------------------------------------------------------------------- IF (Findex /= old_Findex) offset_type_c1 = '' ! Memory should not carry over between fault sections. !--------------------------------------------------------------------------- !Extract offset_type_c1 (D, L, N, P, R, T), according to Peter Bird's NeoKinema/Slippery conventions. !MEASURED_COMPONENT_OF_SLIP = column(U) = column(21) may be one of: ! A = slip (total vector) ! B = throw (vertical component of slip) ! C = strike-slip (trace-parallel component of heave) ! D = opening/closing (trace-perpendicular component of heave) !SENSE_OF_MOTION = column(V) = column(22) may contain either one or two codes. ! If there are two codes, they are separated by hypen '-', and the first code is predominant sense. ! Codes may be one of: ! LL = left-lateral ! N = normal ! R = reverse ! RL = right-lateral ! Note that mixed codes must have exactly one of the dip-slip senses (N or R) ! and exactly one of the strike-slip senses (LL or RL). IF (LEN_TRIM(column(21)) == 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "MEASURED_COMPONENT_OF_SLIP in column 21 is blank." END IF n_hyphen = SCAN(column(22), '-') IF (n_hyphen > 0) THEN dominant_sense_c2 = column(22)(1:(n_hyphen-1)) secondary_sense_c2 = TRIM(column(22)((n_hyphen+1):50)) ELSE ! only dominant indicator given dominant_sense_c2 = TRIM(column(22)) secondary_sense_c2 = '' END IF IF (column(21) == 'A') THEN ! total slip !Note that RAKE in column(W) = column(23) is an integer in degrees: ! 0 = left-lateral, 90 = reverse, 180 = right-lateral, 270 (or -90) = normal. IF (LEN_TRIM(column(23)) == 0) THEN ! RAKE is blank; use dominant sense. IF (LEN_TRIM(column(22)) == 0) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "SENSE_OF_MOTION (col. 22) and RAKE (col. 23) are both blank." END IF IF (dominant_sense_c2 == "LL") THEN offset_type_c1 = 'L' offset_magnifier = one ELSE IF (dominant_sense_c2 == 'N') THEN IF (old_offset_type_c1 == 'D') THEN ! continue prevailing classification; do NOT trigger new_group offset_type_c1 = 'D' offset_magnifier = DCOS(D_dip_radians) ELSE ! establish offset type 'N' offset_type_c1 = 'N' offset_magnifier = DSIN(N_dip_radians) END IF ELSE IF (dominant_sense_c2 == 'R') THEN IF (old_offset_type_c1 == 'P') THEN ! continue prevailing classification; do NOT trigger new_group offset_type_c1 = 'P' offset_magnifier = DCOS(P_dip_radians) ELSE ! establish offset type 'T' offset_type_c1 = 'T' offset_magnifier = DSIN(T_dip_radians) END IF ELSE IF (dominant_sense_c2 == "RL") THEN offset_type_c1 = 'R' offset_magnifier = one ELSE warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Illegal SENSE_OF_MOTION in column 22." END IF ELSE ! make use of RAKE (see conventions in comments above); warning: SENSE_OF_MOTION might be blank! READ (column(23), *) rake_in_degrees rake_in_radians = rake_in_degrees / 57.29577951D0 IF (DSIN(rake_in_radians) >= 0.7071D0) THEN ! reverse faulting IF (old_offset_type_c1 == 'P') THEN ! continue prevailing classification; do NOT trigger new_group offset_type_c1 = 'P' offset_magnifier = DSIN(rake_in_radians) * DCOS(P_dip_radians) ELSE ! establish offset type 'T' offset_type_c1 = 'T' offset_magnifier = DSIN(rake_in_radians) * DSIN(T_dip_radians) END IF ELSE IF (DSIN(rake_in_radians) <= -0.7071D0) THEN ! normal faulting IF (old_offset_type_c1 == 'D') THEN ! continue prevailing classification; do NOT trigger new_group offset_type_c1 = 'D' offset_magnifier = -DSIN(rake_in_radians) * DCOS(D_dip_radians) ELSE ! establish offset type 'N' offset_type_c1 = 'N' offset_magnifier = -DSIN(rake_in_radians) * DSIN(N_dip_radians) END IF ELSE ! strike-slip faulting IF (DCOS(rake_in_radians) > 0.0D0) THEN ! left-lateral offset_type_c1 = 'L' offset_magnifier = DCOS(rake_in_radians) ELSE ! right-lateral offset_type_c1 = 'R' offset_magnifier = -DCOS(rake_in_radians) END IF END IF END IF ELSE IF (column(21) == 'B') THEN ! throw offset_magnifier = one IF (dominant_sense_c2 == 'N') THEN offset_type_c1 = 'N' ELSE IF (dominant_sense_c2 == 'R') THEN offset_type_c1 = 'T' ELSE ! look into secondary sense IF (secondary_sense_c2 == 'N') THEN offset_type_c1 = 'N' ELSE IF (secondary_sense_c2 == 'R') THEN offset_type_c1 = 'T' ELSE warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Inconsistent MEASURED_COMPONENT_OF_SLIP (col.21) vs. SENSE_OF_MOTION (col. 22)." END IF END IF ELSE IF (column(21) == 'C') THEN ! strike-slip offset_magnifier = one IF (dominant_sense_c2 == "LL") THEN offset_type_c1 = 'L' ELSE IF (dominant_sense_c2 == "RL") THEN offset_type_c1 = 'R' ELSE ! look into secondary sense IF (secondary_sense_c2 == "LL") THEN offset_type_c1 = 'L' ELSE IF (secondary_sense_c2 == "RL") THEN offset_type_c1 = 'R' ELSE warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Inconsistent MEASURED_COMPONENT_OF_SLIP (col.21) vs. SENSE_OF_MOTION (col. 22)." END IF END IF ELSE IF (column(21) == 'D') THEN ! opening/closing offset_magnifier = one IF (dominant_sense_c2 == 'N') THEN offset_type_c1 = 'D' ELSE IF (dominant_sense_c2 == 'R') THEN offset_type_c1 = 'P' ELSE ! look into secondary sense IF (secondary_sense_c2 == 'N') THEN offset_type_c1 = 'D' ELSE IF (secondary_sense_c2 == 'R') THEN offset_type_c1 = 'P' ELSE warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Inconsistent MEASURED_COMPONENT_OF_SLIP (col.21) vs. SENSE_OF_MOTION (col. 22)." END IF END IF ELSE warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "Illegal code (not A~D) in MEASURED_COMPONENT_OF_SLIP, in column 21." END IF !--------------------------------------------------------------------------- !Reject this line if there is no information on offset distance: IF ((LEN_TRIM(column(41)) == 0).AND.(LEN_TRIM(column(42)) == 0).AND.(LEN_TRIM(column(43)) == 0)) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "No information on offset distance provided." END IF !--------------------------------------------------------------------------- !Reject this line if there is no information on timing of offset: IF ((LEN_TRIM(column(25)) == 0).AND.(LEN_TRIM(column(26)) == 0).AND.(LEN_TRIM(column(27)) == 0).AND. & & (LEN_TRIM(column(29)) == 0).AND.(LEN_TRIM(column(30)) == 0).AND.(LEN_TRIM(column(31)) == 0)) THEN warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "No information on timing of offset provided." END IF !--------------------------------------------------------------------------- !======================================================================================= !Decide whether this line begins a new group (new fault section, or new offset sense) !by comparing new Findex with old_Findex, and new offset_type_c1 with old_offset_type_c1: new_group = ((Findex .NE. old_Findex).OR.(offset_type_c1 .NE. old_offset_type_c1)) !======================================================================================= !Extract fault_name: fault_name = TRIM(column(5)) !--------------------------------------------------------------------------- !Note: Editor's offset (summary for whole fault trace) is not provided in this dataset. editors_offset_km_c9 = '' !Note: Editor's offset sigma (summary for whole fault trace) is not provided in this dataset. editors_offset_sigma_km_c9 = '' !Note: Editor's began_Ma (summary for whole fault trace) is not provided in this dataset. editors_began_Ma_c9 = '' !Note: Editor's ended_Ma (summary for whole fault trace) is not provided in this dataset. editors_ended_Ma_c9 = '' !Note: Editor's neorate (summary for whole fault trace) is not provided in this dataset. editors_neorate_mmpa_c9 = '' !Note: Editor's neorate sigma (summary for whole fault trace) is not provided in this dataset. editors_neorate_sigma_mmpa_c9 = '' !N.B. All are returned left-justified. !--------------------------------------------------------------------------- !Extract source_letter (P = Peter Bird; Q = QFaults; F = SCEC-FAD; R = researcher contribution) column(10) = ADJUSTL(column(10)) ! because many " SCEC-FAD" entries have leading space IF ((column(10)(1:5) == "Peter").OR.(column(10)(1:1) == '3')) THEN source_letter = 'P' ELSE IF ((column(10)(1:1) == 'Q').OR.(column(10)(1:1) == '2')) THEN source_letter = 'Q' ELSE IF ((column(10)(1:4) == "SCEC").OR.(column(10)(1:1) == '1')) THEN source_letter = 'F' ELSE IF ((column(10)(1:1) == 'R').OR.(column(10)(1:1) == 'r').OR.(column(10)(1:1) == '4')) THEN source_letter = 'R' ELSE source_letter = '?' END IF !N.B. This letter just appears in the output table; it is not used in computations. !--------------------------------------------------------------------------- !Extract grade (A = primary, B = secondary, C = tertiary) IF (source_letter == 'P') THEN grade = column(6)(8:8) ELSE IF (source_letter == 'Q') THEN grade = 'A' ELSE IF (source_letter == 'F') THEN grade = 'A' ELSE IF (source_letter == 'R') THEN grade = 'A' ELSE grade = 'C' END IF !--------------------------------------------------------------------------- !Extract datum_name: datum_name = TRIM(column(18)) ! SHORT_CITATION !--------------------------------------------------------------------------- !======================================================================================= !Test for paleotectonic datum == pre-Holocene ending time, !and re-cast as datum with zero displacement since age of the overlap assemblage! !(However, all _c20 variables used to represent the datum on-screen, and in tables, are processed normally, below.) !Note: column(28) = END_TIME_UNITS; column(29) = PREF_END_TIME; column(30) = MAX_END_TIME(earliest); column(31) = MIN_END_TIME IF (LEN_TRIM(column(31)) > 0) THEN ! we have MIN_END_TIME entry READ (column(31), *) min_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN min_end_time_Ma = min_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN min_end_time_Ma = (2007.D0 - min_end_time) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN min_end_time_Ma = (min_end_time + 2007.D0) * 1.0D-6 END IF paleotectonic = (min_end_time_Ma > 0.010D0) pivot_time_Ma = min_end_time_Ma ELSE IF (LEN_TRIM(column(29)) > 0) THEN ! if no MIN_END_TIME, try for PREF_END_TIME READ (column(29), *) nominal_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN nominal_end_time_Ma = nominal_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN nominal_end_time_Ma = (2007.D0 - nominal_end_time) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN nominal_end_time_Ma = (nominal_end_time + 2007.D0) * 1.0D-6 END IF paleotectonic = (nominal_end_time_Ma > 0.010D0) pivot_time_Ma = nominal_end_time_Ma ELSE ! no information on end time; assume fault is still active paleotectonic = .FALSE. END IF IF (paleotectonic) THEN F_is_zero = .TRUE. !Use pivot_time_Ma as start time for offset of zero: IF ((LEN_TRIM(column(30)) > 0).AND.(LEN_TRIM(column(31)) > 0)) THEN !use two bounds on pivot time if available READ (column(30), *) max_age IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN max_age_Ma = max_age * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN max_age_Ma = MAX((2007. - max_age), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN max_age_Ma = (max_age + 2007.) * 1.0D-6 END IF max_age_Ma_c9 = ADJUSTL(ASCII9(max_age_Ma)) READ (column(31), *) min_age IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN min_age_Ma = min_age * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN min_age_Ma = MAX((2007. - min_age), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN min_age_Ma = (min_age + 2007.) * 1.0D-6 END IF min_age_Ma_c9 = ADJUSTL(ASCII9(min_age_Ma)) age_is_Gaussian = .FALSE. age_has_two_bounds = .TRUE. sigma_Delta_t_min_stated = .FALSE. Delta_t_min_range_stated = .FALSE. Delta_t_min = min_age_Ma Delta_t_min_string = min_age_Ma_c9 sigma_Delta_t_max_stated = .FALSE. Delta_t_max_range_stated = .FALSE. Delta_t_max = max_age_Ma Delta_t_max_string = max_age_Ma_c9 !N.B. Calling program, Slippery.f90, currently makes ! no use of age_from_raw_14C or age_from_raw_nuclides ! in this case, so no values are set. (Default = .FALSE.) !******************************************** ELSE IF (LEN_TRIM(column(29)) > 0) THEN !second-best is for a preferred pivot time to be given: READ (column(29), *) nominal_age IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN nominal_age_Ma = nominal_age * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN nominal_age_Ma = MAX((2007. - nominal_age), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN nominal_age_Ma = (nominal_age + 2007.) * 1.0D-6 END IF nominal_age_Ma_c9 = ADJUSTL(ASCII9(nominal_age_Ma)) age_is_Gaussian = .TRUE. sigma_Delta_t_stated = .FALSE. age_range_stated = .FALSE. Delta_t_nominal = nominal_age_Ma Delta_t_string = nominal_age_Ma_c9 IF (nominal_age_Ma <= 0.03) THEN ! <= 30,000 years; see p. 126 of Yeats, Sieh, & Allen [1997] The Geology of Earthquakes. age_from_raw_14C = .TRUE. ELSE IF (nominal_age_Ma <= 0.6) THEN ! <= 600,000 years {rather arbitrary, but younger than Bishop tuff = 0.768 Ma} age_from_raw_nuclides = .TRUE. END IF !******************************************** ELSE IF (LEN_TRIM(column(31)) > 0) THEN !third choice is to have only a lower limit on the pivot age: READ (column(31), *) min_age IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN min_age_Ma = min_age * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN min_age_Ma = MAX((2007. - min_age), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN min_age_Ma = (min_age + 2007.) * 1.0D-6 END IF min_age_Ma_c9 = ADJUSTL(ASCII9(min_age_Ma)) age_is_Gaussian = .FALSE. age_has_two_bounds = .FALSE. age_is_lower_limit = .TRUE. age_is_upper_limit = .FALSE. sigma_Delta_t_min_stated = .FALSE. Delta_t_min_range_stated = .FALSE. Delta_t_min = min_age_Ma Delta_t_min_string = min_age_Ma_c9 !N.B. Calling program, Slippery.f90, currently makes ! no use of age_from_raw_14C or age_from_raw_nuclides ! in this case, so no values are set. (Default = .FALSE.) !******************************************** ENDIF END IF ! paleotectonic !======================================================================================= neotectonic = (.NOT.paleotectonic) !Extract all parameters related to fault-offset F: ! datum_offset_km_c20, datum_offset_sigma_km_c20, ! F_is_Gaussian, F_is_upper_limit, F_is_lower_limit, F_is_zero, ! sigma_F_stated, F_range_stated, ! Delta_F_nominal, sigma_Delta_F, Delta_F_nominal_high, Delta_F_nominal_low, Delta_F_string, ! Delta_F_max, sigma_Delta_F_max, Delta_F_max_high, Delta_F_max_low, ! Delta_F_min, sigma_Delta_F_min, Delta_F_min_high, Delta_F_min_low. !Notes: PaleoSites has all offsets in m; must be converted to km for Slippery! ! column(41) = PREF_OFFSET_(m); column(42) = MAX_OFFSET; column(43) = MIN_OFFSET IF ((LEN_TRIM(column(42)) > 0).AND.(LEN_TRIM(column(43)) > 0)) THEN ! range of offset is the preferred model, if available. READ (column(42), *) max_offset_m max_offset_km = max_offset_m * 1.0D-3 max_offset_km_c9 = ADJUSTL(ASCII9(max_offset_km)) READ (column(43), *) min_offset_m min_offset_km = min_offset_m * 1.0D-3 min_offset_km_c9 = ADJUSTL(ASCII9(min_offset_km)) datum_offset_km_c20 = TRIM(min_offset_km_c9) // '~' // TRIM(max_offset_km_c9) datum_offset_sigma_km_c20 = '' IF (neotectonic) THEN F_is_Gaussian = .TRUE. sigma_F_stated = .FALSE. F_range_stated = .TRUE. Delta_F_nominal_low = min_offset_km Delta_F_nominal_high = max_offset_km END IF ELSE IF (LEN_TRIM(column(41)) > 0) THEN !second best is to have a preferred offset stated READ (column(41), *) nominal_offset_m nominal_offset_km = nominal_offset_m * 1.0D-3 nominal_offset_km_c9 = ADJUSTL(ASCII9(nominal_offset_km)) datum_offset_km_c20 = nominal_offset_km_c9 datum_offset_sigma_km_c20 = '' IF (neotectonic) THEN F_is_Gaussian = .TRUE. sigma_F_stated = .FALSE. F_range_stated = .FALSE. Delta_F_nominal = nominal_offset_km Delta_F_string = nominal_offset_km_c9 ! c20 END IF ELSE IF (LEN_TRIM(column(42)) > 0) THEN !third choice is a maximum offset READ (column(42), *) max_offset_m max_offset_km = max_offset_m * 1.0D-3 max_offset_km_c9 = ADJUSTL(ASCII9(max_offset_km)) datum_offset_km_c20 = '<' // max_offset_km_c9 datum_offset_sigma_km_c20 = '' IF (neotectonic) THEN F_is_Gaussian = .FALSE. F_is_upper_limit = .TRUE. sigma_F_stated = .FALSE. F_range_stated = .FALSE. Delta_F_max = max_offset_km Delta_F_string = max_offset_km_c9 ! c20 END IF ELSE IF (LEN_TRIM(column(43)) > 0) THEN !third choice is a minimum offset READ (column(43), *) min_offset_m min_offset_km = min_offset_m * 1.0D-3 min_offset_km_c9 = ADJUSTL(ASCII9(min_offset_km)) datum_offset_km_c20 = '>' // min_offset_km_c9 datum_offset_sigma_km_c20 = '' IF (neotectonic) THEN F_is_Gaussian = .FALSE. F_is_upper_limit = .FALSE. F_is_lower_limit = .TRUE. sigma_F_stated = .FALSE. F_range_stated = .FALSE. Delta_F_min = min_offset_km Delta_F_string = min_offset_km_c9 ! c20 END IF END IF ! No case provided for all-blank offset fields; we already set warning for this case, above. !---------------- end of offset Delta_F ------------------------------ !Extract all parameters related to age of offset feature, Delta_t: ! datum_began_Ma_c20, datum_ended_Ma_c20, ! age_is_Gaussian, age_has_two_bounds, age_is_lower_limit, age_is_upper_limit, ! sigma_Delta_t_stated, age_range_stated, age_from_raw_14C, age_from_raw_nuclides, ! Delta_t_nominal, sigma_Delta_t, age_nominal_high, age_nominal_low, Delta_t_string, ! sigma_Delta_t_min_stated, Delta_t_min_range_stated, Delta_t_min, sigma_Delta_t_min, ! Delta_t_min_high, Delta_t_min_low, Delta_t_min_string, ! sigma_Delta_t_max_stated, Delta_t_max_range_stated, Delta_t_max, sigma_Delta_t_max, ! Delta_t_max_high, Delta_t_max_low, Delta_t_max_string !Variables related to start time/age of offset feature: !NOTE: column(24) = START_TIME_UNITS, column(25) = PREF_START_TIME, column(26) = MAX_START_TIME, column(27) = MIN_START_TIME IF ((LEN_TRIM(column(26)) > 0).AND.(LEN_TRIM(column(27)) > 0)) THEN !First preference is for a range of ages: READ (column(26), *) max_age IF ((column(24)(1:2) == "KA").OR.(column(24)(1:2) == "ka")) THEN max_age_Ma = max_age * 1.0D-3 ELSE IF ((column(24)(1:2) == "AD").OR.(column(24)(1:2) == "ad")) THEN max_age_Ma = MAX((2007. - max_age), one) * 1.0D-6 ELSE IF ((column(24)(1:2) == "BC").OR.(column(24)(1:2) == "bc")) THEN max_age_Ma = (max_age + 2007.) * 1.0D-6 END IF max_age_Ma_c9 = ADJUSTL(ASCII9(max_age_Ma)) READ (column(27), *) min_age IF ((column(24)(1:2) == "KA").OR.(column(24)(1:2) == "ka")) THEN min_age_Ma = min_age * 1.0D-3 ELSE IF ((column(24)(1:2) == "AD").OR.(column(24)(1:2) == "ad")) THEN min_age_Ma = MAX((2007. - min_age), one) * 1.0D-6 ELSE IF ((column(24)(1:2) == "BC").OR.(column(24)(1:2) == "bc")) THEN min_age_Ma = (min_age + 2007.) * 1.0D-6 END IF min_age_Ma_c9 = ADJUSTL(ASCII9(min_age_Ma)) datum_began_Ma_c20 = TRIM(max_age_Ma_c9) // '~' // TRIM(min_age_Ma_c9) IF (neotectonic) THEN age_is_Gaussian = .FALSE. age_has_two_bounds = .TRUE. sigma_Delta_t_min_stated = .FALSE. Delta_t_min_range_stated = .FALSE. Delta_t_min = min_age_Ma Delta_t_min_string = min_age_Ma_c9 sigma_Delta_t_max_stated = .FALSE. Delta_t_max_range_stated = .FALSE. Delta_t_max = max_age_Ma Delta_t_max_string = max_age_Ma_c9 !N.B. Calling program, Slippery.f90, currently makes ! no use of age_from_raw_14C or age_from_raw_nuclides ! in this case, so no values are set. (Default = .FALSE.) !******************************************** END IF ELSE IF (LEN_TRIM(column(25)) > 0) THEN !second-best is for a preferred start time to be given: READ (column(25), *) nominal_age IF ((column(24)(1:2) == "KA").OR.(column(24)(1:2) == "ka")) THEN nominal_age_Ma = nominal_age * 1.0D-3 ELSE IF ((column(24)(1:2) == "AD").OR.(column(24)(1:2) == "ad")) THEN nominal_age_Ma = MAX((2007. - nominal_age), one) * 1.0D-6 ELSE IF ((column(24)(1:2) == "BC").OR.(column(24)(1:2) == "bc")) THEN nominal_age_Ma = (nominal_age + 2007.) * 1.0D-6 END IF nominal_age_Ma_c9 = ADJUSTL(ASCII9(nominal_age_Ma)) datum_began_Ma_c20 = nominal_age_Ma_c9 IF (neotectonic) THEN age_is_Gaussian = .TRUE. sigma_Delta_t_stated = .FALSE. age_range_stated = .FALSE. Delta_t_nominal = nominal_age_Ma Delta_t_string = nominal_age_Ma_c9 IF (nominal_age_Ma <= 0.03) THEN ! <= 30,000 years; see p. 126 of Yeats, Sieh, & Allen [1997] The Geology of Earthquakes. age_from_raw_14C = .TRUE. ELSE IF (nominal_age_Ma <= 0.6) THEN ! <= 600,000 years {rather arbitrary, but younger than Bishop tuff = 0.768 Ma} age_from_raw_nuclides = .TRUE. END IF END IF ELSE IF (LEN_TRIM(column(26)) > 0) THEN !third choice is to have only a maximum age: READ (column(26), *) max_age IF ((column(24)(1:2) == "KA").OR.(column(24)(1:2) == "ka")) THEN max_age_Ma = max_age * 1.0D-3 ELSE IF ((column(24)(1:2) == "AD").OR.(column(24)(1:2) == "ad")) THEN max_age_Ma = MAX((2007. - max_age), one) * 1.0D-6 ELSE IF ((column(24)(1:2) == "BC").OR.(column(24)(1:2) == "bc")) THEN max_age_Ma = (max_age + 2007.) * 1.0D-6 END IF max_age_Ma_c9 = ADJUSTL(ASCII9(max_age_Ma)) datum_began_Ma_c20 = '<' // max_age_Ma_c9 IF (neotectonic) THEN age_is_Gaussian = .FALSE. age_has_two_bounds = .FALSE. age_is_lower_limit = .FALSE. age_is_upper_limit = .TRUE. sigma_Delta_t_max_stated = .FALSE. Delta_t_max_range_stated = .FALSE. Delta_t_max = max_age_Ma Delta_t_max_string = max_age_Ma_c9 !N.B. Calling program, Slippery.f90, currently makes ! no use of age_from_raw_14C or age_from_raw_nuclides ! in this case, so no values are set. (Default = .FALSE.) !******************************************** END IF ELSE IF (LEN_TRIM(column(27)) > 0) THEN !fourth choice is to have only a lower limit on the age: READ (column(27), *) min_age IF ((column(24)(1:2) == "KA").OR.(column(24)(1:2) == "ka")) THEN min_age_Ma = min_age * 1.0D-3 ELSE IF ((column(24)(1:2) == "AD").OR.(column(24)(1:2) == "ad")) THEN min_age_Ma = MAX((2007. - min_age), one) * 1.0D-6 ELSE IF ((column(24)(1:2) == "BC").OR.(column(24)(1:2) == "bc")) THEN min_age_Ma = (min_age + 2007.) * 1.0D-6 END IF min_age_Ma_c9 = ADJUSTL(ASCII9(min_age_Ma)) datum_began_Ma_c20 = '>' // min_age_Ma_c9 IF (neotectonic) THEN age_is_Gaussian = .FALSE. age_has_two_bounds = .FALSE. age_is_lower_limit = .TRUE. age_is_upper_limit = .FALSE. sigma_Delta_t_min_stated = .FALSE. Delta_t_min_range_stated = .FALSE. Delta_t_min = min_age_Ma Delta_t_min_string = min_age_Ma_c9 !N.B. Calling program, Slippery.f90, currently makes ! no use of age_from_raw_14C or age_from_raw_nuclides ! in this case, so no values are set. (Default = .FALSE.) !******************************************** END IF ELSE ! all start time fields are blank. [ N.B. The case we already rejected was all start AND end times blank. ] warning_count = MIN(20, warning_count + 1) warning_messages(warning_count) = "No start-time information provided." datum_began_Ma_c20 = '' END IF !--------------------- end of start time/age -------------------------------- !Variables related to end-time of fault slip: !Note: column(28) = END_TIME_UNITS; column(29) = PREF_END_TIME; column(30) = MAX_END_TIME(earliest); column(31) = MIN_END_TIME !GPBhere IF (neotectonic) THEN IF (warning_count == 0) THEN datum_ended_Ma_c20 = 'active' END IF ELSE ! paleotectonic IF ((LEN_TRIM(column(30)) > 0).AND.(LEN_TRIM(column(31)) > 0)) THEN !First preference is for a range of end-times: READ (column(30), *) max_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN max_end_time_Ma = max_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN max_end_time_Ma = MAX((2007. - max_end_time), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN max_end_time_Ma = (max_end_time + 2007.) * 1.0D-6 END IF max_end_time_Ma_c9 = ADJUSTL(ASCII9(max_end_time_Ma)) READ (column(31), *) min_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN min_end_time_Ma = min_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN min_end_time_Ma = MAX((2007. - min_end_time), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN min_end_time_Ma = (min_end_time + 2007.) * 1.0D-6 END IF min_end_time_Ma_c9 = ADJUSTL(ASCII9(min_age_Ma)) datum_ended_Ma_c20 = TRIM(max_end_time_Ma_c9) // '~' // TRIM(min_end_time_Ma_c9) ELSE IF (LEN_TRIM(column(29)) > 0) THEN !second-best is for a preferred end time to be given: READ (column(29), *) nominal_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN nominal_end_time_Ma = nominal_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN nominal_end_time_Ma = MAX((2007. - nominal_end_time), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN nominal_end_time_Ma = (nominal_end_time + 2007.) * 1.0D-6 END IF nominal_end_time_Ma_c9 = ADJUSTL(ASCII9(nominal_end_time_Ma)) datum_ended_Ma_c20 = nominal_end_time_Ma_c9 ELSE IF (LEN_TRIM(column(30)) > 0) THEN !third choice is to have only a maximum end time: READ (column(30), *) max_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN max_end_time_Ma = max_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN max_end_time_Ma = MAX((2007. - max_end_time), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN max_end_time_Ma = (max_end_time + 2007.) * 1.0D-6 END IF max_end_time_Ma_c9 = ADJUSTL(ASCII9(max_end_time_Ma)) datum_ended_Ma_c20 = '<' // max_end_time_Ma_c9 ELSE IF (LEN_TRIM(column(31)) > 0) THEN !fourth choice is to have only a lower limit on the end time: READ (column(31), *) min_end_time IF ((column(28)(1:2) == "KA").OR.(column(28)(1:2) == "ka")) THEN min_end_time_Ma = min_end_time * 1.0D-3 ELSE IF ((column(28)(1:2) == "AD").OR.(column(28)(1:2) == "ad")) THEN min_end_time_Ma = MAX((2007. - min_end_time), one) * 1.0D-6 ELSE IF ((column(28)(1:2) == "BC").OR.(column(28)(1:2) == "bc")) THEN min_end_time_Ma = (min_end_time + 2007.) * 1.0D-6 END IF min_end_time_Ma_c9 = ADJUSTL(ASCII9(min_end_time_Ma)) datum_ended_Ma_c20 = '>' // min_end_time_Ma_c9 ELSE ! all end time fields are blank. [ N.B. The case we already rejected was all start AND end times blank. ] datum_ended_Ma_c20 = '' ! N.B. Program flow should (theoretically) never come through here. END IF END IF ! neotectonic, or paleotectonic !--------------------- end of end-time -------------------------------- !Extract (human PaleoSites editor's) datum_neorate_mmpa_c20: IF ((LEN_TRIM(column(36)) > 0).AND.(LEN_TRIM(column(37)) > 0)) THEN ! MIN_SLIP_RATE & MAX_SLIP_RATE both available. datum_neorate_mmpa_c20 = TRIM(column(37)) // '~' // TRIM(column(36)) ELSE IF (LEN_TRIM(column(35)) > 0) THEN ! PREF_SLIP_RATE_(mm/yr) [ N.B. Should read "PREF_OFFSET_RATE_(mm/yr)". ] datum_neorate_mmpa_c20 = TRIM(column(35)) ELSE IF (LEN_TRIM(column(36)) > 0) THEN ! MAX_SLIP_RATE [ N.B. Should read "MAX_OFFSET_RATE". ] datum_neorate_mmpa_c20 = '<' // TRIM(column(36)) ELSE IF (LEN_TRIM(column(37)) > 0) THEN ! MIN_SLIP_RATE [ N.B. Should read "MIN_OFFSET_RATE". ] datum_neorate_mmpa_c20 = '>' // TRIM(column(37)) ELSE datum_neorate_mmpa_c20 = '' END IF !Note: (Human PaleoSites editor's) datum_neorate_mmpa_c20 is not provided in this dataset. datum_neorate_sigma_mmpa_c20 = '' !--------------------------------------------------------------------------- END SUBROUTINE Read_Line_from_PaleoSites_Table SUBROUTINE Read_1_pdf (pdf1) !Reads a single pdf from an already-open file on device 1. Suggestion: OPEN(...., PAD = "YES") !TYPE pdf ! probability density function ! CHARACTER*40 :: name ! DOUBLE PRECISION :: x_min, x_max ! limits on range of independent variable ! LOGICAL :: logarithmic ! should never be T unless x_min > 0. ! DOUBLE PRECISION, DIMENSION(nbins) :: pd ! probability density (per unit of x) in each bin ! DOUBLE PRECISION, DIMENSION(nbins) :: xlist ! independent variable value at the center of each bin ! DOUBLE PRECISION, DIMENSION(0:nbins) :: rxlist ! independent variable value at the right of each bin ! LOGICAL :: normalized ! only T if integral of this pdf is 1.0000D0 !END TYPE pdf IMPLICIT NONE CHARACTER(132) :: temp_line TYPE (pdf), INTENT(OUT) :: pdf1 INTEGER :: i READ (1, "(A)") pdf1%name IF (pdf1%name(1:4) == "New ") THEN temp_line = pdf1%name pdf1%name = TRIM("Inherited " // temp_line(5:132)) END IF READ (1, "(2ES12.4)") pdf1%x_min, pdf1%x_max READ (1, "(L1)") pdf1%logarithmic pdf1%rxlist(0) = pdf1%x_min DO i = 1, nbins READ (1, "(3ES12.4)") pdf1%xlist(i), pdf1%pd(i), pdf1%rxlist(i) END DO READ (1, "(L1)") pdf1%normalized END SUBROUTINE Read_1_pdf SUBROUTINE Read_6_pdfs (L_D_prior_pdf, L_L_prior_pdf, L_N_prior_pdf, L_P_prior_pdf, L_R_prior_pdf, L_T_prior_pdf) IMPLICIT NONE TYPE (pdf), INTENT(OUT) :: L_D_prior_pdf, L_L_prior_pdf, L_N_prior_pdf, L_P_prior_pdf, L_R_prior_pdf, L_T_prior_pdf CHARACTER(80) :: priors_pathfile INTEGER :: ios 10 WRITE (*, "(' Enter name for (existing) file that contains 6 old L_x_prior pdfs: ')") READ (*, "(A)") priors_pathfile OPEN (UNIT = 1, FILE = TRIM(priors_pathfile), STATUS = "OLD", IOSTAT = ios, PAD = "YES") IF (ios /= 0) THEN WRITE (*, "(' ',A/' could not be read [in current directory].' & & /'Try again [perhaps including path?].')") TRIM(priors_pathfile) CALL Pause() GO TO 10 END IF CALL Read_1_pdf (L_D_prior_pdf) CALL Read_1_pdf (L_L_prior_pdf) CALL Read_1_pdf (L_N_prior_pdf) CALL Read_1_pdf (L_P_prior_pdf) CALL Read_1_pdf (L_R_prior_pdf) CALL Read_1_pdf (L_T_prior_pdf) CLOSE (1) END SUBROUTINE Read_6_pdfs SUBROUTINE Silent_Pause() IMPLICIT NONE READ (*,*) END SUBROUTINE Silent_Pause SUBROUTINE Write_1_pdf (pdf1) !Writes a single pdf into an already-open file on device 1. !TYPE pdf ! probability density function ! CHARACTER*40 :: name ! DOUBLE PRECISION :: x_min, x_max ! limits on range of independent variable ! LOGICAL :: logarithmic ! should never be T unless x_min > 0. ! DOUBLE PRECISION, DIMENSION(nbins) :: pd ! probability density (per unit of x) in each bin ! DOUBLE PRECISION, DIMENSION(nbins) :: xlist ! independent variable value at the center of each bin ! DOUBLE PRECISION, DIMENSION(0:nbins) :: rxlist ! independent variable value at the right of each bin ! LOGICAL :: normalized ! only T if integral of this pdf is 1.0000D0 !END TYPE pdf IMPLICIT NONE TYPE (pdf), INTENT(IN) :: pdf1 INTEGER :: i WRITE (1, "(A)") pdf1%name WRITE (1, "(2ES12.4)") pdf1%x_min, pdf1%x_max WRITE (1, "(L1)") pdf1%logarithmic DO i = 1, nbins WRITE (1, "(3ES12.4)") pdf1%xlist(i), pdf1%pd(i), pdf1%rxlist(i) END DO ! N.B. This format is wasteful of CR/LF bytes, but it permits the left 2 columns to be ! captured in Excel for easy plotting of the pdf as an x/y piecewise-linear graph or bar graph. ! Note that pdf1%rxlist(0) = pdf1%x_min is not ever written out explicitly, but pdf1%x_min IS, ! so the companion routine Read_1_pdf has no trouble inferring its value. WRITE (1, "(L1)") pdf1%normalized END SUBROUTINE Write_1_pdf SUBROUTINE Write_6_pdfs (L_D_prior_pdf, L_L_prior_pdf, L_N_prior_pdf, L_P_prior_pdf, L_R_prior_pdf, L_T_prior_pdf) IMPLICIT NONE TYPE (pdf), INTENT(IN) :: L_D_prior_pdf, L_L_prior_pdf, L_N_prior_pdf, L_P_prior_pdf, L_R_prior_pdf, L_T_prior_pdf CHARACTER(80) :: priors_pathfile INTEGER :: ios 10 WRITE (*, "(' Enter name for NEW .txt file to contain 6 L_X_posterior pdfs/cpfs: ')") READ (*, "(A)") priors_pathfile OPEN (UNIT = 1, FILE = TRIM(priors_pathfile), STATUS = "NEW", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Could not create a NEW file with that name.')") CALL Pause() GO TO 10 END IF CALL Write_1_pdf (L_D_prior_pdf) CALL Write_1_pdf (L_L_prior_pdf) CALL Write_1_pdf (L_N_prior_pdf) CALL Write_1_pdf (L_P_prior_pdf) CALL Write_1_pdf (L_R_prior_pdf) CALL Write_1_pdf (L_T_prior_pdf) CLOSE (1) END SUBROUTINE Write_6_pdfs END PROGRAM Slippery