PROGRAM Assign_LRs !A utility program, belonging to the Shells_v5.0 / ShellSet group !of dynamic F-E codes for modeling neotectonics. !Uses nodal data (e.g., longitude, latitude, elevation, ...) !in an existing F-E grid (.FEG) file to assign Lithospheric Rheology !(LR) integer index numbers to all triangular continuum elements and !all fault elements. The particular assignment strategy is hard-coded, !but it will be relatively easy to change it and recompile this utility. !By Peter Bird, Department of Earth, Planetary, and Space Sciences, !University of California, Los Angeles, CA 90095-1567, !20 December 2022 (for Earth5-2023a.feg) & 31 January 2023 (for Earth5-2023b.feg) ! & 22 February 2023 (for Earth5-2023c.feg). USE DSphere ! In file DSPhere.f90, provided by Peter Bird, UCLA IMPLICIT NONE CHARACTER*132 :: longer_line, shorter_line, title1 CHARACTER*80, DIMENSION(:), ALLOCATABLE :: orogen_title INTEGER :: i, iOrogen, ios, iUnitD, iUnitG, iUnitO, iUnitT, & & j, & & LRi, LRn, & & mxDOF, mxEl, mxFEl, mxNode, maxPerOrogen, & & n, n1000, nFakeN, nRealN, numEl, nFl, numNod, numOrogens INTEGER, DIMENSION(:), ALLOCATABLE :: continuum_LRi ! will get dimensioned (1:mxEl). INTEGER, DIMENSION(:), ALLOCATABLE :: fault_LRi ! will get dimensioned (1:mxFEl). INTEGER, DIMENSION(:), ALLOCATABLE :: orogen_length ! will get dimensioned (1:numOrogens); = maximum # of digitized points in any outline. INTEGER, DIMENSION(:, :), ALLOCATABLE :: nodes ! will get dimensioned (1:3, 1:mxEl). INTEGER, DIMENSION(:, :), ALLOCATABLE :: nodeF ! will get dimensioned (1:4, 1:mxFEl). LOGICAL :: brief, inside, seafloor LOGICAL, DIMENSION(:), ALLOCATABLE :: checkN ! Will get dimensioned to (1:mxNode). LOGICAL, DIMENSION(:), ALLOCATABLE :: checkE ! Will get dimensioned to (1:mxEl). LOGICAL, DIMENSION(:), ALLOCATABLE :: checkF ! Will get dimensioned to (1:mxFEl). REAL*8 :: elevation, ELon, NLat, offMax REAL*8, DIMENSION(3) :: tvec, uvec, uvec1, uvec2, uvec3 REAL*8, DIMENSION(:), ALLOCATABLE :: cooling_curvature, & ! These will all get dimensioned as (1:mxNode). & density_anomaly, & & dQdTdA, & & elev, & & tLNode, & & xNode, yNode, zMNode REAL*8, DIMENSION(:), ALLOCATABLE :: offset ! will get dimensioned as (1:mxFEl). REAL*8, DIMENSION(:,:), ALLOCATABLE :: fDip ! will get dimensioned as (1:2, 1:mxFEl). REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: orogenUvecs ! will get dimensioned as (1:3, 1:numOrogens, 1:maxPerOrogen) ! iUnitG = device number associated with the grid input file: DATA iUnitG / 1 / ! iUnitO = device number associated with the grid output file: DATA iUnitO / 2 / ! iUnitD = device number associated with input of selected_PB2002_orogens.dig: DATA iUnitD / 3 / ! iUnitT = device number associated with the logFile (ASCII text output, including large tables): DATA iUnitT / 21 / ! numOrogens = # of orogens in the "selected_PB2002_orogens.dig" input file: DATA numOrogens / 1 / ! maxPerOrogen = maximum # of digitization points in the outline of any one orogen (see above): DATA maxPerOrogen / 1000 / !------------------ Beginning of executable code: ------------------------------------------- WRITE (*, *) WRITE (*, "(' PROGRAM Assign_LRs:')") WRITE (*, *) WRITE (*, "(' A utility program, belonging to the Shells_v5.0 / ShellSet group')") WRITE (*, "(' of dynamic F-E codes for modeling neotectonics.')") WRITE (*, "(' Uses nodal data (e.g., longitude, latitude, elevation, ...)')") WRITE (*, "(' in an existing F-E grid (.FEG) file to assign Lithospheric Rheology')") WRITE (*, "(' (LR) integer index numbers to all triangular continuum elements and')") WRITE (*, "(' all fault elements. The particular assignment strategy is hard-coded,')") WRITE (*, "(' but it will be relatively easy to change it and recompile this utility.')") WRITE (*, "(' By Peter Bird, Department of Earth, Planetary, and Space Sciences,')") WRITE (*, "(' University of California, Los Angeles, CA 90095-1567,')") WRITE (*, "(' 20 December 2022 (Earth5-2023a.feg) & 31 January 2023 (Earth5-2023b.feg),')") WRITE (*, "(' & 22 February 2023 (Earth5-2023c.feg).')") WRITE (*, *) CALL Pause() !======================== BEGIN lengthy extract of code from Shells_v5.0: =================================== ! Preview .feg file to determine array sizes: WRITE (*, 101) iUnitG 101 FORMAT (/' Attempting to read finite element grid from unit', I3/) READ (iUnitG, * , IOSTAT = ios) IF (ios /= 0) THEN WRITE(*, "(' ERR','OR: File not found, or file is empty,' & & /' or file is too short.')") CALL Pause() STOP END IF READ (iUnitG, * , IOSTAT = ios) numNod IF (ios /= 0) THEN WRITE(*, "(' ERR','OR: File not found, or file is empty,' & & /' or file is too short.')") CALL Pause() STOP END IF mxNode = numNod mxDOF = 2 * mxNode DO 102 i = 1, numNod READ (iUnitG, * , IOSTAT = ios) IF (ios /= 0) THEN WRITE(*, "(' ERR','OR:File not found, or file is empty,' & & /' or file is too short.')") CALL Pause() STOP END IF 102 CONTINUE READ (iUnitG, * , IOSTAT = ios) numEl IF (ios /= 0) THEN WRITE(*, "(' ERR','OR: File not found, or file is empty,' & & /' or file is too short.')") CALL Pause() STOP END IF mxEl = numEl !Initialize survey to find LRn = MAX(continuum_LRi(1:mxEl), fault_LRi(1:MXFel) LRn = 0 ! until incremented below... DO 103 i = 1, numEl READ (iUnitG, "(A)", IOSTAT = ios) longer_line IF (ios /= 0) THEN WRITE(*, "(' ERR','OR:File not found, or file is empty,' & & /' or file is too short.')") CALL Pause() STOP END IF CALL Extract_LRi(longer_line, LRi, shorter_line) LRn = MAX(LRn, LRi) 103 CONTINUE nFl = 0 READ (iUnitG, * , IOSTAT = ios) n IF (ios == 0) nFl = n nFl = MAX(nFl, 0) mxFEl = nFl DO 105 i = 1, nFl READ (iUnitG, "(A)", IOSTAT = ios) longer_line IF (ios /= 0) THEN WRITE(*, "(' ERR','OR:File not found, or file is empty,' & & /' or file is too short.')") CALL Pause() STOP END IF CALL Extract_LRi(longer_line, LRi, shorter_line) LRn = MAX(LRn, LRi) 105 CONTINUE REWIND (UNIT = iUnitG) ! to prepare for CALL GetNet, below... ! ALLOCATE adjustable arrays (except those whose sizes ! are based on results from CALL KSize): ! DIMENSIONs using size variable mxNode: ALLOCATE ( checkN(mxNode), & & cooling_curvature(mxNode), & & density_anomaly(mxNode), & & dQdTdA(mxNode), & & elev(mxNode), & & tLNode(mxNode), & & xNode(mxNode), yNode(mxNode), zMNode(mxNode) ) ! DIMENSIONs using size variable mxEl: ALLOCATE ( checkE(mxEl), continuum_LRi(mxEl), nodes(3, mxEl) ) ! DIMENSIONs using size variable mxFEl: ALLOCATE ( checkF(mxFEl), & & fault_LRi(mxFEl), fDip(2, mxFEl), & & nodeF(4, mxFEl), & & offset(mxFEl) ) ! Input finite element grid and nodal data (up to 6 fields): WRITE(*, "(/' Attempting to create (new) log file on unit ', I3, '. Choose a filename: ' /)") iUnitT brief = .TRUE. ! minimize echoing; could be set .FALSE. for debugging? CALL GetNet (iUnitG, iUnitT, & ! input & mxDOF, mxEl, mxFEl, mxNode, & & brief, continuum_LRi, cooling_curvature, & ! output & density_anomaly, & & dQdTdA, elev, fault_LRi, fDip, & & nFakeN, nFl, nodeF, nodes, nRealN, & & numEl, numNod, n1000, offMax, offset, & & title1, tLNode, xNode, yNode, zMNode, & & checkE, checkF, checkN) ! work WRITE (*, "(' Finite element grid file has been read.')") !=========================== END lengthy extract of code from Shells_v5.0 ==================================== !============================ BEGIN logic of PROGRAM Assign_LRs: ============================================= continuum_LRi = 0 ! whole array (so nothing is undefined) fault_LRi = 0 ! whole array (so nothing is undefined) !Allocate space for orogen outlines to be read: ALLOCATE ( orogen_title(1:numOrogens) ) ALLOCATE ( orogen_length(1:numOrogens) ) ALLOCATE ( orogenUvecs(1:3, 1:numOrogens, 1:maxPerOrogen) ) orogen_title = ' ' ! whole list, to simplify debugging orogen_length = 0 ! whole list, to simplify debugging orogenUvecs = 0.0D0 ! whole array, to simplify any debugging !Get the orogen outlines: WRITE (*, "(/' Ready to read outlines of ', I2, ' selected orogens from .DIG file on unit ', I3 / ' Enter input .DIG file name: '/)") numOrogens, iUnitD DO i = 1, numOrogens READ (iUnitD, "(A)") longer_line orogen_title(i) = TRIM(longer_line) j = 0 recording: DO READ (iUnitD, "(A)") longer_line IF (longer_line(1:3) == "***") EXIT recording READ (longer_line, *) ELon, NLat CALL DLonLat_2_Uvec(Elon, NLat, uvec) j = j + 1 orogen_length(i) = j ! (although this may be increased on next pass) orogenUvecs(1:3, i, j) = uvec(1:3) END DO recording END DO ! i = 1, numOrogens (polygon outlines, in .DIG format) WRITE (*, *) WRITE (*, "(' Successfully read outlines of ', I2, ' orogens.')") numOrogens CALL Pause() WRITE (*, *) WRITE (*, "(' Working......')") !For each continuum element center, check whether it is: (1) seafloor?; and (2) in any orogen? DO i = 1, numEl !center-point of element: CALL DThetaPhi_2_Uvec(xNode(nodes(1, i)), yNode(nodes(1, i)), uvec1) CALL DThetaPhi_2_Uvec(xNode(nodes(2, i)), yNode(nodes(2, i)), uvec2) CALL DThetaPhi_2_Uvec(xNode(nodes(3, i)), yNode(nodes(3, i)), uvec3) tvec(1:3) = uvec1(1:3) + uvec2(1:3) + uvec3(1:3) CALL DMake_Uvec(tvec, uvec) ! uvec marks the center of the triangle !elevation of center-point of element: elevation = (elev(nodes(1, i)) + elev(nodes(2, i)) + elev(nodes(3, i))) / 3.0D0 seafloor = (elevation <= -2000.D0) !Now determine which orogen (if any) encloses this point? iOrogen = 0 ! (meaning "not in any orogen"; but this is just the default in case of failure) testing: DO j = 1, numOrogens !================================================================= inside = Within(uvec, numOrogens, maxPerOrogen, orogenUvecs, j, orogen_length(j)) !================================================================= IF (inside) THEN iOrogen = j EXIT testing END IF END DO testing !================= begin logic for Earth5-2023a, b, & d.feg ====================== IF (iOrogen > 0) THEN ! inside one of the selected orogens: continuum_LRi(i) = 0 ! so its rheology is adjustable by ShellSet ELSE continuum_LRi(i) = 1 ! so its rheology is fixed when running Shellset END IF !================= end logic for Earth5-2023a, b, & d.feg ====================== !!================= begin logic for Earth5-2023c.feg ====================== !IF (seafloor) THEN ! continuum_LRi(i) = 0 ! so its rheology is adjustable by ShellSet !ELSE ! continuum_LRi(i) = 1 ! so its rheology is fixed when running Shellset !END IF !!================= end logic for Earth5-2023c.feg ====================== END DO ! i = 1, numEl !For each fault-element center, check whether it is: (1) seafloor?; and (2) in any orogen? DO i = 1, nFl !center-point of element: CALL DThetaPhi_2_Uvec(xNode(nodeF(1, i)), yNode(nodeF(1, i)), uvec1) CALL DThetaPhi_2_Uvec(xNode(nodeF(2, i)), yNode(nodeF(2, i)), uvec2) tvec(1:3) = uvec1(1:3) + uvec2(1:3) CALL DMake_Uvec(tvec, uvec) ! uvec marks the center of the fault-element trace !elevation of center-point of fault-element: elevation = (elev(nodeF(1, i)) + elev(nodeF(2, i))) / 2.0D0 seafloor = (elevation <= -2000.D0) !Now determine which orogen (if any) encloses this point? iOrogen = 0 ! (meaning "not in any orogen"; but this is just the default in case of failure) retesting: DO j = 1, numOrogens !================================================================= inside = Within(uvec, numOrogens, maxPerOrogen, orogenUvecs, j, orogen_length(j)) !================================================================= IF (inside) THEN iOrogen = j EXIT retesting END IF END DO retesting !================= begin logic for Earth5-2023a, b, & d .feg ====================== IF (iOrogen > 0) THEN ! inside one of the selected orogens: fault_LRi(i) = 0 ! so its rheology is adjustable by ShellSet ELSE fault_LRi(i) = 1 ! so its rheology is fixed when running Shellset END IF !================= end logic for Earth5-2023a, b, & d.feg ====================== !!================= begin logic for Earth5-2023c.feg ====================== !IF (seafloor) THEN ! fault_LRi(i) = 0 ! so its rheology is adjustable by ShellSet !ELSE ! fault_LRi(i) = 1 ! so its rheology is fixed when running Shellset !END IF !!================= end logic for Earth5-2023c.feg ====================== END DO ! i = 1, nFl !============================ END logic of PROGRAM Assign_LRs: =============================================== !============================ OUTPUT the F-E grid (.FEG file), now with LR#s included: ======================= CALL PutNet (iUnitO, & ! INTENT(IN) & brief, & ! INTENT(IN) & continuum_LRi, & ! INTENT(IN) & dQdTdA, elev, & ! INTENT(IN) & fault_LRi, fDip, & ! INTENT(IN) & mxEl, mxFEl, mxNode, n1000, & ! INTENT(IN) & nFakeN, nFl, nodeF, nodes, & ! INTENT(IN) & nRealN, numEl, numNod, offset, & ! INTENT(IN) & title1, tLNode, xNode, yNode, zMNode, & ! INTENT(IN) & density_anomaly, & ! INTENT(IN) & cooling_curvature) ! INTENT(IN) !======================= CONCLUSION of this utility program: ================================================= WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS SUBROUTINE Extract_LRi (longer_line, & ! input & LRi, shorter_line) ! output ! New routine added for Shells_v5.0+ to support multiple !"Lithospheric Rheology" (abbreviated as "LR") integer codes, ! in any line of the input .feg file which define an element !(either a triangular continuum element, or a ! linear fault element). ! CHARACTER*80, INTENT(IN) :: longer_line is the whole ! element-definition line from the .feg file. ! INTEGER, INTENT(OUT) :: LRi is the rheologic code ! (or 0, if no such code was found). ! CHARACTER*80, INTENT(OUT) :: shorter_line has the " LRi" portion removed (if any), ! so it can be interpreted by the same code as in Shells_v4.1-. IMPLICIT NONE CHARACTER*80, INTENT(IN) :: longer_line INTEGER, INTENT(OUT) :: LRi CHARACTER*80, INTENT(OUT) :: shorter_line CHARACTER*80 :: string INTEGER :: longer_length, LR_start_byte longer_length = LEN_TRIM(longer_line) LR_start_byte = INDEX(longer_line, "LR") IF (LR_start_byte > 0) THEN ! the "LR" flag was found IF (longer_length > (LR_start_byte + 1)) THEN ! some byte(s) follow the "LR" string = longer_line((LR_start_byte + 2):longer_length) READ (string, *) LRi shorter_line = longer_line(1:(LR_start_byte - 1)) ELSE ! "LR" is present, but nothing follows it; infer 0. LRi = 0 shorter_line = longer_line(1:(LR_start_byte - 1)) END IF ELSE ! no "LR" flag is present LRi = 0 shorter_line = longer_line END IF END SUBROUTINE Extract_LRi SUBROUTINE GetNet (iUnit7, iUnit8, & ! input & mxDOF, mxEl, mxFEl, mxNode, & & brief, continuum_LRi, cooling_curvature, & ! output & density_anomaly, & & dQdTdA, elev, fault_LRi, fDip, & & nFakeN, nFl, nodeF, nodes, nRealN, & & numEl, numNod, n1000, offMax, offset, & & title1, tLNode, xNode, yNode, zMNode, & & checkE, checkF, checkN) ! work ! Read finite element grid from unit iUnit7 (assumed already OPENed). ! Echoes the important values to unit iUnit8 (assumed already OPENed). IMPLICIT NONE ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER, INTENT(IN) :: iUnit7, iUnit8, mxDOF, mxEl, mxFEl, mxNode ! input LOGICAL, INTENT(OUT) :: brief ! output INTEGER, INTENT(OUT) :: continuum_LRi ! output REAL*8, INTENT(OUT) :: cooling_curvature, density_anomaly, dQdTdA, elev ! output INTEGER, INTENT(OUT) :: fault_LRi ! output REAL*8, INTENT(OUT) :: fDip ! output INTEGER, INTENT(OUT) :: nFakeN, nFl, nodeF, nodes, nRealN, numEl, numNod, n1000 ! output REAL*8, INTENT(OUT) :: offMax, offset ! output CHARACTER*132, INTENT(OUT) :: title1 ! output REAL*8, INTENT(OUT) :: tLNode, xNode, yNode, zMNode ! output LOGICAL checkE, checkF, checkN ! work ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CHARACTER*80 :: longer_line, shorter_line LOGICAL allOK INTEGER i, index, j, k, l, LRi, n, nrt2 REAL*8 cooling_curvature_cpm2, density_anomaly_kgpm3, dips, elevi, & & off, pLat, pLon, qi, tli, vector, xi, yi, zmi DIMENSION checkE(mxEl), checkF(mxFEl), checkN(mxNode), & & continuum_LRi(mxEl), & & cooling_curvature(mxNode), & & density_anomaly(mxNode), & & dQdTdA(mxNode), elev(mxNode), & & fault_LRi(mxFEl), & & fDip(2, mxFEl), nodeF(4, mxFEl), & & nodes(3, mxEl), offset(mxFEl), tLNode(mxNode), & & xNode(mxNode), yNode(mxNode), zMNode(mxNode) DIMENSION dips(3), vector(9) title1 = ' ' READ (iUnit7, 2) title1 2 FORMAT (A) WRITE (iUnit8, 3) TRIM(title1) 3 FORMAT(/' Title of finite-element grid ='/' ', A) ! Read number of nodes, plus out-dated parameters that once ! permitted boundary nodes to be specially numbered as ! "fake" nodes with numbers from n1000+1 ... n1000+nFakeN. ! This option is no longer supported by my programs! ! (Option "brief" suppresses most output.) READ (iUnit7, * ) numNod, nRealN, nFakeN, n1000, brief IF (numNod /= (nRealN + nFakeN)) THEN WRITE (iUnit8, 4) numNod, nRealN, nFakeN 4 FORMAT (/' ERR0R: numNod (',I6,') IS NOT EQUAL TO SUM' & & /' OF nRealN (',I6,') AND nFakeN (',I6,').') CALL Pause() STOP END IF IF (nRealN > n1000) THEN WRITE (iUnit8, 5) nRealN, n1000 5 FORMAT (/' ERR0R: nRealN (',I6,') IS GREATER THAN' & & /' n1000 (',I6,').') CALL Pause() STOP END IF IF (numNod > mxNode) THEN WRITE (iUnit8, 10) numNod 10 FORMAT(/' INCREASE ARRAY-SIZE maxNod TO BE AT LEAST' & & /' THE NUMBER OF NODES (', I6, ') AND RECOMPILE.') CALL Pause() STOP END IF nrt2 = nRealN * 2 IF (nrt2 > mxDOF) THEN WRITE (iUnit8, 12) nrt2 12 FORMAT (/' INCREASE ARRAY-SIZE maxDOF TO ', I6, & & ' AND RECOMPILE.') CALL Pause() STOP END IF IF (brief) THEN WRITE (iUnit8, 35) 35 FORMAT(/' (Since option ""brief"" = .TRUE., grid will not be echoed here.)') ELSE WRITE (iUnit8, 40) numNod 40 FORMAT (/' There are',I5,' nodes in the grid') WRITE (iUnit8, 50) 50 FORMAT (/ & & 77X,' mantle'/ & & 77X,' crustal lithosphere'/ & & ' node E-longitude N-latitude', & & ' theta phi elevation', & & ' heat-flow thickness thickness'/) END IF DO 90 k = 1, numNod checkN(k) = .FALSE. 90 CONTINUE DO 100 k = 1, numNod CALL ReadN (iUnit7, iUnit8, 9, & ! input & vector) ! output index = vector(1) + 0.5D0 IF (index > nRealN) THEN IF ((index <= n1000).OR. & & (index > (n1000 + nFakeN))) THEN WRITE (iUnit8, 91) index 91 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER: ',I6) CALL Pause() STOP END IF END IF pLon = vector(2) pLat = vector(3) IF (ABS(pLat) > 90.01) THEN WRITE (iUnit8, 92) index 92 FORMAT (/' ERR0R: ABS(latitude) > 90 AT NODE ',I6) CALL Pause() STOP END IF IF (ABS(pLat) > 89.99D0) THEN WRITE (iUnit8, 93) index 93 FORMAT (/' ERR0R: NODE ',I6,' LIES ON A POLE.' & & /' THIS IS A SINGULAR POINT OF THE' & & ,' SPHERICAL COORDINATE SYSTEM.' & & /' MOVE THIS NODE, AT LEAST SLIGHTLY.') CALL Pause() STOP END IF xi = (90.0D0 - pLat) * 0.0174532925199433D0 yi = pLon * 0.0174532925199433D0 elevi = vector(4) qi = vector(5) zmi = vector(6) tli = vector(7) density_anomaly_kgpm3 = vector(8) cooling_curvature_cpm2 = vector(9) IF (index <= nRealN) THEN i = index ELSE i = nRealN + index - n1000 END IF checkN(i) = .TRUE. xNode(i) = xi yNode(i) = yi elev(i) = elevi dQdTdA(i) = qi IF (qi < 0.0D0) THEN WRITE (iUnit8, 96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL Pause() STOP END IF IF (zmi < 0.0D0) THEN WRITE (iUnit8, 97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') CALL Pause() STOP END IF zMNode(i) = zmi IF (tli < 0.0D0) THEN WRITE (iUnit8, 98) 98 FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS IS', & & ' NON-PHYSICAL.') CALL Pause() STOP END IF tLNode(i) = tli IF (.NOT.brief) THEN WRITE (iUnit8, 99) INDEX, pLon, pLat, xi, yi, elevi, & & qi, zmi, tli 99 FORMAT (' ',I10,0P,2F12.3,2F11.5,1P,3E10.2,E12.2) END IF density_anomaly(i) = density_anomaly_kgpm3 cooling_curvature(i) = cooling_curvature_cpm2 100 CONTINUE allOK = .TRUE. DO 101 i = 1, numNod allOK = allOK.AND.checkN(i) 101 CONTINUE IF (.NOT.allOK) THEN WRITE (iUnit8, 102) 102 FORMAT(' THE FOLLOWING NODES WERE NEVER READ:') DO 104 i = 1, numNod IF (i <= nRealN) THEN index = i ELSE index = n1000 + i - nRealN END IF IF (.NOT.checkN(i)) WRITE(iUnit8, 103) index 103 FORMAT (' ',36X,I6) 104 CONTINUE CALL Pause() STOP END IF ! Read triangular elements: READ (iUnit7, * ) numEl IF (numEl > mxEl) THEN WRITE (iUnit8, 108) numEl 108 FORMAT(/' INCREASE PARAMETER maxEl TO BE AT LEAST EQUAL' & & /' TO THE NUMBER OF ELEMENTS (',I6,') AND RECOMPILE.') CALL Pause() STOP END IF DO 109 k = 1, numEl checkE(k) = .FALSE. 109 CONTINUE IF (.NOT.brief) THEN WRITE (iUnit8, 110) numEl 110 FORMAT(/' There are ',I6,' triangular continuum elements.'/ & & ' (Node numbers for each are given at corners, counter', & & 'clockwise'/ / & & ' element C1 C2 C3 Lithospheric_Rheology#') END IF DO 200 k = 1, numEl ! (Elements need not be input in order, but must all be present.) READ (iUnit7, "(A)") longer_line CALL Extract_LRi(longer_line, LRi, shorter_line) continuum_LRi(k) = LRi READ (shorter_line, * ) i, (nodes(j, i), j = 1, 3) IF ((i < 1).OR.(i > numEl)) THEN WRITE (iUnit8, 115) i 115 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I6) CALL Pause() STOP END IF checkE(i) = .TRUE. IF (.NOT.brief) THEN IF (LRi == 0) THEN WRITE (iUnit8, 120) i, (nodes(j, i), j = 1, 3) 120 FORMAT (' ', I6, ':', 3I10) ELSE WRITE (iUnit8, 121) i, (nodes(j, i), j = 1, 3), LRi 121 FORMAT (' ', I6, ':', 3I10, ' LR', I8) END IF END IF DO 130 j = 1, 3 n = nodes(j, i) IF (n > nRealN) n = nRealN + (n - n1000) IF ((n <= 0).OR.(n > numNod)) THEN WRITE (iUnit8, 125) nodes(j, i) 125 FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.') CALL Pause() STOP END IF nodes(j, i) = n 130 CONTINUE 200 CONTINUE allOK = .TRUE. DO 201 i = 1, numEl allOK = allOK.AND.checkE(i) 201 CONTINUE IF (.NOT.allOK) THEN WRITE (iUnit8, 202) 202 FORMAT (' THE FOLLOWING ELEMENTS WERE NEVER READ:') DO 204 i = 1, numEl IF (.NOT.checkE(i)) WRITE(iUnit8, 203) i 203 FORMAT (' ',39X,I6) 204 CONTINUE CALL Pause() STOP END IF ! Read fault elements: READ (iUnit7, * ) nFl IF (nFl > mxFEl) THEN WRITE (iUnit8, 220)nFl 220 FORMAT (/' INCREASE PARAMETER maxFEl TO BE AT LEAST EQUAL' & & /' TO THE NUMBER OF FAULTS (',I6,') AND RECOMPILE.') CALL Pause() STOP END IF offMax = 0.0D0 DO 222 i = 1, nFl checkF(i) = .FALSE. 222 CONTINUE IF (.NOT.brief) WRITE(iUnit8, 230) nFl 230 FORMAT(/ /' There are ', I6, ' great-circle fault elements.') IF ((.NOT.brief).AND.(nFl > 0)) WRITE(iUnit8, 231) 231 FORMAT (/' (The 4 node numbers defining each element must be', & & ' in a counterclockwise order:'/ & & ' n1, and n2 are in left-to-right seguence on the', & & ' near side,'/ & & ' then n3 is opposite n2, and n4 is opposite n1.'/, & & ' (Fault dips are given at n1, n2,', & & ' in degrees from horizontal;'/ & & ' positive dips are toward n1 and n2, respectively, '/ & & ' while negative dips are toward n4 and n3.)'/ & & ' Offset is the total past slip of the fault.'/ / & & ' Element n1 n2 n3 n4 dip1 dip2', & & ' offset Lithospheric_Rheology#'/) 240 FORMAT (' ', I6, ':', 4I5, 1X, 2F6.1, 1X, F9.0) DO 300 k = 1, nFl off = 0.0D0 READ(iUnit7, "(A)") longer_line CALL Extract_LRi(longer_line, LRi, shorter_line) fault_LRi(k) = LRi READ(shorter_line, * ) i, (nodeF(j, k), j = 1, 4), (dips(l), l = 1, 2), off IF ((i < 1).OR.(i > nFl)) THEN WRITE (iUnit8, 241) i 241 FORMAT (/' ERR0R: ILLEGAL FAULT NUMBER: ', I6) CALL Pause() STOP END IF checkF(i) = .TRUE. IF (.NOT.brief) THEN IF (LRi == 0) THEN WRITE (iUnit8, 240) i, (nodeF(j, i), j = 1, 4), (dips(l), l = 1, 2), off ELSE WRITE (iUnit8, 242) i, (nodeF(j, i), j = 1, 4), (dips(l), l = 1, 2), off, LRi 242 FORMAT (' ', I6, ':', 4I5, 1X, 2F6.1, 1X, F9.0, " LR", I8) END IF END IF DO 250 j = 1, 4 n = nodeF(j, i) IF (n > nRealN) n = nRealN + (n - n1000) IF ((n <= 0).OR.(n > numNod)) THEN WRITE (iUnit8, 243) nodeF(j, i), i 243 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER (', I6, ') IN FAULT ',I6) CALL Pause() STOP END IF nodeF(j, i) = n 250 CONTINUE DO 260 l = 1, 2 IF (ABS(dips(l)) > 90.0D0) THEN WRITE(iUnit8, 252) dips(l) 252 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN', & & ' RANGE OF -90. TO +90. DEGREES.'/ & & ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE', & & ' HORIZONAL;'/ & & ' A + PREFIX (OR NONE) INDICATES A DIP', & & ' TOWARD THE n1-n2 SIDE;'/ & & ' A - PREFIX INDICATES A DIP TOWARD', & & ' THE n4-n3 SIDE.)') CALL Pause() STOP END IF IF (dips(l) < 0.0D0) dips(l) = 180.0D0 + dips(l) fDip(l, i) = dips(l) * 0.0174532925199433D0 260 CONTINUE IF (off < 0.0D0) THEN WRITE (iUnit8, 280) off, k 280 FORMAT (' ILLEGAL FAULT OFFSET OF ', ES10.2, & & ' FOR FAULT ELEMENT', I6 / & & ' OFFSETS MAY NOT BE NEGATIVE.') CALL Pause() STOP END IF offset(i) = off offMax = MAX(offMax, off) 300 CONTINUE allOK = .TRUE. DO 301 i = 1, nFl allOK = allOK.AND.checkF(i) 301 CONTINUE IF (.NOT.allOK) THEN WRITE(iUnit8, 302) 302 FORMAT(' THE FOLLOWING FAULTS WERE NEVER READ:') DO 304 i = 1, nFl IF (.NOT.checkF(i)) WRITE(iUnit8, 303) i 303 FORMAT(' ',36X,I6) 304 CONTINUE CALL Pause() STOP ELSE IF (offMax > 0.0D0) THEN WRITE (iUnit8, 400) offMax 400 FORMAT (/' Greatest fault offset read was ',1P,D10.2) ELSE WRITE (iUnit8, 401) 401 FORMAT (/' Since fault offsets are all zero,', & & ' input parameter Byerly will have no effect.') END IF END IF IF (.NOT. brief) WRITE (iUnit8, 999) 999 FORMAT (' --------------------------------------------------', & & '-----------------------------') END SUBROUTINE GetNet SUBROUTINE PAUSE() IMPLICIT NONE WRITE(*, "(' Press [Enter]..'\)") READ(*, * ) END SUBROUTINE Pause SUBROUTINE PutNet (iUnitO, & ! INTENT(IN) & brief, & ! INTENT(IN) & continuum_LRi, & ! INTENT(IN) & dQdTdA, elev, & ! INTENT(IN) & fault_LRi, fDip, & ! INTENT(IN) & mxEl, mxFEl, mxNode, n1000, & ! INTENT(IN) & nFakeN, nFl, nodeF, nodes, & ! INTENT(IN) & nRealN, numEl, numNod, offset, & ! INTENT(IN) & title1, tLNode, xNode, yNode, zMNode, & ! INTENT(IN) & chemical_delta_rho_list, & ! INTENT(IN) & cooling_curvature_list) ! INTENT(IN) ! Writes finite-element grid (.FEG file) to unit "iUnitO". IMPLICIT NONE !Arguments: INTEGER, INTENT(IN) :: iUnitO LOGICAL, INTENT(IN) :: brief INTEGER, INTENT(IN) :: continuum_LRi REAL*8, INTENT(IN) :: dQdTdA, elev INTEGER, INTENT(IN) :: fault_LRi REAL*8, INTENT(IN) :: fDip INTEGER, INTENT(IN) :: mxEl, mxFEl, mxNode, n1000, & & nFakeN, nFl, nodeF, nodes, & & nRealN, numEl, numNod REAL*8, INTENT(IN) :: offset CHARACTER*80, INTENT(IN) :: title1 REAL*8, INTENT(IN) :: tLNode, xNode, yNode, zMNode, & & chemical_delta_rho_list, & & cooling_curvature_list DIMENSION chemical_delta_rho_list(mxNode), & & continuum_LRi(mxEl), & & cooling_curvature_list(mxNode), & & dQdTdA(mxNode), elev(mxNode), & & fault_LRi(mxFEl), & & fDip(2, mxFEl), & & nodeF(4, mxFEl), nodes(3, mxEl), & & offset(mxFEl), tLNode(mxNode), & & xNode(mxNode), yNode(mxNode), zMNode(mxNode) !Internal variables: CHARACTER*8 :: LR_c8 INTEGER :: i, iPrint, k, LRi INTEGER, DIMENSION(4) :: nP REAL*8 :: pLat, pLon REAL*8, DIMENSION(2) :: dips WRITE(*, "(/ /' READY TO CREATE OUTPUT .FEG FILE ON UNIT ', I3)") iUnitO WRITE(*, "(' Please enter a (new) filename with "".feg"" ending: ')") WRITE (iUnitO, 1) title1 1 FORMAT (A80) WRITE (iUnitO, 2) numNod, nRealN, nFakeN, n1000, brief 2 FORMAT(4I8, L8, ' (numNod, nRealN, nFakeN, N1000, BRIEF)') DO 100 i = 1, numNod IF (i <= nRealN) THEN iPrint = i ELSE iPrint = n1000 + (i - nRealN) END IF pLat = 90.0D0 - xNode(i) * 57.2957795130823D0 pLon = yNode(i) * 57.2957795130823D0 WRITE (iUnitO, 91) i, pLon, pLat, elev(i), dQdTdA(i), zMNode(i), & & tLNode(i), chemical_delta_rho_list(i), & & cooling_curvature_list(i) 91 FORMAT (I8, 2F11.5, 6ES10.2) 100 CONTINUE WRITE (iUnitO, 110) numEl 110 FORMAT (I10,' (numEl = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') DO 200 i = 1, numEl DO 150 k = 1, 3 IF (nodes(k, i) <= nRealN) THEN nP(k) = nodes(k, i) ELSE nP(k) = n1000 + (nodes(k, i) - nRealN) END IF 150 CONTINUE LRi = continuum_LRi(i) IF (LRi == 0) THEN WRITE (iUnitO, 160) i, (nP(k), k = 1, 3) 160 FORMAT (I8, 3I8) ELSE WRITE (LR_c8, "(I8)") LRi LR_c8 = ADJUSTL(LR_c8) WRITE (iUnitO, "(I8, 3I8, ' LR', A)") i, (nP(k), k = 1, 3), TRIM(LR_c8) END IF 200 CONTINUE WRITE (iUnitO, 210) nFl 210 FORMAT (I10,' (nFl = NUMBER OF CURVILINEAR FAULT ELEMENTS)') DO 300 i = 1, nFl DO 220 k = 1, 4 IF (nodeF(k, i) <= nRealN) THEN nP(k) = nodeF(k, i) ELSE nP(k) = n1000 + (nodeF(k, i) - nRealN) END IF 220 CONTINUE DO 230 k = 1, 2 dips(k) = fDip(k, i) dips(k) = dips(k) * 57.2957795130823D0 IF (dips(k) > 90.01D0) dips(k) = dips(k) - 180.0D0 230 CONTINUE LRi = fault_LRi(i) IF (LRi == 0) THEN WRITE (iUnitO, 250) i, (nP(k), k = 1, 4), (dips(k), k = 1, 2), offset(i) 250 FORMAT (I8, 4I8, 2F6.1, ES10.2) ELSE WRITE (LR_c8, "(I8)") LRi LR_c8 = ADJUSTL(LR_c8) WRITE (iUnitO, "(I8, 4I8, 2F6.1, ES10.2, ' LR', A)") i, (nP(k), k = 1, 4), (dips(k), k = 1, 2), offset(i), TRIM(LR_c8) END IF 300 CONTINUE END SUBROUTINE PutNet SUBROUTINE ReadN (iUnitP, iUnitT, n, & ! input & vector) ! output ! A utility routine designed to permit -Faults- input files ! to also be used by -Plates-, which expects more numbers ! in some records. ! This routine attempts to READ "n" floating-point values ! (using * FORMAT) from the next record on device iUnitP. ! If anything goes wrong, the missing values are set to zero. ! Since June 2005, this feature also allows the reading of ! both old-format (OrbData) .feg files (with 4 nodal data), ! and new-format (OrbData5) .feg files (with 6 nodal data), ! by program -Shells-. IMPLICIT NONE ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER, INTENT(IN) :: iUnitP, iUnitT, n ! input REAL*8, INTENT(OUT) :: vector ! output ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CHARACTER*1 c CHARACTER*132 line INTEGER i, number LOGICAL anyIn, dotted, expon, signed DIMENSION vector(n) line = ' ' READ (iUnitP, 1) line 1 FORMAT (A) number = 0 anyIn = .FALSE. expon = .FALSE. signed = .FALSE. dotted = .FALSE. DO 10 i = 1, 132 c = line(i:i) IF ((c == ' ').OR.(c == ',').OR.(c == '/')) THEN signed = .FALSE. expon = .FALSE. dotted = .FALSE. IF (anyIn) THEN number = number + 1 anyIn = .FALSE. END IF ELSE IF ((c == '+').OR.(c == '-')) THEN IF (signed) THEN GO TO 50 ELSE signed = .TRUE. END IF ELSE IF ((c == 'E').OR.(c == 'D').OR. & & (c == 'e').OR.(c == 'd')) THEN IF (expon) THEN GO TO 50 ELSE expon = .TRUE. signed = .FALSE. dotted = .TRUE. END IF ELSE IF (c == '.') THEN IF (dotted) THEN GO TO 50 ELSE dotted = .TRUE. END IF ELSE IF ((c == '0').OR.(c == '1').OR.(c == '2').OR. & & (c == '3').OR.(c == '4').OR.(c == '5').OR. & & (c == '6').OR.(c == '7').OR.(c == '8').OR. & & (c == '9')) THEN signed = .TRUE. anyIn = .TRUE. ELSE GO TO 50 END IF 10 CONTINUE IF (anyIn) number = number + 1 50 IF (number == 0) THEN WRITE (iUnitT, 91) n, line 91 FORMAT (/' ERR0R: A LINE OF PARAMETER INPUT WHICH', & & ' WAS SUPPOSED TO CONTAIN 1-',I2,' NUMBERS'/ & & ' COULD NOT BE INTERPRETED. LINE FOLLOWS:'/ & & ' ',A80) CALL Pause() STOP ELSE number = MIN(number, n) BACKSPACE iUnitP READ (iUnitP, * ) (vector(i), i = 1, number) IF (number < n) THEN DO 99 i = number + 1, n vector(i) = 0.0D0 99 CONTINUE END IF END IF RETURN END SUBROUTINE ReadN LOGICAL FUNCTION Within(uvec, numPoly, maxPoly, plate_outline_uvecs, iPlate, outline_count) ! Determines whether uvec is inside the circuit of plate_outline_uvecs(1:outline_count), ! where the convention is that plate_outline_uvecs(1) == plate_outline_uvecs(outline_count). USE DSphere ! Fortran MODULE DSphere is in file DSphere.f90, provided by Peter Bird of UCLA. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec INTEGER, INTENT(IN) :: numPoly, maxPoly REAL*8, DIMENSION(3, numPoly, maxPoly), INTENT(IN) :: plate_outline_uvecs INTEGER, INTENT(IN) :: iPlate ! which of the outlines are we using in this test? INTEGER, INTENT(IN) :: outline_count ! number of digitization points in this particular outline INTEGER :: i REAL*8, DIMENSION(3) :: tuvec_0, tuvec_1 REAL*8 :: angle_0, angle_1, angle_sum, d_angle angle_sum = 0.0D0 tuvec_0(1:3) = plate_outline_uvecs(1:3, iPlate, 1) angle_0 = DRelative_Compass(from_uvec = uvec, to_uvec = tuvec_0) !result is azimuth, clockwise from N, in radians DO i = 2, outline_count tuvec_1(1:3) = plate_outline_uvecs(1:3, iPlate, i) angle_1 = DRelative_Compass(from_uvec = uvec, to_uvec = tuvec_1) !If uvec is inside, then typically angle_1 < angle_0 (except for cycle shifts) d_angle = -(angle_1 - angle_0) ! reversing sign, so d_angle will typically be positive if uvec is inside. d_angle = ATAN2(SIN(d_angle), COS(d_angle)) ! getting rid any cycle shifts! angle_sum = angle_sum + d_angle !prepare for next plate-boundary step: tuvec_0 = tuvec_1 angle_0 = angle_1 END DO !If uvec is inside, then angle_sum should be somewhere close to 2*Pi. Within = (angle_sum > 3.0D0) .AND. (angle_sum < 9.0D0) !but angle_sum will be either ~0.0 or around -2.0*Pi, if point is outside. END FUNCTION Within END PROGRAM Assign_LRs