PROGRAM StrainRates2016 !Computes velocities and (2-D) strain-rates on a fine (e.g., 0.01-degree) grid !based on a preferred NeoKinema model whose .feg file covers most/all of this area. !Output is provided for 3 velocity/strain-rate fields on this grid: ! (1) Long-term-average ! (2) Mean-coseismic ! (3) Interseismic = (Long-term-average) - (Mean-coseismic) !For all grid points that fall within a faulted element of the NeoKinema .feg, !the NeoKinema elements are used to interpolate NeoKinema nodal velocities, !so strain-rates are uniform in small triangles, and there are no singularities. !However, for grid points that fall in un-faulted NeoKinema elements, the !NeoKinema interpolation of NK nodal velocities is only used for field (1); !for field (2) the values are computed freshly from dislocation solutions !applied to NeoKinema fault patches, at the grid receiver points. !This better represents the expected arctan-function behavior away from faults. !The difference between StrainRates2016 and the older StrainRates2015 is that !SR2016 adds triangular dislocation patches (Nikkhoo & Walter, 2015) in gores !created where the traces of non-vertical NeoKinema faults have a bend between !two consecutive digitizing steps. (This is like the correction that was added !in NeoKinema_v5.0, so for full self-consistency, StrainRates2016 should only !be applied to solutions coming from NeoKinema_v5.0+.) !For each field, output is provided in each of 2 formats: ! (A) Table file, like that provided to David Sandwell and SCEC CGM. Each line has: ! longitude latitude V_East V_North e^dot_EW e^dot_NS e^dot_NE ! where units are degrees-East, degrees-North, m/s, and /s. ! Grid points that fall outside the NeoKinema .feg area are just omitted. ! (B) Velocities in NeoKinema output format, ! allowing both velocities and strain-rates to be plotted by NeoKineMap. ! Note that a regularly-spaced .feg called "StrainRates2016.feg" is also ! written to make this possible; change the p[token].nki file to refer ! to this new .feg file. Note that strain-rates plotted in this way ! will NOT be exactly equal to, or colocated with, those given in (A). !The user is asked to select the grid limits and the NeoKinema model !to be displayed. The suggested default answers produce a model of strain-rates !in California based on NeoKinema model NSHM-WUS_2013001. !By Peter Bird ! Department of Earth, Planetary, and Space Sciences ! University of California Los Angeles ! October-November 2015 (StrainRates2015) & November 2016 (StrainRates2016) ! & May 2021 (providing for case of NK model with no faults). USE DSphere ! in Peter Bird's file DSphere.f90 USE DDislocation ! in Peter Bird's file DDislocation.f90 USE DTriangular ! in Peter Bird's file DTriangular.f90 IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------------------------------- !Defining the output grids (which should mostly lie within the NK model area!): REAL*8 :: d_lon = 0.01, d_lat = 0.01 REAL*8 :: lon_min = -125.00, lon_max = -114.00 ! Cape Mendocino to the NV-UT border REAL*8 :: lat_min = 31.00, lat_max = 41.00 ! northernmost Mexico to Cape Mendocino ! N.B. The NSHM-WUS_2013001 NeoKinema model would not show accurate ! interseismic strain-rates in Cascadia, because it was computed ! without including the Cascadia Trench interplate megathrust, and after ! Cascadian elastic strain-rates were removed from the GPS field by Rob McCaffrey. INTEGER :: rows, columns ! to be computed based on the above REAL*8, DIMENSION(:,:), ALLOCATABLE :: longTermAverage_V_E, longTermAverage_V_N, longTermAverage_eDot_EW, longTermAverage_eDot_NS, longTermAverage_eDot_NE, & & meanCoseismic_V_E, meanCoseismic_V_N, meanCoseismic_eDot_EW, meanCoseismic_eDot_NS, meanCoseismic_eDot_NE, & & interseismic_V_E, interseismic_V_N, interseismic_eDot_EW, interseismic_eDot_NS, interseismic_eDot_NE !Note that computing and writing the 15 matrices above is the main new task of this utility program. LOGICAL, DIMENSION(:,:), ALLOCATABLE :: withinNKfeg, inNKfaultZone ! <-- important logical flags, determining: whether output is written, and how computed INTEGER, DIMENSION(:,:), ALLOCATABLE :: grid_nodeNumber ! <-- after considering that grid points outside the NK feg won't get any node number. !----------------------------------------------------------------------------------------------------------------------------------------------------------- !NeoKinema variables and arrays: CHARACTER*132 :: path_in = ' ' CHARACTER*132 :: parameter_file = "p_NSHM-WUS_2013001.nki" CHARACTER*132 :: parameter_pathfile CHARACTER*132 :: token INTEGER :: n_refine REAL*8 :: mu_, xi_ REAL*8 :: sigma_offnormal_degrees REAL*8 :: m_per_km = 1000.0D0 REAL*8 :: R ! Earth radius, in m REAL*8 :: locking_depth_m_min, locking_depth_m_max, locking_depth_m_subduction_min, locking_depth_m_subduction_max LOGICAL :: faults_give_sigma_1h, floating_frame, conservative_geodetic_adjustment CHARACTER*132 :: e_file, f_dig, f_nki, f_nko, s_dat, gps_file, gp2_file, x_feg, x_bcs, feg_file, vel_file CHARACTER*132 :: e_pathfile, f_dig_pathfile, f_nki_pathfile, f_nko_pathfile, s_dat_pathfile, gps_pathfile, gp2_pathfile, & & x_feg_pathfile, x_bcs_pathfile, feg_pathfile, vel_pathfile REAL*8 :: geodesy_weight CHARACTER*2 :: reference_plate_c2 REAL*8 :: L0, A0 INTEGER :: stress_interpolation_method INTEGER :: numnod, numel, nfl REAL*8, DIMENSION(3) :: uvec1 REAL*8, DIMENSION(:,:), ALLOCATABLE :: node_uvec INTEGER, DIMENSION(:,:), ALLOCATABLE :: nodes REAL*8, DIMENSION(:), ALLOCATABLE :: vw, vw_interseismic, vw_meanCoseismic ! vw = "vw_longTermAverage" LOGICAL, DIMENSION(:), ALLOCATABLE :: cracking INTEGER :: f_dig_count, f_highest, read_status CHARACTER*50 :: c50 REAL*8, DIMENSION(:), ALLOCATABLE :: f_dip_degrees REAL*8, DIMENSION(:,:), ALLOCATABLE :: trace INTEGER, DIMENSION(:,:), ALLOCATABLE :: trace_loc LOGICAL :: in_trace, got_point, got_index REAL*8 :: t1, t2, dip_degrees INTEGER :: trace_index, start_loc, beyond_loc, end_loc, j1, j2, internal_ios REAL*8, DIMENSION(3) :: uvec CHARACTER*132 :: f_nki_format, f_nko_format INTEGER :: fault_count CHARACTER*6 :: c6 CHARACTER*1 :: first_byte INTEGER, DIMENSION(:), ALLOCATABLE :: f_trace INTEGER :: f_nko_count CHARACTER*1, DIMENSION(:), ALLOCATABLE :: f_c1 REAL*8, DIMENSION(:), ALLOCATABLE :: f_locking_depth_m_max REAL*8, DIMENSION(:), ALLOCATABLE :: f_model_offset_rate LOGICAL :: creeping, faults_were_modeled, faults_were_omitted REAL*8 :: model_offset_rate_mmpa CHARACTER*1 :: c1 CHARACTER*4 :: c4 REAL*8, PARAMETER :: s_per_year = 365.25D0 * 24.0D0 * 60.0D0 * 60.0D0 REAL*8, PARAMETER :: normal_dip_degrees = 55.0D0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity REAL*8, PARAMETER :: thrust_dip_degrees = 20.0D0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity REAL*8, PARAMETER :: subduction_dip_degrees = 14.0D0 ! consistent with Bird & Kagan [2004] Table 5 value used in Long_Term_Seismicity LOGICAL*1, DIMENSION(:), ALLOCATABLE :: f_creeps !----------------------------------------------------------------------------------------------------------------------------------------------------------- !GENERAL UTILITY VARIABLES for this conversion program: CHARACTER*132 :: line, NK_vOut_title1, NK_vOut_title2, NK_vOut_title3 INTEGER :: a, b, element, grid_elements, grid_numnod, i, ios, j, k, l_, m, m1, m2, n1, n2, n3, node INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor LOGICAL :: folding, n1_in, n2_in, n3_in, problem REAL*8 :: argume, & & base_azimuth_radians, bphi, btheta, & & dipf, dot_Phi, dot_Theta, Ds, duphi, duthet, dx, dy, dV_E_over_dy, dV_N_over_dx, & & fphi, ftheta, & & half_R2, heave_azimuth_radians, heave_rate_mps, & & latitude, lf, LLxKm, longitude, & & P1_depth_m, P2_depth_m, P3_depth_m, phi_, prebend_trace_azim, postbend_trace_azim, Poisson, & & radius, rate_mps, & & s1, s2, s3, SIGN, Ss, & & test_depth_m, theta_, Ts, & & ULxKm, uEast, uRadial, uSouth, & & v_E, v_N, v_phi, v_theta, & & wedge, zbot, ztop REAL*8, DIMENSION(3) :: eps_dot, mid_uvec, omega_uvec, P1_uvec, P2_uvec, P3_uvec, Phi, & & r_, sliprate_mps_x1x2x3, test_uvec, Theta, tv, tv1, tv2, tv3, v1 REAL*8, DIMENSION(:), ALLOCATABLE :: a_ REAL*8, DIMENSION(:,:), ALLOCATABLE :: center REAL*8, DIMENSION (3,2,2) :: G REAL*8, DIMENSION (3,2,2,2) :: dG INTEGER :: f_index, percent_done, old_percent_done !GPBtop !----------------------------------------------------------------------------------------------------------------------------------------------------------- WRITE (*, *) WRITE (*, "(' PROGRAM StrainRates2016:')") WRITE (*, *) WRITE (*, "(' Computes velocities and (2-D) strain-rates on a fine (e.g., 0.01-degree) grid')") WRITE (*, "(' based on a NeoKinema model whose .feg file covers most/all of this area.')") WRITE (*, "(' Output is provided for 3 velocity/strain-rate fields on this grid:')") WRITE (*, "(' (1) Long-term-average')") WRITE (*, "(' (2) Mean-coseismic')") WRITE (*, "(' (3) Interseismic = (Long-term-average) - (Mean-coseismic)')") WRITE (*, "(' For all grid points that fall within a faulted element of the NeoKinema .feg,')") WRITE (*, "(' the NeoKinema elements are used to interpolate NeoKinema nodal velocities,')") WRITE (*, "(' so strain-rates are uniform in small triangles, and there are no singularities.')") WRITE (*, "(' However, for grid points that fall in un-faulted NeoKinema elements, the')") WRITE (*, "(' NeoKinema interpolation of NK nodal velocities is only used for field (1);')") WRITE (*, "(' for field (2) the values are computed freshly from dislocation solutions')") WRITE (*, "(' applied to NeoKinema fault patches, at the grid receiver points.')") WRITE (*, "(' This better represents the expected arctan-function behavior away from faults.')") WRITE (*, *) CALL Pause() WRITE (*, *) WRITE (*, "(' The difference between StrainRates2016 and the older StrainRates2015 is that')") WRITE (*, "(' SR2016 adds triangular dislocation patches (Nikkhoo & Walter, 2015) in gores')") WRITE (*, "(' created where the traces of non-vertical NeoKinema faults have a bend between')") WRITE (*, "(' two consecutive digitizing steps. (This is like the correction that was added')") WRITE (*, "(' in NeoKinema_v5.0, so for full self-consistency, StrainRates2016 should only ')") WRITE (*, "(' be applied to solutions coming from NeoKinema_v5.0+.)')") WRITE (*, *) CALL Pause() WRITE (*, *) WRITE (*, "(' For each field, output is provided in each of 2 formats:')") WRITE (*, "(' (A) Table file, like that provided to David Sandwell and SCEC CGM.')") WRITE (*, "(' Each line has:')") WRITE (*, "(' longitude latitude V_East V_North e^dot_EW e^dot_NS e^dot_NE')") WRITE (*, "(' where units are degrees-East, degrees-North, m/s, and /s.')") WRITE (*, "(' Grid points that fall outside the NeoKinema .feg area are just omitted.')") WRITE (*, "(' (B) Velocities in NeoKinema output format, ')") WRITE (*, "(' allowing both velocities and strain-rates to be plotted by NeoKineMap.')") WRITE (*, "(' Note that a regularly-spaced .feg called ""StrainRates2016.feg"" is also')") WRITE (*, "(' written to make this possible; change the p[token].nki file to refer')") WRITE (*, "(' to this new .feg file. Note that strain-rates plotted in this way')") WRITE (*, "(' will NOT be exactly equal to, or colocated with, those given in (A).')") WRITE (*, *) CALL Pause() WRITE (*, *) WRITE (*, "(' The user is asked to select the grid limits and the NeoKinema model')") WRITE (*, "(' to be displayed. The suggested default answers produce a model of strain-rates')") WRITE (*, "(' in California based on NeoKinema model NSHM-WUS_2013001.')") WRITE (*, *) WRITE (*, "(' By Peter Bird')") WRITE (*, "(' Department of Earth, Planetary, and Space Sciences')") WRITE (*, "(' University of California Los Angeles')") WRITE (*, "(' October-November 2015 (StrainRates2015) & November 2016 (StrainRates2016)')") WRITE (*, "(' & May 2021 (providing for case of NK model with no faults)')") WRITE (*, *) CALL Pause() !begin log-file: OPEN (UNIT = 60, FILE = "StrainRates2016-log.txt") ! unconditional OPEN; overwrites any existing WRITE (60, "('PROGRAM StrainRates2016:')") WRITE (60, *) WRITE (60, "('Computes velocities and (2-D) strain-rates on a fine (e.g., 0.01-degree) grid')") WRITE (60, "('based on a preferred NeoKinema model whose .feg file covers most/all of this area.')") WRITE (60, "('Output is provided for 3 velocity/strain-rate fields on this grid:')") WRITE (60, "(' (1) Long-term-average')") WRITE (60, "(' (2) Mean-coseismic')") WRITE (60, "(' (3) Interseismic = (Long-term-average) - (Mean-coseismic)')") WRITE (60, "('For all grid points that fall within a faulted element of the NeoKinema .feg,')") WRITE (60, "('the NeoKinema elements are used to interpolate NeoKinema nodal velocities,')") WRITE (60, "('so strain-rates are uniform in small triangles, and there are no singularities.')") WRITE (60, "('However, for grid points that fall in un-faulted NeoKinema elements, the')") WRITE (60, "('NeoKinema interpolation of NK nodal velocities is only used for field (1);')") WRITE (60, "('for field (2) the values are computed freshly from dislocation solutions')") WRITE (60, "('applied to NeoKinema fault patches, at the grid receiver points.')") WRITE (60, "('This better represents the expected arctan-function behavior away from faults.')") WRITE (60, *) WRITE (60, "('The difference between StrainRates2016 and the older StrainRates2015 is that')") WRITE (60, "('SR2016 adds triangular dislocation patches (Nikkhoo & Walter, 2015) in gores')") WRITE (60, "('created where the traces of non-vertical NeoKinema faults have a bend between')") WRITE (60, "('two consecutive digitizing steps. (This is like the correction that was added')") WRITE (60, "('in NeoKinema_v5.0, so for full self-consistency, StrainRates2016 should only ')") WRITE (60, "('be applied to solutions coming from NeoKinema_v5.0+.)')") WRITE (60, *) WRITE (60, "('For each field, output is provided in each of 2 formats:')") WRITE (60, "(' (A) Table file, like that provided to David Sandwell and SCEC CGM.')") WRITE (60, "(' Each line has:')") WRITE (60, "(' longitude latitude V_East V_North e^dot_EW e^dot_NS e^dot_NE')") WRITE (60, "(' where units are degrees-East, degrees-North, m/s, and /s.')") WRITE (60, "(' Grid points that fall outside the NeoKinema .feg area are just omitted.')") WRITE (60, "(' (B) Velocities in NeoKinema output format, ')") WRITE (60, "(' allowing both velocities and strain-rates to be plotted by NeoKineMap.')") WRITE (60, "(' Note that a regularly-spaced .feg called ""StrainRates2016.feg"" is also')") WRITE (60, "(' written to make this possible; change the p[token].nki file to refer')") WRITE (60, "(' to this new .feg file. Note that strain-rates plotted in this way')") WRITE (60, "(' will NOT be exactly equal to, or colocated with, those given in (A).')") WRITE (60, *) WRITE (60, "('The user is asked to select the grid limits and the NeoKinema model')") WRITE (60, "('to be displayed. The suggested default answers produce a model of strain-rates')") WRITE (60, "('in California based on NeoKinema model NSHM-WUS_2013001.')") WRITE (60, *) WRITE (60, "('By Peter Bird')") WRITE (60, "(' Department of Earth, Planetary, and Space Sciences')") WRITE (60, "(' University of California Los Angeles')") WRITE (60, "(' October-November 2015 (StrainRates2015) & November 2016 (StrainRates2016)E')") WRITE (60, "(' & May 2021 (providing for case of NK model with no faults)')") WRITE (60, *) !Obtain NeoKinema parameter filename and grid definition from user: WRITE (*, *) WRITE (60, *) CALL DPrompt_for_String('Which parameter file (p_.nki) describes the NeoKinema model?', parameter_file, parameter_file) WRITE (60, "('Which parameter file (p_.nki) describes the NeoKinema model?: ',A)") TRIM(parameter_file) CALL DPrompt_for_Real('What is the E-W spacing (d_lon) of the output grid?', d_lon, d_lon) WRITE (60, "('What is the E-W spacing (d_lon) of the output grid?: ', F6.3)") d_lon CALL DPrompt_for_Real('What is the N-S spacing (d_lat) of the output grid?', d_lat, d_lat) WRITE (60, "('What is the N-S spacing (d_lat) of the output grid?: ', F6.3)") d_lat CALL DPrompt_for_Real('What is the longitude of the West side of the grid?', lon_min, lon_min) WRITE (60, "('What is the longitude of the West side of the grid?: ', F8.3)") lon_min CALL DPrompt_for_Real('What is the longitude of the East side of the grid?', lon_max, lon_max) WRITE (60, "('What is the longitude of the East side of the grid?: ', F8.3)") lon_max CALL DPrompt_for_Real('What is the latitude of the South side of the grid?', lat_min, lat_min) WRITE (60, "('What is the latitude of the South side of the grid?: ', F8.3)") lat_min CALL DPrompt_for_Real('What is the latitude of the North side of the grid?', lat_max, lat_max) WRITE (60, "('What is the latitude of the North side of the grid?: ', F8.3)") lat_max !compute grid dimensions and allocate & initialize all grid-type working and output arrays: rows = 1 + NINT((lat_max - lat_min) / d_lat) columns = 1 + NINT((lon_max - lon_min) / d_lon) WRITE (*, *) WRITE (*, "(' Grid spacing is d_lon = ',F6.3,' degrees and d_lat = ',F6.3,' degrees.')") d_lon, d_lat WRITE (*, "(' Grid extends from ',F8.3,'N to ',F8.3,'N and has ',I5,' rows.')") lat_min, lat_max, rows WRITE (*, "(' Grid extends from ',F8.3,'E to ',F8.3,'E and has ',I5,' columns.')") lon_min, lon_max, columns WRITE (60, *) WRITE (60, "('Grid spacing is d_lon = ',F6.3,' degrees and d_lat = ',F6.3,' degrees.')") d_lon, d_lat WRITE (60, "('Grid extends from ',F8.3,'N to ',F8.3,'N and has ',I5,' rows.')") lat_min, lat_max, rows WRITE (60, "('Grid extends from ',F8.3,'E to ',F8.3,'E and has ',I5,' columns.')") lon_min, lon_max, columns ALLOCATE ( withinNKfeg(0:(rows+1), 0:(columns+1)) ) ! N.B. Buffer rows and columns are not used; they are just there to prevent subscript-out-of-range abend. withinNKfeg = .FALSE. ALLOCATE ( inNKfaultZone(rows, columns) ) inNKfaultZone = .FALSE. ALLOCATE ( grid_nodeNumber(rows, columns) ) grid_nodeNumber = 0 ALLOCATE ( longTermAverage_V_E(rows, columns) ) longTermAverage_V_E = 0.0D0 ALLOCATE ( longTermAverage_V_N(rows, columns) ) longTermAverage_V_N = 0.0D0 ALLOCATE ( longTermAverage_eDot_EW(rows, columns) ) longTermAverage_eDot_EW = 0.0D0 ALLOCATE ( longTermAverage_eDot_NS(rows, columns) ) longTermAverage_eDot_NS = 0.0D0 ALLOCATE ( longTermAverage_eDot_NE(rows, columns) ) longTermAverage_eDot_NE = 0.0D0 ALLOCATE ( meanCoseismic_V_E(rows, columns) ) meanCoseismic_V_E = 0.0D0 ALLOCATE ( meanCoseismic_V_N(rows, columns) ) meanCoseismic_V_N = 0.0D0 ALLOCATE ( meanCoseismic_eDot_EW(rows, columns) ) meanCoseismic_eDot_EW = 0.0D0 ALLOCATE ( meanCoseismic_eDot_NS(rows, columns) ) meanCoseismic_eDot_NS = 0.0D0 ALLOCATE ( meanCoseismic_eDot_NE(rows, columns) ) meanCoseismic_eDot_NE = 0.0D0 ALLOCATE ( interseismic_V_E(rows, columns) ) interseismic_V_E = 0.0D0 ALLOCATE ( interseismic_V_N(rows, columns) ) interseismic_V_N = 0.0D0 ALLOCATE ( interseismic_eDot_EW(rows, columns) ) interseismic_eDot_EW = 0.0D0 ALLOCATE ( interseismic_eDot_NS(rows, columns) ) interseismic_eDot_NS = 0.0D0 ALLOCATE ( interseismic_eDot_NE(rows, columns) ) interseismic_eDot_NE = 0.0D0 !=============== BEGIN READING THE NEOKINEMA MODEL: ======================================== CALL Get_Parameters() ! using path_in, parameter_file values defined at top in declaration statements WRITE (*, *) WRITE (*, "(' NeoKinema input parameters were read from ',A)") TRIM(parameter_file) WRITE (60, *) WRITE (60, "('NeoKinema input parameters were read from ',A)") TRIM(parameter_file) CALL Pause() !read and memorize the NeoKinema finite-element grid: feg_file = x_feg feg_pathfile = x_feg_pathfile WRITE (*, *) WRITE (*, "(' Attempting to read finite-element grid file ',A)") TRIM(feg_file) WRITE (*, *) WRITE (60, *) WRITE (60, "('Attempting to read finite-element grid file ',A)") TRIM(feg_file) OPEN (UNIT = 21, FILE = feg_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') problem = (ios /= 0) READ (21, *, IOSTAT = ios) ! title line problem = problem .OR. (ios /= 0) READ (21, *, IOSTAT = ios) numnod problem = problem .OR. (ios /= 0) ALLOCATE ( node_uvec(3,numnod) ) DO i = 1, numnod READ (21, *, IOSTAT = ios) j, longitude, latitude problem = problem .OR. (ios /= 0) .OR. (j /= i) CALL DLonLat_2_Uvec (longitude, latitude, uvec1) node_uvec(1:3, i) = uvec1(1:3) END DO ! i = 1, numnod READ (21, *, IOSTAT = ios) numel problem = problem .OR. (ios /= 0) ALLOCATE ( nodes(3,numel) ) DO i = 1, numel READ (21, *, IOSTAT = ios) j, nodes(1,i), nodes(2,i), nodes(3,i) problem = problem .OR. (ios /= 0) END DO ! i = 1, numel IF (problem) THEN WRITE (*,"(' ERROR: Necessary information absent or defective in this file.')") CLOSE (21) CALL Pause() END IF READ (21, *) nfl IF (nfl > 0) CALL No_Fault_Elements_Allowed() CLOSE (21) WRITE (*, "(' Finite-element grid file ',A,' was successfully read.')") TRIM(feg_file) WRITE (60, "('Finite-element grid file ',A,' was successfully read.')") TRIM(feg_file) !read the primary NeoKinema v[token].out solution file: vel_file = 'v' // TRIM(token) // ".out" vel_pathfile = TRIM(path_in)//TRIM(vel_file) WRITE (*, *) WRITE (*, "(' Attempting to read NeoKinema velocity solution file ',A)") TRIM(vel_file) WRITE (*, *) WRITE (60, *) WRITE (60, "('Attempting to read NeoKinema velocity solution file ',A)") TRIM(vel_file) OPEN(UNIT = 22, FILE = vel_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (22) CALL Pause() END IF READ (22, "(A)") NK_vOut_title1 DO i = 2, 3 READ (22,"(A)") line END DO ! first 3 lines of velocity file ALLOCATE ( vw(2*numnod) ) READ (22,*) (vw(i), i = 1, (2*numnod)) CLOSE (22) WRITE (*, "(' NeoKinema velocity solution file ',A,' was successfully read.')") TRIM(vel_file) WRITE (60, "('NeoKinema velocity solution file ',A,' was successfully read.')") TRIM(vel_file) !also read which NeoKinema finite-elements include active faults: e_file = 'e' // TRIM(token) // ".nko" e_pathfile = TRIM(path_in)//TRIM(e_file) WRITE (*, *) WRITE (*, "(' Attempting to read NeoKinema element strain-rate file ',A)") TRIM(e_file) WRITE (*, *) WRITE (60, *) WRITE (60, "('Attempting to read NeoKinema element strain-rate file ',A)") TRIM(e_file) OPEN(UNIT = 22, FILE = e_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (22) CALL Pause() END IF ALLOCATE ( cracking(numel) ) DO i = 1, numel READ (22, *) j, cracking(j) END DO CLOSE (22) WRITE (*, "(' NeoKinema element strain-rate file ',A,' was successfully read.')") TRIM(e_file) WRITE (60, "('NeoKinema element strain-rate file ',A,' was successfully read.')") TRIM(e_file) !read NeoKinema input file f.dig (if used): IF (faults_were_modeled) THEN ! normal case, where an f_dig file should be read: WRITE (*, *) WRITE (*, "(' Attempting to read NeoKinema fault-traces file ',A)") TRIM(f_dig) WRITE (*, *) WRITE (60, *) WRITE (60, "('Attempting to read NeoKinema fault-traces file ',A)") TRIM(f_dig) f_dig_pathfile = TRIM(path_in)//TRIM(f_dig) OPEN(UNIT = 21, FILE = f_dig_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN ! f*.dig not found WRITE (*,"(' ERROR: File not found.')") CLOSE (21) CALL Pause() STOP ELSE ! Scan for total length and highest trace number: f_dig_count = 0 ! total number of digitized points in the file f_highest = 0 ! highest trace index, e.g. 2096 from F2096N loop_thru: DO READ (21, "(A)", IOSTAT = read_status) c50 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN f_dig_count = f_dig_count + 1 ELSE IF ((c50(1:1) == 'F') .OR. (c50(1:1) == 'f')) THEN READ (c50,"(1X,I4)", IOSTAT = read_status) i IF (read_status == 0) THEN f_highest = MAX (f_highest, i) END IF END IF ELSE; EXIT loop_thru; END IF END DO loop_thru CLOSE (21) ! (will be re-read) !allocate memory for any dip_degrees found (or 0 otherwise): ALLOCATE ( f_dip_degrees(f_highest) ) f_dip_degrees = 0.0D0 ! whole array; this is a flag showing that no dip was found in f*.dig !allocate arrays to hold list of uvecs, and pointers to them: ALLOCATE ( trace(3, f_dig_count) ) ALLOCATE ( trace_loc(2, f_highest) ) trace_loc = 0 ! whole array; if any trace has no points (or only 1) this will be a warning ! re-read f_dig file and memorize contents: OPEN (UNIT = 21, FILE = f_dig_pathfile, STATUS = "OLD", ACTION = "READ", PAD = "YES") in_trace = .FALSE. i = 0 read_dig: DO READ (21, "(A)", IOSTAT = read_status) c50 IF (read_status == 0) THEN ! read was successful IF ((c50(1:2) == ' +') .OR. (c50(1:2) == ' -')) THEN got_point = .TRUE. got_index = .FALSE. READ (c50,*) t1, t2 ! E longitude, N latitude IF (ABS(t2) > 90.001D0) THEN PRINT "(' Bad latitude ',F10.2,' in ',A)", t2, TRIM(f_dig) CALL Pause() STOP END IF ELSE IF ((c50(1:1) == 'F') .OR. (c50(1:1) == 'f')) THEN got_point = .FALSE. READ (c50,"(1X,I4)") j2 ! new trace number (NOTE: Any slip-sense bytes are ignored.) got_index = .TRUE. trace_index = j2 ELSE IF (c50(1:3) == "***") THEN ! '*** end of line segment ***' got_point = .FALSE. got_index = .FALSE. ELSE ! might be a dip_degrees line; try to read it: start_loc = INDEX(c50, "dip_degrees") IF (start_loc > 0) THEN ! found "dip_degrees" in line beyond_loc = start_loc + 11 ! first byte which is not part of "dip_degrees" end_loc = LEN_TRIM(c50) c50 = c50(beyond_loc:end_loc) READ (c50, *, IOSTAT = internal_ios) dip_degrees IF (internal_ios == 0) f_dip_degrees(trace_index) = dip_degrees END IF ! found "dip_degrees" in c50 END IF ! different types of line within f*.dig ELSE; EXIT read_dig; END IF ! terminate file-reading if another line is not available IF (in_trace) THEN IF (got_point) THEN i = i + 1 CALL DLonLat_2_Uvec(t1, t2, uvec) trace(1:3, i) = uvec(1:3) trace_loc(2, j1) = i ELSE ! *** end ... in_trace = .FALSE. END IF ELSE ! (not in_trace) IF (got_index) THEN j1 = j2 ! new index becomes current index ELSE IF (got_point) THEN i = i + 1 CALL DLonLat_2_Uvec(t1, t2, uvec) trace(1:3, i) = uvec(1:3) trace_loc(1, j1) = i in_trace = .TRUE. ELSE ! *** end ... END IF END IF ! (in_trace) END DO read_dig CLOSE (21) ! close f_dig for the last time END IF ! f_dig not found, or WAS found WRITE (*, "(' NeoKinema fault-traces file ',A,' was successfully read.')") TRIM(f_dig) WRITE (60, "('NeoKinema fault-traces file ',A,' was successfully read.')") TRIM(f_dig) ELSE ! faults_were_omitted f_dig_count = 0 ! total number of digitized points in the file f_highest = 0 ! highest trace index, e.g. 2096 from F2096N !allocate memory for any dip_degrees found (or 0 otherwise): ALLOCATE ( f_dip_degrees(99) ) ! arbitrary (small) size; should never be used(?) f_dip_degrees = 0.0D0 ! whole array; this is a flag showing that no dip was found in f*.dig !allocate arrays to hold list of uvecs, and pointers to them: ALLOCATE ( trace(3, 99) ) ! arbitrary (small) size; should never be used(?) ALLOCATE ( trace_loc(2, 99) ) ! arbitrary (small) size; should never be used(?) trace_loc = 0 ! whole array; if any trace has no points (or only 1) this will be a warning WRITE (*, "(' NO NeoKinema fault-traces file was provided.')") WRITE (60, "('NO NeoKinema fault-traces file was provided.')") END IF ! f_dig was used, or not !Read f_nki file (if used): IF (faults_were_modeled) THEN ! normal case, where a f_nki file should be read: ALLOCATE ( f_locking_depth_m_max(f_highest) ) f_locking_depth_m_max = 0.0D0 ! whole vector !read the f_nki input file to get maximum locking depths: WRITE (*, *) WRITE (*, "(' Attempting to read NeoKinema fault data input file ',A)") TRIM(f_nki) WRITE (*, *) WRITE (60, *) WRITE (60, "('Attempting to read NeoKinema fault data input file ',A)") TRIM(f_nki) f_nki_pathfile = TRIM(path_in)//TRIM(f_nki) OPEN(UNIT = 22, FILE = f_nki_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (22) CALL Pause() STOP END IF OPEN(UNIT = 22, FILE = f_nki_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') READ (22, "(A)") f_nki_format READ (22, *) ! f_nki header line read_f_nki: DO ! reading to end of file READ (22, f_nki_format, IOSTAT = ios) c6, c50, t1, t2, creeping, ULxKm, LLxKm ! c50, t1, t2 and creeping are not used IF (ios == -1) EXIT read_f_nki ! EOF c4 = c6(2:5) ! integer part READ (c4,"(I4)") j IF (LLxKm < 0.0D0) LLxKm = locking_depth_m_max / m_per_km f_locking_depth_m_max(j) = LLxKm * m_per_km END DO read_f_nki CLOSE (22) WRITE (*, "(' NeoKinema fault data input file ',A,' was successfully read.')") TRIM(f_nki) WRITE (60, "('NeoKinema fault data input file ',A,' was successfully read.')") TRIM(f_nki) ELSE ! faults_were_omitted, and thus no file should be read ALLOCATE ( f_locking_depth_m_max(99) ) ! arbitrary (small) size; should never be used(?) f_locking_depth_m_max = 0.0D0 ! whole vector WRITE (*, "(' NO NeoKinema fault data input file was provided.')") WRITE (60, "('NO NeoKinema fault data input file was provided.')") END IF ! f_nki was READ, or not !read the (trace-mean) offset-rates from the NeoKinema solution (if faults were used): IF (faults_were_modeled) THEN ! normal case; should read f_nko file: f_nko = 'f' // TRIM(token) // ".nko" WRITE (*, *) WRITE (*, "(' Attempting to read NeoKinema fault offset-rate file ',A)") TRIM(f_nko) WRITE (*, *) WRITE (60, *) WRITE (60, "('Attempting to read NeoKinema fault offset-rate file ',A)") TRIM(f_nko) f_nko_pathfile = TRIM(path_in)//TRIM(f_nko) OPEN(UNIT = 22, FILE = f_nko_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') IF (ios /= 0) THEN WRITE (*,"(' ERROR: File not found.')") CLOSE (22) CALL Pause() STOP END IF !Read through once to find the number of fault offset-rate (component) lines: READ (22, "(A)") f_nko_format READ (22, *) ! f_nko header line f_nko_count = 0 DO READ (22, f_nko_format, IOSTAT = ios) c6 IF (ios == -1) EXIT ! EOF first_byte = c6(1:1) IF ((first_byte == 'F').OR.(first_byte == 'f')) THEN f_nko_count = f_nko_count + 1 END IF END DO CLOSE (22) ALLOCATE ( f_trace(f_nko_count) ) ! integer from F1234D string (e.g., 1234) ALLOCATE ( f_c1(f_nko_count) ) ! offset-type character: R, L, T, P, S, N, D ALLOCATE ( f_model_offset_rate(f_nko_count) ) ! heave or throw component in m/s ALLOCATE ( f_creeps(f_nko_count) ) !Second reading is to remember heave rates (some converted from throw rates); may be negative at this point OPEN(UNIT = 22, FILE = f_nko_pathfile, STATUS = 'OLD', IOSTAT = ios, PAD = 'YES') READ (22, "(A)") f_nko_format READ (22, *) ! f_nko header line read_f_nko: DO i = 1, f_nko_count ! reading to end of file READ (22, f_nko_format, IOSTAT = ios) c6, c50, t1, t2, creeping, model_offset_rate_mmpa ! c50, t1, and t2 are not used IF (ios == -1) EXIT read_f_nko ! EOF c4 = c6(2:5) ! integer part READ (c4,"(I4)") j f_trace(i) = j f_c1(i) = c6(6:6) f_model_offset_rate(i) = model_offset_rate_mmpa * 1.0D-3 / s_per_year ! heave or throw component in m/s f_creeps(i) = creeping END DO read_f_nko CLOSE (22) WRITE (*, "(' NeoKinema fault offset-rate file ',A,' was successfully read.')") TRIM(f_nko) WRITE (60, "('NeoKinema fault offset-rate file ',A,' was successfully read.')") TRIM(f_nko) ELSE ! no faults were used in the NeoKinema model! f_nko_count = 0 ALLOCATE ( f_trace(99) ) ! arbitrary (small) number; should never be used(?) ALLOCATE ( f_c1(99) ) ! ditto ALLOCATE ( f_model_offset_rate(99) ) ! ditto ALLOCATE ( f_creeps(99) ) ! ditto WRITE (*, "(' NO NeoKinema fault offset-rate file was read.')") WRITE (60, "('NO NeoKinema fault offset-rate file was read.')") END IF ! f_nko was read, or not !=============== BEGIN COMPUTATIONS: ======================================================= !Compute vw_meanCoseismic at NeoKinema nodes. ![Note: The NeoKinema output file vw_interseismic_[token].out is NOT used, because that was ! computed using active fault traces broken into segments, which is VERY SLIGHTLY ! DIFFERENT from using active fault traces broken into digitizing steps, as this program does. WRITE (60, *) WRITE (60, "('Computing vw_meanCoseismic at each of ',I8,' NeoKinema nodes...')") numnod ALLOCATE ( vw_meanCoseismic(2*numnod) ) vw_meanCoseismic = 0.0D0 WRITE (*, *) WRITE (*, *) DO i = 1, numnod r_(1:3) = node_uvec(1:3, i) CALL MeanCoseismic_Velocity_at_Point(r_, v_E, v_N) ! included in this file, below "CONTAINS" vw_meanCoseismic(2*i) = +v_E ! v_phi = v_E vw_meanCoseismic(2*i-1) = -v_N ! v_theta = v_S = -v_N WRITE (*, "('+Computing vw_meanCoseismic at each of ',I8,' NeoKinema nodes...',I8)") numnod, i END DO ! i = 1, numnod WRITE (*, "('+Computing vw_meanCoseismic at each of ',I8,' NeoKinema nodes...DONE ')") numnod !Complete characterization of velocities at NeoKinema nodes: WRITE (*, *) WRITE (*, "(' At each NeoKinema node, vw_interseismic = vw - vw_meanCoseismic.')") WRITE (60, *) WRITE (60, "('At each NeoKinema node, vw_interseismic = vw - vw_meanCoseismic.')") ALLOCATE ( vw_interseismic(2*numnod) ) vw_interseismic = vw - vw_meanCoseismic ! whole vector !-------------- For the rectangular output grid: ------------------------------------------- !DO/DO all grid points, and: ! (1) Determine if they fall within the NeoKinema .feg area ! IF SO, ! (2) Determine whether they fall within a faulting element of the NeoKinema .feg ! IF SO, ! (3) Determine all 3 vector velocities and 3 strain-rate tensors from NeoKinema nodal velocities and nodal functions. ! ELSE, ! (4) Compute v_longTermAverage velocity vector & strain-rate tensor from NeoKinema nodal velocities and nodal functions. ! (5) Compute v_meanCoseismic vector from NeoKinema fault traces and offset-rates. ! (6) v_interseismic = v_longTermAverage - v_meanCoseismic ! END IF ! END IF !END DO/END DO !DO/DO all grid points, and: ! (1) Determine if they fall within the NeoKinema .feg area ! IF SO, ! IF (meanCoseismic & interseismic strain-rates need to be determined), ! (7) Compute meanCosesimic strain-rate tensor by F-D on gridded velocities; ! (8) Compute interseismic strain-rate tensor by simple subtraction: interseismic = longTermAverage - meanCoseismic. ! END IF ! END IF !END DO/END DO !prepare (for computing element internal coordinates) by learning the spherical triangle elements: WRITE (*, *) WRITE (*, "(' Learning the NeoKinema spherical-triangle elements...')") WRITE (60, *) WRITE (60, "('Learning the NeoKinema spherical-triangle elements...')") ALLOCATE ( a_(numel) ) ALLOCATE ( neighbor(3, numel) ) ALLOCATE ( center(3, numel) ) half_R2 = R**2 / 2.0D0 CALL Plane_Area(folding) IF (folding) THEN WRITE (*, "(' ERROR: SUBROUTINE Plane_Area reports folding.')") WRITE (60, "('ERROR: SUBROUTINE Plane_Area reports folding.')") CALL Pause() STOP END IF ! find element center points DO l_ = 1, numel v1(1:3) = node_uvec(1:3, nodes(1, l_)) + node_uvec(1:3, nodes(2, l_)) + node_uvec(1:3, nodes(3, l_)) CALL DMake_Uvec(v1, tv) center(1:3, l_) = tv(1:3) END DO ! find neighboring elements on all 3 sides of each neighbor = 0 ! initializing the whole array (3, numel) homes: DO l_ = 1, numel sides: DO j = 1, 3 ! 3 sides k = 1 + MOD (j, 3) a = nodes(k, l_) ! 1st node along side b = nodes(1 + MOD (k, 3), l_) ! 2nd node along side strangers: DO m = 1, numel IF ((nodes(1, m) == b) .AND. (nodes(2, m) == a)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF ((nodes(2, m) == b) .AND. (nodes(3, m) == a)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF ((nodes(3, m) == b) .AND. (nodes(1, m) == a)) THEN neighbor(j, l_) = m EXIT strangers END IF END DO strangers END DO sides END DO homes WRITE (*, *) WRITE (*, "(' Starting main computational double-loop on all grid points...')") WRITE (60, *) WRITE (60, "('Starting main computational double-loop on all grid points...')") WRITE (*, *) old_percent_done = 0 !DO/DO all grid points, and: DO i = 1, rows latitude = lat_max - (i-1) * d_lat DO j = 1, columns longitude = lon_min + (j-1) * d_lon CALL DLonLat_2_Uvec(longitude, latitude, r_) CALL Internal(r_, l_, s1, s2, s3) !(1) Determine if they fall within the NeoKinema .feg area: withinNKfeg(i, j) = (l_ > 0) !IF SO, IF (withinNKfeg(i, j)) THEN !(2) Determine whether they fall within a faulting element of the NeoKinema .feg inNKfaultZone(i, j) = cracking(l_) !IF SO, IF (inNKfaultZone(i, j)) THEN !(3) Determine all 3 vector velocities and 3 strain-rate tensors from NeoKinema nodal velocities and nodal functions. CALL DUvec_2_ThetaPhi(r_, theta_, phi_) CALL Gjxy (l_, r_, G) ! computes nodal function G(j, x, y) at position r_ in element l_. CALL Del_Gjxy_del_thetaphi (l_, r_, dG) ! computes derivatives of nodal function G at position r_ in element l_. n1 = nodes(1, l_); n2 = nodes(2, l_); n3 = nodes(3, l_) !(3.1) longTermAverage velocity vector: v_theta = (vw(2*n1-1) * G(1,1,1)) + (vw(2*n1) * G(1,2,1)) + & & (vw(2*n2-1) * G(2,1,1)) + (vw(2*n2) * G(2,2,1)) + & & (vw(2*n3-1) * G(3,1,1)) + (vw(2*n3) * G(3,2,1)) v_phi = (vw(2*n1-1) * G(1,1,2)) + (vw(2*n1) * G(1,2,2)) + & & (vw(2*n2-1) * G(2,1,2)) + (vw(2*n2) * G(2,2,2)) + & & (vw(2*n3-1) * G(3,1,2)) + (vw(2*n3) * G(3,2,2)) longTermAverage_V_E(i, j) = v_phi longTermAverage_V_N(i, j) = -v_theta !(3.2) meanCoseismic velocity vector: v_theta = (vw_meanCoseismic(2*n1-1) * G(1,1,1)) + (vw_meanCoseismic(2*n1) * G(1,2,1)) + & & (vw_meanCoseismic(2*n2-1) * G(2,1,1)) + (vw_meanCoseismic(2*n2) * G(2,2,1)) + & & (vw_meanCoseismic(2*n3-1) * G(3,1,1)) + (vw_meanCoseismic(2*n3) * G(3,2,1)) v_phi = (vw_meanCoseismic(2*n1-1) * G(1,1,2)) + (vw_meanCoseismic(2*n1) * G(1,2,2)) + & & (vw_meanCoseismic(2*n2-1) * G(2,1,2)) + (vw_meanCoseismic(2*n2) * G(2,2,2)) + & & (vw_meanCoseismic(2*n3-1) * G(3,1,2)) + (vw_meanCoseismic(2*n3) * G(3,2,2)) meanCoseismic_V_E(i, j) = v_phi meanCoseismic_V_N(i, j) = -v_theta !(3.3) interseismic velocity vector: v_theta = (vw_interseismic(2*n1-1) * G(1,1,1)) + (vw_interseismic(2*n1) * G(1,2,1)) + & & (vw_interseismic(2*n2-1) * G(2,1,1)) + (vw_interseismic(2*n2) * G(2,2,1)) + & & (vw_interseismic(2*n3-1) * G(3,1,1)) + (vw_interseismic(2*n3) * G(3,2,1)) v_phi = (vw_interseismic(2*n1-1) * G(1,1,2)) + (vw_interseismic(2*n1) * G(1,2,2)) + & & (vw_interseismic(2*n2-1) * G(2,1,2)) + (vw_interseismic(2*n2) * G(2,2,2)) + & & (vw_interseismic(2*n3-1) * G(3,1,2)) + (vw_interseismic(2*n3) * G(3,2,2)) interseismic_V_E(i, j) = v_phi interseismic_V_N(i, j) = -v_theta !(3.4) longTermAverage strain-rate tensor: CALL E_rate(R, l_, nodes, G, dG, theta_, vw, eps_dot) ! returns values of REAL*8, DIMENSION(3) :: eps_dot ! whose components are: _thetaTheta, _thetaPhi, _phiPhi longTermAverage_eDot_EW(i, j) = +eps_dot(3) longTermAverage_eDot_NS(i, j) = +eps_dot(1) longTermAverage_eDot_NE(i, j) = -eps_dot(2) !(3.5) meanCoseismic strain-rate tensor: CALL E_rate(R, l_, nodes, G, dG, theta_, vw_meanCoseismic, eps_dot) ! returns values of REAL*8, DIMENSION(3) :: eps_dot ! whose components are: _thetaTheta, _thetaPhi, _phiPhi meanCoseismic_eDot_EW(i, j) = +eps_dot(3) meanCoseismic_eDot_NS(i, j) = +eps_dot(1) meanCoseismic_eDot_NE(i, j) = -eps_dot(2) !(3.6) interseismic strain-rate tensor: CALL E_rate(R, l_, nodes, G, dG, theta_, vw_interseismic, eps_dot) ! returns values of REAL*8, DIMENSION(3) :: eps_dot ! whose components are: _thetaTheta, _thetaPhi, _phiPhi interseismic_eDot_EW(i, j) = +eps_dot(3) interseismic_eDot_NS(i, j) = +eps_dot(1) interseismic_eDot_NE(i, j) = -eps_dot(2) ELSE ! .NOT.inNKfaultZone(i, j) !(4) Compute v_longTermAverage velocity vector & strain-rate tensor from NeoKinema nodal velocities and nodal functions. CALL DUvec_2_ThetaPhi(r_, theta_, phi_) CALL Gjxy (l_, r_, G) ! computes nodal function G(j, x, y) at position r_ in element l_. CALL Del_Gjxy_del_thetaphi (l_, r_, dG) ! computes derivatives of nodal function G at position r_ in element l_. n1 = nodes(1, l_); n2 = nodes(2, l_); n3 = nodes(3, l_) !(4.1) longTermAverage velocity vector: v_theta = (vw(2*n1-1) * G(1,1,1)) + (vw(2*n1) * G(1,2,1)) + & & (vw(2*n2-1) * G(2,1,1)) + (vw(2*n2) * G(2,2,1)) + & & (vw(2*n3-1) * G(3,1,1)) + (vw(2*n3) * G(3,2,1)) v_phi = (vw(2*n1-1) * G(1,1,2)) + (vw(2*n1) * G(1,2,2)) + & & (vw(2*n2-1) * G(2,1,2)) + (vw(2*n2) * G(2,2,2)) + & & (vw(2*n3-1) * G(3,1,2)) + (vw(2*n3) * G(3,2,2)) longTermAverage_V_E(i, j) = v_phi longTermAverage_V_N(i, j) = -v_theta !(4.2) longTermAverage strain-rate tensor: CALL E_rate(R, l_, nodes, G, dG, theta_, vw, eps_dot) ! returns values of REAL*8, DIMENSION(3) :: eps_dot ! whose components are: _thetaTheta, _thetaPhi, _phiPhi longTermAverage_eDot_EW(i, j) = +eps_dot(3) longTermAverage_eDot_NS(i, j) = +eps_dot(1) longTermAverage_eDot_NE(i, j) = -eps_dot(2) !(5) Compute v_meanCoseismic vector from NeoKinema fault traces and offset-rates, using DDislocation. CALL MeanCoseismic_Velocity_at_Point(r_, v_E, v_N) ! included in this file, below "CONTAINS" meanCoseismic_V_E(i, j) = v_E meanCoseismic_V_N(i, j) = v_N !(6) v_interseismic = v_longTermAverage - v_meanCoseismic interseismic_V_E(i, j) = longTermAverage_V_E(i, j) - meanCoseismic_V_E(i, j) interseismic_V_N(i, j) = longTermAverage_V_N(i, j) - meanCoseismic_V_N(i, j) END IF ! inNKfaultZone(i, j), or NOT END IF ! withinNKfeg(i, j). N.B. There is no ELSE; I just leave the velocity and strain-rate as initialized to 0.0D0; they will never be written. percent_done = NINT(100.0D0 * (j + columns*(i-1)) / (rows * columns)) IF (percent_done > old_percent_done) THEN WRITE (*, "('+ main computational double-loop is ',I3,'% done ...')") percent_done old_percent_done = percent_done END IF END DO ! j = 1, columns END DO ! i = 1, rows WRITE (*, *) !DO/DO all grid points (2nd huge double-loop): DO i = 1, rows latitude = lat_max - (i-1) * d_lat DO j = 1, columns longitude = lon_min + (j-1) * d_lon !(1) Determine if they fall within the NeoKinema .feg area: IF (withinNKfeg(i, j)) THEN !IF SO, !IF (meanCoseismic & interseismic strain-rates need to be determined), IF ((meanCoseismic_eDot_EW(i, j) == 0.0D0).AND.(meanCoseismic_eDot_NS(i, j) == 0.0D0).AND.(meanCoseismic_eDot_NE(i, j) == 0.0D0)) THEN !(7) Compute meanCosesimic strain-rate tensor by F-D on gridded velocities: !(7.1) meanCoseismic_eDot_EW: IF ((j == 1).OR.(.NOT. withinNKfeg(i, j-1))) THEN ! problem; no column or no velocity to the left! dx = d_lon * radians_per_degree * R * COS(latitude * radians_per_degree) meanCoseismic_eDot_EW(i, j) = (meanCoseismic_V_E(i, j+1) - meanCoseismic_V_E(i, j)) / dx dV_N_over_dx = (meanCoseismic_V_N(i, j+1) - meanCoseismic_V_N(i, j)) / dx ELSE IF ((j == columns).OR.(.NOT. withinNKfeg(i, j+1))) THEN ! problem; no column or no velocity to the right! dx = d_lon * radians_per_degree * R * COS(latitude * radians_per_degree) meanCoseismic_eDot_EW(i, j) = (meanCoseismic_V_E(i, j) - meanCoseismic_V_E(i, j-1)) / dx dV_N_over_dx = (meanCoseismic_V_N(i, j) - meanCoseismic_V_N(i, j-1)) / dx ELSE ! general case; interior of the grid, with valid velocities both W and E dx = 2.0D0 * d_lon * radians_per_degree * R * COS(latitude * radians_per_degree) meanCoseismic_eDot_EW(i, j) = (meanCoseismic_V_E(i, j+1) - meanCoseismic_V_E(i, j-1)) / dx dV_N_over_dx = (meanCoseismic_V_N(i, j+1) - meanCoseismic_V_N(i, j-1)) / dx END IF !(7.2) meanCoseismic_eDot_NS: IF ((i == 1).OR.(.NOT. withinNKfeg(i-1, j))) THEN ! problem; no column or no velocity above! dy = d_lat * radians_per_degree * R meanCoseismic_eDot_NS(i, j) = (meanCoseismic_V_N(i, j) - meanCoseismic_V_N(i+1, j)) / dy dV_E_over_dy = (meanCoseismic_V_E(i, j) - meanCoseismic_V_E(i+1, j)) / dy ELSE IF ((i == rows).OR.(.NOT. withinNKfeg(i+1, j))) THEN ! problem; no column or no velocity below! dy = d_lat * radians_per_degree * R meanCoseismic_eDot_NS(i, j) = (meanCoseismic_V_N(i-1, j) - meanCoseismic_V_N(i, j)) / dy dV_E_over_dy = (meanCoseismic_V_E(i-1, j) - meanCoseismic_V_E(i, j)) / dy ELSE ! general case; interior of the grid, with valid velocities both above and below dy = 2.0D0 * d_lat * radians_per_degree * R meanCoseismic_eDot_NS(i, j) = (meanCoseismic_V_N(i-1, j) - meanCoseismic_V_N(i+1, j)) / dy dV_E_over_dy = (meanCoseismic_V_E(i-1, j) - meanCoseismic_V_E(i+1, j)) / dy END IF !(7.3) meanCoseismic_eDot_NE: meanCoseismic_eDot_NE(i, j) = (dV_N_over_dx + dV_E_over_dy) / 2.0D0 !(8) Compute interseismic strain-rate tensor by simple subtraction: interseismic = longTermAverage - meanCoseismic. interseismic_eDot_EW(i, j) = longTermAverage_eDot_EW(i, j) - meanCoseismic_eDot_EW(i, j) interseismic_eDot_NS(i, j) = longTermAverage_eDot_NS(i, j) - meanCoseismic_eDot_NS(i, j) interseismic_eDot_NE(i, j) = longTermAverage_eDot_NE(i, j) - meanCoseismic_eDot_NE(i, j) END IF ! meanCoseismic & interseismic strain-rates need to be determined END IF ! withinNKfeg(i, j) END DO ! j = 1, columns END DO ! i = 1, rows !=============== BEGIN OUTPUT SECTION: ===================================================== !grid-type .feg file, needed for plotting gridded velocities/strain-rates with NeoKineMap: OPEN (UNIT = 1, FILE = "StrainRates2016.feg") ! unconditional OPEN; overwrites any existing WRITE (1, "('StrainRates2016.feg:')") !count and number the output nodes that are within the NeoKinema feg area: grid_numnod = 0 DO i = 1, rows DO j = 1, columns IF (withinNKfeg(i, j)) THEN grid_numnod = grid_numnod + 1 grid_nodeNumber(i, j) = grid_numnod END IF END DO END DO !write out those nodes that are defined: WRITE (1, "(2I8,' 0 9000000 T (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)')") grid_numnod, grid_numnod DO i = 1, rows latitude = lat_max - (i-1) * d_lat DO j = 1, columns longitude = lon_min + (j-1) * d_lon IF (withinNKfeg(i, j)) THEN WRITE (1, "(I8,1X,F8.3,1X,F7.3,' 0.00E+00 0.00E+00 0.00E+00 0.00E+00')") grid_nodeNumber(i, j), longitude, latitude END IF END DO ! j = 1, columns END DO ! i = 1, rows WRITE (*, *) WRITE (*, "(' StrainRates2016.feg contains ',I8,' nodes.')") grid_numnod WRITE (60, *) WRITE (60, "('StrainRates2016.feg contains ',I8,' nodes.')") grid_numnod !count the output elements that are entirely within the NK feg area: grid_elements = 0 DO i = 1, (rows-1) DO j = 1, (columns-1) !lower-left or SW element, starting with top/N node: n1_in = withinNKfeg(i, j) n2_in = withinNKfeg(i+1, j) n3_in = withinNKfeg(i+1, j+1) IF (n1_in.AND.n2_in.AND.n3_in) grid_elements = grid_elements + 1 !upper-right or NE element, starting with bottom/S node: n1_in = withinNKfeg(i+1, j+1) n2_in = withinNKfeg(i, j+1) n3_in = withinNKfeg(i, j) IF (n1_in.AND.n2_in.AND.n3_in) grid_elements = grid_elements + 1 END DO ! j = 1, (columns-1) END DO ! i = 1, (rows-1_ !write out the non-empty elements: WRITE (1, "(I8)") grid_elements element = 0 DO i = 1, (rows-1) DO j = 1, (columns-1) !lower-left or SW element, starting with top/N node: n1_in = withinNKfeg(i, j) n2_in = withinNKfeg(i+1, j) n3_in = withinNKfeg(i+1, j+1) IF (n1_in.AND.n2_in.AND.n3_in) THEN element = element + 1 WRITE (1, "(4I8)") element, grid_nodeNumber(i, j), grid_nodeNumber(i+1, j), grid_nodeNumber(i+1, j+1) END IF !upper-right or NE element, starting with bottom/S node: n1_in = withinNKfeg(i+1, j+1) n2_in = withinNKfeg(i, j+1) n3_in = withinNKfeg(i, j) IF (n1_in.AND.n2_in.AND.n3_in) THEN element = element + 1 WRITE (1, "(4I8)") element, grid_nodeNumber(i+1, j+1), grid_nodeNumber(i, j+1), grid_nodeNumber(i, j) END IF END DO ! j = 1, (columns-1) END DO ! i = 1, (rows-1_ WRITE (*, "(' StrainRates2016.feg contains ',I8,' elements.')") grid_elements WRITE (60, "('StrainRates2016.feg contains ',I8,' elements.')") grid_elements WRITE (1, "(' 0')") ! NFL CLOSE (UNIT = 1) ! StrainRates2016.feg !write out 3 nodal-velocity files, in NeoKinema v[token].out format, referring to StrainRates2016.feg: CALL Write_Nodal_Velocities("longTermAverage", longTermAverage_V_E, longTermAverage_V_N) CALL Write_Nodal_Velocities("meanCoseismic", meanCoseismic_V_E, meanCoseismic_V_N) CALL Write_Nodal_Velocities("interseismic", interseismic_V_E, interseismic_V_N) !write out 3 point-table files, in format designed for use by David Sandwell and the SCEC Community Geodetic Model: CALL Write_GSM_Table("longTermAverage", longTermAverage_V_E, longTermAverage_V_N, longTermAverage_eDot_EW, longTermAverage_eDot_NS, longTermAverage_eDot_NE) CALL Write_GSM_Table("meanCoseismic", meanCoseismic_V_E, meanCoseismic_V_N, meanCoseismic_eDot_EW, meanCoseismic_eDot_NS, meanCoseismic_eDot_NE) CALL Write_GSM_Table("interseismic", interseismic_V_E, interseismic_V_N, interseismic_eDot_EW, interseismic_eDot_NS, interseismic_eDot_NE) !=============== END: ====================================================================== WRITE (60, *) WRITE (60, "('Job completed.')") CLOSE (UNIT = 60) ! log-file WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS SUBROUTINE Bad_Parameters() !prints admonition screen, calls Pause(), and STOP's. IMPLICIT NONE WRITE (*, *) WRITE (*, "(' ==========================================================================')") WRITE (*, "(' ERROR DURING READING OF NeoKinema PARAMETER FILE')") WRITE (*, "(' ')") WRITE (*, "(' Please review the proper format, which is listed in the introductory')") WRITE (*, "(' comment lines of NeoKinema.f90.')") WRITE (*, "(' ')") WRITE (*, "(' This program is smart enough to read either NeoKinema v1.x, OR')") WRITE (*, "(' NeoKinema v2.x~v3.x, OR NeoKinema v4.x flavors of the')") WRITE (*, "(' input parameter file (and decide which is intended),')") WRITE (*, "(' but it cannot correct for other unexpected errors.')") WRITE (*, "(' ')") WRITE (*, "(' After you have read this message, this program will stop.')") WRITE (*, "(' Please correct the input parameter file and re-start it.')") WRITE (*, "(' ==========================================================================')") Call Pause() STOP END SUBROUTINE Bad_Parameters SUBROUTINE Could_Not_Find_File(pathfilename) !prevents multiple duplications of this very simple code: IMPLICIT NONE CHARACTER*(*), INTENT(IN), OPTIONAL :: pathfilename IF (PRESENT(pathfilename)) THEN WRITE (*, "(' ERROR: Could not find a necessary input file:'/' ',A)") TRIM(pathfilename) ELSE WRITE (*, "(' ERROR: Could not find a necessary input file')") END IF CALL Pause() STOP END SUBROUTINE Could_Not_Find_File CHARACTER*10 FUNCTION DASCII10(x) ! Returns a right-justified 10-byte (or shorter) version of a ! floating-point number, in "human" format, with no more ! than 4 significant digits. IMPLICIT NONE REAL*8, INTENT(IN) :: x CHARACTER*10 :: temp10 CHARACTER*20 :: temp20 INTEGER :: j, k1, k10, zeros LOGICAL :: punt REAL*8 :: x_log DOUBLE PRECISION :: y IF (x == 0.0D0) THEN DASCII10=' 0' RETURN ELSE IF (x > 0.0D0) THEN punt = (x >= 999999999.5D0).OR.(x < 0.0000100D0) ELSE ! x < 0.0 punt = (x <= -99999999.5D0).OR.(x > -0.000100D0) END IF IF (punt) THEN ! need exponential notation; use Fortran utility WRITE (temp10,'(1P,D10.3)') x !consider possible improvements, from left to right: IF (temp10(3:6) == '.000') THEN ! right-shift 4 spaces over it temp20(7:10) = temp10(7:10) temp20(5:6) = temp10(1:2) temp20(1:4) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(5:6) == '00') THEN ! right-shift 2 spaces over it temp20(7:10) = temp10(7:10) temp20(3:6) = temp10(1:4) temp20(1:2) = ' ' temp10 = temp20(1:10) ELSE IF (temp10(6:6) == '0') THEN ! right-shift 1 space over it temp20(7:10) = temp10(7:10) temp20(2:6) = temp10(1:5) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(8:8) == '+') THEN ! right-shift over + sign in exponent temp20(9:10) = temp10(9:10) temp20(2:8) = temp10(1:7) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF IF (temp10(9:9) == '0') THEN ! right-shift over leading 0 in exponent temp20(10:10) = temp10(10:10) temp20(2:9) = temp10(1:8) temp20(1:1) = ' ' temp10 = temp20(1:10) END IF DASCII10 = temp10 ELSE ! can represent without exponential notation x_log = DLOG10(ABS(x)) zeros = DInt_Below(x_log) - 3 y = (10.D0**zeros) * NINT(ABS(x) / (10.D0**zeros)) IF (x < 0.0D0) y = -y WRITE (temp20,"(F20.9)") y ! byte 11 is the '.' !Avoid results like "0.7400001" due to rounding error! IF (temp20(19:20) == '01') temp20(19:20) = '00' !Find first important byte from right; change 0 -> ' ' k10 = 10 ! (if no non-0 found to right of .) right_to_left: DO j = 20, 12, -1 IF (temp20(j:j) == '0') THEN temp20(j:j) = ' ' ELSE k10 = j EXIT right_to_left END IF END DO right_to_left !put leading (-)0 before . of fractions, if it fits IF (x > 0.0D0) THEN IF (temp20(10:11) == ' .') temp20(10:11) = '0.' ELSE ! x < 0.0 IF (k10 <= 18) THEN IF (temp20(9:11) == ' -.') temp20(9:11) = '-0.' END IF END IF k1 = k10 - 9 DASCII10 = temp20(k1:k10) END IF ! punt, or not END FUNCTION DASCII10 SUBROUTINE Del_Gjxy_del_thetaphi (l_, r_, dG) IMPLICIT NONE INTEGER, INTENT(IN) :: l_ ! element number REAL*8, DIMENSION(3), INTENT(IN) :: r_ ! position vector DOUBLE PRECISION, DIMENSION (3,2,2,2), INTENT(OUT) :: dG ! computes array of 2 derivitives of each of the 2 components of ! each of the 6 nodal functions for element l_ at ! position r_ (Cartesian unit vector). ! Results are in 1./radian (dimensionless), NOT 1./m or 1./degree ! ! It is user's responsibility that element l_ contains r_. INTEGER, SAVE :: l_last = 0 ! remembers l_ from previous invocation ! Subscripts of dG(j, x, y, m): INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2)? INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function? INTEGER :: m ! 1:2 = theta (S-ward) or phi (E-ward) derivitive? DOUBLE PRECISION, DIMENSION(3,2) :: del_r_ ! theta- and phi-derivitives of r_ (in 3-D) DOUBLE PRECISION, DIMENSION(3,2) :: local ! local Theta, Phi unit vectors at r_ (xyz, SE) DOUBLE PRECISION, DIMENSION(3,2,2) :: del_local ! theta-, phi- derivitives of local DOUBLE PRECISION, DIMENSION(3,3), SAVE :: corner ! positions vector of corner nodes (xyz, 123) DOUBLE PRECISION, DIMENSION(3,3,2), SAVE :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) DOUBLE PRECISION, DIMENSION(3) :: tr_, tv, tvi, tvo, tvp, tv1, tv2, tv3, vfa, vfb ! temporary vector factors DOUBLE PRECISION :: cos_phi, cos_theta, factor, phi, sin_phi, sin_theta INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ DO j = 1, 3 tvi(1:3) = node_uvec(1:3, nodes(j, l_)) ! promote to R8 {in earlier versions of this routine} CALL DMake_Uvec(tvi, tvo) ! and check length corner(1:3, j) = tvo(1:3) CALL DLocal_Theta(tvo, tvp) post(1:3, j, 1) = tvp(1:3) CALL DLocal_Phi(tvo, tvp) post(1:3, j, 2) = tvp(1:3) END DO END IF ! begin calculations which depend on r_ tvi(1:3) = r_(1:3) ! promote to R8 {in earlier versions of this routine} CALL DMake_Uvec(tvi, tr_) ! and check length CALL DLocal_Theta(tr_, tv) local(1:3, 1) = tv(1:3) CALL DLocal_Phi(tr_, tv) local(1:3, 2) = tv(1:3) ! Note: these functions will catch polar points; don't test again phi = DATAN2(tr_(2), tr_(1)) cos_phi = DCOS(phi) sin_phi = DSIN(phi) cos_theta = tr_(3) sin_theta = DSQRT(tr_(1)**2 + tr_(2)**2) del_r_(1:3,1) = local(1:3,1) ! d.r_/d.theta = Theta del_r_(1:3,2) = local(1:3,2) * sin_theta ! d.r_/d.phi = Phi * SIN(theta) del_local(1:3,1,1) = -tr_(1:3) ! d.Theta/d.theta = - r_ del_local(1:3,1,2) = local(1:3,2) * cos_theta ! d.Theta/d.phi = Phi * COS(theta) del_local(1:3,2,1) = (/ 0.0D0, 0.0D0, 0.0D0 /) ! d.Phi/d.theta = 0 del_local(1:3,2,2) = (/ -cos_phi, -sin_phi, 0.0D0 /) ! d.Phi/d.phi = (-COS(phi),-SIN(phi,0) DO j = 1, 3 ! 3 corner nodes of element i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1(1:3) = corner(1:3, i1) tv2(1:3) = corner(1:3, i2) tv3(1:3) = corner(1:3, i3) CALL DCross(tv2, tv3, vfa) factor = 1.0D0 / DDot(tv1, vfa) vfb(1:3) = vfa(1:3) * factor DO x = 1, 2 ! unit velocity at node is S or E DO y = 1, 2 ! S- or E- component of nodal function tv1(1:3) = post(1:3, j, x) tvi(1:3) = local(1:3, y) DO m = 1, 2 ! theta- or phi-derivitive tv(1:3) = del_r_(1:3, m) tvo(1:3) = del_local(1:3, y, m) dG(j, x, y, m) = & & (DDot(tv, vfb) * DDot(tv1, tvi)) + & & (DDot(tr_, vfb) * DDot(tv1, tvo)) END DO END DO END DO END DO END SUBROUTINE Del_Gjxy_del_thetaphi INTEGER FUNCTION DInt_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL*8, INTENT(IN) :: x INTEGER :: i REAL*8 :: y i = INT(x) IF (x >= 0.0D0) THEN DInt_Below = i ELSE ! x < 0. y = 1.0D0*i IF (y <= x) THEN DInt_Below = i ELSE ! most commonly DInt_Below = i - 1 END IF END IF END FUNCTION DInt_Below SUBROUTINE DPrompt_for_Real (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a real value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! Note that prompt_text should usually end with '?'. ! It can be more than 52 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 REAL*8, INTENT(IN) :: default REAL*8, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, point, written LOGICAL :: finished REAL*8 :: trial !------------------------------------------------------------------------------------ !This code worked (provided 4 significant digits), but left unecessary trailing zeros ("20.00"; "6.000E+07") !IF (((ABS(default) >= 0.1).AND.(ABS(default) < 1000.)).OR.(default == 0.0)) THEN ! ! Provide 4 significant digits by using Gxx.4 (the suffix shows significant digits, NOT digits after the decimal point!) ! WRITE (suggested,"(G11.4)") default !ELSE ! ! Use 1P,E because it avoids wasted and irritating leading 0 ("0.123E+4"). ! WRITE (suggested,"(1P,E11.3)") default !END IF !------------------------------------------------------------------------------------ !So I replaced it with the following: !(1) Use ASCII10 to get 4 significant digits (but no unecessary trailing zeroes): suggested = DASCII10(default) !(2) Be sure that the number contains some sign that it is floating-point, not integer: IF (INDEX(suggested, '.') == 0) THEN IF ((INDEX(suggested, 'E') == 0).AND.(INDEX(suggested, 'e') == 0).AND. & & (INDEX(suggested, 'D') == 0).AND.(INDEX(suggested, 'd') == 0)) THEN suggested = ADJUSTL(suggested) point = LEN_TRIM(suggested) + 1 suggested(point:point) = '.' END IF END IF !------------------------------------------------------------------------------------ suggested = ADJUSTL(suggested) bytes = LEN_TRIM(prompt_text) finished = .FALSE. DO WHILE (.NOT. finished) written = 0 DO WHILE ((bytes - written) > 52) blank_at = written + INDEX(prompt_text((written+1):(written+52)),' ',.TRUE.) ! searching L from end for ' ' IF (blank_at < (written + 2)) blank_at = written + 52 WRITE (*,"(' ',A)") prompt_text((written+1):blank_at) written = blank_at END DO WRITE (*,"(' ',A,' [',A']: '\)") prompt_text((written+1):bytes), TRIM(suggested) finished = .TRUE. ! unless changed below READ (*,"(A)") instring IF (LEN_TRIM(instring) == 0) THEN answer = default ELSE !The following lead to occoasional abends !under Digital Visual Fortran 5.0D !(memory violations caught by WinNT): !READ (instring, *, IOSTAT = ios) trial !The following fix leads to a compiler error: !BACKSPACE (*) !READ (*, *, IOSTAT = ios) trial !and the following fix lead to an immediate abend: !BACKSPACE (5) !READ (*, *, IOSTAT = ios) trial !So, I am creating and then reading a dummy file: OPEN (UNIT = 72, FILE = 'trash') WRITE (72, "(A)") instring CLOSE (72) OPEN (UNIT = 72, FILE = 'trash') READ (72, *, IOSTAT = ios) trial CLOSE (72, STATUS = 'DELETE') IF (ios /= 0) THEN ! bad string WRITE (*,"(' ERROR: Your response (',A,') was not recognized.')") TRIM(instring) WRITE (*,"(' Enter an real number using 11 characters (or less).')") WRITE (*,"(' Please try again:')") finished = .FALSE. ELSE answer = trial END IF ! problem with string, or not? END IF ! some bytes were entered END DO ! until finished END SUBROUTINE DPrompt_for_Real SUBROUTINE DPrompt_for_String (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with a character-string value. ! If [Enter] is pressed with nothing preceding it, ! then "answer" takes the value "default". ! 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 CHARACTER*(*), INTENT(IN) :: default CHARACTER*(*), INTENT(OUT) :: answer CHARACTER*80 trial INTEGER :: blank_at, default_bytes, leftover, & & prompt_bytes, written prompt_bytes = LEN_TRIM(prompt_text) default_bytes = LEN_TRIM(default) written = 0 leftover = 79 - prompt_bytes - 4 ! unless changed below DO WHILE ((prompt_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) leftover = 79 - (blank_at - (written+1) + 1) - 4 written = blank_at END DO IF (leftover >= default_bytes) THEN WRITE (*,"(' ',A,' [',A,']')") prompt_text(written+1:prompt_bytes), TRIM(default) ELSE WRITE (*,"(' ',A)") prompt_text(written+1:prompt_bytes) WRITE (*,"(' [',A,']')") TRIM(default) END IF WRITE (*,"(' ?: '\)") READ (*,"(A)") trial IF (LEN_TRIM(trial) == 0) THEN answer = TRIM(default) ELSE answer = TRIM(trial) END IF END SUBROUTINE DPrompt_for_String SUBROUTINE Dumb_s123 (element, vector, s1, s2, s3) ! Finds s1, s2, s3 coordinates of position vector "in" element ! (whether the point is actually in the element or NOT). INTEGER, INTENT(IN) :: element REAL*8, DIMENSION(3), INTENT(IN) :: vector REAL*8, INTENT(OUT) :: s1, s2, s3 INTEGER :: i1, i2, i3 REAL*8, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL*8 :: d1, dc, t IF (element == 0) THEN CALL Prevent('element = 0', 1, 'Dumb_s123') END IF i1 = nodes(1, element) i2 = nodes(2, element) i3 = nodes(3, element) !shorten(?) vector to just touch plane element -> v1 tv1 = center(1:3, element) dc = DDot(vector, tv1) IF (dc <= 0.0D0) CALL Prevent('"Internal" vector >= 90 deg. from element', 1, 'Dumb_s123') tv2 = node_uvec(1:3, i1) d1 = DDot(tv2, tv1) t = d1 / dc v1 = t * vector tvi = node_uvec(1:3,i3) - node_uvec(1:3,i2) tvo = v1(1:3) - node_uvec(1:3,i3) CALL DCross(tvi, tvo, tv) s1 = DDot(tv1, tv) * half_R2 / a_(element) tvi = node_uvec(1:3,i1) - node_uvec(1:3,i3) tvo = v1(1:3) - node_uvec(1:3,i1) CALL DCross(tvi, tvo, tv) s2 = DDot(tv1, tv) * half_R2 / a_(element) s3 = 1.0D0 - s1 - s2 END SUBROUTINE Dumb_s123 SUBROUTINE E_rate(R, l_, nodes, G, dG, theta_, vw, eps_dot) ! Evaluate strain-rate at one point in one spherical continuum element (# l_); ! the specific point is implied by the values of arrays G and dG supplied, ! but note that the value of theta_ must also be consistent. REAL*8, INTENT(IN) :: R ! radius of planet, in m INTEGER, INTENT(IN) :: l_ ! element number INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes DOUBLE PRECISION, DIMENSION(3,2,2), INTENT(IN) :: G ! nodal functions @ selected point DOUBLE PRECISION, DIMENSION(3,2,2,2), INTENT(IN):: dG ! derivitives of nodal functions @ selected point REAL*8, INTENT(IN) :: theta_ ! colatitude, radians DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: vw REAL*8, DIMENSION(3), INTENT(OUT) :: eps_dot DOUBLE PRECISION, DIMENSION(3) :: sums INTEGER :: iv, iw, j DOUBLE PRECISION :: cott, csct, prefix cott = 1.0D0 / DTAN(1.0D0 * theta_) csct = 1.0D0 / DSIN(1.0D0 * theta_) prefix = 1.0D0 / R sums(1:3) = 0.0D0 DO j = 1, 3 iv = 2 * nodes(j, l_) - 1 iw = iv + 1 ! epsilon_dot_sub_theta_theta sums(1) = sums(1) + & & vw(iv) * prefix * dG(j,1,1,1) + & & vw(iw) * prefix * dG(j,2,1,1) ! epsilon_dot_sub_theta_phi sums(2) = sums(2) + & & vw(iv) * prefix * 0.5D0 * (csct * dG(j,1,1,2) + dG(j,1,2,1) - cott * G(j,1,2)) + & & vw(iw) * prefix * 0.5D0 * (csct * dG(j,2,1,2) + dG(j,2,2,1) - cott * G(j,2,2)) ! epsilon_dot_sub_phi_phi sums(3) = sums(3) + & & vw(iv) * prefix * (csct * dG(j,1,2,2) + cott * G(j,1,1)) + & & vw(iw) * prefix * (csct * dG(j,2,2,2) + cott * G(j,2,1)) END DO ! 3 local nodes eps_dot(1:3) = sums(1:3) END SUBROUTINE E_rate SUBROUTINE Get_filename (unit, filename, error) ! Obtains a filename from the beginning of the next line; ! name may be padded with blanks on both sides, and ! may be followed by comments. Note that ! "unit" should have been opened with PAD = "YES". IMPLICIT NONE INTEGER, INTENT (IN) :: unit ! Fortran device number to READ CHARACTER*132, INTENT(OUT) :: filename ! main result of this SUBR LOGICAL, INTENT(OUT) :: error ! IF (ios /= 0) error = .TRUE. CHARACTER(132) :: buffer INTEGER :: i, ios LOGICAL :: past error = .FALSE. ! unless changed below READ (unit,"(A)", IOSTAT = ios) buffer IF (ios /= 0) THEN error = .TRUE. RETURN END IF buffer = ADJUSTL(buffer) ! left-justify past = .FALSE. ! will be T when past end of filename blank_right: DO i=2,132 IF ((buffer(i:i) == ' ') .OR. & (buffer(i:i) == ',') .OR. & (buffer(i:i) == '=') .OR. & (buffer(i:i) == ':')) past = .TRUE. if (past) buffer(i:i) = ' ' END DO blank_right IF (((buffer(1:1) == 'N') .OR. (buffer(1:1) == 'n')) .AND. & ((buffer(2:2) == 'O') .OR. (buffer(2:2) == 'o')) .AND. & ((buffer(3:3) == 'N') .OR. (buffer(3:3) == 'n')) .AND. & ((buffer(4:4) == 'E') .OR. (buffer(4:4) == 'e')) .AND. & (buffer(5:5) == ' ')) buffer = 'none' filename = buffer(1:80) END SUBROUTINE Get_filename SUBROUTINE Get_Parameters ! Reads input parameter file p*.nki ! and memorizes its contents. ! Values reside in global variables. ! This is just a SUBR to avoid repetition of code. ! ! Revised 2014.10.04 to read EITHER NeoKinema v1.x, or NeoKinema v2.x - ! NeoKinema v3.x, OR NeoKinema v.4.x parameter file. ! input parameter files. IMPLICIT NONE CHARACTER*132 :: temp_path_in = ' ' INTEGER :: gps_index_1, gps_index_2, ios, line LOGICAL :: bad_parameter, any_bad_parameters REAL*8 :: t, t1, t2 !NOTE that this routine also uses (and fills) global variables like mu_, xi_, ! sigma_offnormal_degrees, .... temp_path_in = path_in !N.B. parameter_file is a global character variable. !CALL File_List( file_type = "p*.nki", & ! & suggested_file = parameter_file, & ! & using_path = temp_path_in) !CALL DPrompt_for_String('Which parameter file should be used?',parameter_file,parameter_file) parameter_pathfile = TRIM(temp_path_in)//TRIM(parameter_file) !------- Try reading file by assuming that it is in NeoKinema v1.x format: OPEN (UNIT = 1, FILE = parameter_pathfile, STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) IF (ios /= 0) CALL Could_Not_Find_File(parameter_pathfile) any_bad_parameters = .FALSE. ! until an error occurs (below) CALL Get_Filename(1, token, bad_parameter) any_bad_parameters = any_bad_parameters .OR. bad_parameter ! delaying reaction line = 1 ! number of refinements READ (1, *, IOSTAT = ios) n_refine ; line = line + 1 IF (n_refine < 0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! strain-rate uncertainty for rigid blocks READ (1, *, IOSTAT = ios) mu_ ; line = line + 1 IF (mu_ <= 0.D0) bad_parameter = .TRUE. IF (mu_ < DSQRT(1.1D0 * TINY(mu_))) bad_parameter = .TRUE. IF (mu_ > 1.D-10) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! small strain-rate increment (xi_) for imposing stress-directions READ (1, *, IOSTAT = ios) xi_ ; line = line + 1 IF (xi_ <= 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! sigma, in degrees, of angle between [heave vectors of dip-slip faults] & [trace-normal direction] READ (1, *, IOSTAT = ios) sigma_offnormal_degrees ; line = line + 1 IF (sigma_offnormal_degrees < 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! radius of planet READ (1, *, IOSTAT = ios) t ; line = line + 1 IF (t <= 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter R = t * m_per_km ! minimum and maximum locking depths of intraplate faults, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_min = t1 * 1000.0D0 locking_depth_m_max = t2 * 1000.0D0 ! minimum and maximum locking depths of subduction zones, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_subduction_min = t1 * 1000.0D0 locking_depth_m_subduction_max = t2 * 1000.0D0 ! do new active faults count as sigma_1h data? READ (1, *, IOSTAT = ios) faults_give_sigma_1h ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! names of additional input files (or, "none ") CALL Get_filename (1, f_nki, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter faults_were_omitted = (f_nki(1:5) == 'none ').OR.(f_nki(1:5) == 'NONE ').OR.(f_nki(1:5) == 'None ') faults_were_modeled = .NOT.faults_were_omitted IF (faults_were_modeled) THEN CALL Get_filename (1, f_dig, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ELSE ! ignore whatever may be entered in the f_dig filename line READ (1, *) ! read line with f_dig filename, but do not store it line = line + 1 ! count the line, as usual. f_dig = "none" ! replacing whatever this line might say f_nki = "none" ! (standardise the capitalization(?) of this word, to simplify later code) faults_give_sigma_1h = .FALSE. ! (just to be clear) END IF CALL Get_filename (1, s_dat, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gps_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gp2_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ! is velocity reference frame for geodetic data allowed to float? READ (1, *, IOSTAT =ios) floating_frame ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! conservative_geodetic_adjustment? (using geologic slip rates) READ (1, *, IOSTAT = ios) conservative_geodetic_adjustment ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! weight factor for all geodetic data: READ (1, *, IOSTAT = ios) geodesy_weight ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) CALL Get_filename (1, x_feg, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, x_bcs, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter READ (1,"(A)", IOSTAT = ios) reference_plate_c2; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) CLOSE (UNIT = 1) ! close parameter_file IF (.NOT.any_bad_parameters) THEN f_nko_pathfile = TRIM(path_in) // TRIM(f_nko) f_dig_pathfile = TRIM(path_in) // TRIM(f_dig) s_dat_pathfile = TRIM(path_in) // TRIM(s_dat) gps_pathfile = TRIM(path_in) // TRIM(gps_file) x_feg_pathfile = TRIM(path_in) // TRIM(x_feg) x_bcs_pathfile = TRIM(path_in) // TRIM(x_bcs) RETURN END IF !=========== Note: We only reach here if there was an error during reading ! of parameters under assumption of NeoKinema v1.x format. !------- Try reading file by assuming that it is in NeoKinema v2.x~v3.x format, !------- in which case it has new lines added for parameters L0 and A0. OPEN (UNIT = 1, FILE = parameter_pathfile, STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) IF (ios /= 0) CALL Could_Not_Find_File(parameter_pathfile) any_bad_parameters = .FALSE. ! until an error occurs (below) CALL Get_Filename(1, token, bad_parameter) any_bad_parameters = any_bad_parameters .OR. bad_parameter ! delaying reaction line = 1 ! L0 READ (1, *, IOSTAT = ios) L0 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! A0 READ (1, *, IOSTAT = ios) A0 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! number of refinements READ (1, *, IOSTAT = ios) n_refine ; line = line + 1 IF (n_refine < 0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! strain-rate uncertainty for rigid blocks READ (1, *, IOSTAT = ios) mu_ ; line = line + 1 IF (mu_ <= 0.D0) bad_parameter = .TRUE. IF (mu_ < DSQRT(1.1D0 * TINY(mu_))) bad_parameter = .TRUE. IF (mu_ > 1.D-10) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! small strain-rate increment (xi_) for imposing stress-directions READ (1, *, IOSTAT = ios) xi_ ; line = line + 1 IF (xi_ <= 0.) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! sigma, in degrees, of angle between [heave vectors of dip-slip faults] & [trace-normal direction] READ (1, *, IOSTAT = ios) sigma_offnormal_degrees ; line = line + 1 IF (sigma_offnormal_degrees < 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! radius of planet READ (1, *, IOSTAT = ios) t ; line = line + 1 IF (t <= 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter R = t * m_per_km ! minimum and maximum locking depths of intraplate faults, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_min = t1 * 1000.0D0 locking_depth_m_max = t2 * 1000.0D0 ! minimum and maximum locking depths of subduction zones, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_subduction_min = t1 * 1000.0D0 locking_depth_m_subduction_max = t2 * 1000.0D0 ! do new active faults count as sigma_1h data? READ (1, *, IOSTAT = ios) faults_give_sigma_1h ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! names of additional input files (or, "none ") CALL Get_filename (1, f_nki, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter faults_were_omitted = (f_nki(1:5) == 'none ').OR.(f_nki(1:5) == 'NONE ').OR.(f_nki(1:5) == 'None ') faults_were_modeled = .NOT.faults_were_omitted IF (faults_were_modeled) THEN CALL Get_filename (1, f_dig, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ELSE ! ignore whatever may be entered in the f_dig filename line READ (1, *) ! read line with f_dig filename, but do not store it line = line + 1 ! count the line, as usual. f_dig = "none" ! replacing whatever this line might say f_nki = "none" ! (standardise the capitalization(?) of this word, to simplify later code) faults_give_sigma_1h = .FALSE. ! (just to be clear) END IF CALL Get_filename (1, s_dat, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gps_file, bad_parameter) ; line = line + 1 !Protect against reading "INTEGER :: stress_interpolation_method" of a NKv4 parameter file; !do this by searching for required filename extension ".gps". gps_index_1 = INDEX(gps_file, ".gps") gps_index_2 = INDEX(gps_file, ".GPS") bad_parameter = ((gps_index_1 == 0).AND.(gps_index_2 == 0)) any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gp2_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ! is velocity reference frame for geodetic data allowed to float? READ (1, *, IOSTAT =ios) floating_frame ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! conservative_geodetic_adjustment? (using geologic slip rates) READ (1, *, IOSTAT = ios) conservative_geodetic_adjustment ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! weight factor for all geodetic data: geodesy_weight = 1.0D0 ! by definition in NeoKinema v2.x CALL Get_filename (1, x_feg, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, x_bcs, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter READ (1,"(A)", IOSTAT = ios) reference_plate_c2; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) CLOSE (UNIT = 1) ! close parameter_file IF (.NOT.any_bad_parameters) THEN f_nko_pathfile = TRIM(path_in) // TRIM(f_nko) f_dig_pathfile = TRIM(path_in) // TRIM(f_dig) s_dat_pathfile = TRIM(path_in) // TRIM(s_dat) gps_pathfile = TRIM(path_in) // TRIM(gps_file) x_feg_pathfile = TRIM(path_in) // TRIM(x_feg) x_bcs_pathfile = TRIM(path_in) // TRIM(x_bcs) RETURN END IF !=========== Note: We only reach here if there was an error during reading ! of parameters under assumption of NeoKinema v1.x format, AND ! under assumption of NeoKinema v2.x~3.x format. !------- Try reading file by assuming that it is in NeoKinema v4.x format, !------- in which case it has ANOTHER new line added for parameter "stress_interpolation_method" !------- coming after "s_dat" and before "gps_file". OPEN (UNIT = 1, FILE = parameter_pathfile, STATUS = "OLD", ACTION = "READ", PAD = "YES", IOSTAT = ios) IF (ios /= 0) CALL Could_Not_Find_File(parameter_pathfile) any_bad_parameters = .FALSE. ! until an error occurs (below) CALL Get_Filename(1, token, bad_parameter) any_bad_parameters = any_bad_parameters .OR. bad_parameter ! delaying reaction line = 1 ! L0 READ (1, *, IOSTAT = ios) L0 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! A0 READ (1, *, IOSTAT = ios) A0 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! number of refinements READ (1, *, IOSTAT = ios) n_refine ; line = line + 1 IF (n_refine < 0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! strain-rate uncertainty for rigid blocks READ (1, *, IOSTAT = ios) mu_ ; line = line + 1 IF (mu_ <= 0.D0) bad_parameter = .TRUE. IF (mu_ < DSQRT(1.1D0 * TINY(mu_))) bad_parameter = .TRUE. IF (mu_ > 1.D-10) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! small strain-rate increment (xi_) for imposing stress-directions READ (1, *, IOSTAT = ios) xi_ ; line = line + 1 IF (xi_ <= 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! sigma, in degrees, of angle between [heave vectors of dip-slip faults] & [trace-normal direction] READ (1, *, IOSTAT = ios) sigma_offnormal_degrees ; line = line + 1 IF (sigma_offnormal_degrees < 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter ! radius of planet READ (1, *, IOSTAT = ios) t ; line = line + 1 IF (t <= 0.D0) bad_parameter = .TRUE. any_bad_parameters = any_bad_parameters .OR. (ios /= 0) .OR. bad_parameter R = t * m_per_km ! minimum and maximum locking depths of intraplate faults, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_min = t1 * 1000.0D0 locking_depth_m_max = t2 * 1000.0D0 ! minimum and maximum locking depths of subduction zones, in km READ (1, *, IOSTAT = ios) t1, t2 ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) locking_depth_m_subduction_min = t1 * 1000.0D0 locking_depth_m_subduction_max = t2 * 1000.0D0 ! do new active faults count as sigma_1h data? READ (1, *, IOSTAT = ios) faults_give_sigma_1h ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! names of additional input files (or, "none ") CALL Get_filename (1, f_nki, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter faults_were_omitted = (f_nki(1:5) == 'none ').OR.(f_nki(1:5) == 'NONE ').OR.(f_nki(1:5) == 'None ') faults_were_modeled = .NOT.faults_were_omitted IF (faults_were_modeled) THEN CALL Get_filename (1, f_dig, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ELSE ! ignore whatever may be entered in the f_dig filename line READ (1, *) ! read line with f_dig filename, but do not store it line = line + 1 ! count the line, as usual. f_dig = "none" ! replacing whatever this line might say f_nki = "none" ! (standardise the capitalization(?) of this word, to simplify later code) faults_give_sigma_1h = .FALSE. ! (just to be clear) END IF CALL Get_filename (1, s_dat, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ! stress interpolation method (INTEGER index): line = line + 1 READ (1, *, IOSTAT = ios) stress_interpolation_method bad_parameter = (ios /= 0) any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gps_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, gp2_file, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter ! is velocity reference frame for geodetic data allowed to float? READ (1, *, IOSTAT =ios) floating_frame ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! conservative_geodetic_adjustment? (using geologic slip rates) READ (1, *, IOSTAT = ios) conservative_geodetic_adjustment ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) ! weight factor for all geodetic data: geodesy_weight = 1.0D0 ! by definition in NeoKinema v2.x and beyond CALL Get_filename (1, x_feg, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter CALL Get_filename (1, x_bcs, bad_parameter) ; line = line + 1 any_bad_parameters = any_bad_parameters .OR. bad_parameter READ (1,"(A)", IOSTAT = ios) reference_plate_c2; line = line + 1 any_bad_parameters = any_bad_parameters .OR. (ios /= 0) CLOSE (UNIT = 1) ! close parameter_file IF (.NOT.any_bad_parameters) THEN f_nko_pathfile = TRIM(path_in) // TRIM(f_nko) f_dig_pathfile = TRIM(path_in) // TRIM(f_dig) s_dat_pathfile = TRIM(path_in) // TRIM(s_dat) gps_pathfile = TRIM(path_in) // TRIM(gps_file) x_feg_pathfile = TRIM(path_in) // TRIM(x_feg) x_bcs_pathfile = TRIM(path_in) // TRIM(x_bcs) RETURN END IF !=========== Note: We only reach here if there was an error during reading ! of parameters under assumptions of NeoKinema v1.x format, ! NeoKinema v2.x~3,x, and NeoKinema v4.x format! CALL Bad_Parameters() END SUBROUTINE Get_Parameters SUBROUTINE Gjxy (l_, r_, G) IMPLICIT NONE INTEGER, INTENT(IN) :: l_ ! element number REAL*8, DIMENSION(3), INTENT(IN) :: r_ ! position vector DOUBLE PRECISION, DIMENSION (3,2,2), INTENT(OUT) :: G ! computes matrix of 6 vector nodal functions for element l_ at ! position r_ (Cartesian unit vector). ! It is user's responsibility that element l_ contains r_. INTEGER, SAVE :: l_last = 0 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2) INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function DOUBLE PRECISION, DIMENSION(3,2) :: local ! local unit vectors at r_ (xyz, SE) DOUBLE PRECISION, DIMENSION(3,3), SAVE :: corner ! positions vector of corner nodes (xyz, 123) DOUBLE PRECISION, DIMENSION(3,3,2), SAVE :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) DOUBLE PRECISION, DIMENSION(3) :: tr_, tvi, tvo, tvp, tv1, tv2, tv3, vf ! temporary vector factor DOUBLE PRECISION :: f_sup_j ! as in Kong and Bird (1995) [ j == k ] INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ DO j = 1, 3 tvi(1:3) = node_uvec(1:3, nodes(j, l_)) CALL DMake_Uvec(tvi, tvo) corner(1:3, j) = tvo(1:3) CALL DLocal_Theta(tvo, tvp) post(1:3, j, 1) = tvp(1:3) CALL DLocal_Phi(tvo, tvp) post(1:3, j, 2) = tvp(1:3) END DO END IF ! begin computations which depend on r_ tvi(1:3) = r_(1:3) CALL DMake_Uvec(tvi, tr_) CALL DLocal_Theta(tr_, tvo) local(1:3,1) = tvo(1:3) CALL DLocal_Phi(tr_, tvo) local(1:3,2) = tvo(1:3) DO j = 1, 3 i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1(1:3) = corner(1:3, i1) tv2(1:3) = corner(1:3, i2) tv3(1:3) = corner(1:3, i3) CALL DCross(tv2, tv3, vf) f_sup_j = DDot(tr_, vf) / DDot(tv1, vf) DO x = 1, 2 tv1(1:3) = post(1:3, j, x) DO y = 1, 2 tv2(1:3) = local(1:3, y) G(j, x, y) = f_sup_j * DDot(tv1, tv2) END DO END DO END DO END SUBROUTINE Gjxy SUBROUTINE Internal (b_, iele, s1, s2, s3) REAL*8, DIMENSION(3), INTENT(IN) :: b_ INTEGER, INTENT(OUT) :: iele REAL*8, INTENT(OUT) :: s1, s2, s3 ! Input is a Cartesian unit vector in the unit sphere. ! Output is element number and s1, s2, s3 internal coordinates of the ! plane triangle finite element. ! Returns iele = 0 if no element contains the point b_. INTEGER :: attempts, best_iele_yet, i, iet, k, l_ INTEGER, PARAMETER :: history_length = 100 INTEGER, DIMENSION(history_length) :: iele_history LOGICAL :: in_loop REAL*8 :: best_minimum_yet, r2, r2min, s1t, s2t, s3t REAL*8, DIMENSION(3) :: s_temp, tv ! establish defaults (not found) in case of quick exit iele = 0 ! default s1 = 0.0D0; s2 = 0.0D0; s3 = 0.0D0 ! default !find closest element center to initialize search r2min = 4.01D0 DO l_ = 1, numel r2 = (b_(1) - center(1,l_))**2 +(b_(2) - center(2,l_))**2 +(b_(3) - center(3,l_))**2 IF (r2 < r2min) THEN r2min = r2 iet = l_ END IF END DO ! If closest element center is more than 1 radian away, give up. tv = center(1:3, iet) IF (DDot(b_, tv) < 0.5403D0) RETURN ! initialize search memory attempts = 0 best_iele_yet = 0 best_minimum_yet = -999.9D0 !iele_history = 0 ! zeros the whole long vector; good for debugging, but not necessary for production! is_it_here: DO ! first, check for infinite loop, of any peridocity! in_loop = .FALSE. ! but, this may change below... IF (attempts >= 3) THEN ! need to check for a loop retrospection: DO k = (attempts - 1), 1, -1 IF (iele_history(k) == iele_history(attempts)) THEN in_loop = .TRUE. EXIT retrospection END IF ! found a problem END DO retrospection ! k = (attempts - 1), 1, -1 (looking backward) END IF ! any loop is possible yet IF (in_loop) THEN iet = best_iele_yet CALL Dumb_s123 (iet, b_, s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! not in any loop; ! normal operation; search, and grow the history attempts = attempts + 1 IF (attempts > history_length) THEN WRITE (*, *) WRITE (*, "(' ERROR: PARAMETER history_length in SUBROUTINE Internal is too small.')") WRITE (21, *) WRITE (21, "('ERROR: PARAMETER history_length in SUBROUTINE Internal is too small.')") CALL DTraceback END IF iele_history(attempts) = iet CALL Dumb_s123 (iet, b_, s1t, s2t, s3t) IF ((s1t <= s2t) .AND. (s1t <= s3t)) THEN ! s1t is most negative; therefore, most critical IF (s1t > best_minimum_yet) THEN best_iele_yet = iet best_minimum_yet = s1t END IF ! new best solution found IF (s1t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t <= s1t) .AND. (s2t <= s3t)) THEN ! s2t is most negative; therefore, most critical IF (s2t > best_minimum_yet) THEN best_iele_yet = iet best_minimum_yet = s2t END IF ! new best solution found IF (s2t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3t is most negative; therefore, most critical IF (s3t > best_minimum_yet) THEN best_iele_yet = iet best_minimum_yet = s3t END IF ! new best solution found IF (s3t >= 0.0D0) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF END IF END IF ! in/not in a loop END DO is_it_here ! successful completion iele = iet s1 = s1t s2 = s2t s3 = s3t END SUBROUTINE Internal SUBROUTINE MeanCoseismic_Velocity_at_Point(r_, v_E, v_N) !code from main program made into a SUBR to avoid repetition (and possible inconsistency). IMPLICIT NONE ! but global variables and arrays are freely accessed REAL*8, DIMENSION(3), INTENT(IN) :: r_ ! uvec of location where velocity is needed REAL*8, INTENT(OUT) :: v_E, v_N ! outputs, in m/s CHARACTER*1 :: c1 INTEGER :: f_index, k, m, m1, m2 REAL*8 :: argume, bphi, btheta, dip_degrees, dipf, dot_Phi, dot_Theta, fphi, ftheta, lf, phi_, radius, rate_mps, theta_, wedge, zbot, ztop REAL*8, DIMENSION(3) :: mid_uvec, Phi, sliprate_mps_x1x2x3, Theta, tv !(5) Compute v_meanCoseismic horizontal (2-component) vector from NeoKinema fault traces and offset-rates, ! using MODULEs DDislocation and DTriangular: v_E = 0.0D0 v_N = 0.0D0 ! initializing, before sum DO k = 1, f_nko_count ! considering all offset-rate components allowed in this NK model f_index = f_trace(k) c1 = f_c1(k) IF (.NOT.f_creeps(k)) THEN ! this offset rate is seismogenic m1 = trace_loc(1, f_index) m2 = trace_loc(2, f_index) DO m = m1, (m2-1) ! consider all steps along this fault trace (whether within the NK .feg area, or not) !ignore any zero-length steps, with repeated locations IF ((trace(1, m) /= trace(1, m+1)).OR.(trace(2, m) /= trace(2, m+1)).OR.(trace(3, m) /= trace(3, m+1))) THEN !This step along the fault trace has positive length. !Earth radius: radius = R ! from NeoKinema parameter file !Observer location: CALL DUvec_2_ThetaPhi(r_, btheta, bphi) ! r_ is the uvec of the observer (node, or grid-point) location !Mid-point of trace of dislocation patch: tv(1:3) = trace(1:3, m) + trace(1:3, m+1) CALL DMake_Uvec(tv, mid_uvec) CALL DUvec_2_ThetaPhi(mid_uvec, ftheta, fphi) !Half-length (in m) of the trace of this dislocation patch: tv(1:3) = mid_uvec(1:3) - trace(1:3, m) lf = R * DMagnitude(tv) ! where R is Earth radius in m, and tv is dimensionless ! ARGUME is the argument of the trace of the dislocation patch, ! measured in radians counterclockwise from +theta (from South): CALL DLocal_Theta(r_, Theta) CALL DLocal_Phi(r_, Phi) dot_Theta = DDot(tv, Theta) ! S component of (dimensionless) half-length of trace dot_Phi = DDot(tv, Phi) ! E component of (dimensionless) half-length of trace argume = DATan2F(dot_Phi, dot_Theta) ! DIPF is fault dip in radians measured clockwise ! (initially, down) from horizontal on the right side of the ! fault (when looking in direction ARGUME). Usually >= Pi_over_2. !Note: The forward direction is the digitizing direction, which is L-->R on the footwall side; ! thus dips measured from the RHS are all 90-165 degrees (but, converted to radians for dipf). dip_degrees = f_dip_degrees(f_index) ! value from f*.dig, if any! 1 <= f_index <= f_highest IF (dip_degrees == 0.0D0) THEN ! use NeoKinema default value: IF ((c1 == 'R').OR.(c1 == 'r').OR.(c1 == 'L').OR.(c1 == 'l')) THEN dip_degrees = 90.0D0 ELSE IF ((c1 == 'N').OR.(c1 == 'n').OR.(c1 == 'D').OR.(c1 == 'd')) THEN dip_degrees = normal_dip_degrees ELSE IF ((c1 == 'T').OR.(c1 == 't').OR.(c1 == 'P').OR.(c1 == 'p')) THEN dip_degrees = thrust_dip_degrees ELSE IF ((c1 == 'S').OR.(c1 == 's')) THEN dip_degrees = subduction_dip_degrees ELSE WRITE (*, "(' ERROR: Illegal sense-byte ',A,' for fault offset-rate datum #',I6,' on F',I4)") c1, k, f_index WRITE (60, "('ERROR: Illegal sense-byte ',A,' for fault offset-rate datum #',I6,' on F',I4)") c1, k, f_index CALL Pause() STOP END IF END IF dipf = Pi - (dip_degrees * radians_per_degree)! Pi_over_2 < dipf < Pi !depths of top and bottom of the locked patch: ztop = 0.0D0 ! Was: "ztop = f_locking_depth_m_min(k)" in NeoKinema versions before 2.3 zbot = f_locking_depth_m_max(f_index) ! accessing via the trace index number, e.g., "1234" of "F1234". !Note that these were already converted from km-->m, and replaced with default values if negative, in the input section. IF (zbot > 0.0D0) THEN ! SLIP must be already rotated into fault coordinates: ! SLIP(1) is the component parallel to the fault, ! and it is positive for sinistral sense of slip. ! SLIP(2) is the horizontal component in the direction ! perpendicular to the trace of the fault, ! and it is positive for divergence across the trace. ! SLIP(3) is the relative vertical component, and it ! is positive when the right side of the fault ! (when looking along direction ARGUME) is down. sliprate_mps_x1x2x3 = 0.0D0 ! initializing in case enumeration of cases fails rate_mps = f_model_offset_rate(k) ! This is not locally accurate; it uses the trace-averaged offset rate. In m/s. IF ((c1 == 'L').OR.(c1 == 'l')) THEN sliprate_mps_x1x2x3(1) = rate_mps sliprate_mps_x1x2x3(2) = 0.0D0 sliprate_mps_x1x2x3(3) = 0.0D0 ELSE IF ((c1 == 'R').OR.(c1 == 'r')) THEN sliprate_mps_x1x2x3(1) = -rate_mps sliprate_mps_x1x2x3(2) = 0.0D0 sliprate_mps_x1x2x3(3) = 0.0D0 ELSE IF ((c1 == 'D').OR.(c1 == 'd')) THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = rate_mps sliprate_mps_x1x2x3(3) = rate_mps * DTAN(dipf) ELSE IF ((c1 == 'N').OR.(c1 == 'n')) THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = rate_mps / ABS(DTAN(dipf)) sliprate_mps_x1x2x3(3) = -rate_mps ELSE IF ((c1 == 'T').OR.(c1 == 't')) THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = -rate_mps / ABS(DTAN(dipf)) sliprate_mps_x1x2x3(3) = rate_mps ELSE IF ((c1 == 'P').OR.(c1 == 'p').OR.(c1 == 'S').OR.(c1 == 's')) THEN sliprate_mps_x1x2x3(1) = 0.0D0 sliprate_mps_x1x2x3(2) = -rate_mps sliprate_mps_x1x2x3(3) = -rate_mps * DTAN(dipf) END IF ! WEDGE is a tolerance (in radians) for dip; if DIPF is ! within WEDGE of Pi/2, then fault is considered vertical ! and routine Halo is called; otherwise, Aura is called. wedge = 0.2618D0 ! 15 degrees, expressed in radians (as in NeoKinema) CALL DChange (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & radius, & & sliprate_mps_x1x2x3, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output V_E = V_E + duphi ! East = Phi direction (in m/s) V_N = V_N - duthet ! North = -Theta direction (in m/s) !======= BEGIN new triangular-dislocation-patch code added in StrainRates2016: =========================== !GPBhere ! Add triangular dislocation patches in gores of non-vertical ! faults which have any significant bend(s) within their traces (i.e., bends between ! adjacent digitization steps in the input f_[token].DIG file). IF (ABS(Pi_over_2 - dipf) > 0.2618D0) THEN ! fault is significantly non-vertical (small-angle dip is less than 75 degrees). IF (m < (m2-1)) THEN ! there will be another step, so a bend going into the next step is possible IF ((trace(1, m+1) /= trace(1, m+2)).OR.(trace(2, m+1) /= trace(2, m+2)).OR.(trace(3, m+1) /= trace(3, m+2))) THEN !Next step has postitive length (and defined azimuth). tv1(1:3) = trace(1:3, m) ! start point of first step (before the possible bend) tv2(1:3) = trace(1:3, m+1) ! location of the possible bend tv3(1:3) = trace(1:3, m+2) ! which is legal since m < (m2-1); m+2 <= m2. prebend_trace_azim = DCompass(from_uvec = tv1, to_uvec = tv2) postbend_trace_azim = DCompass(from_uvec = tv2, to_uvec = tv3) IF (ABS(SIN(postbend_trace_azim - prebend_trace_azim)) > 0.05236D0) THEN ! significant bend, > 3 degrees (guarding against cycle-shifts) P1_uvec = tv2 ! shallow vertex of the TD is at the point of the bend P1_depth_m = 100.0D0 ! burying the shallowest vertex by 100 m below the free surface, to avoid singularities and abends IF (SIN(postbend_trace_azim - prebend_trace_azim) > 0.0D0) THEN ! bend is to the right (seen from above, moving in digitization direction) SIGN = 1.0D0 ! positive sign; must add a triangular dislocation to fill an empty gore. CALL DTurn_To (azimuth_radians = postbend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P2_uvec) ! output CALL DTurn_To (azimuth_radians = prebend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P3_uvec) ! output ELSE ! bend is to the left (seen from above, moving in digitization direction) SIGN = -1.0D0 ! negative sign; must subtract a triangular dislocation to cancel part of one redundant rectangular patch CALL DTurn_To (azimuth_radians = prebend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P2_uvec) ! output CALL DTurn_To (azimuth_radians = postbend_trace_azim - Pi_over_2, & & base_uvec = P1_uvec, far_radians = (zBot/TAN(Pi - dipf)/R), & ! input & omega_uvec = omega_uvec, result_uvec = P3_uvec) ! output END IF P2_depth_m = zBot P3_depth_m = zBot base_azimuth_radians = DCompass(from_uvec = P2_uvec, to_uvec = P3_uvec) ! determined in the counterclockwise direction; measured clockwise from N !NOTE that this deep "base" of the triangle is always horizontal, because there is always a uniform zBot locking depth along one NeoKinema fault. IF (c1 == 'L') THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians ! for upper plate, above the triangular dislocation ELSE IF (c1 == 'R') THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians + Pi ! for upper plate, above the triangular dislocation ELSE IF (c1 == 'D') THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians + Pi_over_2 ! for upper plate, above the triangular dislocation ELSE IF (c1 == 'N') THEN heave_rate_mps = rate_mps / ABS(TAN(dipf)) heave_azimuth_radians = base_azimuth_radians + Pi_over_2 ! for upper plate, above the triangular dislocation ELSE IF (c1 == 'T') THEN heave_rate_mps = rate_mps / ABS(TAN(dipf)) heave_azimuth_radians = base_azimuth_radians - Pi_over_2 ! for upper plate, above the triangular dislocation ELSE IF ((c1 == 'P').OR.(c1 == 'S')) THEN heave_rate_mps = rate_mps heave_azimuth_radians = base_azimuth_radians - Pi_over_2 ! for upper plate, above the triangular dislocation END IF !NOTE that the heave(-rate) vector is always horizontal, by definition. !Thus, the angle between "base" and "heave" is always in the horizontal plane. Ds = SIGN * (heave_rate_mps/COS(Pi - dipf)) * SIN(-(heave_azimuth_radians - base_azimuth_radians)) ! where thrusting is defined as positive. Ss = SIGN * heave_rate_mps * COS(heave_azimuth_radians - base_azimuth_radians) ! where left-lateral is defined as positive Ts = 0.0D0 ! no sill or dike intrusion component Poisson = 0.25D0 test_uvec(1:3) = r_(1:3) ! unit vector representing the test point (which may be a node, or a grid-point). test_depth_m = 0.0D0 ! always compute the result on the free surface CALL Gore(R, P1_uvec, P2_uvec, P3_uvec, P1_depth_m, P2_depth_m, P3_depth_m, & ! inputs & Ds, Ss, Ts, Poisson, & ! inputs & test_uvec, test_depth_m, & ! inputs & uSouth, uEast, uRadial) ! outputs V_E = V_E + uEast ! East = Phi direction (in m/s) V_N = V_N - uSouth ! North = -Theta direction (in m/s) END IF ! bend at this point is significant (more than 5 degrees). END IF ! next step has positive length (and thus a defined azimuth). END IF ! (m < (m2-1)); there will be another step, so a bend going into that step is possible. END IF! this fault is non-vertical (small-angle dip is less than 75 degrees). !========= END new triangular-dislocation-patch code added in StrainRates2016: =========================== END IF ! zbot > 0.0D0 END IF ! this step along the fault trace has positive length END DO ! m = m1, (m2-1) (OR, equivalently) m = trace_loc(1, f_index), (trace_loc(2, f_index) - 1) END IF ! .NOT.f_creeps(k) END DO ! k = 1, f_nko_count END SUBROUTINE MeanCoseismic_Velocity_at_Point SUBROUTINE No_Fault_Elements_Allowed() ! Called if nfl > 0 in the .feg file: ! Prints explanatory messages and stops execution. IMPLICIT NONE 101 FORMAT (& &/' This .feg (finite element grid) file contains fault elements!'& &/' Fault elements are not allowed in NeoKinema grids, because:'& &/' I. NeoKinema does not require fault elements.'& &/' 1. NeoKinema has logic to add the compliance of any number of'& &/' faults to the continuum (triangle) elements that contain them.'& &/' 2. NeoKinema has logic to infer the heave-rate and slip-rate of'& &/' such implied fault(s) from the computed strain-rate of the'& &/' triangular continuum element(s).'& &/' 3. Graphics program NeoKineMap has logic to plot the heave-rates'& &/' of these faults, and also velocity fields with fault '& &/' disontinuities, without the use of fault elements.') 102 FORMAT (& &/' II. Fault elements cause bad grid topology.'& &/' 1. Fault elements are not read or stored by NeoKineMap.'& &/' 2. With fault elements ignored, the grids on the two sides of'& &/' each fault are not connected in any way.'& &/' 3. The solution process may fail due to a singular stiffness'& &/' matrix during solution of the linear system.'& &/' 4. Even if the solution does not fail, its physical interpretation'& &/' will be problematical.') 103 FORMAT (& &/' III. Fault elements should be eliminated from the .feg file.'& &/' 1. Re-load this .feg file into OrbWin, select command Faults,'& &/' and use the right mouse button to Heal the fault cuts.'& &/' 2. Use command Adjust to move all nodes off of fault traces.'& &/' 3. Save the edited grid and re-number it with OrbNumber.'& &/' 4. Alternatively, use command Tile in OrbWin to create a'& &/' new .feg grid with no fault elements.') WRITE (*, 101) CALL Pause() WRITE (*, 102) CALL Pause() WRITE (*, 103) CALL Pause() WRITE (21, 101) WRITE (21, 102) WRITE (21, 103) STOP END SUBROUTINE No_Fault_Elements_Allowed SUBROUTINE Pause () IMPLICIT NONE CHARACTER*1 c1 WRITE (*,"(' Press [Enter] to continue...'\)") READ (*,"(A)") c1 END SUBROUTINE Pause SUBROUTINE Plane_area (folding) LOGICAL, INTENT(OUT) :: folding !puts areas of plane triangles (below the spherical surface, except at corners) into array a_; !if any area is zero or negative, reports folding INTEGER :: i1, i2, i3, l_ REAL*8, DIMENSION(3) :: a, b, c, t folding = .FALSE. DO l_ = 1, numel ! global, like most variables i1 = nodes(1, l_) i2 = nodes(2, l_) i3 = nodes(3, l_) a = node_uvec(1:3,i2) - node_uvec(1:3,i1) b = node_uvec(1:3,i3) - node_uvec(1:3,i2) CALL DCross (a, b, c) t = node_uvec(1:3,i1) + node_uvec(1:3,i2) + node_uvec(1:3,i3) IF (DDot(t, c) > 0.0D0) THEN a_(l_) = DMagnitude(c) * half_R2 ELSE folding = .TRUE. RETURN END IF END DO END SUBROUTINE Plane_area SUBROUTINE Prevent (bad_thing, line, filename) INTEGER, INTENT(IN) :: line CHARACTER(*), INTENT(IN) :: bad_thing, filename PRINT "(' Error: ',A,' is illegal in line ',I6/' of ',A)", & TRIM(bad_thing), line, TRIM(filename) WRITE (60,"('Error: ',A,' is illegal in line ',I6/' of ',A)") & TRIM(bad_thing), line, TRIM(filename) CALL Pause() STOP END SUBROUTINE Prevent SUBROUTINE Pull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. REAL*8, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL*8 :: factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.0D0) THEN highest = MAXVAL(s) medium = 1.00D0 - lowest - highest IF (medium > 0.0D0) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0.0D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00D0 / (1.00D0 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00D0 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00D0 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0.0D0 s(sideb) = 0.0D0 END IF END IF END SUBROUTINE Pull_in SUBROUTINE Write_Nodal_Velocities(type_string, V_E, V_N) IMPLICIT NONE ! but global StainRates2016 variables are used freely. CHARACTER*(*), INTENT(IN) :: type_string !Note that the above will be either: longTermAverage, meanCoseismic, or interseismic. REAL*8, DIMENSION(rows, columns), INTENT(IN) :: V_E, V_N CHARACTER*80 :: filename INTEGER :: i, j filename = "v_" // TRIM(type_string) // ".out" OPEN (UNIT = 1, FILE = TRIM(filename)) ! absolute OPEN; overwrites any existing file WRITE (1, "(A,', {for NeoKineMap plotting, together with StrainRates2016.feg}')") TRIM(filename) WRITE (1, "(A)") "Contents based on NeoKinema output file: " // TRIM(NK_vOut_title1) WRITE (1, "('This file created by PROGRAM StrainRates2016.')") DO i = 1, rows DO j = 1, columns IF (withinNKfeg(i, j)) WRITE (1, "(ES13.6,1X,ES13.6)") -V_N(i, j), V_E(i, j) ! V_theta, V_phi = V_S, V_E = -V_N, V_E END DO END DO CLOSE (1) WRITE (*, "(' ',A,' was written.')") TRIM(filename) WRITE (60, "(A,' was written.')") TRIM(filename) END SUBROUTINE Write_Nodal_Velocities SUBROUTINE Write_GSM_Table(type_string, V_E, V_N, eDot_EW, eDot_NS, eDot_NE) IMPLICIT NONE ! but global StainRates2016 variables are used freely. CHARACTER*(*), INTENT(IN) :: type_string !Note that the above will be either: longTermAverage, meanCoseismic, or interseismic. REAL*8, DIMENSION(rows, columns), INTENT(IN) :: V_E, V_N, eDot_EW, eDot_NS, eDot_NE CHARACTER*80 :: filename INTEGER :: i, j REAL*8 :: latitude, longitude filename = TRIM(type_string) // "_table.dat" OPEN (UNIT = 1, FILE = TRIM(filename)) ! absolute OPEN; overwrites any existing file WRITE (1, "('E_lon_dg N_lat_d V_E_m/s V_N_m/s eDot_EW eDot_NS eDot_NE (in /s)')") DO i = 1, rows latitude = lat_max - (i-1) * d_lat DO j = 1, columns longitude = lon_min + (j-1) * d_lon IF (withinNKfeg(i, j)) WRITE (1, "(F8.3,1X,F7.3,1X,ES13.6,1X,ES13.6,1X,ES9.2,1X,ES9.2,1X,ES9.2)") & & longitude, latitude, V_E(i, j), V_N(i, j), eDot_EW(i, j), eDot_NS(i, j), eDot_NE(i, j) END DO END DO CLOSE (1) WRITE (*, "(' ',A,' was written.')") TRIM(filename) WRITE (60, "(A,' was written.')") TRIM(filename) END SUBROUTINE Write_GSM_Table END PROGRAM StrainRates2016