PROGRAM OrbNumber ! (Mature version of ! 13 February 2010: supports OrbData5 format, ! and also provides a list of boundary nodes. ! Note that elapsed-time report now covers ONLY time spent ! within SUBROUTINE Number.) ! by Peter Bird ! Department of Earth and Space Sciences ! University of California ! Los Angeles, CA 90095-1567 ! pbird@ess.ucla.edu ! (C) Copyright 1993, 1994, 1996, 1998, 1999, 2005, 2010, 2013 ! 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 formats created by -OrbWeave- ! (for DOS) or OrbWin (for Windows). ! 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 element data). ! 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. ! *** 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 "UNITO". 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. PARAMETER (maxnod = 50000, maxbn = 4000, maxel = 90000, maxfel = 9000) PARAMETER (maxgsi = 900000000) LOGICAL * 1 checke, checkf, checkn, donext, edgefs, edgets, & & lgraph LOGICAL brief, redeem, sphere, orbdata5 CHARACTER*80 title1 DIMENSION atnode(maxnod), & & checke(maxel), checkf(maxfel), checkn(maxnod), & & cooling_curvature(maxnod), & & density_anomaly(maxnod), & & donext(maxnod), dqdtda(maxnod), & & eledat(3, maxel), elev(maxnod), & & edgefs(2, maxfel), edgets(3, maxel), & & 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 (B.NKI 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, 1) 1 FORMAT (/' Attempting to read existing .feg file from UNIT 1..'/) CALL Getnet (input, iuniti, iunitt, & & mxel, mxfel, mxnode, & & output, brief, & & cooling_curvature, density_anomaly, & & dqdtda, eledat, elev, fdip, & & nfl, nodef, nodes, nfaken, nrealn, n1000, & & numel, numnod, offmax, offset, orbdata5, & & title1, tlnode, xnode, ynode, zmnode, & & work, checke, checkf, checkn) WRITE (iunitt, 151) 151 FORMAT (' ....DONE.') WRITE (iunitt, 2) 2 FORMAT (/' SQUARE IS CHECKING FOR CORRECT TOPOLOGY...') CALL Square (input, brief, fdip, iunitt, & & mxbn, mxel, mxfel, mxnode, & & mxstar, nfl, nodef, nodes, & & numel, numnod, & & modify, xnode, ynode, & & output, edgefs, & & edgets, & & ncond, nodcon, & & work, checkn, list) WRITE (iunitt, 151) 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 (input, iunitt, modify, tinsec) CALL Number (input, iunitt, ncond, nodcon, & & mxel, mxfel, mxgsiz, mxnode, & & nfl, nodef, nodes, ntodo, & & numel, numnod, sphere, xnode, ynode, & & modify, tinsec, & & output, ialias, iuser, maxdif, mindif, & & work, donext, idiag, lgraph, & & maxlst) WRITE (iunitt, 151) CALL Clock (input, iunitt, modify, tinsec) WRITE (*, "(' which only counts time spent in SUBROUTINE Number.'/)") ! REPORT RESULTS TO USER ON DEVICE "IUNITT": !CCC WRITE (IUNITT,12) !CC12 FORMAT (/ /' BEST RENUMBERING SCHEME FOUND IN THIS RUN:'/ !CCC + /' OLD-NUMBER NEW-NUMBER INVERSE') !CCC DO 14 I=1,NUMNOD !CCC WRITE (IUNITT,13) I,IALIAS(I),IUSER(I) !CC13 FORMAT (' ',3I12) !CC14 CONTINUE WRITE (iunitt, 15) maxdif 15 FORMAT (/ /' BANDWIDTH BEFORE RENUMBERING = ',I5) WRITE (iunitt, 20) mindif 20 FORMAT (/ /' BANDWIDTH AFTER RENUMBERING = ',I5) IF (redeem) THEN DO 22 i = 1, numnod WRITE (iunito, 21) i, ialias(i), iuser(i) 21 FORMAT (3I10) 22 CONTINUE ELSE IF (mindif < maxdif) THEN ! REORDER THE VALUES IN THE SIMPLE FLOATING-POINT VECTOR ARRAYS: CALL Mixupx (input, ialias, mxnode, numnod, & & modify, xnode, & & work, atnode) CALL Mixupx (input, ialias, mxnode, numnod, & & modify, ynode, & & work, atnode) CALL Mixupx (input, ialias, mxnode, numnod, & & modify, elev, & & work, atnode) CALL Mixupx (input, ialias, mxnode, numnod, & & modify, dqdtda, & & work, atnode) CALL Mixupx (input, ialias, mxnode, numnod, & & modify, zmnode, & & work, atnode) CALL Mixupx (input, ialias, mxnode, numnod, & & modify, tlnode, & & work, atnode) IF (orbdata5) THEN CALL Mixupx (input, ialias, mxnode, numnod, & & modify, density_anomaly, & & work, atnode) CALL Mixupx (input, ialias, mxnode, numnod, & & modify, cooling_curvature, & & work, atnode) END IF ! MODIFY THE VALUES IN THE N-ENTRY INTEGER ARRAYS: n34 = 3 CALL Mixupi (input, ialias, mxel, mxnode, numel, n34, & & modify, nodes) n34 = 4 CALL Mixupi (input, ialias, mxfel, mxnode, nfl, n34, & & modify, nodef) ! OUTPUT MODIFIED FILE ON DEVICE "IUNITO": WRITE (iunitt, 90) 90 FORMAT (' RENUMBERING COMPLETED.'/ & & ' =========================================='/ & & ' About to write output grid (.feg) file' & & ,' on UNIT 2...'/) CALL Putnet (input, iunito, & & brief, & & cooling_curvature, density_anomaly, & & dqdtda, eledat, elev, fdip, & & mxel, mxfel, mxnode, n1000, & & nfaken, nfl, nodef, nodes, & & nrealn, numel, numnod, offset, orbdata5, & & title1, tlnode, xnode, ynode, zmnode) ELSE WRITE (iunitt, 99) 99 FORMAT (/ /' SINCE BANDWIDTH COULD NOT BE IMPROVED,', & & ' GRID WAS NOT MODIFIED,'/' AND NO OUTPUT .FEG', & & ' FILE WAS WRITTEN.') END IF END IF CALL PAUSE() ! CREATE ORDERED LIST OF BOUNDARY NODES, USING NEW NODE NUMBERS: IF (.NOT.sphere) THEN CALL Square (input, brief, fdip, iunitt, & & mxbn, mxel, mxfel, mxnode, & & mxstar, nfl, nodef, nodes, & & numel, numnod, & & modify, xnode, ynode, & & output, edgefs, & & edgets, & & ncond, nodcon, & & work, checkn, list) WRITE (iunitt, 101) ncond 101 FORMAT (/' There are ',I6,' nodes along the boundary.' & & /' Now, creating a NEW file containing an ordered list,' & & /' 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 Getnet (input, iunit7, iunit8, & & mxel, mxfel, mxnode, & & output, brief, & & cooling_curvature, density_anomaly, & & dqdtda, eledat, elev, fdip, & & nfl, nodef, nodes, nfaken, nrealn, n1000, & & numel, numnod, offmax, offset, orbdata5, & & title1, tlnode, xnode, ynode, zmnode, & & work, checke, checkf, checkn) ! READ FINITE ELEMENT GRID FROM UNIT "IUNIT7". ! ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8". CHARACTER*80 title1 LOGICAL allok, brief, orbdata5 ! NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL * 1 checke, checkf, checkn DIMENSION checke(mxel), checkf(mxfel), checkn(mxnode), & & cooling_curvature(mxnode), density_anomaly(mxnode), & & dqdtda(mxnode), eledat(3, mxel), elev(mxnode), & & 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), ELEVATIONS, HEAT-FLOW, AND ISOSTATIC ! GRAVITY ANOMALIES (ONE RECORD PER NODE). ! (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 (input, iunit7, iunit8, 9, & & output, vector) INDEX = vector(1) + 0.5 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.0 - plat) * 0.017453293 yi = plon * 0.017453293 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.0) .OR. & & (temp_cooling_curvature /= 0.0) IF (qi < 0.) THEN WRITE (iunit8, 96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL PAUSE() STOP END IF IF (zmi < 0.) THEN WRITE (iunit8, 97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') CALL PAUSE() STOP END IF zmnode(i) = zmi IF (tli < 0.) 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 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 CORNERS, COUNTER', & & 'CLOCKWISE'/ / & & ' ELEMENT C1 C2 C3') DO 200 k = 1, numel ! (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) CALL Readn (input, iunit7, iunit8, 7, & & output, vector) 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. IF (.NOT.brief) WRITE (iunit8, 120) i, (nodes(j, i), j = 1, 3), & & (eledat(j, i), j = 1, 3) 120 FORMAT (' ',I6,':',3I10,1P,3E10.2) 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. 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 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. READ(iunit7, * ) 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) WRITE (iunit8, 240) i, (nodef(j, i), j = 1, 4), & & (dips(l), l = 1, 2), off 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. + dips(l) fdip(l, i) = dips(l) * 0.017453293 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) 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.) THEN WRITE (iunit8, 400) offmax 400 FORMAT (/' GREATEST FAULT OFFSET READ WAS ',1P,E10.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 (' --------------------------------------------------', & & '-----------------------------') RETURN END SUBROUTINE GetNet SUBROUTINE Putnet (input, iunito, & & brief, & & cooling_curvature, density_anomaly, & & dqdtda, eledat, elev, fdip, & & mxel, mxfel, mxnode, n1000, & & nfaken, nfl, nodef, nodes, & & nrealn, numel, numnod, offset, orbdata5, & & title1, tlnode, xnode, ynode, zmnode) ! WRITES FINITE ELEMENT GRID TO UNIT "IUNITO". CHARACTER*80 title1 LOGICAL anyemu, brief, orbdata5 DIMENSION cooling_curvature(mxnode), density_anomaly(mxnode), & & dqdtda(mxnode), eledat(3, mxel), elev(mxnode), & & fdip(2, mxfel), np(4), & & nodef(4, mxfel), nodes(3, mxel), & & offset(mxfel), tlnode(mxnode), & & xnode(mxnode), ynode(mxnode), zmnode(mxnode) DIMENSION dips(2) 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. - xnode(i) * 57.29577951 plon = ynode(i) * 57.29577951 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 (I6,0P,2F11.5,1P,6E10.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.0) 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 IF (anyemu) THEN WRITE (iunito, 160) i, (np(k), k = 1, 3), (eledat(k, i), k = 1, 3) 160 FORMAT (I5,3I6,1P,E8.1,0P,F7.1,1P,E8.1) ELSE WRITE (iunito, 170) i, (np(k), k = 1, 3) 170 FORMAT (I5,3I6) 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.2958 IF (dips(k) > 90.01) dips(k) = dips(k) - 180. 230 CONTINUE WRITE (iunito, 250) i, (np(k), k = 1, 4), & & (dips(k), k = 1, 2), offset(i) 250 FORMAT (I5,4I6,2F6.1,1P,E10.2) 300 CONTINUE RETURN END SUBROUTINE PutNet SUBROUTINE Mixupx (input, inew, mxnode, numnod, & & modify, DATA, & & work, atnode) ! SORT FLOATING-POINT SCALAR VALUES IN ARRAY "DATA" SO THAT THE ! ENTRY CURRENTLY IN POSITION I ENDS UP IN POSITION "INEW(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 RETURN END SUBROUTINE Mixupx SUBROUTINE Mixupi (input, inew, mxel, mxnode, numel, n34, & & modify, nodes) ! CHANGE "N34"-ENTRY INTEGER SETS IN ARRAY "NODES" SO THAT THE ! INTEGER WHICH IS CURRENTLY "I" IS CHANGED TO "INEW(I)". 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 RETURN END SUBROUTINE Mixupi SUBROUTINE Next (input, i, j, mxel, mxfel, nfl, nodef, nodes, numel, & & output, kfault, kele) ! 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. 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 FIRST: 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 RETURN END SUBROUTINE Next REAL FUNCTION Atan2f (y, x) ! CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS IF ((y /= 0.).OR.(x /= 0.)) THEN Atan2f = ATAN2(y, x) ELSE Atan2f = 0. END IF RETURN END FUNCTION Atan2f SUBROUTINE Square (input, brief, fdip, iunit8, & & mxbn, mxel, mxfel, mxnode, & & mxstar, nfl, nodef, nodes, & & numel, numnod, & & modify, xnode, ynode, & & output, edgefs, & & edgets, & & ncond, nodcon, & & work, checkn, list) ! CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID LOGICAL agreed, allok, brief ! NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL * 1 checkn, edgefs, edgets DIMENSION corner(3, 3) 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) ! (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.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 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.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 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 toler = short / 10. 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.14) dy = dy - 6.28318 IF (dy < -3.14) dy = dy + 6.28318 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. ysum = 0. 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. 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.0) 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 COUNTERCLOCKWISE ! ORDER. ALSO MAKE A LISTS 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 (input, i, j, mxel, mxfel, nfl, nodef, nodes, numel, & & output, kfault, kele) 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 (' --------------------------------------------------', & & '-----------------------------') RETURN END SUBROUTINE Square SUBROUTINE ReadN (input, iunitp, iunitt, n, & & output, vector) ! 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. CHARACTER*1 c CHARACTER*132 line 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, 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. 99 CONTINUE END IF END IF RETURN END SUBROUTINE ReadN SUBROUTINE Number (input, iunitt, length, list, & & mxel, mxfel, mxgsiz, mxnode, & & nfl, nodef, nodes, nrealn, & & numel, numnod, sphere, xnode, ynode, & & modify, tinsec, & & output, ialias, iuser, maxdif, mindif, & & work, donext, idiag, lgraph, & & maxlst) ! 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 NODE NUMBERS ! 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. SCHWARZ 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 WELL 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, BOTH OF THESE CASES DO NOT 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.) ! 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 LOGICAL found, lessok, moreok, sphere 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) m1 = 1 m9 = length !=================================================================== WRITE (*, "(/' THERE WILL BE ',I6,' RENUMBERING ATTEMPTS:')") length 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. size = SQRT(ax * ax + bx * bx) ax = ax / size bx = bx / size ELSE ax = 1. bx = 0. gx = 0. 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.0 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. lessok = (ABS(xless) + ABS(yless)) > 0. IF (moreok.AND.lessok) THEN argmor = Atan2f(ymore, xmore) argles = Atan2f(yless, xless) ELSE IF (moreok) THEN argmor = Atan2f(ymore, xmore) argles = argmor - 3.14159 IF (argles < -3.14159) argles = argles + 6.28318 ELSE IF (lessok) THEN argles = Atan2f(yless, xless) argmor = argles + 3.14159 IF (argmor > 3.14159) argmor = argmor - 6.28318 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.) THEN argcut = argcut + 3.14159 IF (argcut > 3.14159) argcut = argcut - 6.28318 END IF ! *************************************** END IF ! SIZE OF BANDS TO BE USED AS LEVELS IN SPHERICAL NUMBERING sterad = 4. * 3.14159 / numel darc = 0.30 * 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.0 !----------------- 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 !CC WRITE (IUNITT,308) NTODO,LEVEL !C308 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. 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 ((m1 == m9).AND.(ialdif == minmax)) THEN WRITE (iunitt, 9005) i, j, ialias(i), ialias(j), & & minmax 9005 FORMAT (' OLD NODES ',I5,' AND ',I5,' HAVE', & & ' NEW NUMBERS ',I5,' AND ',I5/ & & ' WHOSE DIFFERENCE OF ',I5,' DEFINES', & & ' THE NEW BANDWIDTH.') END IF END IF 9010 CONTINUE 9020 CONTINUE maxlst(m) = mindif IF (istart == list(m)) THEN WRITE (iunitt, 9031) m, list(m), mindif ELSE WRITE (iunitt, 9032) m, list(m), istart, mindif END IF 9031 FORMAT(' ATTEMPT NUMBER ',I5,' BEGINS WITH NODE',I6,';' & & /' WHEN THIS NODE IS #1, THE BANDWIDTH IS ',I5,'.') 9032 FORMAT(' ATTEMPT NUMBER ',I5,' BEGINS NEAR BOUNDARY NODE',I6,',' & & /' AT NEAREST REAL NODE TO THAT, NUMBER ',I5,'.'/ & & ' WHEN THIS NODE IS #1, THE BANDWIDTH IS ',I5,'.') 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 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 RETURN END SUBROUTINE Number SUBROUTINE Clock (input, iunitt, modify, tinsec) ! REPORT(?) ELAPSED TIME BETWEEN CALLS told = tinsec CALL System_clock(ncount, npers, ncmax) IF (npers > 0) THEN tinsec = ncount / npers ELSE tinsec = 0. 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 RETURN END SUBROUTINE Clock SUBROUTINE Pause() IMPLICIT none WRITE(*, "(' Press [Enter]...')") READ(*, * ) RETURN END SUBROUTINE Pause