PROGRAM c_rst_Summary ! A utility program in the Restore4 family. ! Reads a c.rst file OUTPUT from Restore4 (e.g., cUSM_NI_012.0Ma.rst), ! adds up the model increments (during time-steps so far completed) ! towards meeting goals for lengthening (or shortening) restored cross-sections, ! and writes these as 3 new columns, to the right of the c.rst INPUT table. ! The resulting c_rst_Summary file (e.g., cUSM_NI_012.0Ma_rst_Summary.txt) ! can be opened with Excel or other spreadsheet for plotting. ! Note 2 structural differences from the similar companion program p_rst_Summary: !(1) Paleomagnetic data gives 2 targets per site (latitude anomaly, and ! vertical-axis rotation), but cross-sections give only one: ! the length increase (or "stretch") as geologic time moved forward to the present. ! [Cross-sections that show folding and thrust-faulting will have negative stretch goals.] !(2) Paleomagnetic data always span to the present, so modeling of such data are either ! COMPLETE (if the model report time is older than the age of magnetization), or ! INCOMPLETE (if the model report time is younger). However, for cross-sections, there ! is a third group: WAITING (when the section MinAge is older than the model report time). ! By Peter Bird, UCLA, June 2020 !- - - - - - - - - - - - - - - - - - - IMPLICIT NONE CHARACTER(1) :: tab CHARACTER(5) :: c5, c5_code CHARACTER(6) :: c6 CHARACTER(7) :: c7 CHARACTER(8) :: c8 CHARACTER(10) :: c10 CHARACTER(47) :: c47_citation CHARACTER(80) :: output_filename, c_rst_filename CHARACTER(134) :: c_rst_format ! to read cross-section data CHARACTER(134) :: c_rst_titles ! to write cross-section data CHARACTER(134) :: line INTEGER :: c_rst_count ! number of cross-sections in dataset INTEGER :: delta_t_denominator, i, ios, j, high_time, length REAL*8 :: delta_t_Ma ! inferred timestep REAL*8 :: highest_t_Ma ! greatest geologic age reached so far, in Ma REAL*8 :: delta_t_numerator, duration_in_my, fraction, overlap_in_my, t1, t2, t3, t4 REAL*8, DIMENSION(9) :: vector !-------------------------------------------------------------------------------------------------------- CHARACTER*47, DIMENSION(:), ALLOCATABLE :: c_reference ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:,:),ALLOCATABLE :: c_location ! (1:4 = lon1, lat1, lon2, lat2; 1:c_rst_count = cross-section index) CHARACTER*5, DIMENSION(:), ALLOCATABLE :: c_code ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:), ALLOCATABLE :: c_length_Now_km ! present length of cross-section (in km) ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:), ALLOCATABLE :: c_length_Was_km ! restored (goal; datum) length of cross-section (in km) ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:), ALLOCATABLE :: c_restoration_sigma_km ! uncertainty in length of restored cross-section (in km) ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:), ALLOCATABLE :: c_MaxAge_Ma ! geologic age of restored/undeformed cross-section (in Ma) ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:), ALLOCATABLE :: c_MinAge_Ma ! geologic age of final cross-section when deformation ended; overlap assemblage age (in Ma) ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:), ALLOCATABLE :: c_goal_stretch_km ! target (data) increase in length of cross-section (in km) as geologic time moved forward toward present ! (1:c_rst_count = cross-section index) REAL*8, DIMENSION(:,:),ALLOCATABLE :: c_goal_rate_mmpa ! target (data) rate of increase in length of cross-section, as geologic time moved forward toward present, ! in units of km/m.y. = mm/a. !(1:num_timesteps = timestep index; 1:c_rst_count = cross-section index) REAL*8, DIMENSION(:,:),ALLOCATABLE :: c_model_rate_mmpa ! model (Restore4) rate of increase in length of cross-section, as geologic time moved forward toward present, ! in units of km/m.y. = mm/a. !(1:num_timesteps = timestep index; 1:c_rst_count = cross-section index) REAL*8, DIMENSION(0:2) :: c_errors ! 3 norms of the list of errors in stretch (each stretch normalized by its sigma before combining): !(0:2 = L0,L1,L2 norm). INTEGER, DIMENSION(:), ALLOCATABLE :: summary_class012 ! 0 = WAITING; 1 = INCOMPLETE; 2 = COMPLETE REAL*8, DIMENSION(:), ALLOCATABLE :: summary_goal_stretch_km REAL*8, DIMENSION(:), ALLOCATABLE :: summary_goal_sigma_km REAL*8, DIMENSION(:), ALLOCATABLE :: summary_model_stretch_km !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (*, "(' PROGRAM c_rst_Summary:')") WRITE (*, "(' A utility program in the Restore4 family.')") WRITE (*, "(' Reads a c.rst file OUTPUT from Restore4 (e.g., cUSM_NI_012.0Ma.rst),')") WRITE (*, "(' adds up the model increments (during time-steps so far completed)')") WRITE (*, "(' towards meeting goals for lengthening (or shortening) restored cross-sections,')") ! <==(maximum length) WRITE (*, "(' and writes these as 3 new columns, to the right of the c.rst INPUT table.')") WRITE (*, "(' The resulting c_rst_Summary file (e.g., cUSM_NI_012.0Ma_rst_Summary.txt)')") WRITE (*, "(' can be opened with Excel or other spreadsheet for plotting.')") WRITE (*, "(' Note 2 structural differences from the similar program p_rst_Summary:')") WRITE (*, "(' (1) Paleomagnetic data gives 2 targets per site (latitude anomaly, and')") WRITE (*, "(' vertical-axis rotation), but cross-sections give only one: ')") WRITE (*, "(' the length increase (or ""stretch"") as geologic time moved forward.')") WRITE (*, "(' [Cross-sections with folding and thrusting have negative stretch goals.]')") WRITE (*, "(' (2) Paleomagnetic data always span to present, so models of such data are')") WRITE (*, "(' COMPLETE (if the model report time is older than the magnetization), or')") WRITE (*, "(' INCOMPLETE (if the model report time is younger).')") WRITE (*, "(' However, for cross-sections, there is a third group: ')") WRITE (*, "(' WAITING (when the section MinAge is older than the model report time).')") WRITE (*, "(' By Peter Bird, UCLA, June 2020')") CALL Pause() 10 WRITE (*, *) WRITE (*, "(' Enter filename of c.rst file OUTPUT from Restore: ')") READ (*, "(A)") c_rst_filename OPEN (UNIT = 1, FILE = TRIM(c_rst_filename), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File ', A, ' not found (in current folder). Try again...')") TRIM(c_rst_filename) CALL Pause() GO TO 10 END IF !First time through, just count number of sites (c_rst_count) and greatest timestep achieved (top_time): highest_t_Ma = 0.0D0 ! just initializing before search delta_t_Ma = 0.0D0 c_rst_count = 0 delta_t_numerator = 0.0D0 delta_t_denominator = 0 READ (1, "(A)") c_rst_format READ (1, "(A)") c_rst_titles scanning: DO READ (1, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT scanning IF (line(1:1) == '+') THEN ! present location ! (no action needed in this case) ELSE IF (line(1:1) == '*') THEN ! progress toward stretch in one timestep line(1:1) = ' ' READ (line, *) t1, t2 delta_t_numerator = delta_t_numerator + (t2 - t1) delta_t_denominator = delta_t_denominator + 1 highest_t_Ma = MAX(highest_t_Ma, t2) ELSE ! header line with original INPUT line of the INPUT c_rst file READ (line, c_rst_format, IOSTAT = ios) c47_citation, (vector(j), j = 1, 4), c5_code, (vector(j), j = 5, 9) IF (ios /= 0) EXIT scanning ! perhaps input file has empty lines at end? c_rst_count = c_rst_count + 1 END IF END DO scanning CLOSE (1) ! c_rst input file (which was OUTPUT from Restore) !Infer timestep: IF (delta_t_denominator > 0) THEN delta_t_Ma = delta_t_numerator / delta_t_denominator ELSE WRITE (*, *) WRITE (*, "(' ERROR: This seems to be(?) a c.rst INPUT file, not a c.rst OUTPUT file.')") CALL Pause() STOP END IF !Infer high_time (INTEGER): high_time = NINT(highest_t_Ma / delta_t_Ma) WRITE (*, *) WRITE (*, "(' This file describes ', I5, ' restored cross-sections within the current F-E grid.')") c_rst_count WRITE (*, "(' This file contains results of ', I5, ' timesteps, extending back to ', F8.3, ' Ma.')") high_time, highest_t_Ma CALL Pause() !DIMENSION all needed arrays: ALLOCATE ( c_reference(c_rst_count) ) ALLOCATE ( c_location(4, c_rst_count) ) ALLOCATE ( c_code(c_rst_count) ) ALLOCATE ( c_length_Now_km(c_rst_count) ) ALLOCATE ( c_length_Was_km(c_rst_count) ) ALLOCATE ( c_restoration_sigma_km(c_rst_count) ) ALLOCATE ( c_MaxAge_Ma(c_rst_count) ) ALLOCATE ( c_MinAge_Ma(c_rst_count) ) ALLOCATE ( c_goal_stretch_km(c_rst_count) ) ALLOCATE ( c_goal_rate_mmpa(high_time, c_rst_count) ) ALLOCATE ( c_model_rate_mmpa(high_time, c_rst_count) ) !Initialize all arrays (to simplify debugging): c_reference = ' ' c_location = 0.0D0 c_code = ' ' c_length_Now_km = 0.0D0 c_length_Was_km = 0.0D0 c_restoration_sigma_km = 0.0D0 c_MaxAge_Ma = 0.0D0 c_MinAge_Ma = 0.0D0 c_goal_stretch_km = 0.0D0 c_goal_rate_mmpa = 0.0D0 c_model_rate_mmpa = 0.0D0 !Re-read the input p.rst file, and memorize contents: OPEN (UNIT = 1, FILE = TRIM(c_rst_filename), STATUS = "OLD", IOSTAT = ios) READ (1, "(A)") c_rst_format READ (1, "(A)") c_rst_titles i = 0 ! index of site, 1:c_rst_count. Incremented whenever a header-line is read. memorizing: DO READ (1, "(A)", IOSTAT = ios) line IF (ios /= 0) EXIT memorizing IF (line(1:1) == '+') THEN ! present location !(no action needed in this case) ELSE IF (line(1:1) == '*') THEN ! rate of progress toward stretch in one timestep, in mm/a = km/m.y. line(1:1) = ' ' READ (line, *) t1, t2, t3, t4 j = NINT(t2 / delta_t_Ma) c_model_rate_mmpa(j, i) = t3 c_goal_rate_mmpa(j, i) = t4 ELSE ! header line with original INPUT line of the INPUT c_rst file READ (line, c_rst_format, IOSTAT = ios) c47_citation, (vector(j), j = 1, 4), c5_code, (vector(j), j = 5, 9) IF (ios /= 0) EXIT memorizing ! perhaps input file has empty lines at end? i = i + 1 c_reference(i) = TRIM(c47_citation) c_location(1:4, i) = vector(1:4) c_code(i) = c5_code c_length_Now_km(i) = vector(5) c_length_Was_km(i) = vector(6) c_restoration_sigma_km(i) = vector(7) c_MaxAge_Ma(i) = vector(8) c_MinAge_Ma(i) = vector(9) END IF END DO memorizing CLOSE (1) ! c_rst input file (which was OUTPUT from Restore); now DONE with this. !Create extra arrays to describe Summary results. Note that MANY of these are incomplete/prorated. !Therefore we have to have logically distinct, differently-named arrays, even though the values !may be the same for those sites where the simulation is complete. ALLOCATE ( summary_class012(c_rst_count) ) ALLOCATE ( summary_goal_stretch_km(c_rst_count) ) ALLOCATE ( summary_goal_sigma_km(c_rst_count) ) ALLOCATE ( summary_model_stretch_km(c_rst_count) ) summary_class012 = -1 summary_goal_stretch_km = 0.0D0 summary_goal_sigma_km = 0.0D0 summary_model_stretch_km = 0.0D0 !Fill in values of all output arrays: DO i = 1, c_rst_count IF (c_MaxAge_Ma(i) <= highest_t_Ma) THEN summary_class012(i) = 2 ! COMPLETE ELSE IF (c_MinAge_Ma(i) < highest_t_Ma) THEN summary_class012(i) = 1 ! INCOMPLETE ELSE summary_class012(i) = 0 ! WAITING END IF summary_goal_stretch_km(i) = 0.0D0 ! just initializing, before sum over timesteps... summary_model_stretch_km(i) = 0.0D0 ! just initializing, before sum over timesteps... DO j = 1, high_time summary_goal_stretch_km(i) = summary_goal_stretch_km(i) + delta_t_Ma * c_goal_rate_mmpa(j, i) summary_model_stretch_km(i) = summary_model_stretch_km(i) + delta_t_Ma * c_model_rate_mmpa(j, i) END DO IF (summary_class012(i) == 2) THEN ! COMPLETE; so, use original sigmas: summary_goal_sigma_km(i) = c_restoration_sigma_km(i) ELSE IF (summary_class012(i) == 1) THEN ! INCOMPLETE; use pro-rated sigmas: overlap_in_my = MIN((highest_t_Ma - 0.0D0), & ! This is the time-duration we have modeled (within the cross-section time bracket): & (c_MaxAge_Ma(i) - c_MinAge_Ma(i)), & & (highest_t_Ma - c_MinAge_Ma(i))) duration_in_my = c_MaxAge_Ma(i) - c_MinAge_Ma(i) IF (duration_in_my > 0.0D0) THEN fraction = overlap_in_my / duration_in_my ELSE ! bad data??? fraction = 1.0D0 ENDIF summary_goal_sigma_km(i) = fraction * c_restoration_sigma_km(i) ELSE ! (not really important what we do here) summary_goal_sigma_km(i) = c_restoration_sigma_km(i) END IF END DO !OUTPUT all results, using tabs for ease in opening in Excel, and grouping COMPLETE vs. INCOMPLETE vs. WAITING cross-sections: output_filename = TRIM(c_rst_filename) length = LEN_TRIM(c_rst_filename) output_filename((length-3):(length-3)) = '_' ! changing the '.' of ".rst" to the '_' of "_rst" output_filename = TRIM(output_filename) // "_Summary.txt" OPEN (UNIT = 2, FILE = TRIM(output_filename)) ! unconditional OPEN; overwrites any existing file !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - WRITE (2, "('COMPLETE:')") tab = CHAR(9) ! from ASCII table line = "Reference" // tab // "W.Lon" // tab // "W.Lat" // tab // "E.Lon" // tab // "E.Lat" // tab // "Code" // tab // & & "Now" // tab // "Was" // tab // "Sigma" // tab // "Max.Age" // tab // "Min.Age" // tab // & & tab // & ! (extra tab to separate input from output sections of table(s) & "S_goal" // tab // "S_sigma" // tab // "S_model" WRITE (2, "(A)") TRIM(line) DO i = 1, c_rst_count IF (summary_class012(i) == 2) THEN ! COMPLETE line = TRIM(c_reference(i)) WRITE (c8, "(F8.2)") c_location(1, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c8)) WRITE (c7, "(F7.2)") c_location(2, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c7)) WRITE (c8, "(F8.2)") c_location(3, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c8)) WRITE (c7, "(F7.2)") c_location(4, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c7)) line = TRIM(line) // tab // TRIM(ADJUSTL(c_code(i))) WRITE (c6, "(F6.1)") c_length_Now_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") c_length_Was_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") c_restoration_sigma_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c10, "(F10.1)") c_MaxAge_Ma(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c10)) WRITE (c10, "(F10.1)") c_MinAge_Ma(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c10)) line = TRIM(line) // tab ! extra tab to separate 2 groups of columns (data, vs. summary) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (c6, "(F6.1)") summary_goal_stretch_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") summary_goal_sigma_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") summary_model_stretch_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (2, "(A)") TRIM(line) END IF END DO !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - WRITE (2, "('INCOMPLETE:')") tab = CHAR(9) ! from ASCII table line = "Reference" // tab // "W.Lon" // tab // "W.Lat" // tab // "E.Lon" // tab // "E.Lat" // tab // "Code" // tab // & & "Now" // tab // "Was" // tab // "Sigma" // tab // "Max.Age" // tab // "Min.Age" // tab // & & tab // & ! (extra tab to separate input from output sections of table(s) & "S_goal" // tab // "S_sigma" // tab // "S_model" WRITE (2, "(A)") TRIM(line) DO i = 1, c_rst_count IF (summary_class012(i) == 1) THEN ! INCOMPLETE line = TRIM(c_reference(i)) WRITE (c8, "(F8.2)") c_location(1, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c8)) WRITE (c7, "(F7.2)") c_location(2, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c7)) WRITE (c8, "(F8.2)") c_location(3, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c8)) WRITE (c7, "(F7.2)") c_location(4, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c7)) line = TRIM(line) // tab // TRIM(ADJUSTL(c_code(i))) WRITE (c6, "(F6.1)") c_length_Now_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") c_length_Was_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") c_restoration_sigma_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c10, "(F10.1)") c_MaxAge_Ma(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c10)) WRITE (c10, "(F10.1)") c_MinAge_Ma(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c10)) line = TRIM(line) // tab ! extra tab to separate 2 groups of columns (data, vs. summary) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (c6, "(F6.1)") summary_goal_stretch_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") summary_goal_sigma_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") summary_model_stretch_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (2, "(A)") TRIM(line) END IF END DO !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - WRITE (2, "('WAITING:')") tab = CHAR(9) ! from ASCII table line = "Reference" // tab // "W.Lon" // tab // "W.Lat" // tab // "E.Lon" // tab // "E.Lat" // tab // "Code" // tab // & & "Now" // tab // "Was" // tab // "Sigma" // tab // "Max.Age" // tab // "Min.Age" // tab // & & tab // & ! (extra tab to separate input from output sections of table(s) & "S_goal" // tab // "S_sigma" // tab // "S_model" WRITE (2, "(A)") TRIM(line) DO i = 1, c_rst_count IF (summary_class012(i) == 0) THEN ! WAITING line = TRIM(c_reference(i)) WRITE (c8, "(F8.2)") c_location(1, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c8)) WRITE (c7, "(F7.2)") c_location(2, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c7)) WRITE (c8, "(F8.2)") c_location(3, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c8)) WRITE (c7, "(F7.2)") c_location(4, i) line = TRIM(line) // tab // TRIM(ADJUSTL(c7)) line = TRIM(line) // tab // TRIM(ADJUSTL(c_code(i))) WRITE (c6, "(F6.1)") c_length_Now_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") c_length_Was_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") c_restoration_sigma_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c10, "(F10.1)") c_MaxAge_Ma(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c10)) WRITE (c10, "(F10.1)") c_MinAge_Ma(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c10)) line = TRIM(line) // tab ! extra tab to separate 2 groups of columns (data, vs. summary) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (c6, "(F6.1)") summary_goal_stretch_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") summary_goal_sigma_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) WRITE (c6, "(F6.1)") summary_model_stretch_km(i) line = TRIM(line) // tab // TRIM(ADJUSTL(c6)) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (2, "(A)") TRIM(line) END IF END DO !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - CLOSE (2) WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause END PROGRAM c_rst_Summary