C C PROGRAM OrbNumber C (version of 13 February 2010: supports OrbData5 format, C and also provides a list of boundary nodes) C C by Peter Bird C Department of Earth and Space Sciences C University of California C Los Angeles, CA 90095-1567 C pbird@ess.ucla.edu C (C) Copyright 1993, 1994, 1996, 1998, 1999, 2005, 2010 C by Peter Bird and The C Regents of the University of California. C C Program to read a finite element grid (.feg) file, C renumber the nodes to minimize bandwidth, C and write out the result in the same format. C Accepts input files in the formats created by -OrbWeave- C (for DOS) or OrbWin (for Windows). C These .feg files may be of Shells-type (with faults, C and perhaps with nodal data filled in by C -OrbData- or OrbData5). C Or, they may be of NeoKinema type (without faults, C and with no nodal data except optional mu_ values). C Or, they may be of Restore type (without faults, C but with element data). C In the case of .feg files from OrbData5 (June 2005+) with C 6 nodal variables (instead of the previous 4) it recognizes C the new format silently and passes along all the values. C In the case of .feg files for Restore v.3, with element data, C it recognizes and passes along the element data. C C *** NOTE: SET VARIABLE "REDEEM" TO TRUE OR FALSE. C IF TRUE, IT MEANS THAT FAKE NODES (BOUNDARY NODES WITH NO C MECHANICAL DEGREES OF FREEDOM) WILL BE CONVERTED TO REAL C NODES DURING THE RENUMBERING. THIS OPTION IS USEFUL TO C PREPARE A RENUMBERING FOR USE BY THE GRAPHICS PROGRAMS. C IF "REDEEM" IS FALSE, FAKE NODES ARE NOT RENUMBERED, AND C WILL APPEAR AT THE END OF THE LIST AS BEFORE. C *** NOTE: C "REDEEM" ALSO CONTROLS THE TYPE OF OUTPUT FILE; IF TRUE, C THEN ONLY A LIST OF OLD NODE NUMBERS, NEW NODE NUMBERS, AND C A CROSS-REFERENCE (OLD NUMBERS LISTED IN ORDER OF NEW) IS C OUTPUT ON UNIT "UNITO". IF "REDEEM" IS FALSE, HOWEVER, A C COMPLETE, RESORTED FINITE ELEMENT GRID FILE IS WRITTEN. C THIS REFLECTS THE PRACTICAL REALITY THAT ONE USUALLY C CHOOSES "REDEEM=.TRUE." WHEN PREPARING TO DO GRAPHICS, BUT C CHOOSES "REDEEM=.FALSE." WHEN OPTIMIZING A HAND-EDITED GRID C PRIOR TO FINITE ELEMENT CALCULATIONS. C PARAMETER (MAXNOD=20000,MAXBN=2000,MAXEL=30000,MAXFEL=6000) PARAMETER (MAXGSI=130000000) C LOGICAL*1 CHECKE,CHECKF,CHECKN,DONEXT,EDGEFS,EDGETS, + LGRAPH LOGICAL BRIEF,REDEEM,SPHERE,OrbData5 CHARACTER*80 TITLE1 C 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) C C THE INPUT FINITE-ELEMENT GRID FILE (.FEG FILE) IS READ FROM: DATA IUNITI /1/ C C THE OUTPUT FINITE-ELEMENT GRID FILE (.FEG FILE) IS WRITTEN TO: DATA IUNITO /2/ C C THE OUTPUT LIST OF BOUNDARY NODES (B.NKI FILE) IS WRITTEN TO: DATA IUNITB /3/ C C TEXT MESSAGES DURING EXECUTION ARE SENT TO: DATA IUNITT /6/ C C************************ CHECK THE FOLLOWING LINE! ******************* C C USE .TRUE. TO CREATE ONLY A LIST OF NEW NODE NUMBERS VERSUS OLD. REDEEM=.FALSE. C USE .FALSE. TO CREATE A COMPLETE RENUMBERED .FEG FILE C C************************ CHECK THE PRECEDING LINE! ******************* C MXNODE=MAXNOD MXBN=MAXBN MXEL=MAXEL MXFEL=MAXFEL MXSTAR=MAXNOD MXGSIZ=MAXGSI C CALL SYSTEM_CLOCK(NCOUNT,NPERS,NCMAX) C TO PROVIDE TIME-ELAPSED, WE ARE USING JUST THIS ONE INTRINSIC C SUBROUTINE FROM FORTRAN90, BECAUSE THERE IS NO SUCH STANDARD C ROUTINE IN FORTRAN77. IF YOU DO NOT HAVE ACCESS TO FORTRAN90, C THEN EITHER SUBSTITUTE AN EQUIVALENT ROUTINE, OR DO WITHOUT THE C ELAPSED-TIME INFORMATION. IF (NPERS.GT.0) THEN TINSEC=NCOUNT/NPERS ELSE TINSEC=0. ENDIF WRITE (IUNITT,1) 1 FORMAT (/' Initializing the clock;' + /' Attempting to read existing .feg file from UNIT 1..'/) C 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.') C CALL CLOCK (INPUT,IUNITT,MODIFY,TINSEC) C 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) C CALL CLOCK (INPUT,IUNITT,MODIFY,TINSEC) C IF (REDEEM) THEN NTODO=NUMNOD ELSE NTODO=NREALN ENDIF C IF (NCOND.EQ.0) THEN SPHERE=.TRUE. IF (NFL.EQ.0) THEN C ARBITRARILY CHOOSING TO BEGIN WITH NODE 1: NCOND=1 NODCON(1)=1 ELSE C 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 ENDIF ELSE SPHERE=.FALSE. ENDIF C WRITE (IUNITT,10) 10 FORMAT (/' CALLING NUMBER...') C 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) C CALL CLOCK (INPUT,IUNITT,MODIFY,TINSEC) C C REPORT RESULTS TO USER ON DEVICE "IUNITT": C CCCC WRITE (IUNITT,12) CCC12 FORMAT (/ /' BEST RENUMBERING SCHEME FOUND IN THIS RUN:'/ CCCC + /' OLD-NUMBER NEW-NUMBER INVERSE') CCCC DO 14 I=1,NUMNOD CCCC WRITE (IUNITT,13) I,IALIAS(I),IUSER(I) CCC13 FORMAT (' ',3I12) CCC14 CONTINUE C WRITE (IUNITT,15) MAXDIF 15 FORMAT (/ /' BANDWIDTH BEFORE RENUMBERING = ',I5) WRITE (IUNITT,20) MINDIF 20 FORMAT (/ /' BANDWIDTH AFTER RENUMBERING = ',I5) C IF (REDEEM) THEN DO 22 I=1,NUMNOD WRITE (IUNITO,21) I, IALIAS(I), IUSER(I) 21 FORMAT (3I10) 22 CONTINUE ELSE IF (MINDIF.LT.MAXDIF) THEN C C REORDER THE VALUES IN THE SIMPLE FLOATING-POINT VECTOR ARRAYS: C 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 C C MODIFY THE VALUES IN THE N-ENTRY INTEGER ARRAYS: C N34=3 CALL MIXUPI (INPUT,IALIAS,MXEL,MXNODE,NUMEL,N34, + MODIFY,NODES) N34=4 CALL MIXUPI (INPUT,IALIAS,MXFEL,MXNODE,NFL,N34, + MODIFY,NODEF) C C OUTPUT MODIFIED FILE ON DEVICE "IUNITO": C 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.') ENDIF ENDIF CALL PAUSE() C C CREATE ORDERED LIST OF BOUNDARY NODES, USING NEW NODE NUMBERS: C 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 WRITE (IUNITT, 111) 111 FORMAT (/' Job completed.') CALL PAUSE() STOP END C C C 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) C C READ FINITE ELEMENT GRID FROM UNIT "IUNIT7". C ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8". C CHARACTER*80 TITLE1 LOGICAL ALLOK,BRIEF,OrbData5 C C NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL*1 CHECKE,CHECKF,CHECKN C 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) C 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) C C READ NUMBER OF NODES. C INPUT NODAL LOCATIONS (X,Y), ELEVATIONS, HEAT-FLOW, AND ISOSTATIC C GRAVITY ANOMALIES (ONE RECORD PER NODE). C (OPTION "BRIEF" SUPPRESSES MOST OUTPUT) C READ (IUNIT7,*) NUMNOD,NREALN,NFAKEN,N1000,BRIEF BRIEF=.TRUE. C IF (NUMNOD.NE.(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 ENDIF C IF (NREALN.GT.N1000) THEN WRITE (IUNIT8,5) NREALN,N1000 5 FORMAT (/' ERR0R: NREALN (',I6,') IS GREATER THAN' + /' N1000 (',I6,').') CALL PAUSE() STOP ENDIF C IF (NUMNOD.GT.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 ENDIF C 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'/) ENDIF DO 90 K=1,NUMNOD CHECKN(K)=.FALSE. 90 CONTINUE C OrbData5 = .FALSE. C unless proven otherwise, below, by non-zero input values... C DO 100 K=1,NUMNOD CALL READN (INPUT,IUNIT7,IUNIT8,9, + OUTPUT,VECTOR) INDEX=VECTOR(1)+0.5 IF (INDEX.GT.NREALN) THEN IF ((INDEX.LE.N1000).OR. + (INDEX.GT.(N1000+NFAKEN))) THEN WRITE (IUNIT8,91) INDEX 91 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER: ',I6) CALL PAUSE() STOP ENDIF ENDIF 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) C Note: If READN did not find these last 2 values, it would have set them to 0.0 temp_cooling_curvature = VECTOR(9) IF (INDEX.LE.NREALN) THEN I=INDEX ELSE I=NREALN+INDEX-N1000 ENDIF 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.LT.0.) THEN WRITE (IUNIT8,96) 96 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') CALL PAUSE() STOP ENDIF IF (ZMI.LT.0.) THEN WRITE (IUNIT8,97) 97 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') CALL PAUSE() STOP ENDIF ZMNODE(I)=ZMI IF (TLI.LT.0.) THEN WRITE (IUNIT8,98) 98 FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS IS', + ' NON-PHYSICAL.') CALL PAUSE() STOP ENDIF 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) ENDIF 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.LE.NREALN) THEN INDEX=I ELSE INDEX=N1000+I-NREALN ENDIF IF (.NOT.CHECKN(I)) WRITE(IUNIT8,103)INDEX 103 FORMAT (' ',36X,I6) 104 CONTINUE CALL PAUSE() STOP ENDIF C C READ TRIANGULAR ELEMENTS C READ (IUNIT7,*,IOSTAT=IOS) NUMEL IF (NUMEL.GT.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 ENDIF 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 C (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.LT.1).OR.(I.GT.NUMEL)) THEN WRITE (IUNIT8,111) I 111 FORMAT (/' ERR0R: ILLEGAL ELEMENT NUMBER: ',I6) CALL PAUSE() STOP ENDIF 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.GT.NREALN) N=NREALN+(N-N1000) IF ((N.LE.0).OR.(N.GT.NUMNOD)) THEN WRITE (IUNIT8,125) NODES(J,I) 125 FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.') CALL PAUSE() STOP ENDIF 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 ENDIF C C READ FAULT ELEMENTS C READ (IUNIT7,*) NFL IF (NFL.GT.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 ENDIF 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.GT.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.LT.1).OR.(I.GT.NFL)) THEN WRITE (IUNIT8,241) I 241 FORMAT (/' ERR0R: ILLEGAL FAULT NUMBER: ',I6) CALL PAUSE() STOP ENDIF 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.GT.NREALN) N=NREALN+(N-N1000) IF ((N.LE.0).OR.(N.GT.NUMNOD)) THEN WRITE (IUNIT8,243) NODEF(J,I),I 243 FORMAT (/' ERR0R: ILLEGAL NODE NUMBER (',I6, + ') IN FAULT ',I6) CALL PAUSE() STOP ENDIF NODEF(J,I)=N 250 CONTINUE DO 260 L=1,2 IF (ABS(DIPS(L)).GT.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 ENDIF IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L) FDIP(L,I)=DIPS(L)*0.017453293 260 CONTINUE IF (OFF.LT.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 ENDIF 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.GT.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.') ENDIF ENDIF IF (.NOT. BRIEF) WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C 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) C C WRITES FINITE ELEMENT GRID TO UNIT "IUNITO". C CHARACTER*80 TITLE1 LOGICAL ANYEMU,BRIEF, OrbData5 C 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) C WRITE (IUNITO,1) TITLE1 1 FORMAT (A80) C WRITE (IUNITO,2) NUMNOD, NREALN,NFAKEN,N1000,BRIEF 2 FORMAT(4I8,L8,' (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)') C DO 100 I=1,NUMNOD IF (I.LE.NREALN) THEN IPRINT=I ELSE IPRINT=N1000+(I-NREALN) ENDIF 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 C WRITE (IUNITO,110) NUMEL 110 FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') C ANYEMU=.FALSE. DO 130 I=1,NUMEL DO 120 K=1,3 IF (ELEDAT(K,I).NE.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).LE.NREALN) THEN NP(K)=NODES(K,I) ELSE NP(K)=N1000+(NODES(K,I)-NREALN) ENDIF 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 C WRITE (IUNITO,210) NFL 210 FORMAT (I10,' (NFL = NUMBER OF CURVILINEAR FAULT ELEMENTS)') C DO 300 I=1,NFL DO 220 K=1,4 IF (NODEF(K,I).LE.NREALN) THEN NP(K)=NODEF(K,I) ELSE NP(K)=N1000+(NODEF(K,I)-NREALN) ENDIF 220 CONTINUE DO 230 K=1,2 DIPS(K)=FDIP(K,I) DIPS(K)=DIPS(K)*57.2958 IF (DIPS(K).GT.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 C RETURN END C C C SUBROUTINE MIXUPX (INPUT,INEW,MXNODE,NUMNOD, + MODIFY,DATA, + WORK,ATNODE) C C SORT FLOATING-POINT SCALAR VALUES IN ARRAY "DATA" SO THAT THE C ENTRY CURRENTLY IN POSITION I ENDS UP IN POSITION "INEW(I)". C DIMENSION ATNODE(MXNODE),INEW(MXNODE), + DATA(MXNODE) C DO 10 I=1,NUMNOD ATNODE(I)=DATA(I) 10 CONTINUE C DO 20 I=1,NUMNOD DATA(INEW(I))=ATNODE(I) 20 CONTINUE C RETURN END C C C SUBROUTINE MIXUPI (INPUT,INEW,MXEL,MXNODE,NUMEL,N34, + MODIFY,NODES) C C CHANGE "N34"-ENTRY INTEGER SETS IN ARRAY "NODES" SO THAT THE C INTEGER WHICH IS CURRENTLY "I" IS CHANGED TO "INEW(I)". C DIMENSION INEW(MXNODE), + NODES(N34,MXEL) C DO 20 K=1,N34 DO 10 J=1,NUMEL NODES(K,J)=INEW(NODES(K,J)) 10 CONTINUE 20 CONTINUE C RETURN END C C C SUBROUTINE NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) C C DETERMINE WHETHER THERE ARE MORE ELEMENTS ADJACENT TO SIDE J OF C TRIANGULAR CONTINUUM ELEMENT #I. C J = 1 MEANS THE SIDE OPPOSITE NODE # NODES(1,I). C J = 2 MEANS THE SIDE OPPOSITE NODE # NODES(2,I). C J = 3 MEANS THE SIDE OPPOSITE NODE # NODES(3,I). C IF A FAULT ELEMENT IS ADJACENT, ITS NUMBER IS KFAULT; OTHERWISE, C KFAULT IS SET TO ZERO. C IF ANOTHER TRIANGULAR CONTINUUM ELEMENT IS ADJACENT (EVEN ACROSS C FAULT ELEMENT KFAULT!) THEN ITS NUMBER IS KELE; OTHERWISE, KELE = 0. C LOGICAL FOUNDF DIMENSION NODEF(4,MXFEL),NODES(3,MXEL) C C THREE NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE: N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J+1,3)+1,I) C CHECK FOR ADJACENT FAULT ELEMENT FIRST: FOUNDF=.FALSE. KFAULT=0 IF (NFL.GT.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.EQ.N2).AND.(M2.EQ.N1)).OR. + ((M3.EQ.N2).AND.(M4.EQ.N1))) THEN FOUNDF=.TRUE. KFAULT=K GO TO 11 ENDIF 10 CONTINUE 11 IF (.NOT.FOUNDF) KFAULT=0 C IF THERE WAS A FAULT, REPLACE 2 NODE NUMBERS THAT WE SEARCH FOR: IF (FOUNDF) THEN IF (M2.EQ.N1) THEN N1=M3 N2=M4 ELSE N1=M1 N2=M2 ENDIF ENDIF ENDIF C SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT: KELE=0 KLOW=I KHIGH=I C --- BEGIN IRREGULAR LOOP, TO SEARCH OUT NEAREST ELEMENTS FIRST --- 100 KLOW=KLOW-1 IF (KLOW.GE.1) THEN DO 110 L=1,3 M1=NODES(MOD(L, 3)+1,KLOW) M2=NODES(MOD(L+1,3)+1,KLOW) IF ((M2.EQ.N1).AND.(M1.EQ.N2)) THEN KELE=KLOW RETURN ENDIF 110 CONTINUE ENDIF KHIGH=KHIGH+1 IF (KHIGH.LE.NUMEL) THEN DO 120 L=1,3 M1=NODES(MOD(L, 3)+1,KHIGH) M2=NODES(MOD(L+1,3)+1,KHIGH) IF ((M2.EQ.N1).AND.(M1.EQ.N2)) THEN KELE=KHIGH RETURN ENDIF 120 CONTINUE ENDIF IF ((KLOW.GT.1).OR.(KHIGH.LT.NUMEL)) GO TO 100 RETURN END C C C REAL FUNCTION ATAN2F (Y,X) C C CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS C IF ((Y.NE.0.).OR.(X.NE.0.)) THEN ATAN2F=ATAN2(Y,X) ELSE ATAN2F=0. ENDIF RETURN END C C C 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) C C CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID C LOGICAL AGREED,ALLOK,BRIEF C C NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL*1 CHECKN,EDGEFS,EDGETS C 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) C C (1) CHECK THAT ALL NODES ARE CONNECTED TO AT LEAST ONE C CONTINUUM (TRIANGULAR) ELEMENT OR FAULT ELEMENT; C 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'/ 1 ' BELONG TO ANY TRIANGULAR CONTINUUM ELEMENT'/ 2 ' 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 ENDIF C C (2) AVERAGE TOGETHER THE COORDINATES OF ALL NODES AT ONE "POINT" C WRITE (IUNIT8,2) 2 FORMAT (' IN SQUARE: AVERAGING NODE COORDINATES') DO 410 I=1,NUMNOD CHECKN(I)=.FALSE. C (MEANS "NOT YET INVOLVED IN AVERAGING') 410 CONTINUE DO 490 I=1,NFL DO 480 J1=1,2 NJ1=NODEF(J1,I) C (FAULT ENDS ARE THE ONLY PLACES THAT CAN HAVE PROBLEMS) IF (.NOT.CHECKN(NJ1)) THEN LIST(1)=NJ1 CHECKN(NJ1)=.TRUE. C BEGIN LIST OF NEIGHBORS WITH PAIRED NODE J2=5-J1 NJ2=NODEF(J2,I) LIST(2)=NJ2 CHECKN(NJ2)=.TRUE. NINSUM=2 C FIND SHORTEST FAULT CONNECTED TO EITHER ONE DX=XNODE(NJ1)-XNODE(NJ2) DY=YNODE(NJ1)-YNODE(NJ2) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-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.EQ.NL1).OR.(NJ2.EQ.NL1).OR. + (NJ1.EQ.NL2).OR.(NJ2.EQ.NL2).OR. + (NJ1.EQ.NL3).OR.(NJ2.EQ.NL3).OR. + (NJ1.EQ.NL4).OR.(NJ2.EQ.NL4)) THEN DX=XNODE(NL1)-XNODE(NL2) DY=YNODE(NL1)-YNODE(NL2) IF (DY.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NL1)) TEST=SQRT(DX**2+DY**2) SHORT=MIN(SHORT,TEST) ENDIF 470 CONTINUE C 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.GT.3.14) DY=DY-6.28318 IF (DY.LT.-3.14) DY=DY+6.28318 DY=DY*SIN(XNODE(NJ1)) R2=DX**2+DY**2 IF (R2.LT.T2) THEN NINSUM=NINSUM+1 IF (NINSUM.GT.MXSTAR) THEN WRITE(IUNIT8,421) 421 FORMAT(/' INCREASE VALUE' + ,' OF PARAMETER MAXATP.') CALL PAUSE() STOP ENDIF LIST(NINSUM)=K CHECKN(K)=.TRUE. ENDIF ENDIF 471 CONTINUE C (QUICK EXIT IF ALL NODES IN SAME PLACE) AGREED=.TRUE. DO 472 K=2,NINSUM AGREED=AGREED.AND. + (XNODE(LIST(K)).EQ.XNODE(LIST(1))).AND. + (YNODE(LIST(K)).EQ.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.GT.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!') ENDIF ENDIF ENDIF 480 CONTINUE 490 CONTINUE C C (3) CHECK FOR NEGATIVE AREAS C 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.LE.0.0) THEN WRITE (IUNIT8,509) I 509 FORMAT (/ /' ERR0R: AREA OF ELEMENT ',I5, + ' IS NOT POSITIVE.') CALL PAUSE() STOP ENDIF 600 CONTINUE C C (5) MAKE A LIST OF NODES THAT ARE ON THE BOUNDARY AND REQUIRE C BOUNDARY CONDITIONS (NODCON); THESE ARE IN COUNTERCLOCKWISE C ORDER. ALSO MAKE A LISTS OF ELEMENT SIDES WHICH CONTAIN THESE C NODES: EDGETS AND EDGEFS. C 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.GT.0) THEN DO 802 I=1,NFL EDGEFS(1,I)=.FALSE. EDGEFS(2,I)=.FALSE. 802 CONTINUE ENDIF 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.GT.0) THEN C (ORDINARY INTERIOR SIDE) EDGETS(J,I)=.FALSE. ELSE IF (KFAULT.EQ.0) THEN C (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. ENDIF IF (.NOT.CHECKN(N2)) THEN NCOND=NCOND+1 CHECKN(N2)=.TRUE. ENDIF ELSE C (TRIANGULAR ELEMENT HAS AN EXTERIOR FAULT ELEMENT C ADJACENT TO IT) EDGETS(J,I)=.FALSE. N1=NODES(MOD(J, 3)+1,I) IF (NODEF(2,KFAULT).EQ.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. ENDIF 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. ENDIF 808 CONTINUE ENDIF ENDIF 809 CONTINUE 810 CONTINUE IF (NCOND.GT.MXBN) THEN WRITE(IUNIT8,820) NCOND 820 FORMAT(/' INCREASE PARAMETER MAXBN TO',I6,' AND RECOMPILE.') CALL PAUSE() STOP ENDIF C C STOP WORK IF NO BOUNDARY NODES FOUND (GLOBAL GRID) C IF (NCOND.EQ.0) GO TO 899 C C 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 C BEGINNING OF INDEFINATE LOOP WHICH TRACES AROUND THE PERIMETER. C EACH TIME, IT PROGRESSES BY ONE OF 3 STEPS: C -2 NODES AT A TIME ALONG A TRIANGLE SIDE, OR C -2 NODES AT A TIME ALONG A FAULT ELEMENT SIDE, OR C -BY FINDING ANOTHER (CORNER) NODE WHICH SHARES THE SAME LOCATION. C BEGINNING OF MAIN INDEFINATE LOOP: 840 NODE=NODCON(NDONE) C C IMPORTANT: CHECK THAT WE ARE NOT REVISITING A NODE! C THIS WOULD MEAN THAT THERE ARE TOO MANY BOUNDARY NODES C TO FIT IN THE SIMPLY-CONNECTED LOOP, AND THAT THERE C 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.GT.1) CHECKN(NODE)=.FALSE. C X=XNODE(NODE) Y=YNODE(NODE) C 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.EQ.NODE) GO TO 846 ENDIF 842 CONTINUE 844 CONTINUE GO TO 850 846 N2=NODES(MOD(J+1,3)+1,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 NLEFT=NLEFT-1 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C 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).EQ.NODE) THEN N2=NODEF(2,I) GO TO 856 ENDIF ELSE IF (EDGEFS(2,I)) THEN IF (NODEF(3,I).EQ.NODE) THEN N2=NODEF(4,I) GO TO 856 ENDIF ENDIF 854 CONTINUE GO TO 860 856 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 NLEFT=NLEFT-1 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR ANOTHER EXTERIOR CORNER NODE AT SAME LOCATION 860 DO 865 I=1,NUMNOD IF ((I.NE.NODE).AND.CHECKN(I)) THEN IF ( (ABS(XNODE(I)-X).LT.1.E-6) .AND. + (ABS(YNODE(I)-Y).LT.1.E-6) ) GO TO 867 ENDIF 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.LE.NCOND) NODCON(NDONE)=I NLEFT=NLEFT-1 IF (NLEFT.GT.0) GO TO 840 C 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.)') ENDIF 899 CONTINUE C IF (.NOT. BRIEF) WRITE (IUNIT8,9999) 9999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END C C C SUBROUTINE READN (INPUT,IUNITP,IUNITT,N, + OUTPUT,VECTOR) C C A UTILITY ROUTINE DESIGNED TO PERMIT -FAULTS- INPUT FILES C TO ALSO BE USED BY -PLATES-, WHICH EXPECTS MORE NUMBERS C IN SOME RECORDS. C Also supports reading either OrbData files (with 4 nodal C variables) or OrbData5 files (with 6 nodal variables). C THIS ROUTINE ATTEMPTS TO READ 'N' FLOATING-POINT VALUES C (USING * FORMAT) FROM THE NEXT RECORD ON DEVICE 'IUNITP'. C IF ANYTHING GOES WRONG, THE MISSING VALUES ARE SET TO ZERO. C CHARACTER*1 C CHARACTER*132 LINE LOGICAL ANYIN,DOTTED,EXPON,SIGNED DIMENSION VECTOR(N) C LINE=' '// + ' ' READ (IUNITP,1,IOSTAT=IOS) LINE 1 FORMAT (A) C NUMBER=0 ANYIN=.FALSE. EXPON=.FALSE. SIGNED=.FALSE. DOTTED=.FALSE. DO 10 I=1,132 C=LINE(I:I) IF ((C.EQ.' ').OR.(C.EQ.',').OR.(C.EQ.'/')) THEN SIGNED=.FALSE. EXPON=.FALSE. DOTTED=.FALSE. IF (ANYIN) THEN NUMBER=NUMBER+1 ANYIN=.FALSE. ENDIF ELSE IF ((C.EQ.'+').OR.(C.EQ.'-')) THEN IF (SIGNED) THEN GO TO 50 ELSE SIGNED=.TRUE. ENDIF ELSE IF ((C.EQ.'E').OR.(C.EQ.'D').OR. + (C.EQ.'e').OR.(C.EQ.'d')) THEN IF (EXPON) THEN GO TO 50 ELSE EXPON=.TRUE. SIGNED=.FALSE. DOTTED=.TRUE. ENDIF ELSE IF (C.EQ.'.') THEN IF (DOTTED) THEN GO TO 50 ELSE DOTTED=.TRUE. ENDIF ELSE IF ((C.EQ.'0').OR.(C.EQ.'1').OR.(C.EQ.'2').OR. + (C.EQ.'3').OR.(C.EQ.'4').OR.(C.EQ.'5').OR. + (C.EQ.'6').OR.(C.EQ.'7').OR.(C.EQ.'8').OR. + (C.EQ.'9')) THEN SIGNED=.TRUE. ANYIN=.TRUE. ELSE GO TO 50 ENDIF 10 CONTINUE IF (ANYIN) NUMBER=NUMBER+1 C 50 IF (NUMBER.EQ.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 (NUMBER.GE.N) THEN READ (LINE,*) (VECTOR(I),I=1,N) ELSE READ (LINE,*) (VECTOR(I),I=1,NUMBER) DO 99 I=NUMBER+1,N VECTOR(I)=0. 99 CONTINUE ENDIF ENDIF RETURN END C C C 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) C C FIND A GENERAL RENUMBERING OF THE NODES OF THE FINITE-ELEMENT C GRID WHICH IS (NEARLY) OPTIMAL IN TERMS OF MINIMIZING THE C BANDWIDTH OF THE CONNECTIVITY MATRIX (A SURROGATE C FOR THE STIFFNESS MATRIX, BUT SYMMETRICAL, AND WITH ONLY ONE C ROW AND ONE COLUMN PER NODE). C C IN THE SECTION WHICH FILLS THE CONNECTIVITY MATRIX C "LGRAPH", IT IS ASSUMED THAT THE GRID IS MADE UP OF C A SET OF "NUMEL" 3-NODE TRIANGULAR ELEMENTS (WHOSE OLD C NODE NUMBERS ARE STORED AS "NODES(K=1..3,I=1..NUMEL)"), AND C "NFL" 4-NODE FAULT ELEMENTS (WHOSE OLD NODE NUMBERS ARE C STORED AS "NODEF(K=1..4,I=1..NFL)" ). C NFL MAY BE ZERO. C THIS SECTION CAN BE EASILY MODIFIED FOR OTHER TYPES OF ELEMENTS. C IT DOES NOT MATTER WHETHER OR NOT THE NUMBERS DEFINING C EACH ELEMENT OCCUR IN ANY PARTICULAR ORDER, OR NOT. C DEGENERATE ELEMENTS WITH TWO OR MORE IDENTICAL NODE NUMBERS C ARE ALSO PERMITTED. C C THE TOTAL NUMBER OF NODES IN THE GRID IS INPUT PARAMETER C "NUMNOD". THE FIRST "NREALN" OF THESE WILL BE RENUMBERED. C USUALLY, NUMNOD WILL EQUAL NREALN. OBVIOUSLY, NREALN CANNOT C BE GREATER THAN NUMNOD. THE DISTINCTION IS MADE IN CASE C YOU ARE USING A FINITE ELEMENT PROGRAM LIKE "FAULTS" IN C WHICH FAKE NODES (THOSE WITH NO DEGREES OF FREEDOM) ARE C LISTED LAST (WITH NUMBERS FROM NREALN+1 TO NUMNOD), AND C YOU WANT TO KEEP THIS GROUP AT THE END, AND NOT RENUMBER C THEM. C C THE ALGORITHM USED IS AN ORIGINAL VARIANT OF THE C CUTHILL-MCKEE ALGORITHM C DESCRIBED BY H. R. SCHWARZ IN "FINITE ELEMENT METHODS", C "COMPUTATIONAL MATHEMATICS AND APPLICATIONS" SERIES, C ACADEMIC PRESS/HARCOURT-BRACE-JOVANOVITCH, 1988, PP.167-185. C C INSTEAD OF USING DEGREE TO DECIDE WHICH NODES WILL BE NUMBERED C FIRST, THIS VARIANT NUMBERS NODES WITHIN ONE LEVEL ON THE BASIS C OF THEIR AZIMUTH FROM THE INITIAL NODE. C C BECAUSE OF THIS CHOICE OF ALGORITHM, THIS CODE CAN NOT BE EXPECTED C TO WORK WELL FOR GRIDS WHICH ARE NOT SIMPLY CONNECTED (I.E., C THOSE WITH HOLES IN THEM), AND IT CANNOT BE EXPECTED TO WORK C WELL FOR GRIDS WHICH HAVE EXTREMELY CONCAVE BOUNDARIES. C (FORTUNATELY, BOTH OF THESE CASES DO NOT ARISE IN EARTH SCIENCES, C WHERE THE SIDE BOUNDARIES OF THE GRID ARE ALWAYS ARTIFICIAL, AND C ARE ALMOST ALWAYS CHOSEN AS A CONVEX POLYHEDRON.) C C THE PRODUCTS OF THIS SUBPROGRAM ARE: C *"IALIAS" = AN INTEGER VECTOR OF NEW NODE NUMBERS, ARRANGED C IN THE ORDER OF THE OLD NODE NUMBERS. C *"IUSER" = AN INTEGER VECTOR OF OLD NODE NUMBERS, ARRANGED C IN THE ORDER OF THE NEW NODE NUMBERS. C "IALIAS" AND "IUSER" ARE MUTUAL INVERSES, AS SHOWN BY C THE FOLLOWING THEOROMS: C IUSER(IALIAS(I))=I C IALIAS(IUSER(J))=J C * MAXDIF = HALF BANDWIDTH BEFORE RENUMBERING. C * MINDIF = HALF BANDWIDTH AFTER RENUMBERING. C (THE DEFINITION OF HALF-BANDWIDTH USED HERE C DOES NOT INCLUDE THE DIAGONAL, AND CAN BE DEFINED AS THE C MAXIMUM OF THE ABSOLUTE VALUE OF THE DIFFERENCE OF THE C (NEW) NUMBERS IN ALL POSSIBLE PAIRS OF CONNECTED NODES.) C C NOTE: THE FOLLOWING TYPE DECLARATION IS VERY IMPORTANT FOR C ECONOMIZING STORAGE, AS LGRAPH WILL HAVE A SIZE OF C APPROXIMATELY NREALN*(NREALN-2)/2 LOGICAL ENTRIES. LOGICAL*1 DONEXT,LGRAPH C LOGICAL FOUND,LESSOK,MOREOK,SPHERE C 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) C C STATEMENT FUNCTION, GIVING STORAGE LOCATION (VECTOR INDEX) OF THE C ELEMENT OF THE GRAPH MATRIX AT ROW "IROW" AND COLUMN "JCOL": INDEXN(IROW,JCOL)=IDIAG(MIN(IROW,JCOL))+ABS(JCOL-IROW) C IF (NREALN.GT.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 ENDIF C C COMPUTE THE LOCATIONS OF THE DIAGONAL ELEMENTS OF C THE CONNECTIVITY MATRIX "LGRAPH". C (NOTE THAT THE GRAPH IS SYMMETRIC, SO WE ONLY STORE THE C UPPER TRIANGLE, AND WE ACCESS THE MATRIX C INDIRECTLY USING STATEMENT FUNCTION "INDEXN".) C NUSED=0 DO 10 I=1,NREALN IDIAG(I)=NUSED NINROW=NREALN-I NUSED=NUSED+NINROW 10 CONTINUE IF (NUSED.GT.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 ENDIF C C CREATE THE GRAPH MATRIX C WRITE (IUNITT,12) 12 FORMAT (/' CREATING THE GRAPH MATRIX...') DO 20 I=1,MXGSIZ LGRAPH(I)=.FALSE. 20 CONTINUE C 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.LT.1).OR.(IROW.GT.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 ENDIF IF ((JCOL.LT.1).OR.(JCOL.GT.NUMNOD)) THEN WRITE (IUNITT,21) JCOL,I CALL PAUSE() STOP ENDIF IF ((IROW.LE.NREALN).AND.(JCOL.LE.NREALN) + .AND.(IROW.NE.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.LT.1).OR.(IROW.GT.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 ENDIF IF ((JCOL.LT.1).OR.(JCOL.GT.NUMNOD)) THEN WRITE (IUNITT,51) JCOL,I CALL PAUSE() STOP ENDIF IF ((IROW.LE.NREALN).AND.(JCOL.LE.NREALN) + .AND.(IROW.NE.JCOL)) + LGRAPH(INDEXN(IROW,JCOL))=.TRUE. 130 CONTINUE 140 CONTINUE 150 CONTINUE WRITE (IUNITT,151) 151 FORMAT (' ....DONE.') C CALL CLOCK (INPUT,IUNITT,MODIFY,TINSEC) C C C DETERMINE THE ORIGINAL BANDWIDTH C 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) C CALL CLOCK (INPUT,IUNITT,MODIFY,TINSEC) C M1=1 M9=LENGTH C C=================================================================== C 200 DO 10000 M=M1,M9 C C CREATE A RENUMBERING SCHEME, USING INITIAL NODE "LIST(M)" C C (NECESSARY INITIALIZATION; RESET NUMBERING SCHEME TO NULL) DO 201 I=1,NREALN IALIAS(I)=0 IUSER(I)=0 201 CONTINUE C C (BEGIN REAL WORK) C C LOCATE INITIAL NODE C THETA1=XNODE(LIST(M)) PHI1=YNODE(LIST(M)) A1=COS(PHI1)*SIN(THETA1) B1=SIN(PHI1)*SIN(THETA1) G1=COS(THETA1) C PRIME MERIDIAN VECTOR (X) IS PERPENDICULAR TO INITIAL NODE IF (SIN(THETA1).GT.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. ENDIF C 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 C C DETERMINE ARGUMENT OF CUT IN AZIMUTH FROM INITIAL NODE C IF (SPHERE) THEN C LOCATION OF CUT IS ARBITRARY! ARGCUT=0.0 ELSE IF (M.LT.LENGTH) THEN MMORE=M+1 ELSE MMORE=1 ENDIF IF (M.GT.1) THEN MLESS=M-1 ELSE MLESS=LENGTH ENDIF 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) C "X" AND "Y" ARE IN EQUATORIAL PLANE RELATIVE TO C 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)).GT.0. LESSOK=(ABS(XLESS)+ABS(YLESS)).GT.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.LT.-3.14159) ARGLES=ARGLES+6.28318 ELSE IF (LESSOK) THEN ARGLES=ATAN2F(YLESS,XLESS) ARGMOR=ARGLES+3.14159 IF (ARGMOR.GT.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 ENDIF C ******* 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.LT.0.) THEN ARGCUT=ARGCUT+3.14159 IF (ARGCUT.GT.3.14159) ARGCUT=ARGCUT-6.28318 ENDIF C *************************************** ENDIF C C SIZE OF BANDS TO BE USED AS LEVELS IN SPHERICAL NUMBERING C STERAD=4.*3.14159/NUMEL DARC=0.30*SQRT(STERAD) C C NUMBER THE FIRST NODE (OR A NEARBY REAL SUBSTITUTE) C IF (LIST(M).LE.NREALN) THEN ISTART=LIST(M) ELSE ISTART=1 FOUND=.FALSE. DO 215 I=1,NUMEL DO 214 K=1,3 IF (NODES(K,I).EQ.LIST(M)) THEN DO 213 J=1,3 IF (NODES(J,I).LE.NREALN) THEN FOUND=.TRUE. ISTART=NODES(J,I) ENDIF 213 CONTINUE ENDIF 214 CONTINUE 215 CONTINUE IF (.NOT.FOUND) THEN DO 225 I=1,NFL DO 224 K=1,4 IF (NODEF(K,I).EQ.LIST(M)) THEN DO 223 J=1,4 IF (NODEF(J,I).LE.NREALN) THEN FOUND=.TRUE. ISTART=NODEF(J,I) ENDIF 223 CONTINUE ENDIF 224 CONTINUE 225 CONTINUE ENDIF ENDIF IALIAS(ISTART)=1 IUSER(1)=ISTART NDONE=1 LAST1=1 LAST9=1 LEVEL=0 ARCLIM=0.0 C C----------------- BEGIN INDEFINATE LOOP ON LEVEL -------------------- C C MARK THE NODES TO BE NUMBERED IN THE NEXT LEVEL (UNNUMBERED NODES C CONNECTED TO NODES THAT HAVE ALREADY BEEN NUMBERED) C 240 LEVEL=LEVEL+1 ARCLIM=ARCLIM+DARC DO 250 I=1,NREALN C (NECESSARY INITIALIZATION) DONEXT(I)=.FALSE. 250 CONTINUE IF (SPHERE) THEN C LEVEL IS DEFINED BY DISTANCE FROM INITIAL NODE DO 252 I=1,NUMNOD IF (IALIAS(I).EQ.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.LE.ARCLIM) ENDIF 252 CONTINUE ELSE C LEVEL IS DEFINED BY CONNECTIVITY DO 300 I=LAST1,LAST9 NODE=IUSER(I) IF (NODE.GT.1) THEN DO 260 J=1,NODE-1 DONEXT(J)=DONEXT(J).OR. + ( (IALIAS(J).EQ.0) .AND. LGRAPH(INDEXN(J,NODE)) ) 260 CONTINUE ENDIF IF (NODE.LT.NREALN) THEN DO 270 J=NODE+1,NREALN DONEXT(J)=DONEXT(J).OR. + ( (IALIAS(J).EQ.0) .AND. LGRAPH(INDEXN(J,NODE)) ) 270 CONTINUE ENDIF 300 CONTINUE ENDIF C C HOW MANY ARE TO BE NUMBERED AT THIS LEVEL? C NTODO=0 DO 305 I=1,NREALN IF (DONEXT(I)) NTODO=NTODO+1 305 CONTINUE CCC WRITE (IUNITT,308) NTODO,LEVEL CC308 FORMAT (' ',I6,' NODES IN LEVEL ',I4) IF (NTODO.EQ.0) THEN IF (NDONE.EQ.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 ENDIF ENDIF C$$$ C C LOOP ON NODES IN LIST, TAKING THEM IN ORDER OF THEIR ARGUMENT C ANGLES FROM THE INITIAL NODE. C DO 400 N=1,NTODO C C FIND LOWEST-ARGUMENT NODE LEFT TO BE DONE C 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.LT.ARGCUT) ARG=ARG+6.28318 IF (ARG.LT.ARGMIN) THEN ARGMIN=ARG NODE=I ENDIF ENDIF 350 CONTINUE C C NUMBER IT ! ! ! C NDONE=NDONE+1 IUSER(NDONE)=NODE IALIAS(NODE)=NDONE DONEXT(NODE)=.FALSE. 400 CONTINUE C C SUMMARIZE THE LAST LEVEL DONE C LAST1=LAST9+1 LAST9=LAST1+NTODO-1 C C LOOP TO NEXT LEVEL C IF (NDONE.LT.NREALN) GO TO 240 C C------------------- END INDEFINATE LOOP ON LEVEL ----------------- C C DETERMINE THE NEW HALF-BANDWIDTH C 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.EQ.M9).AND.(IALDIF.EQ.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.') ENDIF ENDIF 9010 CONTINUE 9020 CONTINUE MAXLST(M)=MINDIF IF (ISTART.EQ.LIST(M)) THEN WRITE (IUNITT,9031) M, LIST(M),MINDIF ELSE WRITE (IUNITT,9032) M, LIST(M),ISTART,MINDIF ENDIF 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,'.') C C 10000 CONTINUE C================================================================= C C REPEAT THE NUMBERING FROM THE BEST STARTING POINT C IF (M9.GT.M1) THEN MINMAX=NREALN DO 10001 M=1,LENGTH IF (MAXLST(M).LT.MINMAX) THEN MINMAX=MAXLST(M) MSTART=M ENDIF 10001 CONTINUE M1=MSTART M9=MSTART GO TO 200 ENDIF C IF (NUMNOD.GT.NREALN) THEN DO 90000 I=NREALN+1,NUMNOD IALIAS(I)=I IUSER(I)=I 90000 CONTINUE ENDIF C RETURN END C C C SUBROUTINE CLOCK (INPUT,IUNITT,MODIFY,TINSEC) C C REPORT ELAPSED TIME BETWEEN CALLS C TOLD=TINSEC CALL SYSTEM_CLOCK(NCOUNT,NPERS,NCMAX) IF (NPERS.GT.0) THEN TINSEC=NCOUNT/NPERS ELSE TINSEC=0. END IF ELAPS=TINSEC-TOLD WRITE (IUNITT,1) ELAPS 1 FORMAT (/' ',F10.2,' SECONDS HAVE ELAPSED SINCE LAST CHECK.') RETURN END C C C SUBROUTINE PAUSE() IMPLICIT NONE WRITE(*,"(' Press [Enter]...')") READ(*,*) RETURN END