PROGRAM OrbNumber ! Mature version of 13 February 2019: ! *supports OrbData5 format (up to 6 nodal data); ! *supports Restore3+ format (up to 3 element data); ! *supports Shells_v5.0+ format (with option LR# for each element); ! *also provides a list of boundary nodes; and (most exciting) ! *has a QUICK option for cases where a good choice of node #1 ! has already been made (by a previous complete optimization). ! [Note that this last change occurs entirely within subprogram Number().] ! by Peter Bird ! Department of Earth, Planetary, and Space Sciences ! University of California ! Los Angeles, CA 90095-1567 ! pbird@epss.ucla.edu ! (C)Copyright 1993, 1994, 1996, 1998, 1999, 2005, 2010, 2013, 2018, 2019 ! by Peter Bird and The ! Regents of the University of California. ! Program to read a finite element grid (.feg) file, ! renumber the nodes to minimize bandwidth, ! and write out the result in the same format. ! Accepts input files in the format(s) created by OrbWin. ! These .feg files may be of Shells-type (with faults, ! and perhaps with nodal data filled in by ! OrbData or OrbData5). ! Or, they may be of NeoKinema type (without faults, ! and with no nodal data, except optional mu_ values). ! Or, they may be of Restore type (without faults, ! but with optional element data in the case of Restore3+). ! In the case of .feg files from OrbData5 (June 2005+) with ! 6 nodal variables (instead of the previous 4) it recognizes ! the new format silently and passes along all the values. ! In the case of .feg files for Restore v.3, with element data, ! it recognizes and passes along the element data. ! This utility program also provides an additional output file ! with boundary nodes listed in counterclockwise order ! around the perimeter of the grid; this file can be ! hand-edited to create boundary-conditions files like ! those required by Shells, NeoKinema, & Restore. !======================================================================= ! ! N.B. The following note refers to ANCIENT gridding practices from ! the decade of the 1980's, when computer memory was really ! a serious limit (e.g., in original PC's, running DOS). ! In the modern era, my (newer) programs no longer make ! any distinction between "real" and "fake" nodes; ! instead, all nodes are treated as "real". ! Also, we leave "redeem=.FALSE." to control the output mode. ! ! *** Note: Set variable "redeem" to True or False. ! If True, it means that fake nodes (boundary nodes with no ! mechanical degrees of freedom) will be converted to real ! nodes during the renumbering. This option is useful to ! prepare a renumbering for use by the graphics programs. ! If "redeem" is False, fake nodes are not renumbered, and ! will appear at the end of the list as before. ! *** Note: ! "redeem" also controls the type of output file; if True, ! then only a list of old node numbers, new node numbers, and ! a cross-reference (old numbers listed in order of new) is ! output on unit "iUnitO". If "redeem" is False, however, a ! complete, resorted finite element grid file is written. ! This reflects the practical reality that one usually ! chooses "redeem=.TRUE." when preparing to do graphics, but ! chooses "redeem=.FALSE." when optimizing a hand-edited grid ! prior to finite element calculations. !======================================================================= IMPLICIT NONE INTEGER, PARAMETER :: maxNod = 50000, maxBn = 4000, maxEl = 90000, maxFEl = 9000 INTEGER, PARAMETER :: maxGSi = 1250000000 CHARACTER*80 :: title1 INTEGER :: continuum_LRi, fault_LRi, & & i, iAlias, iDiag, ios, iUnitB, iUnitI, iUnitO, iUnitT, iUser, list, LRi, & & maxDif, maxLst, minDif, mxBn, mxStar, mxNode, mxEl, mxFEl, mxGSiz, & & n34, n1000, nFakeN, nFl, nCond, nodCon, nodeF, nodes, nRealN, nToDo, numNod, numEl LOGICAL*1 :: checkE, checkF, checkN, doNext, edgeFs, edgeTs, lGraph LOGICAL :: brief, redeem, sphere, orbData5 REAL*8 :: atNode, cooling_curvature, density_anomaly, dQdtdA, eleDat, elev, & & fDip, offMax, offset, pLat, pLon, tInSec, tLNode, xNode, yNode, zMNode DIMENSION atNode(maxNod), & & checkE(maxEl), checkF(maxFEl), checkN(maxNod), & & continuum_LRi(maxEl), cooling_curvature(maxNod), & & density_anomaly(maxNod), & & doNext(maxNod), dQdtdA(maxNod), & & eleDat(3, maxEl), elev(maxNod), & & edgeFs(2, maxFEl), edgeTs(3, maxEl), & & fault_LRi(maxFEl), fDip(2, maxFEl), & & iAlias(maxNod), iDiag(maxNod), & & iUser(maxNod), lGraph(maxGSi), list(maxNod), & & maxLst(maxNod), nodCon(maxBn), & & nodeF(4, maxFEl), nodes(3, maxEl), & & offset(maxFEl), tLNode(maxNod), & & xNode(maxNod), yNode(maxNod), zMNode(maxNod) ! The input finite=element grid file (.FEG file) is read from: DATA iUnitI / 1 / ! The output finite-element grid file (.FEG file) is written to: DATA iUnitO / 2 / ! The output list of boundary nodes (.PER file) is written to: DATA iUnitB / 3 / ! Text messages during execution are sent to: DATA iUnitT / 6 / !************************ Check the following line! ******************* ! Use .TRUE. to create only a list of new node numbers versus old. redeem = .FALSE. ! Use .FALSE. to create a complete renumbered .FEG file. !************************ Check the preceding line! ******************* mxNode = maxNod mxBn = maxBn mxEl = maxEl mxFEl = maxFEl mxStar = maxNod mxGSiz = maxGSi WRITE (iUnitT, "(' ')") WRITE (iUnitT, "(' PROGRAM OrbNumber')") WRITE (iUnitT, "(' version of 23 July 2020')") WRITE (iUnitT, "(' by Peter Bird, UCLA.')") WRITE (iUnitT, 1) 1 FORMAT (/' Attempting to read existing .FEG file from UNIT 1..'/) CALL GetNet (iUnitI, iUnitT, & ! inputs & mxEl, mxFEl, mxNode, & ! inputs & brief, & ! outputs & continuum_LRi, & ! outputs & cooling_curvature, density_anomaly, & ! outputs & dQdtdA, eleDat, elev, & ! outputs & fault_LRi, fDip, & ! outputs & nFl, nodeF, nodes, nFakeN, nRealN, n1000, & ! outputs & numEl, numNod, offMax, offset, orbData5, & ! outputs & title1, tLNode, xNode, yNode, zMNode, & ! outputs & checkE, checkF, checkN) ! work-space arrays WRITE (iUnitT, 2) 2 FORMAT(' This .FEG file has been read.') WRITE (iUnitT, 3) 3 FORMAT (/' Square is checking for correct topology...') CALL Square (brief, fDip, iUnitT, & ! inputs & mxBn, mxEl, mxFEl, mxNode, & ! inputs & mxStar, nFl, nodeF, nodes, & ! inputs & numEl, numNod, & ! inputs & xNode, yNode, & ! modify & edgeFs, & ! outputs & edgeTs, & ! outputs & nCond, nodCon, & ! outputs & checkN, list) ! work-space arrays WRITE (iUnitT, 151) 151 FORMAT (' ... done.') IF (redeem) THEN nToDo = numNod ELSE nToDo = nRealN END IF IF (nCond == 0) THEN sphere = .TRUE. IF (nFl == 0) THEN ! Arbitrarily choosing to begin with node 1: nCond = 1 nodCon(1) = 1 ELSE ! Or, arbitrarily choose a set of initial nodes: nCond = MIN(300, nRealN) nodCon(1) = 1 nodCon(nCond) = nRealN DO 5 i = 2, (nCond - 1) nodCon(i) = nRealN * (1. * i) / (1. * nCond) + 0.5 5 CONTINUE END IF ELSE sphere = .FALSE. END IF WRITE (iUnitT, 10) 10 FORMAT (/' Calling Number...') CALL Clock (iUnitT, & ! input & tInSec) ! modify CALL Number (iUnitT, nCond, nodCon, & ! inputs & mxEl, mxFEl, mxGSiz, mxNode, & ! inputs & nFl, nodeF, nodes, nToDo, & ! inputs & numEl, numNod, sphere, xNode, yNode, & ! inputs & tInSec, & ! modify & iAlias, iUser, maxDif, minDif, & ! outputs & doNext, iDiag, lGraph, maxLst) ! work-space arrays WRITE (iUnitT, 151) CALL Clock (iUnitT, & ! input & tInSec) ! modify WRITE (*, "(' which only counts time spent in SUBROUTINE Number.'/)") WRITE (iUnitT, *) WRITE (iUnitT, "(' =====================================')") WRITE (iUnitT, 15) maxDif 15 FORMAT (' BANDWIDTH BEFORE RENUMBERING = ',I6) WRITE (iUnitT, 20) minDif 20 FORMAT (' BANDWIDTH AFTER RENUMBERING = ',I6) WRITE (iUnitT, "(' =====================================')") IF (redeem) THEN DO 22 i = 1, numNod WRITE (iUnitO, 21) i, iAlias(i), iUser(i) 21 FORMAT (3I10) 22 CONTINUE ELSE ! Reorder the values in the simple floating-point vector arrays: CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & xNode, & ! modify & atNode) ! work-space array CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & yNode, & ! modify & atNode) ! work-space array CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & elev, & ! modify & atNode) ! work-space array CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & dQdtdA, & ! modify & atNode) ! work-space array CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & zMNode, & ! modify & atNode) ! work-space array CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & tLNode, & ! modify & atNode) ! work-space array IF (orbData5) THEN CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & density_anomaly, & ! modify & atNode) ! work-space array CALL MixUpX (iAlias, mxNode, numNod, & ! inputs & cooling_curvature, & ! modify & atNode) ! work-space array END IF ! Modify the values in the n-entry integer arrays: n34 = 3 CALL MixUpI (iAlias, mxEl, mxNode, numEl, n34, & ! inputs & nodes) ! modify n34 = 4 CALL MixUpI (iAlias, mxFEl, mxNode, nFl, n34, & ! inputs & nodeF) ! modify ! Output modified file on device "iUnitO": WRITE (iUnitT, 90) iUnitO 90 FORMAT (' Renumbering completed.'/ & & ' =========================================='/ & & ' About to write output grid (.feg) file' & & ,' on unit ', I2, '...'/) CALL PutNet (iUnitO, & ! inputs & brief, & ! inputs & continuum_LRi, & ! inputs & cooling_curvature, density_anomaly, & ! inputs & dQdtdA, eleDat, elev, & ! inputs & fault_LRi, fDip, & ! inputs & mxEl, mxFEl, mxNode, n1000, & ! inputs & nFakeN, nFl, nodeF, nodes, & ! inputs & nRealN, numEl, numNod, offset, orbData5, & ! inputs & title1, tLNode, xNode, yNode, zMNode) ! inputs WRITE (iUnitT, *) WRITE (iUnitT, "(' Renumbered .FEG file has been written.')") END IF CALL PAUSE() ! CREATE ORDERED LIST OF BOUNDARY NODES, USING NEW NODE NUMBERS: IF (.NOT.sphere) THEN CALL Square (brief, fDip, iUnitT, & ! inputs & mxBn, mxEl, mxFEl, mxNode, & ! inputs & mxStar, nFl, nodeF, nodes, & ! inputs & numEl, numNod, & ! inputs & xNode, yNode, & ! modify & edgeFs, & ! outputs & edgeTs, & ! outputs & nCond, nodCon, & ! outputs & checkN, list) ! work-sace arrays WRITE (iUnitT, 101) nCond 101 FORMAT (/' There are ',I6,' nodes along the boundary.' & & /' Now, creating a new perimeter (.PER) file containing an ordered list,' & & /' of boundary nodes which can be used as the basis for a boundary-' & & /' conditions file: '/) WRITE (iUnitB, 103) 103 FORMAT (' BC# Node Latitude Longitude') DO 110 i = 1, nCond pLat = 90. - xNode(nodCon(i)) * 57.29577951 pLon = yNode(nodCon(i)) * 57.29577951 WRITE (iUnitB, 104) i, nodCon(i), pLat, pLon 104 FORMAT (I6, I6, F9.3, F10.3) 110 CONTINUE END IF ! .NOT.sphere WRITE (iUnitT, 111) 111 FORMAT (/' Job completed.') CALL PAUSE() STOP END PROGRAM OrbNumber 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 previous versions of OrbNumber. 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, & ! inputs & mxEl, mxFEl, mxNode, & ! inputs & brief, & ! outputs, & continuum_LRi, & ! outputs, & cooling_curvature, density_anomaly, & ! outputs, & dQdtdA, eleDat, elev, & ! outputs, & fault_LRi, fDip, & ! outputs, & nFl, nodeF, nodes, nFakeN, nRealN, n1000, & ! outputs, & numEl, numNod, offMax, offset, orbData5, & ! outputs, & title1, tLNode, xNode, yNode, zMNode, & ! outputs, & checkE, checkF, checkN) ! work-space arrays ! Read finite element grid from unit "iUnit7". ! Echo the important values to a print dataset on unit "iUnit8". IMPLICIT NONE INTEGER, INTENT(IN) :: iUnit7, iUnit8, mxEl, mxFEl, mxNode LOGICAL, INTENT(OUT) :: brief REAL*8, INTENT(OUT) :: cooling_curvature, density_anomaly, & & dQdtdA, eleDat, elev, fDip, offMax, offset, & & tLNode, xNode, yNode, zMNode INTEGER, INTENT(OUT) :: nFl, nodeF, nodes, nFakeN, nRealN, n1000, & & numEl, numNod INTEGER, INTENT(OUT) :: continuum_LRi, fault_LRi ! DIMENSIONed below... LOGICAL, INTENT(OUT) :: orbData5 CHARACTER*80, INTENT(OUT) :: title1 CHARACTER*80 :: longer_line, shorter_line ! for internal use, in CALL Extract_LRi... INTEGER :: i, index, ios, j, k, l, LRi, n LOGICAL*1 :: checkE, checkF, checkN ! work-space arrays LOGICAL :: allOK REAL*8 :: dips, elevI, off, pLat, pLon, qI, & & temp_density_anomaly, temp_cooling_curvature, & & tLI, vector, xI, yI, zMI DIMENSION checke(mxEl), checkf(mxFEl), checkN(mxNode), & & continuum_LRi(mxEl), & & cooling_curvature(mxNode), density_anomaly(mxNode), & & dQdtdA(mxNode), eleDat(3, mxEl), 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) WRITE (iUnit8, 1) iUnit7 1 FORMAT (' Attempting to read finite element grid from unit',I3) title1 = ' '// & & ' ' READ (iUnit7, 2, IOSTAT = ios) title1 2 FORMAT (A80) WRITE (iUnit8, 3) title1 3 FORMAT(/' Title of finite element grid ='/' ',A) ! Read number of nodes. ! Input nodal locations (x,y), and perhaps elevations, heat-flow, and ! perhaps other nodal data (crustal thickness, mantle-lithosphere thickness, ...). ! (Option "brief" suppresses most output.) READ (iUnit7, * ) numNod, nRealN, nFakeN, n1000, brief brief = .TRUE. 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 PARAMETER MAXNOD TO BE AT LEAST' & & /' THE NUMBER OF NODES (',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. Be careful!!!)') 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 orbData5 = .FALSE. ! unless proven otherwise, below, by non-zero input values... DO 100 k = 1, numNod CALL ReadN (iUnit7, iUnit8, 9, & ! inputs & vector) ! outputs 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) xI = (90.0D0 - pLat) * 0.0174532925199433D0 yI = pLon * 0.0174532925199433D0 elevI = vector(4) qI = vector(5) zMI = vector(6) tLI = vector(7) temp_density_anomaly = vector(8) ! Note: If READN did not find these last 2 values, it would ha temp_cooling_curvature = 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 density_anomaly(i) = temp_density_anomaly cooling_curvature(i) = temp_cooling_curvature orbData5 = orbData5 .OR. (temp_density_anomaly /= 0.0D0) .OR. & & (temp_cooling_curvature /= 0.0D0) 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 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, and perhaps optional element data {(3 numbers, for Restore3+) OR (LR#, for Shells_v5.0+)}? READ (iUnit7, * , iostat = ios) 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) 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') DO 200 k = 1, numEl ! (Elements need not be input in order, but all must be present.) READ (iUnit7, "(A)") longer_line CALL Extract_LRi (longer_line, & ! input & LRi, shorter_line) ! output CALL ReadNFromString (shorter_line, iUnit8, 7, & ! inputs [N.B., "7" allows for up to 3 optional element data.] & vector) ! outputs i = NINT(vector(1)) IF ((i < 1).OR.(i > numEl)) THEN WRITE (iUnit8, 111) i 111 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I6) CALL PAUSE() STOP END IF nodes(1, i) = NINT(vector(2)) nodes(2, i) = NINT(vector(3)) nodes(3, i) = NINT(vector(4)) eleDat(1, i) = vector(5) eleDat(2, i) = vector(6) eleDat(3, i) = vector(7) checke(i) = .TRUE. continuum_LRi(i) = LRi IF (.NOT.brief) THEN IF (LRi == 0) THEN WRITE (iUnit8, "(' ', I6, ':', 3I10, 3ES10.2)") i, (nodes(j, i), j = 1, 3), (eleDat(j, i), j = 1, 3) ELSE WRITE (iUnit8, "(' ', I6, ':', 3I10, 3ES10.2, ' LR', I8)") i, (nodes(j, i), j = 1, 3), (eleDat(j, i), j = 1, 3), LRi 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 sequence on the', & & ' near side,'/ & & ' then N3 is opposite N2, 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.)'/ & & ' (The argument of the fault trace is given at N1 and N2,'/ & & ' in degrees counterclockwise from the theta (S) axis.)'/ & & ' Offset is the total past slip of the fault.'/ / & & ' ELEMENT N1 N2 N3 N4 DIP1 DIP2', & & ' OFFSET'/) 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, & ! input & LRi, shorter_line) ! output 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, "(' ', I6, ':', 4I5, 1X, 2F6.1, 1X, F9.0), ' LR', I8)") i, (nodeF(j, i), j = 1, 4), (dips(l), l = 1, 2), off, LRi 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.) 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.) dips(l) = 180.0D0 + dips(l) fDip(l, i) = dips(l) * 0.0174532925199433D0 260 CONTINUE IF (off < 0.) THEN WRITE (iUnit8, 280) off 280 FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2, & & ' FOR FAULT ELEMENT',I6/ & & ' OFFSETS MAY NOT BE NEGATIVE.') CALL PAUSE() STOP END IF offset(i) = off offmax = MAX(offmax, off) fault_LRi(i) = LRi 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, I8) 304 CONTINUE CALL PAUSE() STOP ELSE IF (nFl > 0) THEN IF (offmax > 0.0D0) THEN WRITE (iUnit8, 400) offmax 400 FORMAT (/' Greatest fault offset read was ', ES10.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 PutNet (iUnitO, & ! inputs & brief, & ! inputs & continuum_LRi, & ! inputs & cooling_curvature, density_anomaly, & ! inputs & dQdtdA, eleDat, elev, & ! inputs & fault_LRi, fDip, & ! inputs & mxEl, mxFEl, mxNode, n1000, & ! inputs & nFakeN, nFl, nodeF, nodes, & ! inputs & nRealN, numEl, numNod, offset, orbData5, & ! inputs & title1, tLNode, xNode, yNode, zMNode) ! inputs ! Writes finite element grid to unit "iUnitO". IMPLICIT NONE INTEGER, INTENT(IN) :: iUnitO, mxEl, mxFEl, mxNode, n1000, & & nFakeN, nFl, nodeF, nodes, nRealN, numEl, numNod LOGICAL, INTENT(IN) :: brief, orbData5 REAL*8, INTENT(IN) :: cooling_curvature, density_anomaly, & & dQdtdA, eleDat, elev, fDip, & & offset, tLNode, xNode, yNode, zMNode INTEGER, INTENT(IN) :: continuum_LRi, fault_LRi ! to be DIMENSIONed below ... CHARACTER*80, INTENT(INOUT) :: title1 ! allowing for ADJUSTL() INTEGER :: i, iPrint, k, LRi, nP LOGICAL :: anyemu REAL*8 :: dips, pLat, pLon DIMENSION continuum_LRi(mxEl), & & cooling_curvature(mxNode), density_anomaly(mxNode), & & dQdtdA(mxNode), eleDat(3, mxEl), elev(mxNode), & & fault_LRi(mxFEl), fDip(2, mxFEl), & & nP(4), nodeF(4, mxFEl), nodes(3, mxEl), & & offset(mxFEl), tLNode(mxNode), & & xNode(mxNode), yNode(mxNode), zMNode(mxNode) DIMENSION dips(2) title1 = ADJUSTL(title1) WRITE (iUnitO, "(A)") TRIM(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 IF (orbData5) THEN WRITE (iUnitO, 91) i, pLon, pLat, elev(i), dQdtdA(i), & & zMNode(i), tLNode(i), & & density_anomaly(i), & & cooling_curvature(i) ELSE WRITE (iUnitO, 91) i, pLon, pLat, elev(i), dQdtdA(i), & & zMNode(i), tLNode(i) END IF 91 FORMAT (I8, 2F11.5, 6ES10.2) 100 CONTINUE WRITE (iUnitO, 110) numEl 110 FORMAT (I10,' (numEl = number of triangular Continuum elements)') anyemu = .FALSE. DO 130 i = 1, numEl DO 120 k = 1, 3 IF (eleDat(k, i) /= 0.0D0) THEN anyemu = .TRUE. GO TO 131 END IF 120 CONTINUE 130 CONTINUE 131 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 (anyemu) THEN IF (LRi == 0) THEN WRITE (iUnitO, 160) i, (np(k), k = 1, 3), (eleDat(k, i), k = 1, 3) 160 FORMAT (I8, 3I8, ES8.1, F7.1, ES8.1) ELSE WRITE (iUnitO, "(I8, 3I8, ES8.1, F7.1, ES8.1, ' LR', I8)") i, (np(k), k = 1, 3), (eleDat(k, i), k = 1, 3), LRi END IF ELSE ! anyemu == .FALSE. IF (LRi == 0) THEN WRITE (iUnitO, 170) i, (np(k), k = 1, 3) 170 FORMAT (I8, 3I8) ELSE WRITE (iUnitO, "(I8, 3I8, ' LR', I8)") i, (np(k), k = 1, 3), LRi END IF 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.01) 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 (iUnitO, "(I8, 4I8, 2F6.1, ES10.2, ' LR', I8)") i, (np(k), k = 1, 4), (dips(k), k = 1, 2), offset(i), LRi END IF 300 CONTINUE END SUBROUTINE PutNet SUBROUTINE MixUpX (iNew, mxNode, numNod, & ! inputs & data, & ! modify & atNode) ! work-space array ! Sort floating-point scalar values in array "data" so that the ! entry currently in position i ends up in position "iNew(i)". IMPLICIT NONE INTEGER, INTENT(IN) :: iNew, mxNode, numNod REAL*8, INTENT(INOUT) :: atNode, data INTEGER :: i DIMENSION atNode(mxNode), iNew(mxNode), & & data(mxNode) DO 10 i = 1, numNod atNode(i) = DATA(i) 10 CONTINUE DO 20 i = 1, numNod DATA(inew(i)) = atNode(i) 20 CONTINUE END SUBROUTINE MixUpX SUBROUTINE MixUpI (iNew, mxEl, mxNode, numEl, n34, & ! inputs & nodes) ! modify ! Change "n34"-entry INTEGER sets in array "nodes" so that the ! integer which is currently "i" is changed to "iNew(i)". IMPLICIT NONE INTEGER, INTENT(IN) :: iNew, mxEl, mxNode, numEl, n34 INTEGER, INTENT(INOUT) :: nodes INTEGER :: j, k DIMENSION inew(mxNode), & & nodes(n34, mxEl) DO 20 k = 1, n34 DO 10 j = 1, numEl nodes(k, j) = inew(nodes(k, j)) 10 CONTINUE 20 CONTINUE END SUBROUTINE MixUpI SUBROUTINE Next (i, j, mxEl, mxFEl, nFl, nodeF, nodes, numEl, & ! inputs & kFault, kEle) ! outputs ! Determine whether there are more elements adjacent to side j of ! triangular continuum element #i. ! j = 1 means the side opposite node # nodes(1,i). ! j = 2 means the side opposite node # nodes(2,i). ! j = 3 means the side opposite node # nodes(3,i). ! If a fault element is adjacent, its number is kFault; otherwise, ! kFault is set to zero. ! If another triangular continuum element is adjacent (even across ! fault element kFault!) then its number is kEle; otherwise, kEle = 0. IMPLICIT NONE INTEGER, INTENT(IN) :: i, j, mxEl, mxFEl, nFl, nodeF, nodes, numEl INTEGER, INTENT(OUT) :: kFault, kEle INTEGER :: k, kHigh, kLow, l, m1, m2, m3, m4, n1, n2 LOGICAL :: foundf DIMENSION nodeF(4, mxFEl), nodes(3, mxEl) ! Three node numbers along the side of interest, counterclockwise: n1 = nodes(MOD(j, 3) + 1, i) n2 = nodes(MOD(j + 1, 3) + 1, i) ! Check for adjacent fault element firsst foundf = .FALSE. kFault = 0 IF (nFl > 0) THEN DO 10 k = 1, nFl m1 = nodeF(1, k) m2 = nodeF(2, k) m3 = nodeF(3, k) m4 = nodeF(4, k) IF (((m1 == n2).AND.(m2 == n1)).OR. & & ((m3 == n2).AND.(m4 == n1))) THEN foundf = .TRUE. kFault = k GO TO 11 END IF 10 CONTINUE 11 IF (.NOT.foundf) kFault = 0 ! If there was a fault, replace 2 node numbers that we search for: IF (foundf) THEN IF (m2 == n1) THEN n1 = m3 n2 = m4 ELSE n1 = m1 n2 = m2 END IF END IF END IF ! Search for adjacent triangular continuum element: kEle = 0 kLow = i kHigh = i ! --- Begin irregular loop, to search out nearest elements first --- 100 kLow = kLow - 1 IF (kLow >= 1) THEN DO 110 l = 1, 3 m1 = nodes(MOD(l, 3) + 1, kLow) m2 = nodes(MOD(l + 1, 3) + 1, kLow) IF ((m2 == n1).AND.(m1 == n2)) THEN kEle = kLow RETURN END IF 110 CONTINUE END IF kHigh = kHigh + 1 IF (kHigh <= numEl) THEN DO 120 l = 1, 3 m1 = nodes(MOD(l, 3) + 1, kHigh) m2 = nodes(MOD(l + 1, 3) + 1, kHigh) IF ((m2 == n1).AND.(m1 == n2)) THEN kEle = kHigh RETURN END IF 120 CONTINUE END IF IF ((kLow > 1).OR.(kHigh < numEl)) GO TO 100 END SUBROUTINE Next REAL*8 FUNCTION ATan2F(y, x) ! Corrects for problem of two zero arguments IMPLICIT NONE REAL*8, INTENT(IN) :: y, x IF ((y /= 0.0D0).OR.(x /= 0.0D0)) THEN ATan2F = ATAN2(y, x) ELSE ATan2F = 0.0D0 END IF END FUNCTION ATan2F SUBROUTINE Square (brief, fDip, iUnit8, & ! inputs & mxBn, mxEl, mxFEl, mxNode, & ! inputs & mxStar, nFl, nodeF, nodes, & ! inputs & numEl, numNod, & ! inputs & xNode, yNode, & ! modify & edgeFs, & ! outputs & edgeTs, & ! outputs & nCond, nodCon, & ! outputs & checkN, list) ! work-space arrays ! Check, correct, and complete the geometry of the grid. IMPLICIT NONE INTEGER, INTENT(IN) :: iUnit8, mxBn, mxEl, mxFEl, mxNode, & & mxStar, nFl, nodeF, nodes, numEl, numNod LOGICAL, INTENT(IN) :: brief REAL*8, INTENT(IN) :: fDip REAL*8, INTENT(INOUT) :: xNode, yNode LOGICAL*1, INTENT(OUT) :: edgeFs, edgeTs INTEGER, INTENT(OUT) :: nCond, nodCon, list LOGICAL*1 :: checkN INTEGER :: i, j, j1, j2, k, kEle, kFault, & & n, nInSum, n1, n2, nDone, nGood, nj1, nj2, nl1, nl2, nl3, nl4, nLeft, node LOGICAL :: agreed, allOK REAL*8 :: c1, c2, c3, corner, dot, dx, dy, & & r, r2, rMax, short, t2, test, toler, & & x, xMean, xSum, y, yMean, ySum DIMENSION checkN(mxNode), & & edgeFs(2, mxFEl), edgeTs(3, mxEl), fDip(2, mxFEl), & & list(mxStar), nodCon(mxBn), & & nodeF(4, mxFEl), nodes(3, mxEl), & & xNode(mxNode), yNode(mxNode) DIMENSION corner(3, 3) ! (1) Check that all nodes are connected to at least one ! continuum (triangular) element or fault element: WRITE (iUnit8, 1) 1 FORMAT (' In Square: Are all nodes connected?') DO 110 i = 1, numNod checkN(i) = .FALSE. 110 CONTINUE DO 130 i = 1, numEl DO 120 j = 1, 3 checkN(nodes(j, i)) = .TRUE. 120 CONTINUE 130 CONTINUE DO 136 i = 1, nFl DO 134 j = 1, 4 checkN(nodeF(j, i)) = .TRUE. 134 CONTINUE 136 CONTINUE allOK = .TRUE. DO 140 i = 1, numNod allOK = allOK.AND.checkN(i) 140 CONTINUE IF (.NOT.allOK) THEN WRITE(iUnit8, 150) 150 FORMAT(' BAD GRID TOPOLOGY: FOLLOWING REAL NODES DO NOT'/ & & ' BELONG TO ANY TRIANGULAR CONTINUUM ELEMENT'/ & & ' OR FAULT ELEMENT:') DO 160 i = 1, numNod IF (.NOT.checkN(i)) WRITE (iUnit8, 155) i 155 FORMAT (' ',43X,I6) 160 CONTINUE CALL PAUSE() STOP END IF ! (2) Average together the coordinates of all nodes at one "point": WRITE (iUnit8, 2) 2 FORMAT (' In Square: Averaging node coordinates') DO 410 i = 1, numNod checkN(i) = .FALSE. ! (Means "not yet involved in averaging') 410 CONTINUE DO 490 i = 1, nFl DO 480 j1 = 1, 2 nj1 = nodeF(j1, i) ! (Fault ends are the only places that can have problems.) IF (.NOT.checkN(nj1)) THEN list(1) = nj1 checkN(nj1) = .TRUE. ! Begin list of neighbors with paired node: j2 = 5 - j1 nj2 = nodeF(j2, i) list(2) = nj2 checkN(nj2) = .TRUE. nInSum = 2 ! Find shortest fault connected to either one: dx = xNode(nj1) - xNode(nj2) dy = yNode(nj1) - yNode(nj2) IF (dy > 3.14159D0) dy = dy - 6.28318D0 IF (dy < -3.14159D0) dy = dy + 6.28318D0 dy = dy * SIN(xNode(nj1)) short = SQRT(dx**2 + dy**2) DO 470 k = 1, nFl nl1 = nodeF(1, k) nl2 = nodeF(2, k) nl3 = nodeF(3, k) nl4 = nodeF(4, k) IF ((nj1 == nl1).OR.(nj2 == nl1).OR. & & (nj1 == nl2).OR.(nj2 == nl2).OR. & & (nj1 == nl3).OR.(nj2 == nl3).OR. & & (nj1 == nl4).OR.(nj2 == nl4)) THEN dx = xNode(nl1) - xNode(nl2) dy = yNode(nl1) - yNode(nl2) IF (dy > 3.14159D0) dy = dy - 6.28318D0 IF (dy < -3.14159D0) dy = dy + 6.28318D0 dy = dy * SIN(xNode(nl1)) test = SQRT(dx**2 + dy**2) short = MIN(short, test) END IF 470 CONTINUE ! Collect all nodes within 10% of this distance: toler = short / 10.0D0 t2 = toler**2 DO 471 k = 1, numNod IF (.NOT.checkN(k)) THEN dx = xNode(nj1) - xNode(k) dy = yNode(nj1) - yNode(k) IF (dy > 3.14159D0) dy = dy - 6.28318D0 IF (dy < -3.14159D0) dy = dy + 6.28318D0 dy = dy * SIN(xNode(nj1)) r2 = dx**2 + dy**2 IF (r2 < t2) THEN nInSum = nInSum + 1 IF (nInSum > mxStar) THEN WRITE(iUnit8, 421) 421 FORMAT(/' INCREASE VALUE' & & ,' OF PARAMETER MAXATP.') CALL PAUSE() STOP END IF list(nInSum) = k checkN(k) = .TRUE. END IF END IF 471 CONTINUE ! (Quick exit if all nodes in same place) agreed = .TRUE. DO 472 k = 2, nInSum agreed = agreed.AND. & & (xNode(list(k)) == xNode(list(1))).AND. & & (yNode(list(k)) == yNode(list(1))) 472 CONTINUE IF (agreed) GO TO 480 xSum = 0.0D0 ySum = 0.0D0 DO 473 k = 1, nInSum xSum = xSum + xNode(list(k)) ySum = ySum + yNode(list(k)) 473 CONTINUE xMean = xSum / nInSum yMean = ySum / nInSum rMax = 0.0D0 DO 474 k = 1, nInSum r = SQRT((xNode(list(k)) - xMean)**2 + & & (yNode(list(k)) - yMean)**2) rMax = MAX(rMax, r) 474 CONTINUE DO 475 k = 1, nInSum xNode(list(k)) = xMean yNode(list(k)) = yMean 475 CONTINUE IF (.NOT.brief) THEN IF (rMax > 0.) THEN WRITE(iUnit8, 476) nInSum, & & (list(n), n = 1, nInSum) 476 FORMAT(/ & & ' AVERAGING TOGETHER THE POSITIONS OF', & & ' THESE ',I6,' NODES:',(/' ',12I6)) WRITE (iUnit8, 477) rMax 477 FORMAT (' MAXIMUM CORRECTION TO ', & & 'ANY POSITION IS',1P,E10.2/ & & ' YOU ARE RESPONSIBLE FOR ', & & ' DECIDING WHETHER THIS IS A', & & ' SERIOUS ERR0R!') END IF END IF END IF 480 CONTINUE 490 CONTINUE ! (3) Check for negative areas: WRITE (iUnit8, 3) 3 FORMAT (' In Square: Checking for negative areas') DO 600 i = 1, numEl DO 510 j = 1, 3 n = nodes(j, i) corner(1, j) = COS(yNode(n)) * SIN(xNode(n)) corner(2, j) = SIN(yNode(n)) * SIN(xNode(n)) corner(3, j) = COS(xNode(n)) 510 CONTINUE c1 = corner(2, 1) * corner(3, 2) - corner(3, 1) * corner(2, 2) c2 = corner(3, 1) * corner(1, 2) - corner(1, 1) * corner(3, 2) c3 = corner(1, 1) * corner(2, 2) - corner(2, 1) * corner(1, 2) dot = c1 * corner(1, 3) + c2 * corner(2, 3) + c3 * corner(3, 3) IF (dot <= 0.0D0) THEN WRITE (iUnit8, 509) i 509 FORMAT (/ /' ERR0R: AREA OF ELEMENT ',I5, & & ' IS NOT POSITIVE.') CALL PAUSE() STOP END IF 600 CONTINUE ! (5) Make a list of nodes that are on the boundary and require ! boundary conditions (nodCon); these are in countercloskwise ! order. Also make a list of element sides which contain these ! nodes: edgeTs AND edgeFs. WRITE (iUnit8, 5) 5 FORMAT (' In Square: Creating list of boundary nodes') nCond = 0 DO 801 i = 1, numNod checkN(i) = .FALSE. 801 CONTINUE IF (nFl > 0) THEN DO 802 i = 1, nFl edgeFs(1, i) = .FALSE. edgeFs(2, i) = .FALSE. 802 CONTINUE END IF DO 810 i = 1, numEl DO 809 j = 1, 3 CALL Next (i, j, mxEl, mxFEl, nFl, nodeF, nodes, numEl, & ! inputs & kFault, kEle) ! outputs IF (kEle > 0) THEN ! (ordinary interior side) edgeTs(j, i) = .FALSE. ELSE IF (kFault == 0) THEN ! (exterior side) edgeTs(j, i) = .TRUE. n1 = nodes(MOD(j, 3) + 1, i) n2 = nodes(MOD(j + 1, 3) + 1, i) IF (.NOT.checkN(n1)) THEN nCond = nCond + 1 checkN(n1) = .TRUE. END IF IF (.NOT.checkN(n2)) THEN nCond = nCond + 1 checkN(n2) = .TRUE. END IF ELSE ! (triangular element has an exterior fault element adjacent to it) edgeTs(j, i) = .FALSE. n1 = nodes(MOD(j, 3) + 1, i) IF (nodeF(2, kFault) == n1) THEN edgeFs(2, kFault) = .TRUE. DO 806 k = 3, 4 n = nodeF(k, kFault) IF (.NOT.checkN(n)) THEN nCond = nCond + 1 checkN(n) = .TRUE. END IF 806 CONTINUE ELSE edgeFs(1, kFault) = .TRUE. DO 808 k = 1, 2 n = nodeF(k, kFault) IF (.NOT.checkN(n)) THEN nCond = nCond + 1 checkN(n) = .TRUE. END IF 808 CONTINUE END IF END IF 809 CONTINUE 810 CONTINUE IF (nCond > mxBn) THEN WRITE(iUnit8, 820) nCond 820 FORMAT(/' INCREASE PARAMETER MAXBN TO',I6,' AND RECOMPILE.') CALL PAUSE() STOP END IF ! Stop work if no boundary nodes found (global grid): IF (nCond == 0) GO TO 899 ! Begin circuit with lowest-numbered boundary node WRITE (iUnit8, 825) 825 FORMAT (' In Square: Ordering list of boundary nodes') DO 830 i = 1, numNod IF (checkN(i)) GO TO 831 830 CONTINUE 831 nodCon(1) = i nDone = 1 nLeft = nCond ! Beginning of indefinate loop which traces around the perimeter. ! Each time, it progresses by one of 3 steps: ! -2 nodes at a time along a triangle side, or ! -2 nodes at a time along a fault element side, or ! -by finding another (corner) node which shares the same location. ! ! Beginning of main indefinate loop: 840 node = nodCon(nDone) ! IMPORTANT: Check that we are not revisiting a node! ! This would mean that there are too many boundary nodes ! to fit in the simply-connected loop, and that there ! are excess boundary nodes somewhere unconnected! IF (.NOT.checkN(node)) THEN nGood = nDone - 2 WRITE (iUnit8, 841) nGood, nCond 841 FORMAT(/' ERROR IN GRID, reported by -SQUARE-:' & & /' BOUNDARY IS NOT SIMPLY-CONNECTED.' & & /' Closed loop of ',I6,' nodes does not' & & /' include all ',I6,' boundary nodes.' & & /' Run command Perimeter in OrbWeaver' & & /' for a map of the bad nodes.'/) CALL PAUSE() STOP END IF IF (nDone > 1) checkN(node) = .FALSE. x = xNode(node) y = yNode(node) ! Look for an adjacent triangular element using this node. DO 844 i = 1, numEl DO 842 j = 1, 3 IF (edgeTs(j, i)) THEN n1 = nodes(MOD(j, 3) + 1, i) IF (n1 == node) GO TO 846 END IF 842 CONTINUE 844 CONTINUE GO TO 850 846 n2 = nodes(MOD(j + 1, 3) + 1, i) nDone = nDone + 1 IF (nDone <= nCond) nodCon(nDone) = n2 nLeft = nLeft - 1 IF (nLeft > 0) THEN GO TO 840 ELSE GO TO 870 END IF ! Else, look for an adjacent fault element using this node: 850 DO 854 i = 1, nFl IF (edgeFs(1, i)) THEN IF (nodeF(1, i) == node) THEN n2 = nodeF(2, i) GO TO 856 END IF ELSE IF (edgeFs(2, i)) THEN IF (nodeF(3, i) == node) THEN n2 = nodeF(4, i) GO TO 856 END IF END IF 854 CONTINUE GO TO 860 856 nDone = nDone + 1 IF (nDone <= nCond) nodCon(nDone) = n2 nLeft = nLeft - 1 IF (nLeft > 0) THEN GO TO 840 ELSE GO TO 870 END IF ! Else, look for another exterior corner node at same location: 860 DO 865 i = 1, numNod IF ((i /= node).AND.checkN(i)) THEN IF ( (ABS(xNode(i) - x) < 1.E-6) .AND. & & (ABS(yNode(i) - y) < 1.E-6) ) GO TO 867 END IF 865 CONTINUE WRITE(iUnit8, 866) node 866 FORMAT(' BAD GRID TOPOLOGY: WHILE TRACING PERIMETER,'/ & & ' COULD NOT FIND ANY WAY TO CONTINUE FROM NODE ',I6/ & & ' EITHER THROUGH SHARED BOUNDARY ELEMENTS, OR'/ & & ' THROUGH OTHER BOUNDARY NODES SHARING THE SAME ', & & 'POSITION.') CALL PAUSE() STOP 867 nDone = nDone + 1 IF (nDone <= nCond) nodCon(nDone) = i nLeft = nLeft - 1 IF (nLeft > 0) GO TO 840 ! ! End of indefinate loop which traces around perimeter. ! 870 IF (.NOT.brief) THEN WRITE(iUnit8, 880) 880 FORMAT(/ /' Here follows a list, in consecutive order,'/ & & ' of the nodes which define the perimeter'/ & & ' of the model; these nodes require boundary', & & ' conditions:'/' BC# NODE') DO 890 i = 1, nCond n = nodCon(i) WRITE(iUnit8, 882) i, n 882 FORMAT(' ',2I6) 890 CONTINUE n = nodCon(1) WRITE (iUnit8, 892) n 892 FORMAT(' (Note: Node ',I6,' completes the loop, but will', & & ' not be listed twice.)') END IF 899 CONTINUE IF (.NOT. brief) WRITE (iUnit8, 9999) 9999 FORMAT (' --------------------------------------------------', & & '-----------------------------') END SUBROUTINE Square SUBROUTINE ReadN (iUnitP, iUnitT, n, & ! inputs & vector) ! outputs ! A utility routine designed to permit -Faults- input files ! to also be used by -Plates-, which expects more numbers ! in some records. ! Also supports reading either OrbData files (with 4 nodal ! variables) or OrbData5 files (with 6 nodal variables). ! 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. IMPLICIT NONE INTEGER, INTENT(IN) :: iUnitP, iUnitT, n REAL*8, INTENT(OUT) :: vector CHARACTER*1 :: c CHARACTER*132 :: line INTEGER :: i, ios, numba LOGICAL :: anyIn, dotted, expon, signed DIMENSION vector(n) line = ' ' READ (iUnitP, 1, iostat = ios) line 1 FORMAT (A) numba = 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 numba = numba + 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) numba = numba + 1 50 IF (numba == 0) THEN WRITE (iUnitT, 91) n, TRIM(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 IF (numba >= n) THEN READ (line, * ) (vector(i), i = 1, n) ELSE READ (line, * ) (vector(i), i = 1, numba) DO 99 i = numba + 1, n vector(i) = 0.0D0 99 CONTINUE END IF END IF END SUBROUTINE ReadN SUBROUTINE ReadNFromString (passed_line, iUnitT, n, & ! inputs & vector) ! outputs ! Just like subprogram ReadN (above), except that it takes ! its input from a one-line text-string, instead of from ! a Fortran input-device designator. ! A utility routine designed to permit -Faults- input files ! to also be used by -Plates-, which expects more numbers ! in some records. ! Also supports reading either OrbData files (with 4 nodal ! variables) or OrbData5 files (with 6 nodal variables). ! 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. IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: passed_line INTEGER, INTENT(IN) :: iUnitT, n REAL*8, INTENT(OUT) :: vector CHARACTER*1 :: c CHARACTER*132 :: line INTEGER :: i, ios, numba LOGICAL :: anyIn, dotted, expon, signed DIMENSION vector(n) line = passed_line ! copying, just to be sure we don't change it. numba = 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 numba = numba + 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) numba = numba + 1 50 IF (numba == 0) THEN WRITE (iUnitT, 91) n, TRIM(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 IF (numba >= n) THEN READ (line, * ) (vector(i), i = 1, n) ELSE READ (line, * ) (vector(i), i = 1, numba) DO 99 i = numba + 1, n vector(i) = 0.0D0 99 CONTINUE END IF END IF END SUBROUTINE ReadNFromString SUBROUTINE Number (iUnitT, length, list, & ! inputs & mxEl, mxFEl, mxGSiz, mxNode, & ! inputs & nFl, nodeF, nodes, nRealN, & ! inputs & numEl, numNod, sphere, xNode, yNode, & ! inputs & tInSec, & ! modify & iAlias, iUser, maxDif, minDif, & ! outputs & doNext, iDiag, lGraph, maxLst) ! work-space arrays ! Find a general renumbering of the nodes of the finite-element ! grid which is (nearly) optimal in terms of minimizing the ! bandwidth of the connectivity matrix (a surrogate ! for the stiffness matrix, but symmetrical and with only one ! row and one column per node). ! In the section which fills the connectivity matrix ! "lGraph", it is assumed that the grid is made up of ! a set of "numEl" 3-node triangular elements (whose old ! node numbers are stored as "nodes(k=1..3,i=1..numEl)"), and ! "nFl" 4-node fault elements (whose old node numbers are ! stored as "nodeF(k=1..4,i=1..nFl)" ). ! nFl may be zero. ! This section can be easily modified for other types of elements. ! It does not matter whether or not the numbers defining ! each element occur in any particular order, or not. ! Degenerate elements with two or more identical nodes numbers ! (e.g., at the end of a fault trace) are also permitted. ! The total number of nodes in the grid is input parameter ! "numNod". The first "nRealN" of these will be renumbered. ! Usually, numNod will equal nRealN. Obviously, nRealN cannot ! be greater than numNod. The distinction is made in case ! you are using a finite element program like "Faults" in ! which fake nodes (those with no degrees of freedom) are ! listed last (with numbers from nRealN+1 to numNod), and ! you want to keep this group at the end, and not renumber ! them. ! The algorithm used is an original variant of the ! Cuthill-McKee algorithm ! described by H. R. Sxhwarz in "Finite Element Methods", ! "Computational Mathematics and Applications" series, ! Academic Press/Harcourt-Brace-Jovanovitch, 1988, pp.167-185. ! Instead of using degree to decide which nodes will be numbered ! first, this variant numbers nodes within one level on the basis ! of their azimuth from the initial node. ! Because of this choice of algorithm, this code CAN NOT be expected ! to work for grids which are not simply-connected (i.e., ! those with holes in them), and it cannot be expected to work ! well for grids which have extremely concave boundaries. ! (Fortunately, these cases rarely arise in Earth sciences, ! where the side boundaries of the grid are always artificial, and ! are almost always chosen as a convex polyhedron.) ! The products of this subprogram are: ! *"iAlias" = an INTEGER vector of new node numbers, arranged ! in the order of the old node numbers. ! *"iUser" = an INTEGER vector of old node numbers, arranged ! in the order of the new node numbers. ! "iAlias" and "iUser" are mutual inverses, as shown by ! the following theoroms: ! iUser(iAlias(i)) == i ! iAlias(iUser(j)) == j ! * maxDif = half-bandwidth before renumbering. ! * minDif = half-bandwidth after renumbering. ! (The definition of half-bandwidth used here ! does not include the diagonal, and can be defined as the ! maximum of the absolute value of the difference of the ! (new) numbers in all possible pairs of connected nodes.) IMPLICIT NONE INTEGER, INTENT(IN) :: iUnitT, length, list, & & mxEl, mxFEl, mxGSiz, mxNode, & & nFl, nodeF, nodes, nRealN, & & numEl, numNod LOGICAL, INTENT(IN) :: sphere REAL*8, INTENT(IN) :: xNode, yNode REAL*8, INTENT(INOUT) :: tInSec INTEGER, INTENT(OUT) :: iAlias, iUser, maxDif, minDif ! Note: The following TYPE declaration is very important for ! economizing storage, as lGraph will have a size of ! approximately nRealN*(nRealN-2)/2 LOGICAL entries. LOGICAL*1 :: doNext, lGraph INTEGER :: iDiag, maxLst INTEGER :: i, iAlDif, indexN, iRow, iStart, j, jCol, k, last1, last9, level, & & m, m1, m9, minMax, mLess, mMore, mStart, & & n, nDone, nInRow, node, nToDo, nUsed, speed_option LOGICAL :: found, lessOK, moreOK, recompute REAL*8 :: a, a1, aLess, aMore, arc, arcLim, arg, argCut, argLes, argMin, argMor, ATan2F, ax, ay, & & b, b1, bLess, bMore, bx, by, cross, dArc, dot, & & g, g1, gLess, gMore, gx, gy, & & phi, phi1, pLess, pMore, size, sterad, theta, theta1, tLess, tMore, & & x, xc, xLess, xm, xMore, y, yc, yLess, ym, yMore DIMENSION doNext(mxNode), iAlias(mxNode), iDiag(mxNode), & & iUser(mxNode), lGraph(mxGSiz), list(mxNode), & & maxLst(mxNode), & & nodeF(4, mxFEl), nodes(3, mxEl), & & xNode(mxNode), yNode(mxNode) ! STATEMENT FUNCTION, giving storage location (vector index) of the ! element of the graph matrix at row "iRow" and column "jCol": indexN(iRow, jCol) = iDiag(MIN(iRow, jCol)) + ABS(jCol - iRow) IF (nRealN > numNod) THEN WRITE (iUnitT, 5) nRealN, numNod 5 FORMAT (/ /' SUBPROGRAM NUMBER WAS ASKED TO RENUMBER', & & ' THE FIRST ',I6,' NODES,'/ & & ' OUT OF A TOTAL OF ',I6,'.'/ & & ' THIS IS IMPOSSIBLE.'/ & & ' Input parameter nRealN must not exceed', & & ' input parameter numNod.') CALL PAUSE() STOP END IF ! Compute the locations of the diagonal elements of ! the connectivity matrix "lGraph". ! (Note that the graph is symmetric, so we only store the ! upper triangle, and we access the matrix ! indirectly using STATEMENT FUNCTION "indexN".) nUsed = 0 DO 10 i = 1, nRealN iDiag(i) = nUsed nInRow = nRealN - i nUsed = nUsed + nInRow 10 CONTINUE IF (nUsed > mxGSiz) THEN WRITE (iUnitT, 11) mxGSiz, nUsed 11 FORMAT (/ /' INCREASE PARAMETER MAXGSI ', & & '(CURRENTLY ',I12,') TO AT LEAST ',I12/ & & ' AND RECOMPILE THIS PROGRAM.') CALL PAUSE() STOP END IF ! Create the graph matrix: WRITE (iUnitT, 12) 12 FORMAT (/' Creating the graph matrix...') DO 20 i = 1, mxGSiz lGraph(i) = .FALSE. 20 CONTINUE DO 50 i = 1, numEl DO 40 j = 1, 2 DO 30 k = j + 1, 3 iRow = nodes(j, i) jCol = nodes(k, i) IF ((iRow < 1).OR.(iRow > numNod)) THEN WRITE (iUnitT, 21) iRow, i 21 FORMAT (/ /' ERR0R IN INPUT DATA:'/ & & ' NODE NUMBER ',I5,' IN DEFINITION OF', & & ' TRIANGULAR ELEMENT ',I5,' IS ILLEGAL.') CALL PAUSE() STOP END IF IF ((jCol < 1).OR.(jCol > numNod)) THEN WRITE (iUnitT, 21) jCol, i CALL PAUSE() STOP END IF IF ((iRow <= nRealN).AND.(jCol <= nRealN) & & .AND.(iRow /= jCol)) & & lGraph(indexN(iRow, jCol)) = .TRUE. 30 CONTINUE 40 CONTINUE 50 CONTINUE DO 150 i = 1, nFl DO 140 j = 1, 3 DO 130 k = j + 1, 4 iRow = nodeF(j, i) jCol = nodeF(k, i) IF ((iRow < 1).OR.(iRow > numNod)) THEN WRITE (iUnitT, 51) iRow, i 51 FORMAT (/ /' ERR0R IN INPUT DATA:'/ & & ' NODE NUMBER ',I5,' IN DEFINITION OF', & & ' FAULT ELEMENT ',I5,' IS ILLEGAL.') CALL PAUSE() STOP END IF IF ((jCol < 1).OR.(jCol > numNod)) THEN WRITE (iUnitT, 51) jCol, i CALL PAUSE() STOP END IF IF ((iRow <= nRealN).AND.(jCol <= nRealN) & & .AND.(iRow /= jCol)) & & lGraph(indexN(iRow, jCol)) = .TRUE. 130 CONTINUE 140 CONTINUE 150 CONTINUE WRITE (iUnitT, 151) 151 FORMAT (' ....done.') ! Determine the original bandwidth: WRITE (iUnitT, 152) 152 FORMAT (/' Determining the original bandwidth...') maxDif = 0 DO 154 i = 1, nRealN - 1 DO 153 j = i + 1, nRealN IF (lGraph(indexN(i, j))) maxDif = MAX(maxDif, ABS(j - i)) 153 CONTINUE 154 CONTINUE WRITE (iUnitT, 151) !=================================================================== ! CHECK whether the user wants a quick-and-dirty renumbering, ! or the most thorough, reliable renumbering possible !(at the cost of much greater execution time)? 199 WRITE (*, *) WRITE (*, "(' ----------------------------------------------------')") WRITE (*, "(' CHOOSE DESIRED RENUMBERING METHOD:')") WRITE (*, "(' 1: QUICK method, assuming node #1 is well-chosen.')") WRITE (*, "(' 2: COMPLETE method, testing all boundary nodes.')") WRITE (*, "(' ----------------------------------------------------')") WRITE (*, "(' ENTER integer to choose: '\)") READ (*, *) speed_option IF ((speed_option < 1).OR.(speed_option > 2)) THEN WRITE (*, *) WRITE (*, "(' Your answer was not understood. Choose 1 or 2.')") WRITE (*, *) GO TO 199 END IF !=================================================================== IF (speed_option == 1) THEN m1 = 1 m9 = 1 WRITE (*, "(/' THERE WILL BE ONLY 1 RENUMBERING ATTEMPT:')") ELSE ! speed_option == 2 m1 = 1 m9 = length WRITE (*, "(/' THERE WILL BE ',I6,' RENUMBERING ATTEMPTS:')") length END IF recompute = .FALSE. 200 DO 10000 m = m1, m9 ! Create a renumbering scheme, using initial node "list(m)" ! (Necessary initialization: Reset numbering scheme to null:) DO 201 i = 1, nRealN iAlias(i) = 0 iUser(i) = 0 201 CONTINUE ! (Begin real work) ! Locate initial node: theta1 = xNode(list(m)) phi1 = yNode(list(m)) a1 = COS(phi1) * SIN(theta1) b1 = SIN(phi1) * SIN(theta1) g1 = COS(theta1) ! Prime meridian vector (x) is perpendicular to initial node: IF (SIN(theta1) > 0) THEN ax = -b1 bx = a1 gx = 0.0D0 size = SQRT(ax * ax + bx * bx) ax = ax / size bx = bx / size ELSE ax = 1.0D0 bx = 0.0D0 gx = 0.0D0 END IF ! Third orthogonal vector (y) ay = b1 * gx - g1 * bx by = g1 * ax - a1 * gx gy = a1 * bx - b1 * ax size = SQRT(ay**2 + by**2 + gy**2) ay = ay / size by = by / size gy = gy / size ! Determine argument of cut in azimuth from initial node: IF (sphere) THEN ! Location of cut is arbitrary! argCut = 0.0D0 ELSE IF (m < length) THEN mMore = m + 1 ELSE mMore = 1 END IF IF (m > 1) THEN mLess = m - 1 ELSE mLess = length END IF tMore = xNode(list(mMore)) pMore = yNode(list(mMore)) aMore = COS(pMore) * SIN(tMore) bMore = SIN(pMore) * SIN(tMore) gMore = COS(tMore) tLess = xNode(list(mLess)) pLess = yNode(list(mLess)) aLess = COS(pLess) * SIN(tLess) bLess = SIN(pLess) * SIN(tLess) gLess = COS(tLess) ! "x" and "y" are in equatorial plane relative to ! the current initial node: xMore = aMore * ax + bMore * bx + gMore * gx yMore = aMore * ay + bMore * by + gMore * gy xLess = aLess * ax + bLess * bx + gLess * gx yLess = aLess * ay + bLess * by + gLess * gy moreOK = (ABS(xMore) + ABS(yMore)) > 0.0D0 lessOK = (ABS(xLess) + ABS(yLess)) > 0.0D0 IF (moreOK.AND.lessOK) THEN argMor = ATan2F(yMore, xMore) argLes = ATan2F(yLess, xLess) ELSE IF (moreOK) THEN argMor = ATan2F(yMore, xMore) argLes = argMor - 3.14159D0 IF (argLes < -3.14159D0) argLes = argLes + 6.28318D0 ELSE IF (lessOK) THEN argLes = ATan2F(yLess, xLess) argMor = argLes + 3.14159D0 IF (argMor > 3.14159D0) argMor = argMor - 6.28318D0 ELSE WRITE (iUnitT, 203) list(m) 203 FORMAT (/ /' PROBLEM IN SUBPROGRAM NUMBER:'/ & & ' BOUNDARY AZIMUTH IS UNDEFINED AT NODE ',I5/ & & ' BECAUSE BOTH NEIGHBORING BOUNDARY NODES ARE LOCATED' & & /' AT EXACTLY THE SAME POINT.') CALL PAUSE() STOP END IF ! ******* The bottom line (argCut) ******* argCut = 0.5 * (argMor + argLes) xm = COS(argMor) ym = SIN(argMor) xc = COS(argCut) yc = SIN(argCut) cross = xc * ym - yc * xm IF (cross < 0.0D0) THEN argCut = argCut + 3.14159D0 IF (argCut > 3.14159D0) argCut = argCut - 6.28318D0 END IF ! *************************************** END IF ! Size of bands to be used as levels in spherical numbering: sterad = 4.0D0 * 3.14159D0 / numEl dArc = 0.30D0 * SQRT(sterad) ! Number the first node (or a nearby real substitute): IF (list(m) <= nRealN) THEN iStart = list(m) ELSE iStart = 1 found = .FALSE. DO 215 i = 1, numEl DO 214 k = 1, 3 IF (nodes(k, i) == list(m)) THEN DO 213 j = 1, 3 IF (nodes(j, i) <= nRealN) THEN found = .TRUE. iStart = nodes(j, i) END IF 213 CONTINUE END IF 214 CONTINUE 215 CONTINUE IF (.NOT.found) THEN DO 225 i = 1, nFl DO 224 k = 1, 4 IF (nodeF(k, i) == list(m)) THEN DO 223 j = 1, 4 IF (nodeF(j, i) <= nRealN) THEN found = .TRUE. iStart = nodeF(j, i) END IF 223 CONTINUE END IF 224 CONTINUE 225 CONTINUE END IF END IF iAlias(iStart) = 1 iUser(1) = iStart nDone = 1 last1 = 1 last9 = 1 level = 0 arcLim = 0.0D0 !----------------- Begin indefinate loop on level -------------------- ! Mark the nodes to be numbered in the next level (unnumbered nodes ! connected to nodes that have already been numbered): 240 level = level + 1 arcLim = arcLim + dArc DO 250 i = 1, nRealN ! (Necessary initialization): doNext(i) = .FALSE. 250 CONTINUE IF (sphere) THEN ! Level is defined by distance from initial node: DO 252 i = 1, numNod IF (iAlias(i) == 0) THEN theta = xNode(i) phi = yNode(i) a = COS(phi) * SIN(theta) b = SIN(phi) * SIN(theta) g = COS(theta) dot = a * a1 + b * b1 + g * g1 dot = MIN(1., MAX(-1., dot)) arc = ACOS(dot) doNext(i) = (arc <= arcLim) END IF 252 CONTINUE ELSE ! Level is defined by connectivity: DO 300 i = last1, last9 node = iUser(i) IF (node > 1) THEN DO 260 j = 1, (node - 1) doNext(j) = doNext(j).OR. & & ( (iAlias(j) == 0) .AND. lGraph(indexN(j, node)) ) 260 CONTINUE END IF IF (node < nRealN) THEN DO 270 j = node + 1, nRealN doNext(j) = doNext(j).OR. & & ( (iAlias(j) == 0) .AND. lGraph(indexN(j, node)) ) 270 CONTINUE END IF 300 CONTINUE END IF ! How many are to be numbered at this level? nToDo = 0 DO 305 i = 1, nRealN IF (doNext(i)) nToDo = nToDo + 1 305 CONTINUE ! WRITE (IUNITT,308) nToDo,level ! 308 FORMAT (' ',I6,' nodes in level ',I4) IF (nToDo == 0) THEN IF (nDone == nRealN) THEN GO TO 8000 ELSE IF (.NOT.sphere) THEN WRITE (iUnitT, 311) iStart, nDone, iUser(nDone), nRealN 311 FORMAT (/ /' ERR0R IN GRID TOPOLOGY:'/ & & ' AFTER BEGINNING NUMBERING AT NODE ',I5/ & & ' AND NUMBERING ',I5,' NODES UP THROUGH ', & & 'NODE ',I5,','/' NO CONNECTED NEIGHBORING', & & ' NODES COULD BE FOUND,'/' SO IT IS NOT', & & ' POSSIBLE TO NUMBER ',I5,' NODES.') CALL PAUSE() STOP END IF END IF !$$$ ! Loop on nodes in list, taking them in order of their argument ! angles from the initial node: DO 400 n = 1, nToDo ! Find lowest-argument node left to be done: argMin = 999.99D0 DO 350 i = 1, nRealN IF (doNext(i)) THEN theta = xNode(i) phi = yNode(i) a = COS(phi) * SIN(theta) b = SIN(phi) * SIN(theta) g = COS(theta) x = a * ax + b * bx + g * gx y = a * ay + b * by + g * gy arg = ATan2F(y, x) IF (arg < argCut) arg = arg + 6.28318 IF (arg < argMin) THEN argMin = arg node = i END IF END IF 350 CONTINUE ! NUMBER IT ! ! ! nDone = nDone + 1 iUser(nDone) = node iAlias(node) = nDone doNext(node) = .FALSE. 400 CONTINUE ! Summarize the last level done: last1 = last9 + 1 last9 = last1 + nToDo - 1 ! Loop to next level IF (nDone < nRealN) GO TO 240 !------------------- END INDEFINATE LOOP ON LEVEL ----------------- ! Determine the new half-bandwidth 8000 minDif = 0 DO 9020 i = 1, nRealN - 1 DO 9010 j = i + 1, nRealN IF (lGraph(indexN(i, j))) THEN iAlDif = ABS(iAlias(i) - iAlias(j)) minDif = MAX(minDif, iAlDif) IF (recompute.AND.(iAlDif == minMax)) THEN WRITE (iUnitT, 9005) i, j, iAlias(i), iAlias(j), & & minMax 9005 FORMAT (' Old nodes ',I7,' and ',I7,' have', & & ' new numbers ',I7,' and ',I7/ & & ' whose difference of ',I7,' defines', & & ' the new bandwidth.') END IF END IF 9010 CONTINUE 9020 CONTINUE maxLst(m) = minDif IF (iStart == list(m)) THEN WRITE (iUnitT, 9031) m, m9, list(m), minDif ELSE WRITE (iUnitT, 9032) m, m9, list(m), iStart, minDif END IF 9031 FORMAT(' Attempt number ',I5,'/',I5,' begins with node ',I7,';' & & /' When this node is #1, the half-bandwidth is ',I7,'.') 9032 FORMAT(' Attempt number ',I5,'/',I5,' begins near boundary node ',I7,',' & & /' at nearest real node to that, number ',I7,'.'/ & & ' When this node is #1, the half-bandwidth is ',I7,'.') 10000 CONTINUE !================================================================= ! Repeat the numbering from the best starting point: IF (m9 > m1) THEN minMax = nRealN DO 10001 m = 1, length IF (maxLst(m) < minMax) THEN minMax = maxLst(m) mStart = m END IF 10001 CONTINUE m1 = mStart m9 = mStart ! ****************** recompute = .TRUE. GO TO 200 ! ****************** END IF IF (numNod > nRealN) THEN DO 90000 i = nRealN + 1, numNod iAlias(i) = i iUser(i) = i 90000 CONTINUE END IF END SUBROUTINE Number SUBROUTINE Clock (iUnitT, & ! input & tInSec) ! modify ! Report(?) elapsed time between CALLs IMPLICIT NONE INTEGER, INTENT(IN) :: iUnitT REAL*8, INTENT(INOUT) :: tInSec INTEGER :: nCount, nPerS, nCMax REAL*8 :: elaps, tOld tOld = tInSec CALL System_clock(nCount, nPerS, nCMax) IF (nPerS > 0) THEN tInSec = nCount / nPerS ELSE tInSec = 0.0D0 END IF elaps = tInSec - tOld IF (elaps <= 604800) THEN ! one week; if more, this is initialization CALL WRITE (iUnitT, 1) elaps 1 FORMAT (/' ',F10.2,' seconds have elapsed since last check.') END IF END SUBROUTINE Clock SUBROUTINE Pause() IMPLICIT none WRITE(*, "(' Press [Enter]...')") READ(*, * ) END SUBROUTINE Pause