PROGRAM Continuum_strainrate_reporter ! Computes mean long-term continuum strainrates in vicinity of template points, ! as predicted by one particular NeoKinema model. !(Continuum strainrates do not include strain due to slip on modeled faults, ! or elastic strain accumulation in interseismic times.) ! Input files are 3: ! (1) an e_[token].out file (continuum strainrates output from NeoKinema); ! (2) the .feg grid file used in that run of NeoKinema; ! (3) a template file listing geographic points where tensors are needed. ! One output file is written: longterm_continuum_strainrates.dat, ! in the same format as the template file (3), but with added columns. ! At each surface point, the 3 components quoted in 3 new columns are: ! e-dot_EW, e-dot_NS, and e-dot_NE. ! All are in SI units of per-second. ! Extensional strain components are considered positive, ! as is elongation of a square into a SW-&-NE-pointing parallelogram. ! For e-dot_NE, I use the "tensor" formula: ! e-dot_NE = (1/2) * ((dV_E/d_N) + (dV_N/d_E)), ! where V is the long-term crustal velocity vector at the surface, ! rather than the "engineering" formula which lacks the factor of (1/2). ! By Peter Bird, UCLA, 2011.10.20, improved 2012.01.31 & 2013.01.10, ! for the UCERF3 fault-based deformation model ! reported to Kaj Johnson, Wayne Thatcher, & Ned Field; ! and later for the NSHM-WUS fault-based deformation model. USE Sphere ! spherical geometry by P. Bird IMPLICIT NONE TYPE :: is123 ! element & internal coordinates of any point on surface INTEGER :: element REAL, DIMENSION(3) :: s END TYPE is123 CHARACTER*79 :: e_token_nko_pathfile, & & feg_pathfile, & & template_pathfile, & & title_line INTEGER :: a, b, i, ii, ios, j, jj, k, l_, latlon_ordering, & & m, n_title_lines, numel, numPoints, numnod, SH_azimuth_int INTEGER, DIMENSION(:,:), ALLOCATABLE :: nodes INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor ! list of neighboring finite elements !(1:3 = side crossed; 1:numel = element index l_) LOGICAL :: folding LOGICAL, DIMENSION(:), ALLOCATABLE :: faulted REAL, PARAMETER :: R = 6371000.0 ! radius of Earth in m REAL, PARAMETER :: Delta_degrees = 0.1 ! descriptive of the RELM_GriddedRegion.txt input file. REAL :: d_degrees, & & eDot_EE, eDot_EW, eDot_NE, eDot_NS, eDot_rr, eDot_SE, eDot_SS, & & eDot_1, eDot_2, eDot_3, & & e1h, e2h, log_sup_eDot_i, & & half_R2, s1, s2, s3, sup_eDot_i, t1, t2, & & u1xh, u1yh, u2xh, u2yh REAL, DIMENSION(3) :: tv, v1 REAL, DIMENSION(:), ALLOCATABLE :: a_ ! areas of elements, in m^2 REAL, DIMENSION(:,:,:), ALLOCATABLE :: benchLat, benchLon ! (11,11,numPoints) REAL, DIMENSION(:), ALLOCATABLE :: nodLat, nodLon REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: benchmark_uvec ! (11,11,3,numPoints) REAL, DIMENSION(:,:), ALLOCATABLE :: center REAL, DIMENSION(:,:), ALLOCATABLE :: e_nko REAL, DIMENSION(:,:), ALLOCATABLE :: xyz_nod TYPE(is123), DIMENSION(:,:,:),ALLOCATABLE :: benchmark_is ! (11,11,numPoints) ! locations of geodetic benchmarks in internal coordinates d_degrees = Delta_degrees/11.0 ! There will be 11 x 11 = 121 integration points per template point. WRITE (*, "(' PROGRAM Continuum_strainrate_reporter')") WRITE (*, "(' ')") WRITE (*, "(' Computes mean long-term continuum strainrates in vicinity of template points,')") WRITE (*, "(' as predicted by one particular NeoKinema model.')") WRITE (*, "(' (Continuum strainrates do not include strain due to slip on modeled faults,')") WRITE (*, "(' or effects of temporary elastic strain accumulation in interseismic time.)')") WRITE (*, "(' ')") WRITE (*, "(' Input files are 3:')") WRITE (*, "(' (1) an e_[token].out file (continuum strainrates output from NeoKinema);')") WRITE (*, "(' (2) the .feg grid file used in that run of NeoKinema;')") WRITE (*, "(' (3) a template file listing geographic points where tensors are needed.')") WRITE (*, "(' ')") WRITE (*, "(' One output file is written: longterm_continuum_strainrates.dat,')") WRITE (*, "(' in the same format as the template file (3), but with added columns.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' At each surface point, the 3 components quoted in 3 new columns are:')") WRITE (*, "(' e-dot_EW, e-dot_NS, and e-dot_NE.')") WRITE (*, "(' All are in SI units of per-second.')") WRITE (*, "(' Extensional strain components are considered positive,')") WRITE (*, "(' as is elongation of a square into a SW-&-NE-pointing parallelogram.')") WRITE (*, "(' For e-dot_NE, I use the ""tensor"" formula:')") WRITE (*, "(' e-dot_NE = (1/2) * ((dV_E/d_N) + (dV_N/d_E)),')") WRITE (*, "(' where V is the long-term crustal velocity vector at the surface,')") WRITE (*, "(' rather than the ""engineering"" formula which lacks the factor of (1/2).')") WRITE (*, "(' ')") WRITE (*, "(' By Peter Bird, UCLA, 2011.10.20, improved 2012.01.31 & 2013.01.10,')") WRITE (*, "(' for the UCERF3 fault-based deformation model')") WRITE (*, "(' reported to Kaj Johnson, Wayne Thatcher, & Ned Field,')") WRITE (*, "(' and later for the NSHM-WUS fault-based models.')") WRITE (*, "(' ')") CALL Pause() WRITE (*, "(' ')") WRITE (*, "(' ')") WRITE (*, "(' ')") WRITE (*, "(' ')") !================== INPUT: ================================================ WRITE (*, *) WRITE (*, "(' Enter [Drive:][\path\]filename for the e_[token].nko file: ')") READ (*, "(A)") e_token_nko_pathfile OPEN (UNIT = 1, FILE = TRIM(e_token_nko_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ', A, ' not found (in current folder).')") TRIM(e_token_nko_pathfile) CALL Pause() STOP END IF WRITE (*, *) WRITE (*, "(' Enter file name for matching F-E grid (.feg) file: ')") READ (*, "(A)") feg_pathfile OPEN (UNIT = 2, FILE = TRIM(feg_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ', A, ' not found (in current folder).')") TRIM(feg_pathfile) CALL Pause() STOP END IF WRITE (*, *) WRITE (*, "(' Title line found in this file:')") READ (2, "(A)") title_line WRITE (*, "(' ',A)") TRIM(title_line) READ (2, *) numnod ALLOCATE ( nodLon(numnod) ) ALLOCATE ( nodLat(numnod) ) DO i = 1, numnod READ (2, *) j, nodLon(j), nodLat(j) END DO READ (2, *) numel ALLOCATE ( nodes(3, numel) ) DO i = 1, numel READ (2, *) j, nodes(1, j), nodes(2, j), nodes(3, j) END DO WRITE (*, *) WRITE (*, "(' ', A, ' has been read and memorized.')") TRIM(feg_pathfile) CLOSE (2) ! feg_pathfile ALLOCATE ( e_nko(3, numel) ) ALLOCATE ( faulted(numel) ) DO i = 1, numel READ (1, *) j, faulted(i), e_nko(1, i), e_nko(2, i), e_nko(3, i) !Note that these components are continuum: eDot_SS, eDot_SE, eDot_EE !because NeoKinema uses an internal (theta, phi) = (South, East) system. END DO WRITE (*, *) WRITE (*, "(' ', A, ' has been read and memorized.')") TRIM(e_token_nko_pathfile) CLOSE (1) ! e_token_nko_pathfile WRITE (*, *) WRITE (*, "(' Enter file name for the point-list (template) file to guide reporting:')") READ (*, "(A)") template_pathfile OPEN (UNIT = 3, FILE = TRIM(template_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: ', A, ' not found (in current folder).')") TRIM(template_pathfile) CALL Pause() STOP END IF WRITE (*, *) WRITE (*, "(' Here are the first 5 lines found in this file:')") DO i = 1, 5 READ (3, "(A)") title_line WRITE (*, "(' ',A)") TRIM(title_line) END DO WRITE (*, *) n_title_lines = 0 CALL Prompt_for_Integer('How many title/header lines are there?',n_title_lines,n_title_lines) CLOSE(3) 10 WRITE (*, *) latlon_ordering = 1 WRITE (*, "(' Which format is in use?')") WRITE (*, "(' 1: {latitude, longitude} (used in UCERF3)')") WRITE (*, "(' 2: {longitude, latitude} (used in NSHM)')") CALL Prompt_for_Integer('Enter 1 or 2 to select format:',latlon_ordering,latlon_ordering) IF ((latlon_ordering < 1).OR.(latlon_ordering > 2)) THEN WRITE (*, "(' ERROR: Illegal choice. Try again.')") CALL Pause() GO TO 10 END IF !re-open same file and count points: OPEN (UNIT = 3, FILE = TRIM(template_pathfile), STATUS = "OLD", IOSTAT = ios) IF (n_title_lines > 0) THEN DO i = 1, n_title_lines READ (3, *) END DO END IF numPoints = 0 siting: DO READ (3, *, IOSTAT = ios) t1, t2 ! either {lat, lon} or {lon, lat} IF (ios /= 0) EXIT siting numPoints = numPoints + 1 END DO siting WRITE (*, *) WRITE (*, "(' ', I6, ' (lon, lat) points were found in blank template file.')") numPoints ALLOCATE ( benchLat(11,11,numPoints) ) ALLOCATE ( benchLon(11,11,numPoints) ) CLOSE (3) ! just so as to reset the file !re-open and third time and memorize contents: OPEN (UNIT = 3, FILE = TRIM(template_pathfile), STATUS = "OLD", IOSTAT = ios) IF (n_title_lines > 0) THEN DO i = 1, n_title_lines READ (3, *) END DO END IF DO i = 1, numPoints IF (latlon_ordering == 1) THEN READ (3, *) benchLat(6,6,i), benchLon(6,6,i) ELSE IF (latlon_ordering == 2) THEN READ (3, *) benchLon(6,6,i), benchLat(6,6,i) END IF DO ii = 1, 11 DO jj = 1, 11 benchLat(ii,jj,i) = benchLat(6,6,i) + d_degrees * (6-ii) benchLon(ii,jj,i) = benchLon(6,6,i) + d_degrees * (jj-6) END DO END DO END DO WRITE (*, "(' ', I6, ' (lon, lat) points have been read, multiplied, & memorized.')") numPoints CLOSE (3) ! template_pathfile !================== COMPUTATION: =========================================== WRITE (*, *) half_R2 = (R**2)/2. ! convert node locations to unit vectors WRITE (*, "(' Converting node locations to unit vectors...')") ALLOCATE ( xyz_nod(3, numnod) ) DO i = 1, numnod CALL LonLat_2_Uvec(nodLon(i), nodLat(i), tv) xyz_nod(1:3, i) = tv(1:3) END DO ! find areas of elements WRITE (*, "(' Finding areas of finite elements...')") ALLOCATE ( a_(numel) ) CALL Plane_Area(folding) !? ! find element center points WRITE (*, "(' Finding center points of finite elements...')") ALLOCATE ( center(3, numel) ) DO l_ = 1, numel v1 = xyz_nod(1:3,nodes(1,l_)) + xyz_nod(1:3,nodes(2,l_)) + xyz_nod(1:3,nodes(3,l_)) CALL Unitise(v1, tv) center(1:3, l_) = tv END DO ! find neighboring elements on all 3 sides of each WRITE (*, "(' Finding neighboring elements on all 3 sides...')") ALLOCATE ( neighbor(3, numel) ) neighbor = 0 ! whole array 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 ! convert integration points surrounding template lat/lon locations to unit vectors WRITE (*, "(/' Converting lat/lon coordinates to unit vectors at all integration points...')") ALLOCATE ( benchmark_uvec(11,11,3,numPoints) ) DO i = 1, numPoints DO ii = 1, 11 DO jj = 1, 11 CALL LonLat_2_Uvec(benchLon(ii,jj,i), benchLat(ii,jj,i), tv) benchmark_uvec(ii,jj,1:3,i) = tv(1:3) END DO END DO END DO !determine internal coordinates of each integration point: WRITE (*, "(/' Determining internal coordinates of each integration point...')") ALLOCATE ( benchmark_is(11,11,numPoints) ) DO i = 1, numPoints DO ii = 1, 11 DO jj = 1, 11 tv = benchmark_uvec(ii,jj,1:3,i) CALL Internal (tv, l_, s1, s2, s3) benchmark_is(ii,jj,i)%element = l_ benchmark_is(ii,jj,i)%s(1) = s1 benchmark_is(ii,jj,i)%s(2) = s2 benchmark_is(ii,jj,i)%s(3) = s3 END DO END DO IF (MOD(i,100)==0) WRITE (*,"(' Finished ',I6,' grid points out of ',I6)") i, numPoints END DO !================== OUTPUT: ================================================ WRITE (*, "(/' Computing mean strain-rate tensors at template points...')") OPEN (UNIT = 4, FILE = "longterm_continuum_strainrates.dat") IF (latlon_ordering == 1) THEN !Standard UCERF3-format horizontal-plane tensor header line: WRITE (4, "(' N_lat_deg E_lon_deg e-dot_EW e-dot_NS e-dot_NE')") ELSE IF (latlon_ordering == 2) THEN !Standard NSHM-format horizontal-plane tensor header line: WRITE (4, "(' E_lon_deg N_lat_deg e-dot_EW e-dot_NS e-dot_NE')") END IF !Header lines for creating a stress-direction-pseudodata file (for test plotting): !WRITE (4, "('(A7, 1X,A17,1X, F9.2, F8.2, I10,1X,A1)')") !WRITE (4, "('Site Type Regime DepthLongitudeLatitudeSH_azimuth Quality')") DO i = 1, numPoints eDot_SS = 0.0 ! initializing before sum over 121 integration points eDot_SE = 0.0 eDot_EE = 0.0 DO ii = 1, 11 DO jj = 1, 11 l_ = benchmark_is(ii,jj,i)%element IF (l_ > 0) THEN eDot_SS = eDot_SS + e_nko(1, l_) / 121. eDot_SE = eDot_SE + e_nko(2, l_) / 121. eDot_EE = eDot_EE + e_nko(3, l_) / 121. ELSE ! template lat/lon point lies outside the .FEG area: !(no additions to the 3 sums) END IF END DO END DO eDot_EW = eDot_EE eDot_NS = eDot_SS eDot_NE = -eDot_SE !Additional analysis is conducted for use in confirmation plots: CALL Principal_Axes_22 (e11=eDot_EW, e12=eDot_NE, e22=eDot_NS, & & e1=e1h, e2=e2h, & & u1x=u1xh,u1y=u1yh, u2x=u2xh,u2y=u2yh) ! Find principal values (e1,e2) of the symmetric 2x2 tensor ! e11 e12 ! e12 e22 ! and also the associated unit eigenvectors: ! #1 = (u1x, u1y); #2 = (u2x, u2y). ! The convention is that e1 <= e2. eDot_rr = -(eDot_EW + eDot_NS) eDot_1 = MIN(eDot_rr, e1h, e2h) eDot_3 = MAX(eDot_rr, e1h, e2h) sup_eDot_i = MAX(ABS(eDot_1), ABS(eDot_3)) IF (sup_eDot_i > 0) THEN log_sup_eDot_i = ALOG10(sup_eDot_i) ELSE log_sup_eDot_i = -17.0 END IF SH_azimuth_int = 57.296 * ATAN2F(u1xh, u1yh) !Use line below to create a test-plot of strain-rate intensity: !WRITE (4, "(2F10.3, F10.3)") benchLon(6,6,i), benchLat(6,6,i), log_sup_eDot_i !Use line below to plot e1h azimuths as pseudo-stress-direction-data, for test-plotting by NeoKineMap: !WRITE (4, "('pseudo-data! ',F9.2,F8.2,I10,1X,'A')") benchLon(6,6,i), benchLat(6,6,i), SH_azimuth_int !Standard UCERF3-format 2x2 horizontal-plane tensor output, after program has been tested & confirmed: IF (latlon_ordering == 1) THEN WRITE (4, "(2F10.3, 3ES10.2)") benchLat(6,6,i), benchLon(6,6,i), eDot_EW, eDot_NS, eDot_NE ELSE IF (latlon_ordering == 2) THEN WRITE (4, "(2F10.3, 3ES10.2)") benchLon(6,6,i), benchLat(6,6,i), eDot_EW, eDot_NS, eDot_NE END IF IF (MOD(i,100)==0) WRITE (*,"(' Finished ',I6,' grid points out of ',I6)") i, numPoints END DO CLOSE (4) ! output file !================== (end): ================================================ WRITE (*, *) WRITE (*, "(' Job completed. See new file longterm_continuum_strainrates.dat (and rename it?).')") CALL Pause() CONTAINS REAL FUNCTION Atan2F (y, x) ! corrects for possible abend in case of (0., 0.) IMPLICIT NONE REAL, INTENT(IN) :: x, y IF ((x /= 0.).OR.(y /= 0.)) THEN Atan2F = ATAN2(y,x) ELSE Atan2F = 0. END IF END FUNCTION Atan2F 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, DIMENSION(3), INTENT(IN) :: vector REAL, INTENT(OUT) :: s1, s2, s3 INTEGER :: i1, i2, i3 REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL :: 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 = Dot_3D(vector, tv1) IF (dc <= 0.) CALL Prevent('"Internal" vector >= 90 deg. from element', 1, 'Dumb_s123') tv2 = xyz_nod(1:3, i1) d1 = Dot_3D(tv2, tv1) t = d1 / dc v1 = t * vector tvi = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) tvo = v1(1:3) - xyz_nod(1:3,i3) CALL Cross(tvi, tvo, tv) s1 = Dot_3D(tv1, tv) * half_R2 / a_(element) tvi = xyz_nod(1:3,i1) - xyz_nod(1:3,i3) tvo = v1(1:3) - xyz_nod(1:3,i1) CALL Cross(tvi, tvo, tv) s2 = Dot_3D(tv1, tv) * half_R2 / a_(element) s3 = 1.00 - s1 - s2 END SUBROUTINE Dumb_s123 SUBROUTINE Internal (b_, iele, s1, s2, s3) REAL, DIMENSION(3), INTENT(IN) :: b_ INTEGER, INTENT(OUT) :: iele REAL, 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 :: back1, back2, i, iet, l_ REAL :: r2, r2min, s1t, s2t, s3t REAL, DIMENSION(3) :: s_temp, tv ! establish defaults (not found) in case of quick exit iele = 0 ! default s1 = 0.0; s2 = 0.0; s3 = 0.0 ! default !find closest element center to initialize search r2min = 4.01 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 (Dot_3D(b_, tv) < 0.540) RETURN ! initialize search memory (with impossible numbers) back1 = -1 back2 = -2 is_it_here: DO ! first, check for infinite loop between 2 elements! IF (iet == back2) THEN ! in loop; force location in one or the other! 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 ! normal operation CALL Dumb_s123 (iet, b_, s1t, s2t, s3t) IF ((s1t < s2t) .AND. (s1t < s3t)) THEN ! s1 is most negative; most critical IF (s1t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t < s1t) .AND. (s2t < s3t)) THEN ! s2 is most negative; most critical IF (s2t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3 is most negative; most critical IF (s3t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN back2 = back1 back1 = iet 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 Pause() !The purpose of this routine is to give the user time to read an error message !before the program stops (and the operating system closes its window). !The programmer should WRITE the error, CALL Pause(), and then STOP. IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause SUBROUTINE Plane_area (folding) LOGICAL, INTENT(OUT) :: folding !puts areas of plane triangles (below surface) into array a_; !if any is zero or negative, reports folding INTEGER :: i1, i2, i3, l_ REAL, 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 = xyz_nod(1:3,i2) - xyz_nod(1:3,i1) b = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) CALL Cross (a, b, c) t = xyz_nod(1:3,i1) + xyz_nod(1:3,i2) + xyz_nod(1:3,i3) IF (Dot_3D(t, c) > 0.) THEN a_(l_) = Magnitude(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 (21,"('Error: ',A,' is illegal in line ',I6/' of ',A)") & TRIM(bad_thing), line, TRIM(filename) CALL Pause() !CALL Bad_Parameters() STOP END SUBROUTINE Prevent SUBROUTINE Principal_Axes_22 (e11, e12, e22, & & e1, e2, u1x,u1y, u2x,u2y) ! Find principal values (e1,e2) of the symmetric 2x2 tensor ! e11 e12 ! e12 e22 ! and also the associated unit eigenvectors: ! #1 = (u1x, u1y); #2 = (u2x, u2y). ! The convention is that e1 <= e2. IMPLICIT NONE REAL, INTENT(IN) :: e11, e12, e22 REAL, INTENT(OUT) :: e1, e2, u1x, u1y, u2x, u2y REAL :: c, f1, f11, f12, f2, f22, r, scale, smallest, test, theta ! Smallest number that can be squared without underflowing: smallest = 1.1 * SQRT(TINY(f1)) ! First, check for trivial isotropic case: IF ((e11 == e22).AND.(e12 == 0.)) THEN ! In this case, directions are arbitrary: e1 = e11 u1x = 1. u1y = 0. e2 = e22 u2x = 0. u2y = 1. RETURN END IF ! Else, re-scale matrix to prevent under/overflows: scale = MAX(ABS(e11), ABS(e12), ABS(e22)) f11 = e11 / scale IF (ABS(f11) < smallest) f11 = 0. f12 = e12 / scale IF (ABS(f12) < smallest) f12 = 0. f22 = e22 / scale IF (ABS(f22) < smallest) f22 = 0. ! Find eigenvalues and eigenvectors of scaled matrix: r = SQRT(((f11 - f22)*0.5)**2 + f12**2) c = (f11 + f22)*0.5 f1 = c - r f2 = c + r test = 0.01 * MAX (ABS(f1), ABS(f2)) IF ((ABS(f12) > test).OR.(ABS(f11 - f1) > test)) THEN theta = Atan2F((f11 - f1), -f12) ELSE theta = Atan2F(f12, (f1 - f22)) END IF u1x = COS(theta) u1y = SIN(theta) u2x = u1y u2y = -u1x ! Undo the scaling e1 = scale * f1 e2 = scale * f2 END SUBROUTINE Principal_Axes_22 SUBROUTINE Prompt_for_Integer (prompt_text, default, answer) ! Writes a line to the default (*) unit, with: ! "prompt_text" ["default"]: ! and accepts an answer with an integer 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 INTEGER, INTENT(IN) :: default INTEGER, INTENT(OUT) :: answer CHARACTER*11 :: instring, suggested INTEGER :: blank_at, bytes, ios, trial, written LOGICAL :: finished WRITE (suggested,"(I11)") default 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 integer of 9 digits 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 Prompt_for_Integer SUBROUTINE Pull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. REAL, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.) THEN highest = MAXVAL(s) medium = 1.00 - lowest - highest IF (medium > 0.) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0. sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00 / (1.00 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0. s(sideb) = 0. END IF END IF END SUBROUTINE Pull_in SUBROUTINE Unitise (b_, unit_vec) REAL, DIMENSION(3), INTENT(IN) :: b_ REAL, DIMENSION(3), INTENT(OUT) :: unit_vec REAL :: length length = Magnitude (b_) IF (length /= 0.) THEN unit_vec(1) = b_(1)/length unit_vec(2) = b_(2)/length unit_vec(3) = b_(3)/length ELSE unit_vec(1) = 1. unit_vec(2) = 0. unit_vec(3) = 0. END IF END SUBROUTINE Unitise END PROGRAM Continuum_strainrate_reporter