//IRAAGPBF JOB TIME=(0,10) COMPILE GRID-EDITING PROGRAM INTO LOAD MOD. /*SCHEDULE PRIORITY=0.8 //CLEAR EXEC PGM=IEFBR14 FIRST DELETE ANY OLD MODULES OF THIS //DD1 DD DSN=EFF9GPB.FIXERMOD,DISP=(OLD,DELETE) NAME. //FORT EXEC PGM=FORTVS,REGION=3072K, RUN COMPILER: // PARM='MAP,NODECK,NOLIST' VS-FORTRAN 2.4.0+ //STEPLIB DD DISP=(SHR,PASS),DSN=APP1.FORTVS.COMPILER //SYSIN DD * C PROGRAM FIX C (UCLA EDITION OF MAY, 1990) C ***** COMPATIBLE WITH MAY, 1990 VERSION OF LARAMY ************ C TAKES OUTPUT FROM A FINITE ELEMENT SIMULATION OF CONTINENTAL C DEFORMATION PERFORMED BY "LARAMY" AND PLOTS CONTOUR DIAGRAMS C IN COLOR ON A 3279 TERMINAL. ALLOWS EDITING OF VALUES AT NODES, C AND AUTOMATICALLY INTERPOLATES CHANGES TO INTEGRATION POINTS. C FINALLY, REWRITES A NEW "TAPE" DATASET FOR USE IN FURTHER RUNS. C C PARALLELS PROGRAM "GDDMCOMP", WHICH IS FOR BATCH USE W/O EDITING. C C USES STRATEGIC AND TACTICAL INPUT PARAMETERS IN C CARD FORMAT FROM DEVICE 4; SHOULD CONFORM TO DATA USED C IN THE ORIGINAL RUN OF "LARAMY"; PLOT CONTROLS APPENDED C AT THE END OF THIS DATASET ARE NOT USED BY THIS PROGRAM. C READS OLD OUTPUT "TAPE" AS SOURCE OF DETAILED DATA FROM DEVICE 8. C PRODUCES ONE OR MORE CORRECTED TIMESTEP FILES ON DEVICE 9. C READS STATE OUTLINES FROM UNIT 11 AND OPTIONALLY INCLUDES IN PLOTS. C C NOTICE: THIS PROGRAM AND ASSOCIATED SUBPROGRAMS WERE CREATED BY C PETER BIRD, DEPARTMENT OF EARTH & SPACE SCIENCES, C UNIVERSITY OF CALIFORNIA, LOS ANGELES. C FIVE YEARS OF SUPPORT FROM THE CRUSTAL STRUCTURE AND TECTONICS C PROGRAM OF THE NATIONAL SCIENCE FOUNDATION ARE GRATEFULLY C ACKNOWLEDGED. C THIS PROGRAM IS PUBLIC PROPERTY AND MAY BE REPRODUCED AND RUN C WITHOUT WRITTEN PERMISSION. HOWEVER, PROPER CREDIT SHOULD C BE GIVEN TO THE AUTHOR IN ANY RESULTING PUBLICATIONS. C USERS ARE ENCOURAGED TO CONTACT THE AUTHOR FOR ADVICE, UPDATES, C AND TECHNICAL SUPPORT, AT (213) 825-1126. C C EXTERNAL ROUTINES USED: C -ESSL (IBM PRODUCT) ROUTINES: C DGBF AND DGBS (DOUBLE PRECISION FACTOR AND C SOLVE A GENERAL BANDED LINEAR SYSTEM) C -FORTRAN INTRINSIC FUNCTIONS CALLED BY GENERIC NAMES: C ABS,ASIN,ATAN2,COS,EXP,LOG,MAX,MIN,MOD,SIN,SQRT,TAN C -IBM GRAPHICAL DATA DISPLAY MANAGER (GDDM) ROUTINES: C FSINIT,FSEXIT,GSCLP,GSSATI,GSSSEG,GSLSS,GSUWIN, C GSSCLS,GSCOL,GSLW,GSLT,GSLINE,GSMOVE,GSAREA,GSENDA, C GSARC,GSCHAR,GSCHAP,GSCM,GSCS,GSQCB,GSCB,GSCA, C GSPAT,GSSAVE,FSPCLR,FSTERM C CHARACTER*80 TITLE CHARACTER*8 ASTER,BLANKS CHARACTER*16 ALPHA,BETA CHARACTER*1 Z INTEGER TSHIFT,VSHIFT DOUBLE PRECISION CODE,FLOWIN LOGICAL ALDONE,ALLREP,BOXIT,DIMERR,DOFEM,DOREP, + DRAWST,FAILUR,LISTOP,RESTRT,RETRO, + STATES,TAPE9 LOGICAL BAR,FAILED,IMAGE,INVENT,MENU,TITLES,SHONUM LOGICAL ISTHER,ANYDIS,ANYEDS,NEWVAR,NEWTIM,CHANGE,ODDCOL,ODDROW COMMON /VTMAP/ ISTHER COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE C C FOLLOWING LINE SETS MAXIMUM NUMBER OF ELEMENTS: PARAMETER (N50=280) C FOLLOWING LINE SETS MAXIMUM NUMBER OF NODES: PARAMETER (N121=609) C FOLLOWING LINE SETS SIZE OF STORAGE RESERVED FOR CENTRAL BAND C OF A MATRIX NOMINALLY N121 BY N121: PARAMETER (N9922=86478) C FOLLOWING LINE RESERVES EXTRA STORAGE FOR LINEAR SYSTEM SOLUTIONS PARAMETER (NEXTRA=609) C FOLLOWING LINE SETS MAXIMUM NUMBER OF POINTS IN STATELINE MAPS: PARAMETER (NSTATE=2000) C FOLLOWING LINE SETS NUMBER OF COLORS IN PALETTE: PARAMETER (NCOLOR=13) C C VARIABLE DIMENSIONS CONTAINING VALUE OF N50: DIMENSION 2 AREAC(N50),AREAM(N50),CONINT(7,N50), 3 DETJC(7,N50), 4 DETJM(7,N50), 5 DXSC(6,7,N50), 6 DXSM(6,7,N50),DYSC(6,7,N50),DYSM(6,7,N50), 7 ESUMC(2,2,7,N50),ESUMM(2,2,7,N50),GEOTHA(4,7,N50), 8 GEOTHC(4,7,N50),GEOTHM(4,7,N50),LISTOP(N50), 9 NODES(6,0:N50),OUTSCA(7,N50),OUTVEC(2,7,N50), A OUTV2(2,7,N50), 1 THIKC(7,N50),THIKM(7,N50), 2 XIPC(7,N50),XIPM(7,N50),YIPC(7,N50),YIPM(7,N50) C VARIABLE DIMENSIONS CONTAINING VALUE OF N121: DIMENSION CHANGE(N121),CONDNS(N121),CONNOD(N121),FLOWIN(N121), + PHINOD(N121),XNODC(N121),XNODM(N121),YNODC(N121), + YNODM(N121),VC(2,N121),VM(2,N121), + THNKC(N121),THNKM(N121),VNODE(2,N121), + WC(N121),WM(N121),XTEMP(N121),YTEMP(N121) C VARIABLE DIMENSION CONTAINING VALUE OF N9922: DIMENSION CODE(N9922) C VARIABLE DIMENSION CONTAINING VALUE OF NEXTRA: DIMENSION LWORK(NEXTRA) C VARIABLE DIMENSION CONTAINING VALUE OF NSTATE: DIMENSION DRAWST(NSTATE),STLINK(3,NSTATE), + XST(NSTATE),YST(NSTATE) C C FIXED-DIMENSION ARRAYS OF GENERAL USE AND VARIABLE VALUE: DIMENSION ACREEP(3),ALPHAT(2),BCREEP(3),CCREEP(3),CONDUC(2), + CINT(24),DCREEP(3), + DIFFUS(2),DVPBYE(2,2),DVPDT(2),ECREEP(3),FBLAND(24), + FRIC(2),HMAX(2),HMIN(2), + LOWBLU(24),RADIO(2),RHOBAR(2),TEMLIM(2), + THICKN(2),VPMEAN(2) DIMENSION ISTHER(24,20) C (FIXED-DIMENSION ARRAYS OF LOCAL USE ARE BUILT INTO SUBROUTINES; C FIXED-DIMENSION ARRAYS OF GENERAL USE BUT CONSTANT VALUE ARE C CONTAINED IN COMMON BLOCKS.) DATA (NODES(J,0),J=1,6)/1,1,1,1,1,1/ DATA RINKM /6371./ DATA CINT /24*0.0/, LOWBLU /24*1/, FBLAND /24*0.0/ DATA RETRO /.FALSE./ DATA NUMVAR /10/ C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C DO 20 I=1,24 DO 10 J=1,20 ISTHER(I,J)=.FALSE. 10 CONTINUE 20 CONTINUE C C INITIALIZE GDDM C WRITE(6,30) 30 FORMAT(' INITIALIZING GDDM....') C CALL FSINIT C C DEFAULT INITIALIZATION SECTION C XCENTR=0.0 YCENTR=0.0 WIDE=6.00E8 IMAGE=.TRUE. INVENT=.TRUE. MENU=.TRUE. TITLES=.TRUE. BAR=.TRUE. STATES=.TRUE. SHONUM=.TRUE. IVAR=1 IVLAST=0 ITLAST=0 VSHIFT=0 TSHIFT=0 ANYDIS=.TRUE. ANYEDS=.FALSE. NCONTR=10 RMSVEC=0.36 IBELOW=0 IFT=0 C C MODEL PARAMETER INPUT SECTION C WRITE(6,50) 50 FORMAT(' READING FT04F001=INFILE....') CALL READIN (TITLE ,FRIC ,ACREEP,ECREEP,BCREEP, + CCREEP,DCREEP,CONDUC,DIFFUS, + RADIO ,THICKN,TEMLIM,RHOBAR, + ALPHAT,VPMEAN,DVPDT ,DVPBYE, + RHOAST,RHOH2O,BIOT ,G ,RADIUS, + X0ELON,Y0NLAT,CPNLAT,IBELOW, + TSLAB0,SIGBOT,PUSHHO,ECLOG , + SLABSZ,PUSHUP,NELROW,NELCOL, + BEGAGE,DELTAT,ENDAGE,DXMAX ,DTHMAX, + RAMP ,NDIFUS,MAXITR,OKTOQT, + VISMAX,ETAMAX,HMIN ,HMAX , + ALLREP,MIDREP,TAPE9 ,RESTRT, + KTIME ,CONRAD,DQDTDA,TSURF , + APLANO,WANDES,VDECOL,OLDGRD, + GWIDE ,GHIGH ,GANGLE) WRITE(6,70) 70 FORMAT(' CHECKING ADEQUACY OF DIMENSIONED ARRAYS...') CALL SETDIM (N50,N121,N9922,NEXTRA,NXL, + NCDIM,NDIFF,NELROW,NELCOL,NUMNOD, + NUMEL,DIMERR) IF (DIMERR) STOP C C PREPARATION OF NECESSARY ARRAYS C WRITE(6,80) 80 FORMAT(' READING IN FT11F001=STATE OUTLINES...') ONEKM=RADIUS/RINKM RTAN=RADIUS*TANDEG(90.-CPNLAT) YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) NXYST=0 DO 100 I=1,NSTATE READ(11,*,END=101) PLAT,PLON,DRAWST(I) CALL LLTOXY (INPUT,CPNLAT, + PLAT,PLON, + RADIUS,RTAN,X0ELON,YPOLE, + OUTPUT,X,Y) XST(I)=X YST(I)=Y NXYST=NXYST+1 100 CONTINUE 101 NXYST=NXYST-MOD(NXYST,7) DRAWST(1)=.FALSE. CALL GRIDDR(INPUT,NELROW,NELCOL,NUMEL, + MODIFY,NODES) NCOLN=2*NELCOL+1 NROWN=NUMNOD/NCOLN IR=NROWN/2 JC=1 ILIGHT=(IR-1)*NCOLN+JC C C FIND NUMBER OF TIMESTEPS ON INPUT TAPE C WRITE(6,120) 120 FORMAT(' READING IN FT08F001=TAPE DATASET...') ITIME=999 CALL GOON (ITIME,NTIME, + XNODC,XNODM,YNODC,YNODM,THNKC,CONNOD,THNKM, + GEOTHC,GEOTHM,GEOTHA,VC,VM,WC,WM,ESUMC,ESUMM, + NUMNOD,NUMEL,TIME,ALDONE) ITIME=NTIME DO 200 IR=1,NUMVAR DO 190 JC=1,NTIME ISTHER(IR,JC)=.TRUE. 190 CONTINUE 200 CONTINUE C C C**** BEGINNING OF ENDLESS LOOP ON IMAGES (UNTIL 'Q'=EXIT) ********* 300 CONTINUE C******************************************************************** C C BEGIN A NEW PAGE C CALL FSPCLR C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C FAILED=((ITIME+TSHIFT).GT.NTIME).OR. + ((ITIME+TSHIFT).LT.1) .OR. + ((IVAR+VSHIFT).GT.NUMVAR).OR. + ((IVAR+VSHIFT).LT.1) IF (FAILED) THEN CALL FSALRM INVENT=.TRUE. CALL STOCK(IVAR,ITIME) CALL ASDFLD(13,3,1,1,29,2) CALL ASCPUT(13,29,'YOU REQUESTED COMMAND-STRING:') CALL ASDFLD(14,3,31,1,15,2) CALL ASCPUT(14,15,ALPHA) CALL ASDFLD(15,4,1,1,50,2) CALL ASCPUT(15,50,'BUT THE REQUESTED TIME OR VARIABLE DOES NO &T EXIST.') CALL ASDFLD(16,5,1,1,26,2) CALL ASCPUT(16,26,'TRY A DIFFERENT COMMAND...') ELSE IVAR=IVAR+VSHIFT ITIME=ITIME+TSHIFT IF (INVENT) THEN CALL STOCK(IVAR,ITIME) ELSE NEWTIM=ITIME.NE.ITLAST IF (NEWTIM) THEN CALL PAST(ACREEP,ALPHAT,AREAC,AREAM, 2 BCREEP,CCREEP,CODE,CONDNS,CONINT,CONNOD,CONRAD, 3 DCREEP,DETJC,DETJM, 4 DXSC,DXSM,DYSC,DYSM,ECLOG,ECREEP, 5 ESUMC,ESUMM,ETAMAX,FLOWIN, 6 FRIC,G,GEOTHA,GEOTHC,GEOTHM, 7 ITIME,LISTOP,LWORK,NCDIM,NDIFF,NELCOL, 8 NODES,NUMEL,NUMNOD,NXL,ONEKM,OUTSCA,OUTVEC, 9 OUTV2,PUSHHO,PUSHUP,RADIUS, A RAMP,RHOAST,RHOH2O,RHOBAR,SIGBOT, 1 SLABSZ,CPNLAT,IBELOW,X0ELON,Y0NLAT, 2 TASTH,TEMLIM, 3 THIKC,THIKM,TIME,THNKC,THNKM, 4 TSLAB0,TSURF,VC,VM,VISMAX, 5 WC,WM,XIPC,XIPM, 6 XNODC,XNODM,YIPC,YIPM,YNODC,YNODM, 7 ALDONE,NELROW,WANDES, 8 HMAX,HMIN) ITLAST=ITIME ENDIF NEWVAR=NEWTIM.OR.(IVAR.NE.IVLAST) IF (NEWVAR) THEN DO 310 I=1,NUMNOD CHANGE(I)=.FALSE. 310 CONTINUE ENDIF CALL REPORT(ITIME,XIPC,XIPM,YIPC,YIPM, 2 XNODC,XNODM,YNODC,YNODM,TITLE,VM,NODES, 3 OUTSCA,OUTVEC,VC, 4 THIKM,THIKC,GEOTHA,GEOTHC,CONDUC, 5 GEOTHM,VPMEAN,DVPBYE,DVPDT, 6 TIME ,NUMNOD,NUMEL, 7 G,HMAX,HMIN,RHOAST,RHOH2O, 8 TEMLIM,ONEKM,ESUMC,ESUMM,AREAC,AREAM, 9 CODE,CONDNS,DETJC,DETJM,FAILUR,FLOWIN, A NCDIM,NDIFF,NXL,LWORK,WC,WM, 1 ECLOG,RADIUS,SLABSZ,CPNLAT,IBELOW, 2 X0ELON,Y0NLAT,OUTV2,RAMP, 3 THNKC,TASTH,TSLAB0, 4 ILIGHT,IVAR,NCONTR, 5 RMSVEC,NELCOL,PHINOD,DRAWST, 6 NXYST,XST,YST,NELROW,THNKM,FBLAND,LOWBLU, 7 CINT,WANDES,CONINT,CONNOD,CONRAD, 8 TSURF,PUSHUP,VNODE,NEWVAR,DFCON,XTEMP,YTEMP) IVLAST=IVAR ENDIF ENDIF IF (MENU) CALL CHOICE C C DEFINE DISPLAY-CONTROL-CHARACTER AREA AT LOWER RIGHT C WITH EDIT-CONTROL-CHARACTER AREA BELOW IT. C CALL ASDFLD(1,31,55,1,9,2) CALL ASCPUT(1,9,'DISPLAY: ') CALL ASDFLD(2,31,64,1,17,0) CALL ASCPUT(2,0,' ') CALL ASDFLD(3,32,55,1,9,2) CALL ASCPUT(3,9,' EDIT: ') CALL ASDFLD(4,32,64,1,17,0) CALL ASCPUT(4,0,' ') CALL ASFCUR(2,1,1) IF (ANYEDS) CALL ASFCUR(4,1,1) CALL ASREAD(ITYPE,IVALUE,ICOUNT) C ----------------------------------------------------------- C THIS IS WHERE THE DISPLAY IS INSPECTED AND CHANGES CAN BE C ORDERED IN THE DISPLAY FIELD (WHICH HAS PRIORITY), OR C IN THE EDIT FIELD (IF THE DISPLAY FIELD IS LEFT BLANK). C ----------------------------------------------------------- C C CLEAR MESSAGE FIELDS (IF ANY) C CALL ASDFLD(13,0,0,0,0,0) CALL ASDFLD(14,0,0,0,0,0) CALL ASDFLD(15,0,0,0,0,0) CALL ASDFLD(16,0,0,0,0,0) C C DISPLAY-COMMAND INTERPRETER SECTION C TSHIFT=0 VSHIFT=0 JTENS=0 ANYDIS=.FALSE. CALL ASCGET(2,16,ALPHA) DO 5000 ICHR=1,15 Z=ALPHA(ICHR:ICHR) IF (Z.EQ.'M'.OR.Z.EQ.'m') THEN MENU=.NOT.MENU JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'I'.OR.Z.EQ.'i') THEN INVENT=.NOT.INVENT JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'L'.OR.Z.EQ.'l') THEN STATES=.NOT.STATES JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'G'.OR.Z.EQ.'g') THEN BAR=.NOT.BAR JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'T'.OR.Z.EQ.'t') THEN TITLES=.NOT.TITLES JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'X'.OR.Z.EQ.'x') THEN WIDE=WIDE*0.8 JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'C'.OR.Z.EQ.'c') THEN WIDE=WIDE/0.8 JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'N'.OR.Z.EQ.'n') THEN YCENTR=YCENTR+WIDE/11. JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'S'.OR.Z.EQ.'s') THEN YCENTR=YCENTR-WIDE/11. JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'E'.OR.Z.EQ.'e') THEN XCENTR=XCENTR+WIDE/11. JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'W'.OR.Z.EQ.'w') THEN XCENTR=XCENTR-WIDE/11. JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'F'.OR.Z.EQ.'f') THEN TSHIFT=TSHIFT+1 JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'B'.OR.Z.EQ.'b') THEN TSHIFT=TSHIFT-1 JTENS=0 ANYDIS=.TRUE. ELSE IF (ICHAR(Z).GE.240.AND.ICHAR(Z).LE.249) THEN J=ICHAR(Z)-240+JTENS*10 TSHIFT=J-ITIME JTENS=J ANYDIS=.TRUE. ELSE IF (Z.EQ.'V'.OR.Z.EQ.'v') THEN VSHIFT=VSHIFT+1 JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'U'.OR.Z.EQ.'u') THEN VSHIFT=VSHIFT-1 JTENS=0 ANYDIS=.TRUE. ELSE IF (Z.EQ.'Q'.OR.Z.EQ.'q') THEN ANYDIS=.TRUE. GO TO 9999 ELSE JTENS=0 GO TO 5001 ENDIF 5000 CONTINUE C C EDIT-COMMAND INTERPRETER SECTION C 5001 ANYEDS=.FALSE. CALL ASCGET(4,16,BETA) IPOINT=1 5002 Z=BETA(IPOINT:IPOINT) IPOINT=IPOINT+1 IF (Z.EQ.'I'.OR.Z.EQ.'i') THEN ANYEDS=.TRUE. SHONUM=.NOT.SHONUM IF (IPOINT.LE.15) GO TO 5002 ELSE IF (Z.EQ.'H'.OR.Z.EQ.'h') THEN ANYEDS=.TRUE. IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) IF (X.GT.0.) THEN ILIGHT=X+0.1 ELSE Z=BETA(IPOINT:IPOINT) IPOINT=IPOINT+1 IF (Z.EQ.'E'.OR.Z.EQ.'e'.OR. + Z.EQ.'R'.OR.Z.EQ.'r') THEN IR=(ILIGHT-1)/NCOLN+1 JC=MOD((ILIGHT-1),NCOLN)+1 JC=MIN(JC+1,NCOLN) ILIGHT=(IR-1)*NCOLN+JC ELSE IF (Z.EQ.'W'.OR.Z.EQ.'w'.OR. + Z.EQ.'L'.OR.Z.EQ.'l') THEN IR=(ILIGHT-1)/NCOLN+1 JC=MOD((ILIGHT-1),NCOLN)+1 JC=MAX(JC-1,1) ILIGHT=(IR-1)*NCOLN+JC ELSE IF (Z.EQ.'N'.OR.Z.EQ.'n'.OR. + Z.EQ.'U'.OR.Z.EQ.'u') THEN IR=(ILIGHT-1)/NCOLN+1 JC=MOD((ILIGHT-1),NCOLN)+1 IR=MAX(IR-1,1) ILIGHT=(IR-1)*NCOLN+JC ELSE IF (Z.EQ.'S'.OR.Z.EQ.'s'.OR. + Z.EQ.'D'.OR.Z.EQ.'d') THEN IR=(ILIGHT-1)/NCOLN+1 JC=MOD((ILIGHT-1),NCOLN)+1 IR=MIN(IR+1,NROWN) ILIGHT=(IR-1)*NCOLN+JC ENDIF ENDIF IF (IPOINT.LE.15) GO TO 5002 ENDIF ELSE IF (Z.EQ.'Z'.OR.Z.EQ.'z') THEN IR=(ILIGHT-1)/NCOLN+1 JC=MOD((ILIGHT-1),NCOLN)+1 ODDROW=MOD(IR,2).EQ.1 ODDCOL=MOD(JC,2).EQ.1 IF (ODDROW) THEN IF (ODDCOL) THEN C -CORNER NODES CANNOT BE INTERPOLATED- ELSE VNODE(1,ILIGHT)=0.5*(VNODE(1,ILIGHT-1)+ + VNODE(1,ILIGHT+1)) VNODE(2,ILIGHT)=0.5*(VNODE(2,ILIGHT-1)+ + VNODE(2,ILIGHT+1)) CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES,NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ENDIF ELSE IF (ODDCOL) THEN VNODE(1,ILIGHT)=0.5*(VNODE(1,ILIGHT-NCOLN)+ + VNODE(1,ILIGHT+NCOLN)) VNODE(2,ILIGHT)=0.5*(VNODE(2,ILIGHT-NCOLN)+ + VNODE(2,ILIGHT+NCOLN)) CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES,NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ELSE VNODE(1,ILIGHT)=0.5*(VNODE(1,ILIGHT-NCOLN+1)+ + VNODE(1,ILIGHT+NCOLN-1)) VNODE(2,ILIGHT)=0.5*(VNODE(2,ILIGHT-NCOLN+1)+ + VNODE(2,ILIGHT+NCOLN-1)) CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES,NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ENDIF ENDIF IF (IPOINT.LE.15) GO TO 5002 ELSE IF (Z.EQ.'S'.OR.Z.EQ.'s') THEN ANYEDS=.TRUE. IF (IVAR.EQ.3.OR.IVAR.EQ.4.OR.IVAR.EQ.5.OR.IVAR.EQ.6.OR. + IVAR.EQ.9.OR.IVAR.EQ.10) THEN IF (IPOINT.LE.15) THEN Z=BETA(IPOINT:IPOINT) IPOINT=IPOINT+1 IF (Z.EQ.'I'.OR.Z.EQ.'i'.OR. + Z.EQ.'U'.OR.Z.EQ.'u') THEN IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) DF=X*DFCON CONDNS(ILIGHT)=CONDNS(ILIGHT)+DF CHANGE(ILIGHT)=.TRUE. ENDIF ELSE IF (Z.EQ.'D'.OR.Z.EQ.'d') THEN IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) DF=X*DFCON CONDNS(ILIGHT)=CONDNS(ILIGHT)-DF CHANGE(ILIGHT)=.TRUE. ENDIF ELSE IF (IPOINT.EQ.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) DF=X*DFCON CONDNS(ILIGHT)=CONDNS(ILIGHT)+DF CHANGE(ILIGHT)=.TRUE. ENDIF ENDIF ENDIF ELSE CALL FSALRM ENDIF ELSE IF (Z.EQ.'V'.OR.Z.EQ.'v') THEN ANYEDS=.TRUE. IF (IVAR.EQ.1.OR.IVAR.EQ.2.OR.IVAR.EQ.7.OR.IVAR.EQ.8) THEN IF (IPOINT.LE.15) THEN Z=BETA(IPOINT:IPOINT) IPOINT=IPOINT+1 IF (Z.EQ.'n'.OR.Z.EQ.'N'.OR. + Z.EQ.'u'.OR.Z.EQ.'U') THEN IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) IF (IVAR.LE.2) THEN DF=X*WIDE/11. ELSE DF=X/VFACT ENDIF VNODE(2,ILIGHT)=VNODE(2,ILIGHT)+DF CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES, + NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ENDIF ELSE IF (Z.EQ.'S'.OR.Z.EQ.'s'.OR. + Z.EQ.'D'.OR.Z.EQ.'d') THEN IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) IF (IVAR.LE.2) THEN DF=X*WIDE/11. ELSE DF=X/VFACT ENDIF VNODE(2,ILIGHT)=VNODE(2,ILIGHT)-DF CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES, + NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ENDIF ELSE IF (Z.EQ.'E'.OR.Z.EQ.'e'.OR. + Z.EQ.'R'.OR.Z.EQ.'r') THEN IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) IF (IVAR.LE.2) THEN DF=X*WIDE/11. ELSE DF=X/VFACT ENDIF VNODE(1,ILIGHT)=VNODE(1,ILIGHT)+DF CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES, + NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ENDIF ELSE IF (Z.EQ.'W'.OR.Z.EQ.'w'.OR. + Z.EQ.'L'.OR.Z.EQ.'l') THEN IF (IPOINT.LE.15) THEN CALL NUMBER (INDATA,15,BETA, + MODIFY,IPOINT, + OUTPUT,X) IF (IVAR.LE.2) THEN DF=X*WIDE/11. ELSE DF=X/VFACT ENDIF VNODE(1,ILIGHT)=VNODE(1,ILIGHT)-DF CHANGE(ILIGHT)=.TRUE. CALL FIXDOT (INPUT,ILIGHT,NODES, + NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) ENDIF ENDIF IF (IVAR.LE.2) THEN CALL FLOW (VNODE,NUMNOD,NODES,NUMEL,OUTVEC) CALL DETJUP (INDATA,NUMEL,NUMNOD,NODES,VNODE, + OUTPUT,OUTSCA) ENDIF IF (IPOINT.LE.15) GO TO 5002 ENDIF ELSE CALL FSALRM ENDIF ELSE IF (Z.EQ.'F'.OR.Z.EQ.'f') THEN ANYEDS=.TRUE. CALL FIXER(INDATA,CHANGE,CONDNS,CONDUC, + HMAX,HMIN,IVAR,NODES, + NUMEL,NUMNOD,ONEKM,OUTSCA, + OUTVEC,VNODE, + MODIFY,AREAC,AREAM, + DETJC,DETJM,DXSC,DXSM,DYSC,DYSM, + ESUMC,ESUMM,GEOTHC,GEOTHM, + LISTOP,NUMBAD, + THIKC,THIKM,THNKC,THNKM,VC,VM, + XIPC,XIPM,XNODC,XNODM,YIPC,YIPM, + YNODC,YNODM) ELSE IF (Z.EQ.'W'.OR.Z.EQ.'w') THEN ANYEDS=.TRUE. CALL TAPE (TITLE,TIME, + XNODC,XNODM,YNODC,YNODM,THNKC,CONNOD,THNKM, + GEOTHC,GEOTHM,GEOTHA,VC,VM,WC,WM,ESUMC,ESUMM, + NUMNOD,NUMEL) IFT=IFT+1 WRITE(6,5990)IFT 5990 FORMAT(' OUTPUT DATASET WRITTEN INTO FT09F00',I1) ENDIF C C**** CONCLUSION OF ENDLESS LOOP ON IMAGES (UNTIL 'Q'=EXIT) ********* 6001 GO TO 300 C******************************************************************** C 9999 CALL FSTERM STOP END C C C SUBROUTINE LLTOXY (INPUT,CPNLAT, + PLAT,PLON, + RADIUS,RTAN,X0ELON,YPOLE, + OUTPUT,X,Y) C C CONVERT A (NORTH LATITUDE=PLAT, EAST LONGITUDE=PLON) POSITION C INTO AN (X,Y) POSITION ON A CONIC PROJECTION WITH TANGENT C LATITUDE CPNLAT, WHEN THE (X,Y) ORIGIN IS AT C (NORTH LATITUDE=Y0NLAT, EAST LONGITUDE=X0ELON). C THE CUT NECESSARY IN THIS PROJECTION IS FROM THE POLE NEAREST C TO THE TANGENT LATITUDE (CPNLAT), ALONG A MERIDIAN WHICH C IS ON THE OPPOSITE SIDE OF THE EARTH FROM X0ELON. C IF PLAT IS MORE THAN 90 DEGREES DIFFERENT FROM CPNLAT, THE C POINT DOES NOT FALL ONTO THE PROJECTION AT ALL. TO PREVENT C CRASHES, WE MERELY PLACE IT VERY FAR OUT ON THE PROJECTION. C C NOTE: FOLLOWING TWO LINES ARE PRECOMPUTED AND PASSED TO SAVE TIME: C RTAN=RADIUS*TANDEG(90.-CPNLAT) C YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) C C STATEMENT FUNCTIONS: SINDEG(S)=SIN(S*0.0174533) COSDEG(S)=COS(S*0.0174533) TANDEG(S)=TAN(S*0.0174533) C IF (ABS(PLAT-CPNLAT).GE.90.) PLAT=CPNLAT+89.*(PLAT-CPNLAT)/ + ABS(PLAT-CPNLAT) R=RTAN-RADIUS*TANDEG(PLAT-CPNLAT) DPLON=PLON-X0ELON IF (DPLON.LT.-180.) DPLON=DPLON+360. IF (DPLON.GT. 180.) DPLON=DPLON-360. ANGLE=DPLON*SINDEG(CPNLAT) X=R*SINDEG(ANGLE) Y=YPOLE-R*COSDEG(ANGLE) RETURN END C C C SUBROUTINE STOCK (IVAR,ITIME) C C USE FREE SEGMENTS #1-4 (AT MOST) TO CREATE A GRAPH OF THE C INVENTORY OF AVAILABLE IMAGES, BY VARIABLE AND TIMESTEP C NUMBERS. LEAVE SPACE AT TOP AND BOTTOM FOR OVERLAY OF C TITLES FROM CURRENT IMAGE. C NOTE THAT ENTIRE IMAGE MAY BE OVERLAIN BY MENU. C CHARACTER*5 ASCII,LABEL CHARACTER*42 TEXT EXTERNAL ASCII LOGICAL ISTHER COMMON /VTMAP/ ISTHER DIMENSION ISTHER(24,20),NVCHAR(24),TEXT(24) C DATA NUMVAR /10/ DATA TEXT(1)/'MANTLE:GRID OF ELEMENTS '/ DATA NVCHAR(1)/23/ DATA TEXT(2)/'CRUST: GRID OF ELEMENTS '/ DATA NVCHAR(2)/23/ DATA TEXT(3)/'MANTLE: THICKNESS '/ DATA NVCHAR(3)/17/ DATA TEXT(4)/'CRUST: THICKNESS '/ DATA NVCHAR(4)/17/ DATA TEXT(5)/'MANTLE: HEAT-FLOW '/ DATA NVCHAR(5)/17/ DATA TEXT(6)/'CRUST: HEAT-FLOW '/ DATA NVCHAR(6)/17/ DATA TEXT(7)/'MANTLE: VELOCITY '/ DATA NVCHAR(7)/16/ DATA TEXT(8)/'CRUST: VELOCITY '/ DATA NVCHAR(8)/31/ DATA TEXT(9)/'MANTLE: DILATATION '/ DATA NVCHAR(9)/18/ DATA TEXT(10)/'CRUST: DILATATION '/ DATA NVCHAR(10)/18/ C C C********************************************************************** C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C INITIALIZE SEGMENT 1 (RED BACKGROUND RECTANGLE) C CALL GSSEG(1) C C SET COLOR TO RED C CALL GSCOL(2) C C OPEN AREA WITH INVISIBLE BOUNDARY C CALL GSAREA(0) C C DRAW OUTLINE OF AREA C CALL GSMOVE(5.0,7.0) CALL GSLINE(10.0,7.0) CALL GSLINE(10.0,1.0) CALL GSLINE(5.0,1.0) CALL GSLINE(5.0,7.0) C C SIGNAL END OF AREA C CALL GSENDA C C CLOSE SEGMENT C CALL GSSCLS C C************************************************************* C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT 2 (TEXT LABELS AND RULED LINES) C CALL GSSEG(2) C C SELECT AND LOAD CHARACTER SET C C CALL GSLSS(2,'ADMUWCRP',199) C C SELECT CHARACTER MODE (2 OR 3) C 2 = DEFAULT TYPE, RECTILINEAR, CONSTANT-SPACED C 3 = VECTOR TYPE, ANY ANGLE, PROPORTIONAL SPACING C CALL GSCM(3) C C COMPLETE CONNECTION OF CHARACTER SET C C CALL GSCS(199) C C SET SIZE OF CHARACTER BOXES RELATIVE TO DEFAULT C CALL GSQCB(WIDTH,HEIGHT) WIDTH=1.0*WIDTH HEIGHT=1.0*HEIGHT CALL GSCB(WIDTH,HEIGHT) C C SET LABEL ANGLE TO ZER0 (LIKE THIS) C CALL GSCA(COS(0.0),SIN(0.0)) C C SET COLOR TO CYAN C CALL GSCOL(5) CALL GSLW(2) C CALL GSMOVE(1.0,7.5) CALL GSCHAP(19,'INVENTORY OF IMAGES') CALL GSMOVE(6.5,7.5) CALL GSCHAP(15,'TIMESTEP NUMBER') C C WRITE TIMESTEP NUMBERS AND VERTICAL RULES C DO 210 I=1,20 LABEL=ASCII(1.00*I) XL=5.00+(I-1)*0.25-WIDTH CALL GSMOVE(XL,7.0) CALL GSCHAP(3,LABEL) XR=5.00+I*0.25 CALL GSMOVE(XR,7.00) CALL GSLINE(XR,1.00) 210 CONTINUE C C WRITE VARIABLE NAMES AND HORIZONTAL RULES C DO 220 I=1,NUMVAR YB=7.00-I*0.25 CALL GSMOVE(1.50,YB) CALL GSCHAP(NVCHAR(I),TEXT(I)) CALL GSMOVE(4.9,YB) CALL GSLINE(10.0,YB) 220 CONTINUE C C CLOSE SEGMENT C CALL GSSCLS C C******************************************************************* C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT (3) OF GREEN BOXES SHOWING AVAILABLE FIGURES C CALL GSSEG(3) C C SET COLOR GREEN C CALL GSCOL(4) C C DRAW A BOX FOR EACH EXISTING FIGURE C DO 390 IR=1,24 DO 380 JC=1,20 IF (ISTHER(IR,JC)) THEN XL=5.00+(JC-1)*0.25 XR=XL+0.25 YB=7.00-IR*0.25 YT=YB+0.25 CALL GSAREA(0) CALL GSMOVE(XL,YB) CALL GSLINE(XR,YB) CALL GSLINE(XR,YT) CALL GSLINE(XL,YT) CALL GSLINE(XL,YB) CALL GSENDA ENDIF 380 CONTINUE 390 CONTINUE C C CLOSE SEGMENT C CALL GSSCLS C C**************************************************************** C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT 4 (MARK CURRENT IMAGE ON MAP) C CALL GSSEG(4) C C USE YELLOW PEN C CALL GSCOL(6) C C DRAW A DOT C IR=IVAR JC=ITIME XC=5.00+0.25*(JC-1)+0.125 YC=7.00-0.25*(IR-1)-0.125 XP=XC+0.100 CALL GSAREA(0) CALL GSMOVE(XP,YC) CALL GSARC(XC,YC,360.) CALL GSENDA C C CLOSE SEGMENT C CALL GSSCLS C C******************************************************************* C RETURN END C C C CHARACTER*5 FUNCTION ASCII (T) C C CONVERTS REAL*4 FLOATING-POINT VARIABLE TO "-12.3" C OR "123.4" OR " 0.9" C REAL*4 S,T S=ABS(T) IF (T.LE.-100.0) THEN ASCII='-**.*' ELSE IF (T.GE.1000.0) THEN ASCII='***.*' ELSE I1=INT(S/100.) ASCII(1:1)=CHAR(240+I1) I2=INT((S-100.*I1)/10.) ASCII(2:2)=CHAR(240+I2) I3=INT(S-100.*I1-10.*I2) ASCII(3:3)=CHAR(240+I3) ASCII(4:4)='.' I4=INT(10.*(S-100.*I1-10.*I2-I3)) ASCII(5:5)=CHAR(240+I4) IF(ASCII(1:1).EQ.'0')ASCII(1:1)=' ' IF(ASCII(2:2).EQ.'0'.AND.ASCII(1:1).EQ.' ')ASCII(2:2)=' ' IF (T.LT.0.0) THEN IF (ASCII(2:2).EQ.' ') THEN ASCII(2:2)='-' ELSE ASCII(1:1)='-' ENDIF ENDIF DO 10 I=2,5 IF (ASCII(I:I).EQ.'0') ASCII(I:I)='O' 10 CONTINUE ENDIF RETURN END C C C SUBROUTINE CHOICE C C WRITE MENU ON TOP OF ANY PRE-EXISTING SEGMENTS C C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C CALL GSSEG(7) C C SET COLOR TO YELLOW C CALL GSCOL(6) CALL GSLW(2) C C OPEN AREA WITH VISIBLE BOUNDARY C CALL GSAREA(1) C C DRAW OUTLINE OF AREA WITH RED PEN C CALL GSCOL(2) CALL GSMOVE(0.1,7.9) CALL GSLINE(5.9,7.9) CALL GSLINE(5.9,1.1) CALL GSLINE(0.1,1.1) CALL GSLINE(0.1,7.9) C C SIGNAL END OF AREA C CALL GSENDA C C SET COLOR TO YELLOW C CALL GSCOL(6) CALL GSLW(2) C C OPEN AREA WITH VISIBLE BOUNDARY C CALL GSAREA(1) C C DRAW OUTLINE OF AREA WITH RED PEN C CALL GSCOL(2) CALL GSMOVE(6.1,7.9) CALL GSLINE(11.1,7.9) CALL GSLINE(11.1,1.1) CALL GSLINE(6.1,1.1) CALL GSLINE(6.1,7.9) C C SIGNAL END OF AREA C CALL GSENDA C C CLOSE FIRST MENU SEGMENT (#6) AND BEGIN SECOND (#7) C CALL GSSCLS C C*************************************** C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C CALL GSSEG(8) C C WRITE OUT MENU OPTIONS C C SELECT AND LOAD CHARACTER SET C C CALL GSLSS(2,'ADMUWCRP',199) C C SELECT CHARACTER MODE (2 OR 3) C 2 = DEFAULT TYPE, RECTILINEAR, CONSTANT-SPACED C 3 = VECTOR TYPE, ANY ANGLE, PROPORTIONAL SPACING C CALL GSCM(3) C C COMPLETE CONNECTION OF CHARACTER SET C C CALL GSCS(199) C C SET SIZE OF CHARACTER BOXES RELATIVE TO DEFAULT C CALL GSQCB(WIDTH,HEIGHT) WIDTH=1.0*WIDTH HEIGHT=1.0*HEIGHT CALL GSCB(WIDTH,HEIGHT) C C SET LABEL ANGLE TO ZER0 (LIKE THIS) C CALL GSCA(COS(0.0),SIN(0.0)) C C USE BLACK PEN C CALL GSCOL(-1) CALL GSLW(2) C C TEXT OF DISPLAY-COMMAND MENU: CALL GSMOVE(0.6,7.5) CALL GSCHAP(20,'DISPLAY-CONTROL MENU') CALL GSMOVE(0.5,7.0) CALL GSCHAP(31,'M toggle this Menu on/off ') CALL GSMOVE(0.5,6.6) CALL GSCHAP(32,'I toggle Inventory screen on/off') CALL GSMOVE(0.5,6.3) CALL GSCHAP(31,'L toggle state Lines on/off ') CALL GSMOVE(0.5,6.0) CALL GSCHAP(31,'G toggle Graphical scale on/off') CALL GSMOVE(0.5,5.7) CALL GSCHAP(31,'T toggle Titles on/off ') CALL GSMOVE(0.5,5.4) CALL GSCHAP(31,'X eXpand map around center ') CALL GSMOVE(0.5,5.1) CALL GSCHAP(31,'C Contract map around center ') CALL GSMOVE(0.5,4.8) CALL GSCHAP(31,'N move window North by 10% ') CALL GSMOVE(0.5,4.5) CALL GSCHAP(31,'S move window South by 10% ') CALL GSMOVE(0.5,4.2) CALL GSCHAP(31,'E move window East by 10% ') CALL GSMOVE(0.5,3.9) CALL GSCHAP(31,'W move window West by 10% ') CALL GSMOVE(0.5,3.6) CALL GSCHAP(31,'F go Forward in time ') CALL GSMOVE(0.5,3.3) CALL GSCHAP(31,'B go Backwards in time ') CALL GSMOVE(0.1,3.0) CALL GSCHAP(35,'1..99 jump to time window #1,2,..99') CALL GSMOVE(0.5,2.7) CALL GSCHAP(33,'V change to next Variable on file') CALL GSMOVE(0.5,2.4) CALL GSCHAP(31,'U Undo or back up variable ') CALL GSMOVE(0.5,2.1) CALL GSCHAP(31,'Q Quit and exit this program ') CALL GSMOVE(0.2,1.3) CALL GSCHAP(35,'YOU MAY COMBINE UP TO 15 CHARACTERS') C C TEXT OF EDIT-COMMAND MENU: CALL GSMOVE(6.6,7.5) CALL GSCHAP(17,'EDIT-CONTROL MENU') CALL GSMOVE(6.5,7.1) CALL GSCHAP(31,'I toggle node Indeces on/off ') CALL GSMOVE(6.5,6.8) CALL GSCHAP(31,'H Highlight node ___ ') CALL GSMOVE(6.5,6.5) CALL GSCHAP(31,' 609 (give node number) or ') CALL GSMOVE(6.5,6.2) CALL GSCHAP(31,' U Up one from current node,or') CALL GSMOVE(6.5,5.9) CALL GSCHAP(31,' D Down one from current, or ') CALL GSMOVE(6.5,5.6) CALL GSCHAP(31,' R Right one from current, or ') CALL GSMOVE(6.5,5.3) CALL GSCHAP(31,' L Left one from current. ') CALL GSMOVE(6.5,5.0) CALL GSCHAP(31,'Z align midpoint node (ONLY) ') CALL GSMOVE(6.5,4.7) CALL GSCHAP(31,'S Scalar value adjustment: ') CALL GSMOVE(6.5,4.4) CALL GSCHAP(31,' I Increase, or ') CALL GSMOVE(6.5,4.1) CALL GSCHAP(31,' D Decrease ') CALL GSMOVE(6.5,3.8) CALL GSCHAP(31,' 1.5 (number of contours) ') CALL GSMOVE(6.5,3.5) CALL GSCHAP(31,'V Vector value adjustment: ') CALL GSMOVE(6.5,3.2) CALL GSCHAP(31,' N North, or ') CALL GSMOVE(6.5,2.9) CALL GSCHAP(31,' S South, or ') CALL GSMOVE(6.5,2.6) CALL GSCHAP(31,' E East, or ') CALL GSMOVE(6.5,2.3) CALL GSCHAP(31,' W West: ') CALL GSMOVE(6.5,2.05) CALL GSCHAP(31,' .15 (number of inches) ') CALL GSMOVE(6.5,1.8) CALL GSCHAP(31,'F FIX CHANGES INTO ARRAYS ') CALL GSMOVE(6.5,1.55) CALL GSCHAP(25,'W WRITE CORRECTED DATASET ') CALL GSMOVE(6.4,1.30) CALL GSCHAP(31,'YOU MAY COMBINE UP TO 15 CHAR.s') C C CLOSE SECOND MENU SEGMENT C CALL GSSCLS RETURN END C C C SUBROUTINE READIN(TITLE ,FRIC ,ACREEP,ECREEP,BCREEP, + CCREEP,DCREEP,CONDUC,DIFFUS, + RADIO ,THICKN,TEMLIM,RHOBAR, + ALPHAT,VPMEAN,DVPDT ,DVPBYE, + RHOAST,RHOH2O,BIOT ,G ,RADIUS, + X0ELON,Y0NLAT,CPNLAT,IBELOW, + TSLAB0,SIGBOT,PUSHHO,ECLOG , + SLABSZ,PUSHUP,NELROW,NELCOL, + BEGAGE,DELTAT,ENDAGE,DXMAX ,DTHMAX, + RAMP ,NDIFUS,MAXITR,OKTOQT, + VISMAX,ETAMAX,HMIN ,HMAX , + ALLREP,MIDREP,TAPE9 ,RESTRT, + KTIME ,CONRAD,DQDTDA,TSURF , + APLANO,WANDES,VDECOL,OLDGRD, + GWIDE ,GHIGH ,GANGLE) C C READS STRATEGIC AND TACTICAL INPUT PARAMETERS FROM DEVICE 5 C CHARACTER*80 TITLE LOGICAL ALLREP,OLDGRD,RESTRT,TAPE9 DIMENSION ACREEP(3),ALPHAT(2),BCREEP(3),CCREEP(3), + CONDUC(2),DCREEP(3),DIFFUS(2),DVPBYE(2,2), + DVPDT(2),ECREEP(3),FRIC(2),HMAX(2),HMIN(2), + RADIO(2),RHOBAR(2),TEMLIM(2),THICKN(2),VPMEAN(2) 1 FORMAT(A80) READ(4,*) READ(4,1) TITLE READ(4,*) READ(4,*) FRIC(1),FRIC(2) READ(4,*) ACREEP(1),ACREEP(3) READ(4,*) ACREEP(2) READ(4,*) ECREEP(1),ECREEP(3) READ(4,*) ECREEP(2) IF (ECREEP(2).NE.ECREEP(1)) THEN ECREEP(2)=ECREEP(1) ENDIF READ(4,*) BCREEP(1),BCREEP(3) READ(4,*) BCREEP(2) READ(4,*) CCREEP(1),CCREEP(3) READ(4,*) CCREEP(2) READ(4,*) DCREEP(1),DCREEP(3) READ(4,*) DCREEP(2) READ(4,*) CONDUC(1),CONDUC(2) READ(4,*) DIFFUS(1),DIFFUS(2) READ(4,*) RADIO(1),RADIO(2) READ(4,*) THICKN(1),THICKN(2) READ(4,*) TEMLIM(1),TEMLIM(2) READ(4,*)(RHOBAR(I),I=1,2) READ(4,*) ALPHAT(1),ALPHAT(2) READ(4,*) VPMEAN(1),VPMEAN(2) READ(4,*) DVPDT(1),DVPDT(2) READ(4,*) DVPBYE(1,1),DVPBYE(1,2) READ(4,*) DVPBYE(2,1),DVPBYE(2,2) READ(4,*) RHOAST READ(4,*) RHOH2O READ(4,*) BIOT BIOT=MAX(0.0,MIN(1.0,BIOT)) READ(4,*) G READ(4,*) RADIUS READ(4,*) X0ELON READ(4,*) Y0NLAT READ(4,*) CPNLAT IF (ABS(CPNLAT).LT.0.01) CPNLAT=0.01 READ(4,*) IBELOW READ(4,*) READ(4,*) TSURF READ(4,*) TSLAB0 READ(4,*) SIGBOT READ(4,*) WANDES READ(4,*) PUSHHO READ(4,*) ECLOG READ(4,*) SLABSZ READ(4,*) PUSHUP READ(4,*) READ(4,*) NELROW READ(4,*) NELCOL READ(4,*) BEGAGE READ(4,*) DELTAT READ(4,*) ENDAGE READ(4,*) DXMAX READ(4,*) DTHMAX READ(4,*) RAMP READ(4,*) NDIFUS READ(4,*) MAXITR READ(4,*) OKTOQT READ(4,*) VISMAX READ(4,*) ETAMAX READ(4,*) READ(4,*) HMIN(1),HMIN(2) READ(4,*) HMAX(1),HMAX(2) READ(4,*) ALLREP READ(4,*) MIDREP READ(4,*) TAPE9 READ(4,*) READ(4,*) RESTRT READ(4,*) KTIME READ(4,*) READ(4,*) CONRAD READ(4,*) DQDTDA READ(4,*) APLANO IF ((APLANO.LE.0.0).AND.(.NOT.RESTRT)) WANDES=0.0 READ(4,*) VDECOL READ(4,*) OLDGRD IF (RESTRT) OLDGRD=.FALSE. READ(4,*) READ(4,*) READ(4,*) GWIDE READ(4,*) GHIGH READ(4,*) GANGLE RETURN END C C C SUBROUTINE SETDIM (N50,N121,N9922,NEXTRA,NXL, + NCDIM,NDIFF,NELROW,NELCOL,NUMNOD, + NUMEL,DIMERR) C C CALCULATES AMOUNTS OF VARIABLE STORAGE SPACE NEEDED VS. AVAILABLE C LOGICAL DIMERR DATA NXEL/275/,NXNOD/21/,NXCDIM/2/,NXXTR/1/ NUMEL=2*NELROW*NELCOL NUMNOD=(2*NELROW+1)*(2*NELCOL+1) NDIFF=2*(2*NELCOL+1) NCDIM=(3*NDIFF+16)*NUMNOD C NCDIM IS PER CONVENTIONS OF SUBROUTINE LIBRARY ESSL, BANDED MATRIX NXL=NUMNOD I11=NUMEL*NXEL I12=N50*NXEL I21=NUMNOD*NXNOD I22=N121*NXNOD I41=NCDIM*NXCDIM I42=N9922*NXCDIM I51=NXL*NXXTR I52=NEXTRA*NXXTR NSUM1=I11+I21+I41+I51 NSUM2=I12+I22+I42+I52 DIMERR=(NUMEL.GT.N50).OR. + (NUMNOD.GT.N121).OR. + (NCDIM.GT.N9922).OR. + (NXL.GT.NEXTRA) IF(DIMERR) WRITE(6,2) 2 FORMAT('0ONE OR MORE VARIABLES ABOVE ARE TOO LARGE.'/ + '0ALL OF THEM MUST BE WITHIN LIMITS BEFORE EXECUTION', + ' CAN PROCEED.'/ + '0INCREASE VALUES OF N50, N121,', + ' NEXTRA, AND/OR N9922', + ' IN PROGRAM FIX PARAMETER STATEMENT.') RETURN END C C C SUBROUTINE GRIDDR(INPUT, NELROW,NELCOL,NUMEL, + MODIFY,NODES) C C CREATES TOPOLOGY OF FINITE ELEMENT GRIDS FOR CRUST AND MANTLE C DIMENSION NODES(6,0:NUMEL) C C NOTE:NUMEL=2*NELROW*NELCOL NROWN=2*NELROW+1 NCOLN=2*NELCOL+1 C NOTE:NUMNOD=NROWN*NCOLN DO 30 I=1,NELROW DO 20 J=1,NELCOL K1=2*NELCOL*(I-1)+2*(J-1)+1 NODES(1,K1)=2*NCOLN*(I-1)+2*(J-1)+1 NODES(2,K1)=NODES(1,K1)+2*NCOLN NODES(3,K1)=NODES(1,K1)+2 NODES(4,K1)=NODES(1,K1)+NCOLN NODES(5,K1)=NODES(4,K1)+1 NODES(6,K1)=NODES(1,K1)+1 K2=K1+1 NODES(1,K2)=NODES(2,K1)+2 NODES(2,K2)=NODES(3,K1) NODES(3,K2)=NODES(2,K1) NODES(4,K2)=NODES(5,K1)+1 NODES(5,K2)=NODES(5,K1) NODES(6,K2)=NODES(2,K1)+1 20 CONTINUE 30 CONTINUE RETURN END C C C SUBROUTINE PAST(ACREEP,ALPHAT,AREAC,AREAM, 2 BCREEP,CCREEP,CODE,CONDNS,CONINT,CONNOD,CONRAD, 3 DCREEP,DETJC,DETJM, 4 DXSC,DXSM,DYSC,DYSM,ECLOG,ECREEP, 5 ESUMC,ESUMM,ETAMAX,FLOWIN, 6 FRIC,G,GEOTHA,GEOTHC,GEOTHM, 7 ITIME,LISTOP,LWORK,NCDIM,NDIFF,NELCOL, 8 NODES,NUMEL,NUMNOD,NXL,ONEKM,OUTSCA,OUTVEC, 9 OUTV2,PUSHHO,PUSHUP,RADIUS, A RAMP,RHOAST,RHOH2O,RHOBAR,SIGBOT, 1 SLABSZ,CPNLAT,IBELOW,X0ELON,Y0NLAT, 2 TASTH,TEMLIM, 3 THIKC,THIKM,TIME,THNKC,THNKM, 4 TSLAB0,TSURF,VC,VM,VISMAX, 5 WC,WM,XIPC,XIPM, 6 XNODC,XNODM,YIPC,YIPM,YNODC,YNODM, 7 ALDONE,NELROW,WANDES, 8 HMAX,HMIN) C C RE-ESTABLISHES ALL IN-PROGRESS ARRAYS FROM REPORT ON "TAPE" C PREVIOUSLY WRITTEN BY SUBROUTINE "TAPE". C DOUBLE PRECISION CODE,FLOWIN LOGICAL ALDONE,CRUST,FAILUR,LISTOP,LOCKIN,LOCKWC,MANTLE DIMENSION ACREEP(2), 2 ALPHAT(2),AREAC(NUMEL),AREAM(NUMEL),BCREEP(2), 3 CCREEP(2),CODE(NCDIM),CONDNS(NUMNOD), 4 CONINT(7,NUMEL),CONNOD(NUMNOD),DCREEP(2), 5 DETJC(7,NUMEL), 6 DETJM(7,NUMEL), 7 DXSC(6,7,NUMEL),DXSM(6,7,NUMEL), 8 DYSC(6,7,NUMEL),DYSM(6,7,NUMEL),ECREEP(2), A ESUMC(2,2,7,NUMEL),ESUMM(2,2,7,NUMEL),FLOWIN(NUMNOD), 1 FRIC(2),GEOTHA(4,7,NUMEL), 3 GEOTHC(4,7,NUMEL),GEOTHM(4,7,NUMEL), 4 HMAX(2),HMIN(2),LISTOP(NUMEL),LWORK(NXL), 5 NODES(6,0:NUMEL),OUTSCA(7,NUMEL),OUTVEC(2,7,NUMEL), 6 OUTV2(2,7,NUMEL), 7 RHOBAR(2) DIMENSION 2 TEMLIM(2), 3 THIKC(7,NUMEL),THIKM(7,NUMEL),THNKC(NUMNOD), 4 THNKM(NUMNOD), 6 VC(2,NUMNOD),VM(2,NUMNOD), 7 WC(NUMNOD),WM(NUMNOD),XIPC(7,NUMEL), 8 XIPM(7,NUMEL),XNODC(NUMNOD),XNODM(NUMNOD), 9 YIPC(7,NUMEL),YIPM(7,NUMEL),YNODC(NUMNOD), A YNODM(NUMNOD) TEMPM(Z,M,I)=MIN(TEMLIM(2),GEOTHM(1,M,I) + +GEOTHM(2,M,I)*Z + +GEOTHM(3,M,I)*Z**2 + +GEOTHM(4,M,I)*Z**3) CALL GOON (ITIME,NTIME, + XNODC,XNODM,YNODC,YNODM,THNKC,CONNOD,THNKM, + GEOTHC,GEOTHM,GEOTHA,VC,VM,WC,WM,ESUMC,ESUMM, + NUMNOD,NUMEL,TIME,ALDONE) CALL INTERP (XNODC,NODES,NUMEL,NUMNOD,XIPC) CALL INTERP (YNODC,NODES,NUMEL,NUMNOD,YIPC) CALL INTERP (XNODM,NODES,NUMEL,NUMNOD,XIPM) CALL INTERP (YNODM,NODES,NUMEL,NUMNOD,YIPM) CALL INTERP (THNKC,NODES,NUMEL,NUMNOD,THIKC) CALL INTERP (THNKM,NODES,NUMEL,NUMNOD,THIKM) DO 10 M=1,7 DO 9 I=1,NUMEL THIKC(M,I)=MAX(THIKC(M,I),HMIN(1)) THIKM(M,I)=MAX(THIKM(M,I),HMIN(2)) THIKC(M,I)=MIN(THIKC(M,I),HMAX(1)) THIKM(M,I)=MIN(THIKM(M,I),HMAX(2)) 9 CONTINUE 10 CONTINUE TASTH=TEMPM(THIKM(5,NUMEL),5,NUMEL) CALL AREAS (NODES,AREAC,XNODC,YNODC,NUMNOD,NUMEL) CALL AREAS (NODES,AREAM,XNODM,YNODM,NUMNOD,NUMEL) CALL DERIV (NUMEL,NUMNOD,NODES,XNODC,YNODC,AREAC, + DETJC,DXSC,DYSC,NUMBAD,LISTOP) CALL DERIV (NUMEL,NUMNOD,NODES,XNODM,YNODM,AREAM, + DETJM,DXSM,DYSM,NUMBAD,LISTOP) LOCKIN=.FALSE. LOCKWC=.FALSE. CALL INTERP (CONNOD,NODES,NUMEL,NUMNOD,CONINT) RETURN END C C C SUBROUTINE GOON (ITIME,NTIME, + XNODC,XNODM,YNODC,YNODM,THNKC,CONNOD,THNKM, + GEOTHC,GEOTHM,GEOTHA,VC,VM,WC,WM,ESUMC,ESUMM, + NUMNOD,NUMEL,TIME,ALDONE) C C READS 'TAPE' WITH THE ARRAYS NEEDED IN ORDER TO C RESTART PROGRAM OR COMPUTE A SET OF PLOTS; C ONLY ESSENTIAL INTEGRATED VARIABLES ARE READ; C PARAMETERS MUST BE RE-INPUT BY "INPUT", AND ALL C RECONSTRUCTABLE ARRAYS MUST BE RECOMPUTED. C CHARACTER*80 TITLE LOGICAL ALDONE DIMENSION CONNOD(NUMNOD),ESUMC(2,2,7,NUMEL),ESUMM(2,2,7,NUMEL), + GEOTHA(4,7,NUMEL),GEOTHC(4,7,NUMEL),GEOTHM(4,7,NUMEL), + THNKC(NUMNOD),THNKM(NUMNOD), + VC(2,NUMNOD),VM(2,NUMNOD), + WC(NUMNOD),WM(NUMNOD), + XNODC(NUMNOD),XNODM(NUMNOD), + YNODC(NUMNOD),YNODM(NUMNOD) ALDONE=.FALSE. NTIME=0 1001 FORMAT(A80) 1002 FORMAT(1P,8E9.2) 1003 FORMAT(0P,8F9.5) 1004 FORMAT(10X,E10.4) 1005 FORMAT(1P,6E13.6) 1006 FORMAT(1P,8E10.3) 1007 FORMAT(0P,F10.3,1P,3E10.3,0P,F10.3,1P,3E10.3) REWIND(8) DO 2000 IT=1,ITIME READ (8,1001,END=8001) TITLE READ (8,1004) TIME TMY=TIME/3.15576E13 NTIME=NTIME+1 READ (8,1001) READ (8,1005) (XNODC(I),I=1,NUMNOD) READ (8,1001) READ (8,1005) (XNODM(I),I=1,NUMNOD) READ (8,1001) READ (8,1005) (YNODC(I),I=1,NUMNOD) READ (8,1001) READ (8,1005) (YNODM(I),I=1,NUMNOD) READ (8,1001) READ (8,1006) (THNKC(I),I=1,NUMNOD) READ (8,1001) READ (8,1006) (CONNOD(I),I=1,NUMNOD) READ (8,1001) READ (8,1006) (THNKM(I),I=1,NUMNOD) READ (8,1001) READ (8,1007) (((GEOTHC(I,J,K),I=1,4),J=1,7),K=1,NUMEL) READ (8,1001) READ (8,1007) (((GEOTHM(I,J,K),I=1,4),J=1,7),K=1,NUMEL) READ (8,1001) READ (8,1007) (((GEOTHA(I,J,K),I=1,4),J=1,7),K=1,NUMEL) READ (8,1001) READ (8,1005) ((VC(I,J),I=1,2),J=1,NUMNOD) READ (8,1001) READ (8,1005) ((VM(I,J),I=1,2),J=1,NUMNOD) READ (8,1001) READ (8,1002) (WC(I),I=1,NUMNOD) READ (8,1001) READ (8,1002) (WM(I),I=1,NUMNOD) READ (8,1001) READ (8,1003) ((ESUMC(1,1,J,K),ESUMC(1,2,J,K), + ESUMC(2,1,J,K),ESUMC(2,2,J,K),J=1,7), + K=1,NUMEL) READ (8,1001) READ (8,1003,END=9001) ((ESUMM(1,1,J,K),ESUMM(1,2,J,K), + ESUMM(2,1,J,K),ESUMM(2,2,J,K),J=1,7), + K=1,NUMEL) 2000 CONTINUE RETURN C C CODE FOR CASE OF MISSING NEXT REPORT C 8001 ALDONE=.TRUE. RETURN C C CODE FOR CASE OF INCOMPLETE LAST REPORT C 9001 DO 9005 I=1,2 DO 9004 J=1,2 DO 9003 K=1,7 DO 9002 L=1,NUMEL ESUMM(I,J,K,L)=0. 9002 CONTINUE 9003 CONTINUE 9004 CONTINUE 9005 CONTINUE RETURN END C C C SUBROUTINE AREAS (NODES,AREA,XNOD,YNOD,NUMNOD,NUMEL) C C COMPUTE AREAS OF ELEMENTS IN GRID AS IF THEY HAD STRAIGHT C SIDES. EFFECT OF SIDE CURVATURE WILL BE HANDLED LATER BY C MULTIPLYING BY DETERMINANT OF JACOBIAN MATRIX FOR THE SIDE- C BENDING MAPPING. NOTE THAT AREA MAY BE NEGATIVE, BUT ELEMENT C IS OK IF DETERMINANT IN DERIV IS ALSO NEGATIVE. C DIMENSION AREA(NUMEL),NODES(6,0:NUMEL),XNOD(NUMNOD),YNOD(NUMNOD) DO 100 INDEX=1,NUMEL I1=NODES(1,INDEX) I2=NODES(2,INDEX) I3=NODES(3,INDEX) AREA(INDEX)= 0.5*(XNOD(I1)*YNOD(I2)-XNOD(I2)*YNOD(I1) + +XNOD(I2)*YNOD(I3)-XNOD(I3)*YNOD(I2) + +XNOD(I3)*YNOD(I1)-XNOD(I1)*YNOD(I3)) 100 CONTINUE RETURN END C C C SUBROUTINE DERIV (NUMEL,NUMNOD,NODES,XNOD,YNOD,AREA, + DETJ,DXS,DYS,NUMBAD,LISTOP) C C CALCULATES DXS AND DYS, THE X-DERIVITIVE AND Y-DERIVITIVE C OF EACH OF THE 6 NODAL FUNCTIONS OF A DEFORMED-TRIANGLE C FINITE ELEMENT, AT EACH OF THE 7 INTEGRATION POINTS IN C THAT ELEMENT. ALSO PROVIDES DETJ, THE DETERMINANT OF THE C JACOBIAN MATRIX FOR THE TRANSFORMATION IN WHICH INTERNAL C POINTS OF A TRIANGLE WITH STRAIGHT SIDES ARE MAPPED INTO C NEW LOCATIONS AS SIDES BEND (BUT CORNERS STAY FIXED). C DOUBLE PRECISION POINTS LOGICAL LISTOP DIMENSION AREA(NUMEL),B(4),C(4),DETJ(7,NUMEL),DN(6,2), + DXS(6,7,NUMEL),DYS(6,7,NUMEL), + LISTOP(NUMEL),NODES(6,0:NUMEL),POINTS(5,7), + X(6),XNOD(NUMNOD),Y(6),YNOD(NUMNOD) COMMON /L1L2L3/ POINTS NUMBAD=0 DO 500 I=1,NUMEL LISTOP(I)=.FALSE. DO 100 J=1,6 NODE=NODES(J,I) X(J)=XNOD(NODE) Y(J)=YNOD(NODE) 100 CONTINUE B(1)=Y(2)-Y(3) B(2)=Y(3)-Y(1) B(3)=Y(1)-Y(2) B(4)=B(1) C(1)=X(3)-X(2) C(2)=X(1)-X(3) C(3)=X(2)-X(1) C(4)=C(1) AI2=1./(2.*AREA(I)) DO 400 M=1,7 DO 200 J=1,3 DN(J,1)=AI2*B(J)*(4.*POINTS(J,M)-1.) DN(J+3,1)=AI2*4.*(B(J)*POINTS(J+1,M) + +B(J+1)*POINTS(J,M)) DN(J,2)=AI2*C(J)*(4.*POINTS(J,M)-1.) DN(J+3,2)=AI2*4.*(C(J)*POINTS(J+1,M) + +C(J+1)*POINTS(J,M)) 200 CONTINUE AJ11=0. AJ12=0. AJ21=0. AJ22=0. DO 300 J=1,6 AJ11=AJ11+DN(J,1)*X(J) AJ12=AJ12+DN(J,1)*Y(J) AJ21=AJ21+DN(J,2)*X(J) AJ22=AJ22+DN(J,2)*Y(J) 300 CONTINUE DETJAC=AJ11*AJ22-AJ12*AJ21 DETJ(M,I)=DETJAC IF ((AREA(I)*DETJAC).LT.0.) LISTOP(I)=.TRUE. AJ11S=AJ11 AJ11=AJ22/DETJAC AJ12=-AJ12/DETJAC AJ21=-AJ21/DETJAC AJ22=AJ11S/DETJAC DO 350 J=1,6 DXS(J,M,I)=AJ11*DN(J,1)+AJ12*DN(J,2) DYS(J,M,I)=AJ21*DN(J,1)+AJ22*DN(J,2) 350 CONTINUE 400 CONTINUE IF (LISTOP(I)) NUMBAD=NUMBAD+1 500 CONTINUE RETURN END C C C SUBROUTINE EXTRAP (AREA,CODE,DETJ,FAILUR, + FLOWIN,FPOLES,NCDIM,NDIFF, + NODES,NUMEL,NUMNOD,NXL,VALUES,LWORK, + LOCKIN,LOCKWC,PHINOD,NELCOL) C C SMOOTHS VALUES OF A SCALAR FIELD KNOWN AT THE INTEGRATION C POINTS (VALUES) TO PRODUCE VALUES AT THE NODES (FPOLES). C INCLUDES OPTION (LOCKIN) TO SET VALUES TO ZERO AT INLAND EDGES, C AND AN OPTION TO SET NODE VALUES TO PHINOD(I) AT WEST EDGE. C LOGICAL FAILUR,LOCKIN,LOCKWC DOUBLE PRECISION CODE,FLOWIN,PHI,WEIGHT DIMENSION AREA(NUMEL),CODE(NCDIM),DETJ(7,NUMEL), + FLOWIN(NUMNOD),FPOLES(NUMNOD),NODES(6,0:NUMEL), + PHI(6,7),PHINOD(NUMNOD), + WEIGHT(7),VALUES(7,NUMEL),LWORK(NXL) COMMON /WGTVEC/ WEIGHT COMMON /PHITAB/ PHI CALL BUILDC (AREA,CODE,DETJ, + NCDIM,NDIFF,NODES,NUMEL,NUMNOD) DO 200 I=1,NUMNOD FLOWIN(I)=0. 200 CONTINUE DO 800 M=1,7 WT=WEIGHT(M) DO 700 I=1,NUMEL VALDA=VALUES(M,I)*AREA(I)*DETJ(M,I)*WT DO 600 J=1,6 K=NODES(J,I) FLOWIN(K)=FLOWIN(K)+PHI(J,M)*VALDA 600 CONTINUE 700 CONTINUE 800 CONTINUE IF (LOCKIN) CALL EBCS (NELCOL,FLOWIN,NUMNOD,NDIFF, + CODE,NCDIM) IF (LOCKWC) THEN NACROS=2*NELCOL+1 NROW=NUMNOD/NACROS DO 850 IR=1,NROW IN=(IR-1)*NACROS+1 VALUE=PHINOD(IN) CALL FIXVAL (IN,FLOWIN,NUMNOD,NDIFF,CODE,NCDIM,VALUE) 850 CONTINUE ENDIF CALL SOLVER (CODE,NCDIM,FLOWIN,NUMNOD,NDIFF,LWORK,NXL,FAILUR) DO 900 I=1,NUMNOD FPOLES(I)=FLOWIN(I) 900 CONTINUE RETURN END C C C SUBROUTINE INTERP (FPOLES,NODES,NUMEL,NUMNOD,VALUES) + C C INTERPOLATES SCALAR FROM NODES TO INTEGRATION POINTS C DOUBLE PRECISION PHI DIMENSION FPOLES(NUMNOD),NODES(6,0:NUMEL), + PHI(6,7),VALUES(7,NUMEL) COMMON /PHITAB/ PHI DO 100 M=1,7 DO 10 I=1,NUMEL VALUES(M,I)=0. 10 CONTINUE 100 CONTINUE DO 200 K=1,6 DO 190 M=1,7 DO 180 I=1,NUMEL VALUES(M,I)=VALUES(M,I)+FPOLES(NODES(K,I))* + PHI(K,M) 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END @PROCESS NOVECTOR C C SUBROUTINE BUILDC (AREA,CODE,DETJ, + NCDIM,NDIFF,NODES,NUMEL,NUMNOD) C C WARNING!!! THE "@PROCESS NOVECTOR" STATEMENT C ABOVE IS NEEDED FOR BUILDC UNDER VS FORTRAN 2.4.0 C BECAUSE OF A COMPILER BUG. IF THIS ROUTINE IS COMPILED C WITH THE VECTOR (DEFAULT) OPTION, IT WILL BE INCORRECT C AND WILL GIVE VERY ODD RESULTS THAT ARE HARD TO TRACE. C C CREATES SMOOTHING MATRIX CODE (CROSS-PRODUCTS OF PHI) C DOUBLE PRECISION CODE,PHI,WEIGHT DIMENSION AREA(NUMEL),CODE(NCDIM),DETJ(7,NUMEL), + NODES(6,0:NUMEL),PHI(6,7),WEIGHT(7) COMMON /PHITAB/ PHI COMMON /WGTVEC/ WEIGHT INDEXK(IR,JC,NBAND)=2*NBAND+1+IR-JC+(JC-1)*(3*NBAND+16) C USE: INDEXK(IR,JC,NDIFF) C BECAUSE THESE TENSORS HAVE HALF THE RANK OF STIFF DO 10 I=1,NCDIM CODE(I)=0. 10 CONTINUE DO 100 I=1,NUMEL DO 90 I6=1,6 DO 80 J6=1,6 IR=NODES(I6,I) JC=NODES(J6,I) SUM=0. DO 70 M=1,7 SUM=SUM+WEIGHT(M)*DETJ(M,I)* + PHI(I6,M)*PHI(J6,M) 70 CONTINUE K=INDEXK(IR,JC,NDIFF) CODE(K)=CODE(K)+SUM*AREA(I) 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END C C C SUBROUTINE SOLVER (ABD,NKDIM,BX,NTNM,NBAND,IPVT,NXL,FAILUR) C C SETS UP FOR CALL TO THE LIBRARY ROUTINE WHICH ACTUALLY C SOLVES THE LINEAR EQUATION SYSTEM C C CURRENT VERSION IS PER CONVENTIONS OF IBM'S ESSL LIBRARY, C DOUBLE PRECISION VERSION. C DOUBLE PRECISION ABD,BX LOGICAL FAILUR DIMENSION ABD(NKDIM),BX(NTNM),IPVT(NXL) C C ****************** IMPORTANT ************************************ C ALL INDEXK DEFINITION STATEMENTS IN PACKAGE MUST AGREE ! INDEXK(IR,JC,NBAND)=2*NBAND+1+IR-JC+(JC-1)*(3*NBAND+16) C****************************************************************** C N=NTNM ML=NBAND MU=NBAND LDA=2*ML+MU+16 CALL DGBF(ABD,LDA,N,ML,MU,IPVT) CALL DGBS(ABD,LDA,N,ML,MU,IPVT,BX) C C PREVENT LATER UNDERFLOWS DUE TO SLOPPY ENFORCEMENT OF 0 VALUES C DO 10 I=1,N SIZE=ABS(BX(I)) IF (SIZE.LE.1.E-35) BX(I)=0. 10 CONTINUE RETURN END C C C SUBROUTINE EBCS (NELCOL,FLOWIN,NUMNOD,NDIFF, + CODE,NCDIM) C C ADDS EDGE BOUNDARY CONDITIONS OF ZERO SCALAR VALUE ON ALL C INLAND BOUNDARY NODES BY OPERATIONS ON MATRIX CODE C AND VECTOR FLOWIN C DOUBLE PRECISION CODE,FLOWIN DIMENSION CODE(NCDIM),FLOWIN(NUMNOD) NACROS=NELCOL*2+1 DO 10 J=1,NACROS CALL FIXVAL (J,FLOWIN,NUMNOD,NDIFF,CODE,NCDIM,0.) 10 CONTINUE JL=NUMNOD JF=JL-NACROS+1 DO 20 J=JF,JL CALL FIXVAL (J,FLOWIN,NUMNOD,NDIFF,CODE,NCDIM,0.) 20 CONTINUE IL=NUMNOD/NACROS-1 DO 30 I=2,IL J=NACROS*I CALL FIXVAL (J,FLOWIN,NUMNOD,NDIFF,CODE,NCDIM,0.) 30 CONTINUE RETURN END C C C SUBROUTINE FIXVAL (I,FLOWIN,NUMNOD,NDIFF,CODE,NCDIM,VALUE) C C SETS SCALAR TO VALUE AT NODE I C CAUTION: C NOTE THAT THIS ROUTINE MAY DESTROY SYMMETRY OF MATRIX CODE C THAT WAS PREVIOUSLY SYMMETRICAL !!! C DOUBLE PRECISION CODE,FLOWIN DIMENSION CODE(NCDIM),FLOWIN(NUMNOD) INDEXK(IR,JC,NBAND)=2*NBAND+1+IR-JC+(JC-1)*(3*NBAND+16) C USE: INDEXK(IR,JC,NDIFF) KD=INDEXK(I,I,NDIFF) PIVOT=CODE(KD) JF=MAX(1,I-NDIFF) JL=MIN(NUMNOD,I+NDIFF) DO 10 J=JF,JL K=INDEXK(I,J,NDIFF) CODE(K)=0. 10 CONTINUE CODE(KD)=MAX(1.,PIVOT) FLOWIN(I)=CODE(KD)*VALUE RETURN END C C C SUBROUTINE REPORT (ITIME,XIPC,XIPM,YIPC,YIPM, + XNODC,XNODM,YNODC,YNODM,TITLE,VM,NODES, + OUTSCA,OUTVEC,VC, + THIKM,THIKC,GEOTHA,GEOTHC,CONDUC, + GEOTHM,VPMEAN,DVPBYE,DVPDT, + TIME2,NUMNOD,NUMEL, + G,HMAX,HMIN,RHOAST,RHOH2O, + TEMLIM,ONEKM,ESUMC,ESUMM,AREAC,AREAM, A CODE,CONDNS,DETJC,DETJM,FAILUR,FLOWIN, + NCDIM,NDIFF,NXL,LWORK,WC,WM, + ECLOG,RADIUS,SLABSZ,CPNLAT,IBELOW, + X0ELON,Y0NLAT,OUTV2,RAMP, + THNKC,TASTH,TSLAB0, + ILIGHT,IVAR,NCONTR, + RMSVEC,NELCOL,PHINOD,DRAWST, + NXYST,XST,YST,NELROW,THNKM,FBLAND,LOWBLU, + CINT,WANDES,CONINT,CONNOD,CONRAD, B TSURF,PUSHUP,VNODE,NEWVAR,DFCON,XTEMP,YTEMP) C C CREATES MAPS OF IMPORTANT QUANTITIES. C ACTUAL DISPLAY OF THE DATA IS HANDLED IN "PAINT" OR "ETCH". C THIS ROUTINE JUST SETS UP THE VALUES AT NODES AND INTEGRATION C POINTS IN TEMPORARY ARRAYS: C C SCALAR VECTOR C ------ ------ C NODES: CONDNS VNODE C INTEGRATION POINTS: OUTSCA OUTVEC C C AND THEN CALLS THE GRAPHICS ROUTINES. C CHARACTER*80 TITLE CHARACTER*42 TEXT,VUNITS LOGICAL ALLPOS,DOAROW,DRAWST, + FAILUR,LOCKIN,LOCKWC,NEWVAR,SOMNEG DOUBLE PRECISION CODE,FLOWIN DIMENSION AREAC(NUMEL),AREAM(NUMEL), 2 CODE(NCDIM),CONDUC(2),CONDNS(NUMNOD), 3 CINT(24),CONINT(7,NUMEL),CONNOD(NUMNOD), 4 DETJC(7,NUMEL),DETJM(7,NUMEL), 5 DRAWST(NXYST),DVPBYE(2,2), 6 DVPDT(2), 7 ESUMC(2,2,7,NUMEL),ESUMM(2,2,7,NUMEL),FBLAND(24), 8 FLOWIN(NUMNOD),GEOTHA(4,7,NUMEL), 9 GEOTHC(4,7,NUMEL),GEOTHM(4,7,NUMEL), A HMAX(2),HMIN(2),LOWBLU(24),LWORK(NXL),NODES(6,0:NUMEL), 1 NVCHAR(24),NVUCHR(24), 2 OUTSCA(7,NUMEL),OUTVEC(2,7,NUMEL),OUTV2(2,7,NUMEL), 3 PHINOD(NUMNOD) DIMENSION TEMLIM(2), 2 TEXT(24),THIKC(7,NUMEL),THIKM(7,NUMEL), 3 THNKC(NUMNOD),THNKM(NUMNOD), 5 VPMEAN(2),VC(2,NUMNOD),VM(2,NUMNOD),VNODE(2,NUMNOD), 7 VUNITS(24),WC(NUMNOD),WM(NUMNOD), 8 XIPC(7,NUMEL),XIPM(7,NUMEL), 9 XNODC(NUMNOD),XNODM(NUMNOD), A XST(NXYST),YST(NXYST),YIPC(7,NUMEL),YIPM(7,NUMEL), 1 YNODC(NUMNOD),YNODM(NUMNOD),XTEMP(NUMNOD), 2 YTEMP(NUMNOD) C DATA TEXT(1)/'MANTLE: GRID OF ELEMENTS '/ DATA NVCHAR(1)/24/ DATA VUNITS(1)/' '/ DATA NVUCHR(1)/0/ DATA TEXT(2)/'CRUST: GRID OF ELEMENTS '/ DATA NVCHAR(2)/24/ DATA VUNITS(2)/' '/ DATA NVUCHR(2)/0/ DATA TEXT(3)/'MANTLE: THICKNESS '/ DATA NVCHAR(3)/17/ DATA VUNITS(3)/'CM '/ DATA NVUCHR(3)/2/ DATA TEXT(4)/'CRUST: THICKNESS '/ DATA NVCHAR(4)/17/ DATA VUNITS(4)/'CM '/ DATA NVUCHR(4)/2/ DATA TEXT(5)/'MANTLE: HEAT-FLOW '/ DATA NVCHAR(5)/17/ DATA VUNITS(5)/'ERG/CM**2-SEC '/ DATA NVUCHR(5)/13/ DATA TEXT(6)/'CRUST: HEAT-FLOW '/ DATA NVCHAR(6)/35/ DATA VUNITS(6)/'ERG/CM**2-SEC '/ DATA NVUCHR(6)/13/ DATA TEXT(7)/'MANTLE: VELOCITY '/ DATA NVCHAR(7)/16/ DATA VUNITS(7)/'CM/SEC '/ DATA NVUCHR(7)/6/ DATA TEXT(8)/'CRUST: VELOCITY '/ DATA NVCHAR(8)/16/ DATA VUNITS(8)/'CM/SEC '/ DATA NVUCHR(8)/6/ DATA TEXT(9)/'MANTLE: DILATATION '/ DATA NVCHAR(9)/18/ DATA VUNITS(9)/' '/ DATA NVUCHR(9)/0/ DATA TEXT(10)/'CRUST: DILATATION '/ DATA NVCHAR(10)/18/ DATA VUNITS(10)/' '/ DATA NVUCHR(10)/0/ DATA ICALL/0/ C LOCKIN=.FALSE. LOCKWC=.FALSE. T2MA=TIME2/(1.E6*365.25*24.*60.*60.) C C IF (IVAR.EQ.1) THEN IF (NEWVAR) THEN DO 101 I=1,NUMNOD VNODE(1,I)=XNODM(I) VNODE(2,I)=YNODM(I) 101 CONTINUE DO 103 M=1,7 DO 102 I=1,NUMEL OUTSCA(M,I)=DETJM(M,I) OUTVEC(1,M,I)=XIPM(M,I) OUTVEC(2,M,I)=YIPM(M,I) 102 CONTINUE 103 CONTINUE ENDIF CALL INLAND (INPUT,VNODE,NUMEL,NUMNOD, + NELROW,XNODC,YNODC, + OUTPUT,CONDNS) CALL ETCH (OUTSCA,DRAWST,ITIME,IVAR, + NODES,NUMEL,NUMNOD,NVCHAR,NXYST,OUTVEC, + TEXT,TITLE,T2MA,VNODE, + XST,YST,CONDNS, + ILIGHT,IVAR,NELROW,XNODC,XNODM,YNODC,YNODM) C C ELSE IF (IVAR.EQ.2) THEN IF (NEWVAR) THEN DO 201 I=1,NUMNOD VNODE(1,I)=XNODC(I) VNODE(2,I)=YNODC(I) 201 CONTINUE DO 203 M=1,7 DO 202 I=1,NUMEL OUTSCA(M,I)=DETJC(M,I) OUTVEC(1,M,I)=XIPC(M,I) OUTVEC(2,M,I)=YIPC(M,I) 202 CONTINUE 203 CONTINUE ENDIF DO 204 I=1,NUMNOD XTEMP(I)=VNODE(1,I) YTEMP(I)=VNODE(2,I) 204 CONTINUE CALL INLAND (INPUT,VNODE,NUMEL,NUMNOD, + NELROW,XTEMP,YTEMP, + OUTPUT,CONDNS) CALL ETCH (OUTSCA,DRAWST,ITIME,IVAR, + NODES,NUMEL,NUMNOD,NVCHAR,NXYST,OUTVEC, + TEXT,TITLE,T2MA,VNODE, + XST,YST,CONDNS, + ILIGHT,IVAR,NELROW,XNODC,XNODM,YNODC,YNODM) C C ELSE IF (IVAR.EQ.3) THEN IF (NEWVAR) THEN DO 301 I=1,NUMNOD CONDNS(I)=THNKM(I) 301 CONTINUE DO 303 M=1,7 DO 302 I=1,NUMEL OUTSCA(M,I)=THIKM(M,I) 302 CONTINUE 303 CONTINUE ENDIF ALLPOS=.FALSE. CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPM,XNODM,YIPM,YNODM) DOAROW=.FALSE. CALL PAINT (NODES,XNODM,YNODM,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.4) THEN IF (NEWVAR) THEN DO 401 I=1,NUMNOD CONDNS(I)=THNKC(I) 401 CONTINUE DO 403 M=1,7 DO 402 I=1,NUMEL OUTSCA(M,I)=THIKC(M,I) 402 CONTINUE 403 CONTINUE ENDIF ALLPOS=.FALSE. CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPC,XNODC,YIPC,YNODC) DOAROW=.FALSE. CALL PAINT (NODES,XNODC,YNODC,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.5) THEN IF (NEWVAR) THEN CALL HEAT(GEOTHM,NUMEL,CONDUC,OUTSCA) ALLPOS=.FALSE. CALL EXTRAP(AREAM,CODE,DETJM,FAILUR, + FLOWIN,CONDNS,NCDIM,NDIFF, + NODES,NUMEL,NUMNOD,NXL,OUTSCA,LWORK, + LOCKIN,LOCKWC,PHINOD,NELCOL) ENDIF CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPM,XNODM,YIPM,YNODM) DOAROW=.FALSE. CALL PAINT (NODES,XNODM,YNODM,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.6) THEN IF (NEWVAR) THEN CALL HEAT(GEOTHC,NUMEL,CONDUC,OUTSCA) ALLPOS=.FALSE. CALL EXTRAP(AREAC,CODE,DETJC,FAILUR, + FLOWIN,CONDNS,NCDIM,NDIFF, + NODES,NUMEL,NUMNOD,NXL,OUTSCA,LWORK, + LOCKIN,LOCKWC,PHINOD,NELCOL) ENDIF CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPC,XNODC,YIPC,YNODC) DOAROW=.FALSE. CALL PAINT (NODES,XNODC,YNODC,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.7) THEN IF (NEWVAR) THEN DO 701 I=1,NUMNOD VNODE(1,I)=VM(1,I) VNODE(2,I)=VM(2,I) 701 CONTINUE ENDIF CALL FLOW (VNODE,NUMNOD,NODES,NUMEL,OUTVEC) ALLPOS=.TRUE. SOMNEG=.NOT.ALLPOS CALL MAGNIN (VNODE,NUMNOD,CONDNS) CALL MAGNIT (OUTVEC,NUMEL,OUTSCA,SOMNEG) CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPM,XNODM,YIPM,YNODM) DOAROW=.TRUE. CALL PAINT (NODES,XNODM,YNODM,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.8) THEN IF (NEWVAR) THEN DO 801 I=1,NUMNOD VNODE(1,I)=VC(1,I) VNODE(2,I)=VC(2,I) 801 CONTINUE ENDIF CALL FLOW (VNODE,NUMNOD,NODES,NUMEL,OUTVEC) ALLPOS=.TRUE. SOMNEG=.NOT.ALLPOS CALL MAGNIN (VNODE,NUMNOD,CONDNS) CALL MAGNIT (OUTVEC,NUMEL,OUTSCA,SOMNEG) CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPC,XNODC,YIPC,YNODC) DOAROW=.TRUE. CALL PAINT (NODES,XNODC,YNODC,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.9) THEN IF (NEWVAR) THEN DO 920 M=1,7 DO 910 I=1,NUMEL RELARE=ESUMM(1,1,M,I)*ESUMM(2,2,M,I)- + ESUMM(1,2,M,I)*ESUMM(2,1,M,I) RELARE=MIN(RELARE,5.0) RELARE=MAX(RELARE,0.3) OUTSCA(M,I)=RELARE 910 CONTINUE 920 CONTINUE CALL EXTRAP (AREAM,CODE,DETJM,FAILUR, + FLOWIN,CONDNS,NCDIM,NDIFF, + NODES,NUMEL,NUMNOD,NXL,OUTSCA,LWORK, + LOCKIN,LOCKWC) ENDIF ALLPOS=.FALSE. CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPM,XNODM,YIPM,YNODM) DOAROW=.FALSE. CALL PAINT (NODES,XNODM,YNODM,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C ELSE IF (IVAR.EQ.10) THEN IF (NEWVAR) THEN DO 1020 M=1,7 DO 1010 I=1,NUMEL RELARE=ESUMC(1,1,M,I)*ESUMC(2,2,M,I)- + ESUMC(1,2,M,I)*ESUMC(2,1,M,I) RELARE=MIN(RELARE,5.0) RELARE=MAX(RELARE,0.3) OUTSCA(M,I)=RELARE 1010 CONTINUE 1020 CONTINUE CALL EXTRAP (AREAC,CODE,DETJC,FAILUR, + FLOWIN,CONDNS,NCDIM,NDIFF, + NODES,NUMEL,NUMNOD,NXL,OUTSCA,LWORK, + LOCKIN,LOCKWC) ENDIF ALLPOS=.FALSE. CALL INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD, + DFCON,NCONTR,XIPC,XNODC,YIPC,YNODC) DOAROW=.FALSE. CALL PAINT (NODES,XNODC,YNODC,TITLE,TEXT,IVAR,T2MA, + CONDNS,DFCON,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) ENDIF RETURN END C C C SUBROUTINE FLOW (V,NUMNOD,NODES,NUMEL,OUTVEC) C C CALCULATES VELOCITY VECTORS AT INTEGRATION POINTS, FROM NODAL VALUES C DOUBLE PRECISION PHI COMMON /PHITAB/ PHI DIMENSION NODES(6,0:NUMEL),OUTVEC(2,7,NUMEL), + PHI(6,7),V(2,NUMNOD) DO 50 M=1,7 DO 40 I=1,NUMEL OUTVEC(1,M,I)=0. OUTVEC(2,M,I)=0. 40 CONTINUE 50 CONTINUE DO 100 J=1,6 DO 90 M=1,7 DO 80 I=1,NUMEL OUTVEC(1,M,I)=OUTVEC(1,M,I)+V(1,NODES(J,I)) + *PHI(J,M) OUTVEC(2,M,I)=OUTVEC(2,M,I)+V(2,NODES(J,I)) + *PHI(J,M) 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END C C C SUBROUTINE HEAT (GEOTH,NUMEL,CONDUC,OUTSCA) + C C CALCULATES SURFACE HEAT-FLOW C DIMENSION CONDUC(2),GEOTH(4,7,NUMEL),OUTSCA(7,NUMEL) DO 100 M=1,7 DO 90 I=1,NUMEL OUTSCA(M,I)=GEOTH(2,M,I)*CONDUC(1) 90 CONTINUE 100 CONTINUE 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 BLOCK DATA BD1 C C DEFINE PHI (NODAL FUNCTIONS) AND WEIGHT (GAUSSIAN INTEGRATION C WEIGHTS) OF 6-NODED TRIANGULAR FINITE ELEMENT FOR THE C SEVEN INTEGRATION POINTS IN EACH ELEMENT, DEFINED BY INTERNAL C COORDINATES POINTS(5,7), WHERE POINTS(1-3,M)=L1-3 OF C INTEGRATION POINT NUMBER M. C DOUBLE PRECISION PHI,POINTS,WEIGHT COMMON /L1L2L3/ POINTS COMMON /PHITAB/ PHI COMMON /WGTVEC/ WEIGHT DIMENSION PHI(6,7),POINTS(5,7),WEIGHT(7) DATA PHI / + -0.1111111111111111,-0.1111111111111111,-0.1111111111111111, + 0.4444444444444444, 0.4444444444444444, 0.4444444444444444, + -0.0525839022774079,-0.0280749439026853,-0.0280749439026853, + 0.1122997756107412, 0.8841342388612960, 0.1122997756107412, + -0.0280749439026853,-0.0525839022774079,-0.0280749439026853, + 0.1122997756107412, 0.1122997756107412, 0.8841342388612960, + -0.0280749439026853,-0.0280749439026853,-0.0525839022774079, + 0.8841342388612960, 0.1122997756107412, 0.1122997756107412, + 0.4743526114618935,-0.0807685938011933,-0.0807685938011933, + 0.3230743752047730, 0.0410358257309469, 0.3230743752047730, + -0.0807685938011933, 0.4743526114618935,-0.0807685938011933, + 0.3230743752047730, 0.3230743752047730, 0.0410358257309469, + -0.0807685938011933,-0.0807685938011933, 0.4743526114618935, + 0.0410358257309469, 0.3230743752047730, 0.3230743752047730/ DATA POINTS / + 0.3333333333333333, 0.3333333333333333, 0.3333333333333333, + 0.3333333333333333, 0.3333333333333333, + 0.0597158733333333, 0.4701420633333333, 0.4701420633333333, + 0.0597158733333333, 0.4701420633333333, + 0.4701420633333333, 0.0597158733333333, 0.4701420633333333, + 0.4701420633333333, 0.0597158733333333, + 0.4701420633333333, 0.4701420633333333, 0.0597158733333333, + 0.4701420633333333, 0.4701420633333333, + 0.7974269866666667, 0.1012865066666667, 0.1012865066666667, + 0.7974269866666667, 0.1012865066666667, + 0.1012865066666667, 0.7974269866666667, 0.1012865066666667, + 0.1012865066666667, 0.7974269866666667, + 0.1012865066666667, 0.1012865066666667, 0.7974269866666667, + 0.1012865066666667, 0.1012865066666667/ DATA WEIGHT / 0.2250000000000000, + 0.1323941500000000, 0.1323941500000000, 0.1323941500000000, + 0.1259391833333333, 0.1259391833333333, 0.1259391833333333/ END C C C SUBROUTINE PAINT (NODES,XNOD,YNOD,TITLE,TEXT,JV,T2MA, + FUNC,CINT,NUMNOD,NUMEL,ALLPOS, + SCALE, + DRAWST,NXYST,XST,YST, + NVCHAR,NVUCHR,VUNITS, + FBLAND,LOWBLU,ITIME, + DOAROW,VNODE,RMSVEC, + ILIGHT,IVAR,NELROW) C C PLOTS CONTOUR DIAGRAMS AND A GRID OF STATE OUTLINES. C LABELS WITH VARIABLE AND TIME ABOVE, MODEL TITLE BELOW. C PLACES COLORBAR WITH CONTOUR VALUES AND UNITS ON RIGHT. C PARAMETER (NCOLOR=12) CHARACTER*80 TITLE,TITLE2 CHARACTER*42 TEXT,TTEXT,VUNITS CHARACTER*5 TMYCHR,CLCHR,ASCII EXTERNAL ASCII LOGICAL ALLPOS,BAR,DOAROW,DRAWST,IMAGE,INSHOW, + STATES,TITLES DIMENSION DRAWST(NXYST),FBLAND(24), + FUNC(NUMNOD),ICOLOR(NCOLOR),LOWBLU(24), + NODES(6,0:NUMEL),NVCHAR(24),NVUCHR(24), + VNODE(2,NUMNOD), + TEXT(24),VUNITS(24), + XNOD(NUMNOD),YNOD(NUMNOD), + XST(NXYST),YST(NXYST) COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE C C SELECT COLORS:128 = WHITE, 12 C 80 = PINK, 11 C 77 = RED, 10 C 69 = DARK RED, 9 C 93 = ORANGE, 8 C 125 = YELLOW, 7 C 117 = YELLOW/GREEN, 6 C 113 = GREEN 5 C 116 = TURQUOISE, 4 C 68 = BLUE, 3 ___________ C 66 = DARK BLUE, 2 __CINT___ C 65 = BLACK. 1 C DATA ICOLOR/65,66,68,116,113,117,125,93,69,77,80,128/ C C STATEMENT FUNCTION: INSHOW(X,Y)=(ABS(X-XCENTR).LE.0.5*WIDE).AND. + (ABS(Y-YCENTR).LE.0.5*WIDE) C C DETERMINE FUNCTION VALUE THAT WILL HAVE MEDIAN COLOR C FTOPS= -9.9E59 FLOOR= +9.9E59 DO 5 I=1,NUMNOD IF (INSHOW(XNOD(I),YNOD(I))) THEN FTOPS=MAX(FTOPS,FUNC(I)) FLOOR=MIN(FLOOR,FUNC(I)) ENDIF 5 CONTINUE FMIDLE=(FTOPS+FLOOR)/2. IFLIP=1 C C CLEAR CURRENT PAGE C CALL FSPCLR C C********************************************************************** IF (IMAGE) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C INITIALIZE SEGMENT 1 (CONTOURED ELEMENTS) C CALL GSSEG(1) C C LOAD 64-COLOR PALLETTE C CALL GSLSS(3,'ADMCOLSD',0) C CALL CONTEL (NODES,XNOD,YNOD,FUNC,CINT,NUMNOD,NUMEL, + FMAX,FMIN,NCOLOR,ICOLOR,FMIDLE,IFLIP, + NBLUE,NYELOW,ALLPOS) C C CLOSE SEGMENT OF COLORED ELEMENTS C CALL GSSCLS C C************************************************************* C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT 2 (VECTORS OR TENSOR SYMBOLS, OR DUMMY) C TENSORS INCLUDED WITH OTHER VARIABLES) C CALL GSSEG(2) C IF (DOAROW) THEN CALL ARROW (NUMNOD,VNODE,RMSVEC, + XNOD,YNOD) ELSE C C DUMMY LINE SO SEGMENT WON'T BE EMPTY C CALL GSCOL(8) CALL GSMOVE(10.5,8.2) CALL GSLINE(10.6,8.2) ENDIF C C CLOSE SEGMENT C CALL GSSCLS ENDIF C C******************************************************************* IF (STATES) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT FOR STATE LINES (3) C CALL GSSEG(3) C C USE BLUE OR YELLOW PEN FOR MAX. CONTRAST C IF (NBLUE.LE.NYELOW) THEN CALL GSCOL(1) ELSE CALL GSCOL(6) ENDIF C C USE HEAVY LINE C CALL GSLW(2) C CALL USMAP (INPUT,DRAWST,NXYST,XST,YST) C C CLOSE SEGMENT WITH STATE LINES C CALL GSSCLS C ENDIF C**************************************************************** C C SEGMENT 4 (HIGHLIGHT AND SOMETIMES NODE NUMBERS) C CALL HILITE(INDATA,ILIGHT,IVAR,NELROW,NUMNOD, + VNODE,XNOD,XNOD,YNOD,YNOD) C C**************************************************************** IF (BAR) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT 5 (COLOR BAR AND CONTOUR INTERVALS) C CALL GSSEG(5) C C MASK OUT AREA WITH BLACK C CALL GSCOL(-1) CALL GSMOVE(9.75,-5.) CALL GSAREA(0) CALL GSLINE(9.75,13.5) CALL GSLINE(16.,13.5) CALL GSLINE(16.,-5.) CALL GSLINE(9.75,-5.) CALL GSENDA C C DETERMINE SCALE FACTORS FOR COLOR BAR C IF (ALLPOS) FMIN=MAX(FMIN,0.) RANGE=FMAX-FMIN STEPS=RANGE/CINT YPERST=MIN(0.5,6.0/MAX(STEPS,0.01)) YTOP=4.25+YPERST*STEPS/2. YBOT=4.25-YPERST*STEPS/2. ORIGIN=4.25-((FMAX+FMIN)/(2.*CINT))*YPERST NSTEPT=IUNDER(FMAX/CINT) NSTEPB=IUNDER(FMIN/CINT) C C LOAD 64-COLOR PALLETTE C CALL GSLSS(3,'ADMCOLSD',0) C C SELECT AND LOAD CHARACTER SET C C CALL GSLSS(2,'ADMUWCRP',199) C C SELECT CHARACTER MODE (2 OR 3) C 2 = DEFAULT TYPE, RECTILINEAR, CONSTANT-SPACED C 3 = VECTOR TYPE, ANY ANGLE, PROPORTIONAL SPACING C CALL GSCM(3) C C COMPLETE CONNECTION OF CHARACTER SET C C CALL GSCS(199) C C SET SIZE OF CHARACTER BOXES RELATIVE TO DEFAULT C CALL GSQCB(WIDTH,HEIGHT) WIDTH=WIDTH*1.00 HEIGHT=HEIGHT*1.00 CALL GSCB(WIDTH,HEIGHT) C C SET LABEL ANGLE TO ZER0 (LIKE THIS) C CALL GSCA(COS(0.0),SIN(0.0)) C C ADD UNITS C TWIDE=WIDTH*NVUCHR(JV) IF (TWIDE.LE.0.25) THEN X=10.75-TWIDE/2. ELSE X=11.0-TWIDE ENDIF Y=YTOP+0.7*HEIGHT YOLD=Y YNEXT=Y-1.1*HEIGHT C C USE WHITE FOR CONTOUR LEVEL LABELS C CALL GSCOL(7) CALL GSLW(2) C CALL GSCHAR(X,Y,NVUCHR(JV),VUNITS(JV)) C C DRAW BOXES AND CONTOUR LABELS C IPOW=IUNDER(ALOG10(ABS(CINT))+0.00001) IF (IPOW.EQ.1) IPOW=0 IF (IPOW.EQ.-1)IPOW=0 CIPOW=IPOW CIINT=ABS(CINT)/10.**CIPOW C DO 1050 I=NSTEPT,NSTEPB,-1 FTOP=(I+1)*CINT IF (I.EQ.NSTEPT) FTOP=FMAX FBOT=I*CINT IF (I.EQ.NSTEPB) FBOT=FMIN YTOP=FTOP*YPERST/CINT+ORIGIN YBOT=FBOT*YPERST/CINT+ORIGIN F=(FTOP+FBOT)/2. N=IHUE (NCOLOR,CINT,FMIDLE,IFLIP,F) CALL GSPAT(ICOLOR(N)) IF (N.EQ.1) THEN ICONT=-2 CALL GSLW(1) ELSE IF (N.EQ.NCOLOR) THEN ICONT=-1 CALL GSLW(2) ELSE ICONT=8 CALL GSLW(2) ENDIF CALL GSCOL(7) CALL GSAREA(1) CALL GSMOVE(11.0,YBOT) CALL GSCOL(ICONT) CALL GSLINE(11.0,YTOP) CALL GSLINE(10.5,YTOP) CALL GSLINE(10.5,YBOT) CALL GSLINE(11.0,YBOT) CALL GSENDA C C USE WHITE FOR CONTOUR LEVEL LABELS C CALL GSCOL(7) CALL GSLW(2) C ARG=1.001*FTOP/10.**CIPOW CLCHR=ASCII(1.001*FTOP/10.**CIPOW) X=10.5-5.5*WIDTH Y=YTOP-0.5*HEIGHT NPRINT=4 IF (CIINT.LT.1.0) NPRINT=5 IF (I.EQ.NSTEPT) NPRINT=5 IF (Y.LE.YNEXT) THEN YOLD=Y YNEXT=Y-1.1*HEIGHT CALL GSCHAR(X,Y,NPRINT,CLCHR) ENDIF IF (I.EQ.NSTEPB) THEN CLCHR=ASCII(1.001*FBOT/10.**CIPOW) X=10.5-5.5*WIDTH Y=YBOT-0.5*HEIGHT IF (Y.LE.YNEXT) CALL GSCHAR(X,Y,5,CLCHR) ENDIF 1050 CONTINUE C C ADD 10**N MULTIPLIER C IF (ABS(CIPOW).GT.0.1) THEN X=11.0-7.0*WIDTH Y=YBOT-2.0*HEIGHT CALL GSCHAR(X,Y,4,'X 1O') IF (ABS(CIPOW-1.).GT.0.1) THEN Y=YBOT-1.3*HEIGHT X=11.0-3.0*WIDTH CLCHR=ASCII(CIPOW) IF (CLCHR(1:1).EQ.' ') X=X-WIDTH IF (CLCHR(2:2).EQ.' ') X=X-WIDTH CALL GSCHAR(X,Y,3,CLCHR) ENDIF ENDIF C C CLOSE SEGMENT OF COLOR BAR C CALL GSSCLS ENDIF C C******************************************************************* IF (TITLES) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT 6 (TITLE + VARIABLE + TIME) C CALL GSSEG(6) C C MASK OUT AREAS WITH BLACK C CALL GSCOL(8) CALL GSMOVE(-5.,-5.) CALL GSAREA(0) CALL GSLINE(-5.,0.5) CALL GSLINE(16.,0.5) CALL GSLINE(16.,-5.) CALL GSLINE(-5.,-5.) CALL GSENDA CALL GSMOVE(-5.,8.08) CALL GSAREA(0) CALL GSLINE(-5.,13.) CALL GSLINE(16.,13.) CALL GSLINE(16.,8.08) CALL GSLINE(-5.,8.08) CALL GSENDA C C SELECT AND LOAD CHARACTER SET C C CALL GSLSS(2,'ADMUWCRP',199) C C SELECT CHARACTER MODE (2 OR 3) C 2 = DEFAULT TYPE, RECTILINEAR, CONSTANT-SPACED C 3 = VECTOR TYPE, ANY ANGLE, PROPORTIONAL SPACING C CALL GSCM(3) C C COMPLETE CONNECTION OF CHARACTER SET C C CALL GSCS(199) C C SET SIZE OF CHARACTER BOXES RELATIVE TO DEFAULT C CALL GSQCB(WIDTH,HEIGHT) WIDTH=1.00*WIDTH HEIGHT=1.25*HEIGHT CALL GSCB(WIDTH,HEIGHT) C C SET LABEL ANGLE TO ZER0 (LIKE THIS) C CALL GSCA(COS(0.0),SIN(0.0)) C C USE WHITE FOR TEXT C CALL GSCOL(7) C C WRITE MODEL TITLE C DO 4010 I=1,80 IF (TITLE(I:I).EQ.'0') THEN TITLE2(I:I)='O' ELSE TITLE2(I:I)=TITLE(I:I) ENDIF 4010 CONTINUE CALL GSCHAR(-0.4,0.1,80,TITLE2) C C WRITE VARIABLE AND TIME IDENTIFIERS C CALL GSCHAR(-0.4,8.1,NVCHAR(JV),TEXT(JV)) CALL GSCHAP(4,' at ') TMYCHR=ASCII(T2MA) CALL GSCHAP(5,TMYCHR) CALL GSCHAP(5,' Ma (') CALL EPOCH (IMPUT,T2MA,OUTPUT,NECHAR,TTEXT) CALL GSCHAP(NECHAR,TTEXT) CALL GSCHAP(1,')') C C CLOSE SEGMENT WITH TEXT LABELS C CALL GSSCLS ENDIF C C**************************************************************** C RETURN END C C C SUBROUTINE ETCH (DETJ,DRAWST,ITIME,JV, + NODES,NUMEL,NUMNOD,NVCHAR,NXYST,OUTVEC, + TEXT,TITLE,T2MA,VNODE, + XST,YST,CONDNS, + ILIGHT,IVAR,NELROW) C C PLOTS THE FINITE ELEMENT GRID AND STATE OUTLINES. C LABELS WITH GRID LEVEL AND TIME ABOVE, MODEL TITLE BELOW. C NOTE: INPUT "CONDNS" IS USED ONLY IN THIS PARTICULAR PROGRAM, C (AND ONLY IF .NOT.SHONUM ) C TO HOLD DISTANCE FROM WESTERN EDGE OF CRUSTAL GRID. C IT IS IN CM, AND THEN IS CONVERTED TO KILOMETERS. C CHARACTER*80 TITLE,TITLE2 CHARACTER*42 TEXT,TTEXT CHARACTER*5 TMYCHR,ASCII,LABEL EXTERNAL ASCII LOGICAL BAR,DRAWST,IMAGE,NOMAP,S4,S5,S6,STATES,TITLES,SHONUM DIMENSION DETJ(7,NUMEL),DRAWST(NXYST),IDSEG(1), + NODES(6,0:NUMEL),NVCHAR(24),OUTVEC(2,7,NUMEL), + TEXT(24),VNODE(2,NUMNOD), + XST(NXYST), + YST(NXYST) DIMENSION CONDNS(NUMNOD) COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE C C STATEMENT FUNCTIONS: XPLT(X)=5.5+5.5*(X-XCENTR)/(0.5*WIDE) YPLT(Y)=4.25+5.5*(Y-YCENTR)/(0.5*WIDE) C C C CLEAR CURRENT PAGE C CALL FSPCLR C C********************************************************************** IF (IMAGE) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C INITIALIZE SEGMENT 1 (FINITE ELEMENT GRID) C CALL GSSEG(1) C C PLOT ALL ELEMENT SIDES (MANY ARE DRAWN TWICE) C CALL GSCOL(4) CALL GSLW(2) DO 10 I=1,NUMEL S4=.TRUE. S5=.TRUE. S6=.TRUE. CALL AROUND (I,S4,S5,S6,NODES,VNODE,NUMNOD,NUMEL) 10 CONTINUE C C PLOT ALL NODES, AND (?) LABEL DISTANCE FROM EDGE OF PLATE C CALL GSCOL(6) CALL GSLW(1) CALL GSMS(6) CALL GSMSC(0.08) C C SELECT AND LOAD CHARACTER SET C C CALL GSLSS(2,'ADMUWCRP',199) C C SELECT CHARACTER MODE (2 OR 3) C 2 = DEFAULT TYPE, RECTILINEAR, CONSTANT-SPACED C 3 = VECTOR TYPE, ANY ANGLE, PROPORTIONAL SPACING C CALL GSCM(3) C C COMPLETE CONNECTION OF CHARACTER SET C C CALL GSCS(199) C C SET SIZE OF CHARACTER BOXES RELATIVE TO DEFAULT C HOWHI=YPLT(50.E5)-YPLT(0.0) NOMAP=(.NOT.SHONUM).AND.(HOWHI.GE.0.10) HOWHI=MIN(HOWHI,0.25) CALL GSQCB(W,H) CALL GSCB(HOWHI*W/H,HOWHI) CALL GSCA(1.,0.) DO 20 I=1,NUMNOD XP=XPLT(VNODE(1,I)) YP=YPLT(VNODE(2,I)) CALL GSMARK(XP,YP) IF (NOMAP) THEN XP=XP-2.6*HOWHI*W/H YP=YP+0.2*HOWHI WRITE(LABEL,19)CONDNS(I) 19 FORMAT(-5P,F5.0) CALL GSCHAR(XP,YP,4,LABEL) ENDIF 20 CONTINUE C C PLOT SYMBOLS AT INTEGRATION POINTS: C GREEN IF DETJ.GT.0.5 C YELLOW IF DETJ BETWEEN 0. AND 0.5 C RED IF DETJ .LT. 0. C R=4.0E5 DO 30 I=1,NUMEL DO 25 M=1,7 XC=OUTVEC(1,M,I) YC=OUTVEC(2,M,I) XP=XC+R IF (DETJ(M,I).GT.0.5) THEN CALL GSCOL(4) ELSE IF (DETJ(M,I).GT.0.0) THEN CALL GSCOL(6) ELSE CALL GSCOL(2) ENDIF CALL GSMOVE(XPLT(XP),YPLT(YC)) CALL GSAREA(1) CALL GSARC(XPLT(XC),YPLT(YC),360.) CALL GSENDA 25 CONTINUE 30 CONTINUE C C CLOSE SEGMENT OF FINITE ELEMENT GRID C CALL GSSCLS ENDIF C C**************************************************************** IF (IMAGE) THEN C C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C OPEN SEGMENT 2 (DUMMY SEGMENT, CORRESPONDS TO VECTORS AND C TENSORS INCLUDED WITH OTHER VARIABLES) C CALL GSSEG(2) C C DUMMY LINE SO SEGMENT WON'T BE EMPTY C CALL GSCOL(8) CALL GSMOVE(10.5,8.2) CALL GSLINE(10.6,8.2) C C CLOSE DUMMY SEGMENT C CALL GSSCLS ENDIF C C******************************************************************* IF (STATES) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT FOR STATE LINES (3) C CALL GSSEG(3) C C USE BLUE PEN C CALL GSCOL(1) C C USE HEAVY LINE C CALL GSLW(2) C CALL USMAP (INPUT,DRAWST,NXYST,XST,YST) C C CLOSE SEGMENT WITH STATE LINES C CALL GSSCLS C ENDIF C**************************************************************** C C SEGMENT 4 (HIGHLIGHT AND SOMETIMES NODE NUMBERS) C CALL HILITE(INDATA,ILIGHT,IVAR,NELROW,NUMNOD, + VNODE,XNODC,XNODM,YNODC,YNODM) C C**************************************************************** IF (BAR) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C OPEN SEGMENT 5 (DUMMY SEGMENT, CORRESPONDS TO COLOR BAR C INCLUDED WITH OTHER VARIABLES) C CALL GSSEG(5) C C DUMMY LINE SO SEGMENT WON'T BE EMPTY C CALL GSCOL(8) CALL GSMOVE(10.5,8.2) CALL GSLINE(10.6,8.2) C C CLOSE DUMMY SEGMENT C CALL GSSCLS ENDIF C C******************************************************************* IF (TITLES) THEN C C ESTABLISH COORDINATE SYSTEM FOR ENTIRE PROGRAM C CALL GSUWIN(0.0,11.0,0.0,8.5) C C C OPEN SEGMENT 6 (TITLE + VARIABLE + TIME) C CALL GSSEG(6) C C MASK OUT AREAS WITH BLACK C CALL GSCOL(8) CALL GSMOVE(-5.,-5.) CALL GSAREA(0) CALL GSLINE(-5.,0.4) CALL GSLINE(16.,0.4) CALL GSLINE(16.,-5.) CALL GSLINE(-5.,-5.) CALL GSENDA CALL GSMOVE(-5.,8.0) CALL GSAREA(0) CALL GSLINE(-5.,13.) CALL GSLINE(16.,13.) CALL GSLINE(16.,8.0) CALL GSLINE(-5.,8.0) CALL GSENDA C C SELECT AND LOAD CHARACTER SET C C CALL GSLSS(2,'ADMUWCRP',199) C C SELECT CHARACTER MODE (2 OR 3) C 2 = DEFAULT TYPE, RECTILINEAR, CONSTANT-SPACED C 3 = VECTOR TYPE, ANY ANGLE, PROPORTIONAL SPACING C CALL GSCM(3) C C COMPLETE CONNECTION OF CHARACTER SET C C CALL GSCS(199) C C SET SIZE OF CHARACTER BOXES RELATIVE TO DEFAULT C CALL GSQCB(WIDTH,HEIGHT) WIDTH=1.00*WIDTH HEIGHT=1.25*HEIGHT CALL GSCB(WIDTH,HEIGHT) C C SET LABEL ANGLE TO ZER0 (LIKE THIS) C CALL GSCA(COS(0.0),SIN(0.0)) C C USE WHITE FOR TEXT C CALL GSCOL(7) C C WRITE MODEL TITLE C DO 410 I=1,80 IF (TITLE(I:I).EQ.'0') THEN TITLE2(I:I)='O' ELSE TITLE2(I:I)=TITLE(I:I) ENDIF 410 CONTINUE CALL GSCHAR(-0.4,0.0,80,TITLE2) C C WRITE VARIABLE AND TIME IDENTIFIERS C CALL GSCHAR(-0.4,8.1,NVCHAR(JV),TEXT(JV)) CALL GSCHAP(4,' at ') TMYCHR=ASCII(T2MA) CALL GSCHAP(5,TMYCHR) CALL GSCHAP(5,' Ma (') CALL EPOCH (IMPUT,T2MA,OUTPUT,NECHAR,TTEXT) CALL GSCHAP(NECHAR,TTEXT) CALL GSCHAP(1,')') C C CLOSE SEGMENT WITH TEXT LABELS C CALL GSSCLS ENDIF C C**************************************************************** C RETURN END C C C SUBROUTINE CONTEL (NODES,XNOD,YNOD,FUNC,DFCON,NUMNOD,NUMEL, + FGMAX,FGMIN,NCOLOR,ICOLOR,FMIDLE,IFLIP, + NBLUE,NYELOW,ALLPOS) C C CONTOURS AND COLORS A SCALAR FIELD ON THE FINITE ELEMENT GRID. C INSTEAD OF FOLLOWING CONTOURS ACROSS ELEMENT BOUNDARIES, IT C CONTOURS EACH ELEMENT SEPARATELY. C HOWEVER, NOTE THAT ELEMENTS ARE SKIPPED WHEN ENTIRELY OUTSIDE C THE FRAME OF THE PICTURE. C PARAMETER(NINLIN=400,NWORK=4000) LOGICAL ALLPOS,ANEDGE,BEGCON,BEGNXT,BITSEG,CENTER,CIRCLE, + COLOR,DASHED,DONE,ENDCON,FINISH,GONOUT, + HITLIM,INSIDE, + NEIGHB,SIDMAX,SIDMIN, + THRU,Z LOGICAL ANYIN,BAR,IMAGE,INSHOW,STATES,TITLES REAL LOWEST COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE DIMENSION NODES(6,0:NUMEL),ICOLOR(NCOLOR), + XNOD(NUMNOD),YNOD(NUMNOD),FUNC(NUMNOD) C C LOCAL STORAGE FOR PROPERTIES OF ONE ELEMENT: DIMENSION IN(6),XN(6),YN(6),FN(6),DS(3),SSE(3),DSIN(3) C C LOCAL STORAGE FOR INITIAL POINTS OF CONTOUR SEGMENTS: DIMENSION PS(5,NINLIN),PS2(5,NINLIN),DONE(NINLIN) C C LOCAL STORAGE FOR SHAPES OF CONTOUR AND EDGE SEGMENTS: DIMENSION SPACE(2,NWORK),ISPPNT(0:NINLIN),ISPLEN(0:NINLIN), & FOFSEG(NINLIN),ANEDGE(NINLIN),MENU(NINLIN), & NTOGO(NINLIN) C DATA DSTEP/0.05/ C C STATEMENT FUNCTIONS: C PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6)= + F1*(-S1+2.*S1**2)+ + F2*(-S2+2.*S2**2)+ + F3*(-S3+2.*S3**2)+ + F4*(4.*S1*S2)+ + F5*(4.*S2*S3)+ + F6*(4.*S3*S1) XPLT(X)=5.5+5.5*(X-XCENTR)/(0.5*WIDE) YPLT(Y)=4.25+5.5*(Y-YCENTR)/(0.5*WIDE) INSHOW(X,Y)=(ABS(X-XCENTR).LE.0.5*WIDE).AND. + (ABS(Y-YCENTR).LE.0.5*WIDE) C C GLOBAL INITIALIZATION (WHOLE GRID) C LIMINT=4./DSTEP FGMIN= 1.E60 FGMAX=-1.E60 NBLUE=0 NYELOW=0 DO 9999 IEL=1,NUMEL C C LOCAL INITIALIZATION (ONE ELEMENT) C NPS=0 ISPNUM=0 ISPLEN(0)=0 ISPPNT(0)=1 CENTER=.FALSE. TSIDE=1.E60 DO 5 J=1,6 K=NODES(J,IEL) IN(J)=K XN(J)=XNOD(K) YN(J)=YNOD(K) FN(J)=FUNC(K) 5 CONTINUE I1=IN(1) I2=IN(2) I3=IN(3) I4=IN(4) I5=IN(5) I6=IN(6) X1=XN(1) X2=XN(2) X3=XN(3) X4=XN(4) X5=XN(5) X6=XN(6) Y1=YN(1) Y2=YN(2) Y3=YN(3) Y4=YN(4) Y5=YN(5) Y6=YN(6) C C IS ANY OF IT VISIBLE? C ANYIN=INSHOW(X1,Y1).OR.INSHOW(X2,Y2).OR.INSHOW(X3,Y3).OR. + INSHOW(X4,Y4).OR.INSHOW(X5,Y5).OR.INSHOW(X6,Y6) IF (.NOT.ANYIN) GO TO 9999 C FMAX=MAX(FN(1),FN(2),FN(3),FN(4),FN(5),FN(6)) FMIN=MIN(FN(1),FN(2),FN(3),FN(4),FN(5),FN(6)) RANGE=MAX((FMAX-FMIN),DFCON) C C PREVENT DEGENERATE CASES WHERE NODES FALL EXACTLY ON CONTOURS C DO 10 J=1,6 I=FN(J)/DFCON IF ((I*DFCON).EQ.FN(J)) FN(J)=FN(J)+0.01*RANGE 10 CONTINUE F1=FN(1) F2=FN(2) F3=FN(3) F4=FN(4) F5=FN(5) F6=FN(6) IHIC= -9999 ILOC= +9999 C C*************************************************************** C C EXAMINE SIDES FOR EXTREMA AND MARK CONTOUR INTERSECTIONS C DO 100 MSIDE=1,3 N1=MSIDE N2=MOD(MSIDE,3)+1 SIDE=SQRT((XN(N1)-XN(N2))**2+(YN(N1)-YN(N2))**2) TSIDE=MIN(TSIDE,SIDE) NM=MSIDE+3 DFDS1= -3.*FN(N1)+4.*FN(NM)- FN(N2) DFDS2= FN(N1)-4.*FN(NM)+3.*FN(N2) D2FDS=4.*FN(N1)-8.*FN(NM)+4.*FN(N2) IF ((DFDS1*DFDS2.GE.0.).OR.(D2FDS.EQ.0.0)) THEN FMX=AMAX1(FN(N1),FN(N2)) FMN=AMIN1(FN(N1),FN(N2)) CALL DOSIDE (FMX,FMN,DFCON,FN,N1,N2,NM,PS,NPS) ELSE SEXT= -DFDS1/D2FDS FEXT=FN(N1)+ + DFDS1*SEXT+ + 0.5*D2FDS*SEXT**2 FMAX=MAX(FMAX,FEXT) FMIN=MIN(FMIN,FEXT) C C FIND INTERSECTIONS OF CONTOURS WITH SIDE CONTAINING EXTREMUM C FMX=AMAX1(FN(N1),FEXT) FMN=AMIN1(FN(N1),FEXT) CALL DOPART (FEXT,FMX,FMN,DFCON,FN, + N1,N2,NM,0.,SEXT,PS,NPS) FMX=AMAX1(FEXT,FN(N2)) FMN=AMIN1(FEXT,FN(N2)) CALL DOPART (FEXT,FMX,FMN,DFCON,FN, + N1,N2,NM,SEXT,1.,PS,NPS) ENDIF 100 CONTINUE RTESTR=TSIDE*DSTEP C C*************************************************************** C C SORT THE POINTS FOUND BY CLOCKWISE PARAMETER S C DO 200 INPS=1,NPS S1=PS(1,INPS) S2=PS(2,INPS) S3=PS(3,INPS) F=PS(4,INPS) IF (F.GE.0.) THEN IC=F/DFCON+0.1 ELSE IT= -F/DFCON+0.1 IC= -IT ENDIF IHIC=MAX(IHIC,IC) ILOC=MIN(ILOC,IC) IF (S3.EQ.0.) THEN SN=S2 ELSE IF (S1.EQ.0.) THEN SN=1.+S3 ELSE SN=2.+S1 ENDIF PS(5,INPS)=SN 200 CONTINUE SNOW= -0.1 DO 300 INPS=1,NPS LOWEST=3.1 JMOVE=INPS DO 250 JNPS=1,NPS IF (PS(5,JNPS).GT.SNOW) THEN IF (PS(5,JNPS).LT.LOWEST) THEN LOWEST=PS(5,JNPS) JMOVE=JNPS ENDIF ENDIF 250 CONTINUE DO 270 K=1,5 PS2(K,INPS)=PS(K,JMOVE) 270 CONTINUE SNOW=LOWEST 300 CONTINUE DO 320 I=1,5 DO 310 J=1,NPS PS(I,J)=PS2(I,J) 310 CONTINUE 320 CONTINUE C C CREATE TABLE OF ELEMENT-SIDE SEGMENTS C S=0. NPSD=0 BEGNXT=.FALSE. C BEGIN NEW SEGMENT 400 ISPNUM=ISPNUM+1 ISPPNT(ISPNUM)=ISPPNT(ISPNUM-1)+ISPLEN(ISPNUM-1) IF (ISPPNT(ISPNUM).GT.NWORK) THEN WRITE(6,401) IEL GO TO 9999 ENDIF 401 FORMAT(' INSUFFICIENT WORKSPACE IN CONTEL. ELEMENT ', & I5,' WILL NOT BE SHOWN.') ISPLEN(ISPNUM)=1 NTOGO(ISPNUM)=1 ANEDGE(ISPNUM)=.TRUE. BEGCON=BEGNXT NINSEG=1 IF(S.LE.1.) THEN S1=1.-S S2=S S3=0. ELSE IF (S.LE.2.) THEN S1=0. S2=2.-S S3=S-1. ELSE S1=S-2. S2=0. S3=3.-S ENDIF X=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) Y=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) F=PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6) SUMSEG=F SPACE(1,ISPPNT(ISPNUM))=X SPACE(2,ISPPNT(ISPNUM))=Y C FIND NEXT POINT 500 IS=IABOVE(S/DSTEP+0.05) IF (S.LT.1.) THEN SLIM=1. ELSE IF (S.LT.2.) THEN SLIM=2. ELSE SLIM=3. ENDIF ST=MIN(SLIM,IS*DSTEP) THRU=.FALSE. ENDCON=.FALSE. BEGNXT=.FALSE. IF (NPSD.LT.NPS) THEN IF (PS(5,NPSD+1).LE.ST) THEN NPSD=NPSD+1 IF ((.NOT.ALLPOS).OR.(PS(4,NPSD).GT.0.0)) THEN THRU=.TRUE. ENDCON=.TRUE. BEGNXT=.TRUE. ST=PS(5,NPSD) ENDIF ENDIF ENDIF IF (ST.EQ.SLIM) THRU=.TRUE. C UPDATE REPRESENTATIVE FUNCTION VALUE FOR SEGMENT NINSEG=NINSEG+1 IF(ST.LE.1.) THEN S1=1.-ST S2=ST S3=0. ELSE IF (ST.LE.2.) THEN S1=0. S2=2.-ST S3=ST-1. ELSE S1=ST-2. S2=0. S3=3.-ST ENDIF F=PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6) BITSEG=THRU.AND.(NINSEG.EQ.2).AND.BEGCON.AND.ENDCON + .AND.(ABS(F-SUMSEG).LT.(0.1*DFCON)) IF (BITSEG) THEN IF (F.GT.FOFSEG(ISPNUM-1)) THEN FOFSEG(ISPNUM)=F+0.5*ABS(DFCON) ELSE FOFSEG(ISPNUM)=F-0.5*ABS(DFCON) ENDIF ELSE SUMSEG=SUMSEG+F FOFSEG(ISPNUM)=SUMSEG/NINSEG ENDIF C RECORD NEXT POINT IN LIST X=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) Y=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) ISPLEN(ISPNUM)=ISPLEN(ISPNUM)+1 IF ((ISPPNT(ISPNUM)+ISPLEN(ISPNUM)-1).GT.NWORK) THEN WRITE(6,401) IEL GO TO 9999 ENDIF SPACE(1,ISPPNT(ISPNUM)+ISPLEN(ISPNUM)-1)=X SPACE(2,ISPPNT(ISPNUM)+ISPLEN(ISPNUM)-1)=Y S=ST IF (S.LT.3.0) THEN IF (THRU) THEN GO TO 400 ELSE GO TO 500 ENDIF ENDIF C C CLEAN-UP THE SEGMENT F VALUES TO MID-RANGE NUMBERS C DO 550 I=1,ISPNUM T=FOFSEG(I)/DFCON IF (ALLPOS) T=MAX(T,0.0) T=IUNDER(T)+0.5 FOFSEG(I)=T*DFCON 550 CONTINUE C C*************************************************************** C C SEARCH FOR EXTREMUM WITHIN DOMAIN OF ELEMENT C CDET=16.*(-F6**2+2.*F6*F5+2.*F6*F4-2.*F6*F2-F5**2+2.*F5*F4 + -2.*F5*F1-F4**2-2.*F4*F3+F3*F2+F3*F1+F2*F1) IF (ABS(CDET).LT. 1.E-40) GO TO 1000 S2EXT=(4.*F6**2-4.*F6*F5-4.*F6*F4-F6*F3+2.*F6*F2-F6*F1+F5* + F3+3.*F5*F1+3.*F4*F3+F4*F1-F3*F2-2.*F3*F1-F2*F1)/(4.*(F6 + **2-2.*F6*F5-2.*F6*F4+2.*F6*F2+F5**2-2.*F5*F4+2.*F5*F1+F4 + **2+2.*F4*F3-F3*F2-F3*F1-F2*F1)) S3EXT=(-4.*F6*F4+3.*F6*F2+F6*F1-4.*F5*F4+F5*F2+3.*F5*F1+4. + *F4**2+2.*F4*F3-F4*F2-F4*F1-F3*F2-F3*F1-2.*F2*F1)/(4.*(F6 + **2-2.*F6*F5-2.*F6*F4+2.*F6*F2+F5**2-2.*F5*F4+2.*F5*F1+F4 + **2+2.*F4*F3-F3*F2-F3*F1-F2*F1)) S1EXT=1.0-S2EXT-S3EXT IF (S1EXT.GT.0.99999.OR.S1EXT.LT.0.00001) GO TO 1000 IF (S2EXT.GT.0.99999.OR.S2EXT.LT.0.00001) GO TO 1000 IF (S3EXT.GT.0.99999.OR.S3EXT.LT.0.00001) GO TO 1000 C C REJECT SADDLE POINTS C DISCA=F1-2.*F4+F2 DISCB=F2-2.*F5+F3 DISCC=F3-2.*F6+F1 CENTER=((DISCA.GT.0.).AND.(DISCB.GT.0.).AND.(DISCC.GT.0.)) + .OR.((DISCA.LT.0.).AND.(DISCB.LT.0.).AND.(DISCC.LT.0.)) IF (.NOT.CENTER) GO TO 1000 XEXT=PHIVAL(S1EXT,S2EXT,S3EXT,X1,X2,X3,X4,X5,X6) YEXT=PHIVAL(S1EXT,S2EXT,S3EXT,Y1,Y2,Y3,Y4,Y5,Y6) FEXT=PHIVAL(S1EXT,S2EXT,S3EXT,F1,F2,F3,F4,F5,F6) FMAX=MAX(FMAX,FEXT) FMIN=MIN(FMIN,FEXT) C C FIND CONTOUR STARTING/STOPPING POINT ALONG CHORD FROM A NODE TO EXT. C NCL=1 DIFF=ABS(F1-FEXT) DS(1)=S1EXT-1. DS(2)=S2EXT DS(3)=S3EXT DO 600 J=2,6 DFF=ABS(FN(J)-FEXT) IF (DFF.LT.DIFF) THEN NCL=J DIFF=DFF DS(1)=S1EXT DS(2)=S2EXT DS(3)=S3EXT IF (J.EQ.2) DS(2)=S2EXT-1. IF (J.EQ.3) DS(3)=S3EXT-1. IF (J.EQ.4.OR.J.EQ.6) DS(1)=S1EXT-0.5 IF (J.EQ.4.OR.J.EQ.5) DS(2)=S2EXT-0.5 IF (J.EQ.5.OR.J.EQ.6) DS(3)=S3EXT-0.5 ENDIF 600 CONTINUE CALL DOLINE (FEXT,DFCON,FN,NCL,DS,S1EXT,S2EXT,S3EXT, + IHIC,ILOC,PS,NPS) C C END OF CODE RELATED TO CASE OF AN INTERNAL EXTREMEM C 1000 IF (NPS.EQ.0) GO TO 9001 C C************************************************************* C C INTEGRATE ALL CONTOUR SEGMENTS C DO 1150 K=1,NPS DONE(K)=.FALSE. 1150 CONTINUE DO 9000 N=1,NPS C C INTEGRATE ONE CONTOUR SEGMENT C IF (.NOT.DONE(N)) THEN C C INITIALIZE INTEGRATION OF CONTOUR C DONE(N)=.TRUE. FVALUE=PS(4,N) IF (ALLPOS.AND.(FVALUE.LE.0.0)) GO TO 9000 ISPNUM=ISPNUM+1 FOFSEG(ISPNUM)=FVALUE ISPPNT(ISPNUM)=ISPPNT(ISPNUM-1)+ISPLEN(ISPNUM-1) IF (ISPPNT(ISPNUM).GT.NWORK) THEN WRITE(6,401) IEL GO TO 9999 ENDIF ISPLEN(ISPNUM)=1 NTOGO(ISPNUM)=2 ANEDGE(ISPNUM)=.FALSE. S1=PS(1,N) S2=PS(2,N) S3=PS(3,N) INSIDE=(S1*S2*S3).GT.0.0 S1OLD=S1 S2OLD=S2 S3OLD=S3 X=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) Y=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) SPACE(1,ISPPNT(ISPNUM))=X SPACE(2,ISPPNT(ISPNUM))=Y ANGLE=0. IF (CENTER) ANGLE=ATAN2((Y-YEXT),(X-XEXT)) ANGLEP=ANGLE ROT=0. DFDS2=-4.*S3*F6+4.*S3*F5-4.*S3*F4+4.*S3*F1-8.*S2*F4+4.*S2* + F2+4.*S2*F1+4.*F4-F2-3.*F1 DFDS3=-8.*S3*F6+4.*S3*F3+4.*S3*F1-4.*S2*F6+4.*S2*F5-4.*S2* + F4+4.*S2*F1+4.*F6-F3-3.*F1 GSIZE=MAX(ABS(DFDS2),ABS(DFDS3)) IF (GSIZE.EQ.0.0) GO TO 8999 DFDS2=DFDS2/GSIZE DFDS3=DFDS3/GSIZE GRADF=SQRT(DFDS2**2+DFDS3**2) GRADFX=DFDS2/GRADF GRADFY=DFDS3/GRADF ROUNDX= +GRADFY ROUNDY= -GRADFX DS2=ROUNDX*DSTEP*0.1 DS3=ROUNDY*DSTEP*0.1 C C REVERSE INTEGRATION STEP DIRECTION IF CONTOUR POINTS OUTWARD C S2P=S2+DS2 S3P=S3+DS3 S1P=1.00-S2P-S3P COUNTR=1. IF ( (S1P.LT.0..OR.S1P.GT.1.) + .OR.(S2P.LT.0..OR.S2P.GT.1.) + .OR.(S3P.LT.0..OR.S3P.GT.1.)) COUNTR= -1. NSEG=0 C C BEGIN LOOP OF INTEGRATION OF CONTOUR LINE C-------------------------------------------- C 3000 NSEG=NSEG+1 C EXTRAPOLATE TO NEXT POINT BY FORWARD METHOD DS2=ROUNDX*COUNTR*DSTEP DS3=ROUNDY*COUNTR*DSTEP S2P=S2+DS2 S3P=S3+DS3 S1P=1.00-S2P-S3P C RECOMPUTE SAME STEP BY BACKWARD METHOD DFDS2=-4.*S3P*F6+4.*S3P*F5-4.*S3P*F4 + +4.*S3P*F1-8.*S2P*F4+4.*S2P* + F2+4.*S2P*F1+4.*F4-F2-3.*F1 DFDS3=-8.*S3P*F6+4.*S3P*F3+4.*S3P*F1 + -4.*S2P*F6+4.*S2P*F5-4.*S2P* + F4+4.*S2P*F1+4.*F6-F3-3.*F1 GSIZE=MAX(ABS(DFDS2),ABS(DFDS3)) IF (GSIZE.EQ.0.0) GO TO 8999 DFDS2=DFDS2/GSIZE DFDS3=DFDS3/GSIZE GRADF=SQRT(DFDS2**2+DFDS3**2) GRADFX=DFDS2/GRADF GRADFY=DFDS3/GRADF ROUNDX= +GRADFY ROUNDY= -GRADFX DS2P=ROUNDX*DSTEP*COUNTR DS3P=ROUNDY*DSTEP*COUNTR C ACTUAL INTEGRATION STEP BY TRAPEZOIDAL METHOD DS2=0.5*(DS2+DS2P) DS3=0.5*(DS3+DS3P) DSLEN=SQRT(DS2**2+DS3**2) IF((DSLEN/DSTEP).LT.0.10) GO TO 8999 S2=S2+DS2 S3=S3+DS3 S1=1.00-S2-S3 C CORRECT CONTOUR TO ACTUAL VALUE DESIRED TRIAL=PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6) ERR=TRIAL-FVALUE IF (ABS(ERR).GE.DFCON) GO TO 8999 DFDS2=-4.*S3 *F6+4.*S3 *F5-4.*S3 *F4 + +4.*S3 *F1-8.*S2 *F4+4.*S2 * + F2+4.*S2 *F1+4.*F4-F2-3.*F1 DFDS3=-8.*S3 *F6+4.*S3 *F3+4.*S3 *F1 + -4.*S2 *F6+4.*S2 *F5-4.*S2 * + F4+4.*S2 *F1+4.*F6-F3-3.*F1 GSIZE=MAX(ABS(DFDS2),ABS(DFDS3)) IF (GSIZE.EQ.0.0) GO TO 8999 DFDS2=DFDS2/GSIZE DFDS3=DFDS3/GSIZE GRADF=SQRT(DFDS2**2+DFDS3**2) GRADFX=DFDS2/GRADF GRADFY=DFDS3/GRADF DISTNC= -ERR/(GRADF*GSIZE) IF (ABS(DISTNC).GT.DSTEP) DISTNC= + DISTNC*DSTEP/ABS(DISTNC) S2=S2+DISTNC*GRADFX S3=S3+DISTNC*GRADFY S1=1.00-S2-S3 C DECIDE WHETHER CONTOUR IS FINISHED OR NOT HITLIM=NSEG.GE.LIMINT IF (HITLIM) WRITE(6,3501)FVALUE,I 3501 FORMAT(' ',1PE10.2,' CONTOUR IN ELEMENT ',I3, + ' SEEMS TO BE IN LOOP. TERMINATED.') GONOUT=(S1.LT.0..OR.S1.GT.1.).OR. + (S2.LT.0..OR.S2.GT.1.).OR. + (S3.LT.0..OR.S3.GT.1.) FINISH=GONOUT.OR.HITLIM IF (CENTER) THEN XT=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) YT=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) ANGLEP=ATAN2((YT-YEXT),(XT-XEXT)) DROT=MIN(ABS(ANGLEP-ANGLE), & 6.2832-ABS(ANGLEP-ANGLE)) ROT=ROT+DROT CIRCLE=ROT.GE.6.2832 FINISH=FINISH.OR.CIRCLE IF (CIRCLE.AND.INSIDE) THEN S1=PS(1,N) S2=PS(2,N) S3=PS(3,N) ENDIF ENDIF C IF VECTOR EXTENDS OUTSIDE OF THE ELEMENT, SHORTEN IT ....... IF (GONOUT) THEN RAT=1.0 IF(S1.GT.1.)RAT=AMIN1(RAT,((1.-S1OLD)/(S1-S1OLD))) IF(S2.GT.1.)RAT=AMIN1(RAT,((1.-S2OLD)/(S2-S2OLD))) IF(S3.GT.1.)RAT=AMIN1(RAT,((1.-S3OLD)/(S3-S3OLD))) IF(S1.LT.0.)RAT=AMIN1(RAT,((0.-S1OLD)/(S1-S1OLD))) IF(S2.LT.0.)RAT=AMIN1(RAT,((0.-S2OLD)/(S2-S2OLD))) IF(S3.LT.0.)RAT=AMIN1(RAT,((0.-S3OLD)/(S3-S3OLD))) RAT=AMAX1(RAT,0.0) S2=S2OLD+(S2-S2OLD)*RAT S3=S3OLD+(S3-S3OLD)*RAT S1=1.00-S2-S3 C .... AND CROSS OFF THE CORRESPONDING SIDE-CROSSING POINT IF ((N.LT.NPS).AND.(.NOT.INSIDE)) THEN XE=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) YE=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) MATE=N R2MIN=9.9E59 NP1=N+1 DO 4000 M=NP1,NPS TEST=PS(1,M)*PS(2,M)*PS(3,M) IF ((.NOT.DONE(M)).AND. + (PS(4,M).EQ.FVALUE).AND. + (TEST.EQ.0.0) ) THEN XT=PHIVAL(PS(1,M),PS(2,M),PS(3,M), + X1,X2,X3,X4,X5,X6) YT=PHIVAL(PS(1,M),PS(2,M),PS(3,M), + Y1,Y2,Y3,Y4,Y5,Y6) R2=(XT-XE)**2+(YT-YE)**2 IF(R2.LT.R2MIN) THEN MATE=M R2MIN=R2 ENDIF ENDIF 4000 CONTINUE DONE(MATE)=.TRUE. S1=PS(1,MATE) S2=PS(2,MATE) S3=PS(3,MATE) ENDIF ENDIF C LOCATE (X,Y) COORDINATES OF DESTINATION POINT X=PHIVAL(S1,S2,S3,X1,X2,X3,X4,X5,X6) Y=PHIVAL(S1,S2,S3,Y1,Y2,Y3,Y4,Y5,Y6) C STORE (X,Y) COORDINATES FOR LATER ISPLEN(ISPNUM)=ISPLEN(ISPNUM)+1 IF ((ISPPNT(ISPNUM)+ISPLEN(ISPNUM)-1).GT.NWORK) THEN WRITE(6,401) IEL GO TO 9999 ENDIF SPACE(1,ISPPNT(ISPNUM)+ISPLEN(ISPNUM)-1)=X SPACE(2,ISPPNT(ISPNUM)+ISPLEN(ISPNUM)-1)=Y C PREPARE FOR NEXT ITERATION IF ONE IS NEEDED S1OLD=S1 S2OLD=S2 S3OLD=S3 IF (.NOT.FINISH) THEN DFDS2=-4.*S3 *F6+4.*S3 *F5-4.*S3 *F4 + +4.*S3 *F1-8.*S2 *F4+4.*S2 * + F2+4.*S2 *F1+4.*F4-F2-3.*F1 DFDS3=-8.*S3 *F6+4.*S3 *F3+4.*S3 *F1 + -4.*S2 *F6+4.*S2 *F5-4.*S2 * + F4+4.*S2 *F1+4.*F6-F3-3.*F1 GSIZE=MAX(ABS(DFDS2),ABS(DFDS3)) IF (GSIZE.EQ.0.0) GO TO 8999 DFDS2=DFDS2/GSIZE DFDS3=DFDS3/GSIZE GRADF=SQRT(DFDS2**2+DFDS3**2) GRADFX=DFDS2/GRADF GRADFY=DFDS3/GRADF ROUNDX= +GRADFY ROUNDY= -GRADFX ANGLE=ANGLEP C C END LOOP OF PUSHING FORWARD ONE CONTOUR SEGMENT C ------------------------------------------ C GO TO 3000 ENDIF C PROVIDE EMERGENCY TERMINATION POINT FOR BEWILDERED CONTOURS 8999 CONTINUE C END OF CODE EXECUTED IF (SEGMENT NOT ALREADY INTEGRATED) ENDIF C CLOSE LOOP ON ALL CONTOUR SEGMENTS 9000 CONTINUE C C**************************************************************** C C BEGIN C CONNECTION OF CONTOUR SEGMENTS AND EDGE SEGMENTS TO CLOSE AREAS C 9001 LEVEL1=IUNDER(FMIN/DFCON) LEVEL2=IUNDER(FMAX/DFCON) IF (ALLPOS) LEVEL1=MAX(LEVEL1,0) DO 9900 IC=LEVEL1,LEVEL2 FCENTR=(IC+0.5)*DFCON IF (ALLPOS) THEN FCENTC=MAX(FCENTR,0.5*DFCON) ELSE FCENTC=FCENTR ENDIF N=IHUE(NCOLOR,DFCON,FMIDLE,IFLIP,FCENTC) CALL GSPAT(ICOLOR(N)) IF (N.EQ.1) THEN ICONT= -2 IAREA=1 ELSE IF (N.EQ.NCOLOR) THEN ICONT= -1 IAREA=1 ELSE ICONT=8 IAREA=0 IF (N.GE.2.AND.N.LE.4) NBLUE=NBLUE+1 IF (N.GE.6.AND.N.LE.8) NYELOW=NYELOW+1 ENDIF IF (ALLPOS.AND.FCENTR.LT.0.0) THEN ICONT=8 IAREA=0 ENDIF CALL GSCOL(7) CALL GSLW(1) CALL GSAREA(IAREA) C BUILD MENU OF RELEVANT SEGMENTS (SO THAT NONE IS USED TWICE) NMENU=0 DO 9020 I=1,ISPNUM IF(ABS((FOFSEG(I)-FCENTR)/DFCON).LE.0.75) THEN IF (NTOGO(I).GT.0) THEN NMENU=MIN(NMENU+1,40) NTOGO(I)=NTOGO(I)-1 MENU(NMENU)=I ENDIF ENDIF 9020 CONTINUE C DRAW SET OF CLOSED AREAS FROM MENU OF RELEVANT SEGMENTS 9050 IF (NMENU.LE.0) GO TO 9899 C BEGIN EACH CLOSED AREA WITH THE TOP SEGMENT IN THE MENU IDISH=1 NADD=+1 I1=ISPPNT(MENU(1)) XORIGN=SPACE(1,I1) YORIGN=SPACE(2,I1) XOR=XPLT(XORIGN) YOR=YPLT(YORIGN) CALL GSMOVE(XOR,YOR) C BEGIN INDEFINATE LOOP ON SEGMENTS IN ONE AREA C ------------------------------------- 9100 NAME=MENU(IDISH) IF (ANEDGE(NAME)) THEN IF (N.EQ.1) THEN IPEN= -1 ELSE IF (N.EQ.NCOLOR) THEN IPEN= -2 ELSE IPEN=8 ENDIF ELSE IPEN=ICONT ENDIF CALL GSCOL(IPEN) IF (NADD.EQ.1) THEN I1=ISPPNT(NAME) I2=I1+ISPLEN(NAME)-1 ELSE I2=ISPPNT(NAME) I1=I2+ISPLEN(NAME)-1 ENDIF DO 9500 I=I1,I2,NADD X=SPACE(1,I) Y=SPACE(2,I) X2=XPLT(X) Y2=YPLT(Y) CALL GSLINE(X2,Y2) 9500 CONTINUE C DROP SEGMENT FROM MENU IMMEDIATELY AFTER DRAWING NMENU=NMENU-1 DO 9510 J=IDISH,NMENU MENU(J)=MENU(J+1) 9510 CONTINUE C IF NECESSARY, STRING TOGETHER OTHER SEGMENTS TO COMPLETE AREA IF (NMENU.LE.0) GO TO 9899 RMIN=1.E60 DO 9600 J=1,NMENU XB=SPACE(1,ISPPNT(MENU(J))) YB=SPACE(2,ISPPNT(MENU(J))) RB=SQRT((X-XB)**2+(Y-YB)**2) IF (RB.LT.RMIN) THEN RMIN=RB IDISH=J NADD=+1 ENDIF XE=SPACE(1,ISPPNT(MENU(J))+ISPLEN(MENU(J))-1) YE=SPACE(2,ISPPNT(MENU(J))+ISPLEN(MENU(J))-1) RE=SQRT((X-XE)**2+(Y-YE)**2) IF (RE.LT.RMIN) THEN RMIN=RE IDISH=J NADD= -1 ENDIF 9600 CONTINUE R=SQRT((X-XORIGN)**2+(Y-YORIGN)**2) IF (R.GT.RMIN) THEN C LOOP NOT FINISHED; GET MORE SEGMENTS GO TO 9100 ELSE C LOOP ESSENTIALLY FINISHED; TIDY UP IF (R.GT.0.) THEN CALL GSLINE(XOR,YOR) ENDIF C ARE THERE ANY MORE CLOSED CURVES TO BE MADE? IF (NMENU.GT.0) GO TO 9050 ENDIF C C CLOSE LOOP ON CONTOUR LEVELS TO BE DRAWN IN ONE ELEMENT 9899 CALL GSENDA 9900 CONTINUE C C *************************************************************** C C CLOSE LOOP ON ALL ELEMENTS FGMAX=MAX(FGMAX,FMAX) FGMIN=MIN(FGMIN,FMIN) 9999 CONTINUE RETURN END C C C SUBROUTINE DOSIDE (FMAX,FMIN,DFCON,FN,N1,N2,NM,PS,NPS) C C FIND BEGINNING/END POINTS OF CONTOURS ALONG SIDE OF ELEMENT C DIMENSION FN(6),PS(5,500) ILOW=IUNDER(FMIN/DFCON) IF (FMIN.EQ.(DFCON*ILOW)) ILOW=ILOW-1 IHI=IABOVE(FMAX/DFCON) IF (FMAX.EQ.(DFCON*IHI)) IHI=IHI+1 NBTWEN=IHI-ILOW-1 IF (NBTWEN.GE.1) THEN NBASE=ILOW A=2.*FN(N1)-4.*FN(NM)+2.*FN(N2) B= -3.*FN(N1)+4.*FN(NM)-FN(N2) DO 10 K=1,NBTWEN N=K+NBASE F=N*DFCON C= FN(N1)-F IF (A.NE.0.) THEN DISC=SQRT(MAX(B**2-4.*A*C,0.)) ROOT1=(-B+DISC)/(2.*A) ROOT2=(-B-DISC)/(2.*A) OUT1=MAX(MAX(0.,-ROOT1),MAX(0.,ROOT1-1.)) OUT2=MAX(MAX(0.,-ROOT2),MAX(0.,ROOT2-1.)) IF (OUT1.LE.OUT2) THEN IF (OUT1.GT.0.10)GO TO 10 S=ROOT1 ELSE IF (OUT2.GT.0.10)GO TO 10 S=ROOT2 ENDIF ELSE IF (FN(N2).NE.FN(N1)) THEN S=(F-FN(N1))/(FN(N2)-FN(N1)) ELSE S=0. ENDIF NPS=NPS+1 PS(N1,NPS)=1.-S PS(N2,NPS)=S N3=6-N1-N2 PS(N3,NPS)=0. PS(4,NPS)=F 10 CONTINUE ENDIF RETURN END C C C SUBROUTINE DOPART (FEXT,FMAX,FMIN,DFCON,FN, + N1,N2,NM,S1,S2,PS,NPS) C C FIND CONTOUR END POINTS ALONG PART OF AN ELEMENT SIDE C DIMENSION FN(6),PS(5,500) ILOW=IUNDER(FMIN/DFCON) IF (FMIN.EQ.(DFCON*ILOW)) ILOW=ILOW-1 IHI=IABOVE(FMAX/DFCON) IF (FMAX.EQ.(DFCON*IHI)) IHI=IHI+1 NBTWEN=IHI-ILOW-1 IF (NBTWEN.GE.1) THEN NBASE=ILOW A=2.*FN(N1)-4.*FN(NM)+2.*FN(N2) B= -3.*FN(N1)+4.*FN(NM)-FN(N2) DO 20 K=1,NBTWEN N=K+NBASE F=N*DFCON C= FN(N1)-F IF (A.NE.0.) THEN DISC=SQRT(MAX(B**2-4.*A*C,0.)) ROOT1=(-B+DISC)/(2.*A) ROOT2=(-B-DISC)/(2.*A) OUT1=MAX(MAX(0.,S1-ROOT1),MAX(0.,ROOT1-S2)) OUT2=MAX(MAX(0.,S1-ROOT2),MAX(0.,ROOT2-S2)) IF (OUT1.LE.OUT2) THEN IF (OUT1.GT.0.10) GO TO 20 S=ROOT1 ELSE IF (OUT2.GT.0.10) GO TO 20 S=ROOT2 ENDIF ELSE IF (FN(N1).NE.FN(N2)) THEN S=(F-FN(N1))/(FN(N2)-FN(N1)) ELSE S=0. ENDIF NPS=NPS+1 PS(N1,NPS)=1.-S PS(N2,NPS)=S N3=6-N1-N2 PS(N3,NPS)=0. PS(4,NPS)=F 20 CONTINUE ENDIF RETURN END C C C SUBROUTINE DOLINE (FEXT,DFCON,FN,NCL,DS,S1EXT,S2EXT,S3EXT, + IHIC,ILOC,PS,NPS) C C FINDS CONTOUR START/END POINT ALONG LINE FROM A NODE TO AN EXTREMUM C DIMENSION FN(6),DS(3),PS(5,500) PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6)= + F1*(-S1+2.*S1**2)+ + F2*(-S2+2.*S2**2)+ + F3*(-S3+2.*S3**2)+ + F4*(4.*S1*S2)+ + F5*(4.*S2*S3)+ + F6*(4.*S3*S1) FMAX=AMAX1(FN(NCL),FEXT) FMIN=AMIN1(FN(NCL),FEXT) ILOW=IUNDER(FMIN/DFCON) IF (FMIN.EQ.(DFCON*ILOW)) ILOW=ILOW-1 IHI=IABOVE(FMAX/DFCON) IF (FMAX.EQ.(DFCON*IHI)) IHI=IHI+1 NBTWEN=IHI-ILOW-1 IF (NBTWEN.GE.1) THEN NBASE=ILOW FMID=PHIVAL((S1EXT-0.5*DS(1)),(S2EXT-0.5*DS(2)), + (S3EXT-0.5*DS(3)),FN(1),FN(2),FN(3), + FN(4),FN(5),FN(6)) A=2.*FN(NCL)-4.*FMID+2.*FEXT B= -3.*FN(NCL)+4.*FMID-FEXT DO 20 K=1,NBTWEN N=K+NBASE F=N*DFCON C= FN(NCL)-F IF (A.NE.0.) THEN DISC=SQRT(MAX(B**2-4.*A*C,0.)) ROOT1=(-B+DISC)/(2.*A) ROOT2=(-B-DISC)/(2.*A) OUT1=MAX(MAX(0.,-ROOT1),MAX(0.,ROOT1-1.)) OUT2=MAX(MAX(0.,-ROOT2),MAX(0.,ROOT2-1.)) IF (OUT1.LE.OUT2) THEN IF (OUT1.GT.0.10) GO TO 20 S=ROOT1 ELSE IF (OUT2.GT.0.10) GO TO 20 S=ROOT2 ENDIF ELSE IF (FEXT.NE.FN(NCL)) THEN S=(F-FN(NCL))/(FEXT-FN(NCL)) ELSE S=0. ENDIF IF (F.GE.0.) THEN IC=F/DFCON+0.1 ELSE IT= -F/DFCON+0.1 IC= -IT ENDIF IF ((IC.GT.IHIC).OR.(IC.LT.ILOC)) THEN NPS=NPS+1 PS(1,NPS)=S1EXT-(1.-S)*DS(1) PS(2,NPS)=S2EXT-(1.-S)*DS(2) PS(3,NPS)=S3EXT-(1.-S)*DS(3) PS(4,NPS)=F ENDIF 20 CONTINUE ENDIF RETURN END C C C SUBROUTINE AROUND (I,S4,S5,S6,NODES,VNODE,NUMNOD,NUMEL) C C DRAW ONE OR MORE SIDE OF AN ELEMENT C LOGICAL BAR,IMAGE,STATES,S4,S5,S6,TITLES COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE DIMENSION NODES(6,0:NUMEL),VNODE(2,NUMNOD) DIMENSION S(3),DS(3) DATA STEP/0.10/, ISTEP/10/ C C STATEMENT FUNCTIONS: PHIVAL(S1,S2,S3,F1,F2,F3,F4,F5,F6)= + F1*(-S1+2.*S1**2)+ + F2*(-S2+2.*S2**2)+ + F3*(-S3+2.*S3**2)+ + F4*(4.*S1*S2)+ + F5*(4.*S2*S3)+ + F6*(4.*S3*S1) XPLT(X)=5.5+5.5*(X-XCENTR)/(0.5*WIDE) YPLT(Y)=4.25+5.5*(Y-YCENTR)/(0.5*WIDE) C I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) I4=NODES(4,I) I5=NODES(5,I) I6=NODES(6,I) X1=VNODE(1,I1) X2=VNODE(1,I2) X3=VNODE(1,I3) X4=VNODE(1,I4) X5=VNODE(1,I5) X6=VNODE(1,I6) Y1=VNODE(2,I1) Y2=VNODE(2,I2) Y3=VNODE(2,I3) Y4=VNODE(2,I4) Y5=VNODE(2,I5) Y6=VNODE(2,I6) DO 100 ISIDE=1,3 IF (ISIDE.EQ.1.AND..NOT.S4) GO TO 100 IF (ISIDE.EQ.2.AND..NOT.S5) GO TO 100 IF (ISIDE.EQ.3.AND..NOT.S6) GO TO 100 J1=ISIDE J2=MOD(ISIDE,3)+1 DO 10 K=1,3 S(K)=0. DS(K)=0. 10 CONTINUE S(J1)=1.00 DS(J1)= -STEP DS(J2)= STEP X=PHIVAL(S(1),S(2),S(3),X1,X2,X3,X4,X5,X6) Y=PHIVAL(S(1),S(2),S(3),Y1,Y2,Y3,Y4,Y5,Y6) XP=XPLT(X) YP=YPLT(Y) CALL GSMOVE(XP,YP) DO 20 K=1,ISTEP DO 15 L=1,3 S(L)=S(L)+DS(L) 15 CONTINUE X=PHIVAL(S(1),S(2),S(3),X1,X2,X3,X4,X5,X6) Y=PHIVAL(S(1),S(2),S(3),Y1,Y2,Y3,Y4,Y5,Y6) XP=XPLT(X) YP=YPLT(Y) CALL GSLINE(XP,YP) 20 CONTINUE 100 CONTINUE RETURN END C C C SUBROUTINE MAGNIT (OUTVEC,NUMEL,OUTSCA,SOMNEG) C C CONVERT VECTOR TO SCALAR MAGNITUDE AT INTEGRATION POINTS C INCLUDES OPTION TO MAKE MAGNITUDES OF RIGHT-POINTING C VECTORS BE NEGATIVE, "UNDOING" THE EFFECT OF VPLOT ON C PRINCIPAL-AXIS "VECTORS". C LOGICAL SOMNEG DIMENSION OUTSCA(7,NUMEL),OUTVEC(2,7,NUMEL) DO 10 M=1,7 DO 9 I=1,NUMEL OUTSCA(M,I)=SQRT(OUTVEC(1,M,I)**2+ + OUTVEC(2,M,I)**2) 9 CONTINUE 10 CONTINUE IF (SOMNEG) THEN DO 20 M=1,7 DO 19 I=1,NUMEL IF(OUTVEC(1,M,I).GT.0.) OUTSCA(M,I)= + -OUTSCA(M,I) 19 CONTINUE 20 CONTINUE ENDIF RETURN END C C C SUBROUTINE MAGNIN (V,NUMNOD,CONDNS) C C CONVERT VECTOR TO SCALAR MAGNITUDE AT NODES C DIMENSION CONDNS(NUMNOD),V(2,NUMNOD) DO 10 I=1,NUMNOD CONDNS(I)=SQRT(V(1,I)**2+V(2,I)**2) 10 CONTINUE RETURN END C C C SUBROUTINE INTRVL (OUTSCA,NUMEL,CONDNS,NUMNOD,DFCON,NCONTR, + XIP,XNOD,YIP,YNOD) C C COMPUTE CONTOUR INTERVAL ROUNDED TO NEAREST 1,2,3,4,5, X 10**P C C IN THIS VERSION, ONLY PARTS OF THE PLOT THAT WILL BE IN VIEW C ARE CONSIDERED WHEN SETTING THE INTERVAL! C LOGICAL BAR,IMAGE,INSHOW,STATES,TITLES,SHONUM COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE DIMENSION CONDNS(NUMNOD),OUTSCA(7,NUMEL), + XIP(7,NUMEL),XNOD(NUMNOD), + YIP(7,NUMEL),YNOD(NUMNOD) C C STATEMENT FUNCTION: INSHOW(X,Y)=(ABS(X-XCENTR).LE.0.5*WIDE).AND. + (ABS(Y-YCENTR).LE.0.5*WIDE) C RLOW=9.9E59 RHI=-9.9E59 DO 20 M=1,7 DO 10 I=1,NUMEL IF (INSHOW(XIP(M,I),YIP(M,I))) THEN RLOW=MIN(RLOW,OUTSCA(M,I)) RHI =MAX(RHI ,OUTSCA(M,I)) ENDIF 10 CONTINUE 20 CONTINUE DO 30 I=1,NUMNOD IF (INSHOW(XNOD(I),YNOD(I))) THEN RLOW=MIN(RLOW,CONDNS(I)) RHI =MAX(RHI ,CONDNS(I)) ENDIF 30 CONTINUE GUESS=(RHI-RLOW)/NCONTR IF (GUESS.GT.0.0) THEN IZERO=IUNDER(ALOG10(GUESS)) FACTOR=GUESS/10.**IZERO IFACTR=FACTOR+0.5 IFACTR=MIN(5,IFACTR) IF (FACTOR.GT.7.) IFACTR=10 DFCON=IFACTR*10.**IZERO ELSE DFCON=1.00 ENDIF RETURN END C C C INTEGER FUNCTION IUNDER (X) C C RETURNS INTEGER .LE. X, UNLIKE INT FUNCTION C IUNDER=INT(X) IF (X.LT.(1.*IUNDER)) IUNDER=IUNDER-1 RETURN END C C C INTEGER FUNCTION IABOVE (X) C C RETURNS INTEGER .GE. X, UNLIKE INT FUNCTION C IF (X.LE.0.) THEN IABOVE=INT(X) ELSE IABOVE=INT(X) IF (X.GT.IABOVE) IABOVE=IABOVE+1 ENDIF RETURN END C C C SUBROUTINE USMAP (INPUT,DRAWST,NXYST,XST,YST) C C PLOTS OUTLINE OF STATES FROM DIGITIZED DATASET. C LOGICAL BAR,DRAW,DRAWST,IMAGE,STATES,TITLES COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE DIMENSION DRAWST(NXYST),XST(NXYST),YST(NXYST) C C STATEMENT FUNCTIONS: XPLT(X)=5.5+5.5*(X-XCENTR)/(0.5*WIDE) YPLT(Y)=4.25+5.5*(Y-YCENTR)/(0.5*WIDE) C DO 100 I=1,NXYST XP=XPLT(XST(I)) YP=YPLT(YST(I)) DRAW=DRAWST(I) IF (DRAW) THEN CALL GSLINE(XP,YP) ELSE CALL GSMOVE(XP,YP) ENDIF 100 CONTINUE RETURN END C C C SUBROUTINE ARROW (NUMNOD,VNODE,RMSVEC, + XNOD,YNOD) C C DRAWS VECTORS WITH RMS LENGTH RMSVEC INCHES FROM NODES C LOGICAL BAR,IMAGE,STATES,TITLES COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE DIMENSION VNODE(2,NUMNOD), + XNOD(NUMNOD),YNOD(NUMNOD) C STATEMENT FUNCTIONS: XPLT(X)=5.5+5.5*(X-XCENTR)/(0.5*WIDE) YPLT(Y)=4.25+5.5*(Y-YCENTR)/(0.5*WIDE) C CALL GSLW(2) CALL GSLT(0) CALL GSCOL(6) SUM=0. DO 100 I=1,NUMNOD SUM=SUM+VNODE(1,I)**2+VNODE(2,I)**2 100 CONTINUE IF (SUM.LE.0.) RETURN VFACT=RMSVEC/SQRT(SUM/NUMNOD) DO 200 I=1,NUMNOD X=XPLT(XNOD(I)) Y=YPLT(YNOD(I)) CALL GSMOVE(X,Y) DX=VFACT*VNODE(1,I) DY=VFACT*VNODE(2,I) XP=X+DX YP=Y+DY CALL GSLINE(XP,YP) AX=DX*(-.217)+DY*(-.125) AY=DX*(+.125)+DY*(-.217) CALL GSLINE(XP+AX,YP+AY) CALL GSMOVE(XP,YP) BX=DX*(-.217)+DY*(+.125) BY=DX*(-.125)+DY*(-.217) CALL GSLINE(XP+BX,YP+BY) 200 CONTINUE RETURN END C C C C SUBROUTINE EPOCH (INPUT,T,OUTPUT,NCHAR,TEXT) C C SELECT NAME OF EPOCH CONTAINING TIME = "T" MY BEFORE PRESENT C PER GEOLOGICAL SOCIETY OF AMERICA'S DNAG 1983 GEOLOGIC TIME SCALE C PARAMETER (NTIME=20) REAL*4 T,TTOP INTEGER NCHAR,NC CHARACTER*40 TEXT,LABELS DIMENSION LABELS(NTIME),NC(NTIME),TTOP(NTIME) DATA LABELS/'Holocene ', & 'Pleistocene ', & 'Late Pliocene ', & 'Early Pliocene ', & 'Late Miocene ', & 'Middle Miocene ', & 'Early Miocene ', & 'Late Oligocene ', & 'Early Oligocene ', & 'Late Eocene ', & 'Middle Eocene ', & 'Early Eocene ', & 'Late Paleocene ', & 'Early Paleocene ', & 'Late Cretaceous-Maastrich. ', & 'Late Cretaceous-Campanian ', & 'Late Cretaceous-Santonian ', & 'Late Cretaceous-Coniacian ', & 'Late Cretaceous-Turonian ', & 'Late Cretaceous-Cenomanian '/ DATA NC/8,11,13,14,12,14,13,14,15,11,13,12,14,15, & 28,27,27,27,26,28/ DATA TTOP/0.01, 1.6, 3.4, 5.3, 11.2, 16.6, & 23.7, 30.0, 36.6, 40.0, 52.0, 57.8, & 63.6, 66.4, 74.5, 84.0, 87.5, 88.5, 91.0, 97.5/ IF (ABS(T).LT.0.001) THEN TEXT='Present' NCHAR=7 ELSE IF (T.LT.0.0) THEN TEXT='Future' NCHAR=6 ELSE IF (T.GT.TTOP(NTIME)) THEN TEXT='?' NCHAR=1 ELSE DO 10 I=1,NTIME IF (T.LE.TTOP(I)) THEN TEXT=LABELS(I) NCHAR=NC(I) RETURN ENDIF 10 CONTINUE ENDIF RETURN END C C C INTEGER FUNCTION IHUE (NCOLOR,CINT,FMIDLE,IFLIP,F) C C RETURNS ORDINAL NUMBER OF COLOR ASSOCIATED WITH FUNCTION VALUE 'F' C WHEN CONTOURED WITH INTERVAL 'CINT', AND WHEN VALUE 'FMIDLE' IS C ROUNDED TO A CONTOUR LEVEL IN THE CENTER OF THE SPECTRUM. C IF IFLIP=+1, RED GOES WITH LOW VALUES AND BLUE WITH HIGH; C IF IFLIP=-1, THE SPECTRUM IS REVERSED. C OUTPUT VALUE SHOULD BE USED AS INDEX IN "ICOLOR" TO SELECT CODE C NUMBER OF THE SHADING PATTERN APPLIED. C THAT IS, ARRAY "ICOLOR" CONTAINS THE SPECTRUM DEFINITION. C C FMP=IUNDER((FMIDLE/CINT)+0.5)*CINT STEPS=IFLIP*(F-FMP)/CINT IHUE=STEPS+(NCOLOR/2.)+1.0 IHUE=MAX(IHUE,1) IHUE=MIN(IHUE,NCOLOR) RETURN END C C C SUBROUTINE TAPE (TITLE,TIME, + XNODC,XNODM,YNODC,YNODM,THNKC,CONNOD,THNKM, + GEOTHC,GEOTHM,GEOTHA,VC,VM,WC,WM,ESUMC,ESUMM, + NUMNOD,NUMEL) C C WRITES TAPE WITH ALL ARRAYS NEEDED IN ORDER TO C RESTART PROGRAM AND CONTINUE IN TIME; C ONLY ESSENTIAL INTEGRATED VARIABLES ARE WRITTEN; C PARAMETERS MUST BE RE-INPUT BY "INPUT", AND ALL C RECONSTRUCTABLE ARRAYS MUST BE RECOMPUTED. C CHARACTER*80 TITLE DIMENSION CONNOD(NUMNOD),ESUMC(2,2,7,NUMEL),ESUMM(2,2,7,NUMEL), + GEOTHA(4,7,NUMEL),GEOTHC(4,7,NUMEL),GEOTHM(4,7,NUMEL), + THNKC(NUMNOD),THNKM(NUMNOD), + VC(2,NUMNOD),VM(2,NUMNOD), + WC(NUMNOD),WM(NUMNOD), + XNODC(NUMNOD),XNODM(NUMNOD), + YNODC(NUMNOD),YNODM(NUMNOD) C 1001 FORMAT(A80) 1002 FORMAT(1P,8E9.2) 1003 FORMAT(0P,8F9.5) 1004 FORMAT(' TIME = ',1P,E10.4,' (',0P,F10.4,')') 1005 FORMAT(1P,6E13.6) 1006 FORMAT(1P,8E10.3) 1007 FORMAT(0P,F10.3,1P,3E10.3,0P,F10.3,1P,3E10.3) C WRITE(9,1001) TITLE TMA=TIME/(1.E6*365.25*24.*60.*60.) WRITE(9,1004) TIME,TMA WRITE(9,1) 1 FORMAT(10X, '(XNODC(I),I=1,NUMNOD)') WRITE(9,1005) (XNODC(I),I=1,NUMNOD) WRITE(9,2) 2 FORMAT(10X, '(XNODM(I),I=1,NUMNOD)') WRITE(9,1005) (XNODM(I),I=1,NUMNOD) WRITE(9,3) 3 FORMAT(10X, '(YNODC(I),I=1,NUMNOD)') WRITE(9,1005) (YNODC(I),I=1,NUMNOD) WRITE(9,4) 4 FORMAT(10X, '(YNODM(I),I=1,NUMNOD)') WRITE(9,1005) (YNODM(I),I=1,NUMNOD) WRITE(9,5) 5 FORMAT(10X, '(THNKC(I),I=1,NUMNOD)') WRITE(9,1006) (THNKC(I),I=1,NUMNOD) WRITE(9,6) 6 FORMAT(10X, '(CONNOD(I),I=1,NUMNOD)') WRITE(9,1006) (CONNOD(I),I=1,NUMNOD) WRITE(9,7) 7 FORMAT(10X, '(THNKM(I),I=1,NUMNOD)') WRITE(9,1006) (THNKM(I),I=1,NUMNOD) WRITE(9,8) 8 FORMAT(10X, '(((GEOTHC(I,J,K),I=1,4),J=1,7),K=1,NUMEL)') WRITE(9,1007) (((GEOTHC(I,J,K),I=1,4),J=1,7),K=1,NUMEL) WRITE(9,9) 9 FORMAT(10X, '(((GEOTHM(I,J,K),I=1,4),J=1,7),K=1,NUMEL)') WRITE(9,1007) (((GEOTHM(I,J,K),I=1,4),J=1,7),K=1,NUMEL) WRITE(9,10) 10 FORMAT(10X, '(((GEOTHA(I,J,K),I=1,4),J=1,7),K=1,NUMEL)') WRITE(9,1007) (((GEOTHA(I,J,K),I=1,4),J=1,7),K=1,NUMEL) WRITE(9,11) 11 FORMAT(10X, '((VC(I,J),I=1,2),J=1,NUMNOD)') WRITE(9,1005) ((VC(I,J),I=1,2),J=1,NUMNOD) WRITE(9,12) 12 FORMAT(10X, '((VM(I,J),I=1,2),J=1,NUMNOD)') WRITE(9,1005) ((VM(I,J),I=1,2),J=1,NUMNOD) WRITE(9,13) 13 FORMAT(10X, '(WC(I),I=1,NUMNOD)') WRITE(9,1002) (WC(I),I=1,NUMNOD) WRITE(9,14) 14 FORMAT(10X, '(WM(I),I=1,NUMNOD)') WRITE(9,1002) (WM(I),I=1,NUMNOD) WRITE(9,15) 15 FORMAT(10X,'((((ESUMC(I,J,K,L),J=1,2),I=1,2),K=1,7),L=1,NUMEL)') WRITE(9,1003)((((ESUMC(I,J,K,L),J=1,2),I=1,2),K=1,7),L=1,NUMEL) WRITE(9,16) 16 FORMAT(10X,'((((ESUMM(I,J,K,L),J=1,2),I=1,2),K=1,7),L=1,NUMEL)') WRITE(9,1003)((((ESUMM(I,J,K,L),J=1,2),I=1,2),K=1,7),L=1,NUMEL) RETURN END C C C SUBROUTINE HILITE(INDATA,ILIGHT,IVAR,NELROW,NUMNOD, + VNODE,XNODC,XNODM,YNODC,YNODM) C C USE GRAPHICS SEGMENT 4 TO HIGHLIGHT NODE #ILIGHT C ALSO, IF(SHONUM), PLOTS NODE LOCATIONS AND NUMBERS FOR OTHERS TOO. C CHARACTER*5 ASCII,HTEXT,TEMP CHARACTER*4 LABEL LOGICAL DONUMR,SHONUM EXTERNAL ASCII COMMON /PLTPRM/ BAR,IMAGE,STATES,TITLES,VFACT,XCENTR,YCENTR, + SHONUM,WIDE DIMENSION VNODE(2,NUMNOD),XNODC(NUMNOD),XNODM(NUMNOD), + YNODC(NUMNOD),YNODM(NUMNOD) C C STATEMENT FUNCTIONS: XPLT(X)=5.5+5.5*(X-XCENTR)/(0.5*WIDE) YPLT(Y)=4.25+5.5*(Y-YCENTR)/(0.5*WIDE) C CALL GSUWIN(0.0,11.0,0.0,8.5) CALL GSSEG(4) IF (SHONUM) THEN CALL GSCOL(7) CALL GSLW(1) CALL GSMS(6) CALL GSMSC(0.08) C CALL GSLSS(2,'ADMUWCRP',199) CALL GSCM(3) C CALL GSCS(199) HOWHI=YPLT(40.E5)-YPLT(0.0) DONUMR=HOWHI.GE.0.10 HOWHI=MIN(HOWHI,0.25) CALL GSQCB(W,H) CALL GSCB(HOWHI*W/H,HOWHI) CALL GSCA(1.,0.) DO 100 I=1,NUMNOD IF (IVAR.LE.2) THEN X=VNODE(1,I) Y=VNODE(2,I) ELSE IF (MOD(IVAR,2).EQ.1) THEN X=XNODM(I) Y=YNODM(I) ELSE X=XNODC(I) Y=YNODC(I) ENDIF XC=XPLT(X) YC=YPLT(Y) CALL GSMARK(XC,YC) IF (DONUMR) THEN XP=XC-2.6*HOWHI*W/H YP=YC+0.2*HOWHI WRITE(LABEL,19)I 19 FORMAT(I4) CALL GSCHAR(XP,YP,4,LABEL) ENDIF 100 CONTINUE ENDIF IF (IVAR.LE.2) THEN X=VNODE(1,ILIGHT) Y=VNODE(2,ILIGHT) ELSE IF (MOD(IVAR,2).EQ.1) THEN X=XNODM(ILIGHT) Y=YNODM(ILIGHT) ELSE X=XNODC(ILIGHT) Y=YNODC(ILIGHT) ENDIF XC=XPLT(X) YC=YPLT(Y) C QUADRANT 1 CALL GSCOL(2) CALL GSAREA(0) CALL GSMOVE(XC,YC) CALL GSLINE(XC+0.1,YC) CALL GSARC(XC,YC,90.) CALL GSLINE(XC,YC) CALL GSENDA C QUADRANT 2 CALL GSCOL(6) CALL GSAREA(0) CALL GSMOVE(XC,YC) CALL GSLINE(XC,YC+0.1) CALL GSARC(XC,YC,90.) CALL GSLINE(XC,YC) CALL GSENDA C QUADRANT 3 CALL GSCOL(4) CALL GSAREA(0) CALL GSMOVE(XC,YC) CALL GSLINE(XC-0.1,YC) CALL GSARC(XC,YC,90.) CALL GSLINE(XC,YC) CALL GSENDA C QUADRANT 4 CALL GSCOL(1) CALL GSAREA(0) CALL GSMOVE(XC,YC) CALL GSLINE(XC,YC-0.1) CALL GSARC(XC,YC,90.) CALL GSLINE(XC,YC) CALL GSENDA CALL ASDFLD(6,16,1,1,5,2) F=ILIGHT+0.1 TEMP=ASCII(F) HTEXT='H ' HTEXT(2:4)=TEMP(1:3) CALL ASCPUT(6,5,HTEXT) CALL GSSCLS RETURN END C C C SUBROUTINE DETJUP (INPUT, NUMEL,NUMNOD,NODES,VNOD, + OUTPUT,DETJ) C C PROVIDES DETJ, THE DETERMINANT OF THE C JACOBIAN MATRIX FOR THE TRANSFORMATION IN WHICH INTERNAL C POINTS OF A TRIANGLE WITH STRAIGHT SIDES ARE MAPPED INTO C NEW LOCATIONS AS SIDES BEND (BUT CORNERS STAY FIXED). C C CUSTOMIZED TO WORK WITH NODE LOCATIONS IN VECTOR STORAGE FORM, C AND WITHOUT REQUIRING SEPARATES "AREAS" ROUTINE. C DOUBLE PRECISION POINTS DIMENSION B(4),C(4),DETJ(7,NUMEL),DN(6,2), + NODES(6,0:NUMEL),POINTS(5,7), + VNOD(2,NUMNOD),X(6),Y(6) COMMON /L1L2L3/ POINTS DO 500 I=1,NUMEL DO 100 J=1,6 NODE=NODES(J,I) X(J)=VNOD(1,NODE) Y(J)=VNOD(2,NODE) 100 CONTINUE I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) B(1)=Y(2)-Y(3) B(2)=Y(3)-Y(1) B(3)=Y(1)-Y(2) B(4)=B(1) C(1)=X(3)-X(2) C(2)=X(1)-X(3) C(3)=X(2)-X(1) C(4)=C(1) AREA= 0.5*(VNOD(1,I1)*VNOD(2,I2)-VNOD(1,I2)*VNOD(2,I1) + +VNOD(1,I2)*VNOD(2,I3)-VNOD(1,I3)*VNOD(2,I2) + +VNOD(1,I3)*VNOD(2,I1)-VNOD(1,I1)*VNOD(2,I3)) AI2=1./(2.*AREA) DO 400 M=1,7 DO 200 J=1,3 DN(J,1)=AI2*B(J)*(4.*POINTS(J,M)-1.) DN(J+3,1)=AI2*4.*(B(J)*POINTS(J+1,M) + +B(J+1)*POINTS(J,M)) DN(J,2)=AI2*C(J)*(4.*POINTS(J,M)-1.) DN(J+3,2)=AI2*4.*(C(J)*POINTS(J+1,M) + +C(J+1)*POINTS(J,M)) 200 CONTINUE AJ11=0. AJ12=0. AJ21=0. AJ22=0. DO 300 J=1,6 AJ11=AJ11+DN(J,1)*X(J) AJ12=AJ12+DN(J,1)*Y(J) AJ21=AJ21+DN(J,2)*X(J) AJ22=AJ22+DN(J,2)*Y(J) 300 CONTINUE DETJAC=AJ11*AJ22-AJ12*AJ21 DETJ(M,I)=DETJAC 400 CONTINUE 500 CONTINUE RETURN END C C C SUBROUTINE FIXER(INPUT,CHANGE,CONDNS,CONDUC, + HMAX,HMIN,IVAR,NODES, + NUMEL,NUMNOD,ONEKM,OUTSCA, + OUTVEC,VNODE, + MODIFY,AREAC,AREAM, + DETJC,DETJM,DXSC,DXSM,DYSC,DYSM, + ESUMC,ESUMM,GEOTHC,GEOTHM, + LISTOP,NUMBAD, + THIKC,THIKM,THNKC,THNKM,VC,VM, + XIPC,XIPM,XNODC,XNODM,YIPC,YIPM, + YNODC,YNODM) C C MAKES CHANGES APPLIED TO THE TEMPORARY ARRAYS PERMANENT C LOGICAL CHANGE,LISTOP DOUBLE PRECISION PHI COMMON /PHITAB/ PHI DIMENSION AREAC(NUMEL),AREAM(NUMEL), + CHANGE(NUMNOD),CONDNS(NUMNOD),CONDUC(2), + DETJC(7,NUMEL),DETJM(7,NUMEL), + DXSC(6,7,NUMEL),DXSM(6,7,NUMEL), + DYSC(6,7,NUMEL),DYSM(6,7,NUMEL), + ESUMC(2,2,7,NUMEL),ESUMM(2,2,7,NUMEL), + GEOTHC(4,7,NUMEL),GEOTHM(4,7,NUMEL), + HMAX(2),HMIN(2),LISTOP(NUMEL), + NODES(6,0:NUMEL), + OUTSCA(7,NUMEL),OUTVEC(2,7,NUMEL),PHI(6,7), + THIKC(7,NUMEL),THIKM(7,NUMEL), + THNKC(NUMNOD),THNKM(NUMNOD), + VC(2,NUMNOD),VM(2,NUMNOD),VNODE(2,NUMNOD), + XIPC(7,NUMEL),XIPM(7,NUMEL), + XNODC(NUMNOD),XNODM(NUMNOD), + YIPC(7,NUMEL),YIPM(7,NUMEL), + YNODC(NUMNOD),YNODM(NUMNOD) C C 1=MANTLE GRID C IF (IVAR.EQ.1) THEN DO 101 I=1,NUMNOD XNODM(I)=VNODE(1,I) YNODM(I)=VNODE(2,I) 101 CONTINUE CALL INTERP (XNODM,NODES,NUMEL,NUMNOD,XIPM) CALL INTERP (YNODM,NODES,NUMEL,NUMNOD,YIPM) CALL AREAS (NODES,AREAM,XNODM,YNODM,NUMNOD,NUMEL) CALL DERIV (NUMEL,NUMNOD,NODES,XNODM,YNODM,AREAM, + DETJM,DXSM,DYSM,NUMBAD,LISTOP) C C 2=CRUSTAL GRID C ELSE IF (IVAR.EQ.2) THEN DO 201 I=1,NUMNOD XNODC(I)=VNODE(1,I) YNODC(I)=VNODE(2,I) 201 CONTINUE CALL INTERP (XNODC,NODES,NUMEL,NUMNOD,XIPC) CALL INTERP (YNODC,NODES,NUMEL,NUMNOD,YIPC) CALL AREAS (NODES,AREAC,XNODC,YNODC,NUMNOD,NUMEL) CALL DERIV (NUMEL,NUMNOD,NODES,XNODC,YNODC,AREAC, + DETJC,DXSC,DYSC,NUMBAD,LISTOP) C C 3=MANTLE THICKNESS C ELSE IF (IVAR.EQ.3) THEN DO 301 I=1,NUMNOD THNKM(I)=CONDNS(I) 301 CONTINUE CALL INTERP (THNKM,NODES,NUMEL,NUMNOD,THIKM) DO 303 M=1,7 DO 302 I=1,NUMEL THIKM(M,I)=MAX(THIKM(M,I),HMIN(2)) THIKM(M,I)=MIN(THIKM(M,I),HMAX(2)) 302 CONTINUE 303 CONTINUE C C 4=CRUSTAL THICKNESS C ELSE IF (IVAR.EQ.4) THEN DO 401 I=1,NUMNOD THNKC(I)=CONDNS(I) 401 CONTINUE CALL INTERP (THNKC,NODES,NUMEL,NUMNOD,THIKC) DO 403 M=1,7 DO 402 I=1,NUMEL THIKC(M,I)=MAX(THIKC(M,I),HMIN(1)) THIKC(M,I)=MIN(THIKC(M,I),HMAX(1)) 402 CONTINUE 403 CONTINUE C C 5=MANTLE HEAT-FLOW C ELSE IF (IVAR.EQ.5) THEN DO 550 I=1,NUMEL I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) I4=NODES(4,I) I5=NODES(5,I) I6=NODES(6,I) IF (CHANGE(I1).OR.CHANGE(I2).OR.CHANGE(I3).OR. + CHANGE(I4).OR.CHANGE(I5).OR.CHANGE(I6)) THEN DO 540 M=1,7 FLUX=0.0 DO 510 K=1,6 FLUX=FLUX+CONDNS(NODES(K,I))*PHI(K,M) 510 CONTINUE DTDZ=FLUX/CONDUC(1) DELDER=DTDZ-GEOTHM(2,M,I) GEOTHM(2,M,I)=DTDZ DELTB=DELDER*THIKM(M,I) DELQUA=DELTB/(THIKM(M,I)**2) GEOTHM(3,M,I)=GEOTHM(3,M,I)-DELQUA 540 CONTINUE ENDIF 550 CONTINUE C C 6=CRUSTAL HEAT-FLOW C ELSE IF (IVAR.EQ.6) THEN DO 650 I=1,NUMEL I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) I4=NODES(4,I) I5=NODES(5,I) I6=NODES(6,I) IF (CHANGE(I1).OR.CHANGE(I2).OR.CHANGE(I3).OR. + CHANGE(I4).OR.CHANGE(I5).OR.CHANGE(I6)) THEN DO 640 M=1,7 FLUX=0.0 DO 610 K=1,6 FLUX=FLUX+CONDNS(NODES(K,I))*PHI(K,M) 610 CONTINUE DTDZ=FLUX/CONDUC(1) DELDER=DTDZ-GEOTHC(2,M,I) GEOTHC(2,M,I)=DTDZ DELTB=DELDER*THIKC(M,I) DELQUA=DELTB/(THIKC(M,I)**2) GEOTHC(3,M,I)=GEOTHC(3,M,I)-DELQUA 640 CONTINUE ENDIF 650 CONTINUE C C 7=MANTLE VELOCITY C ELSE IF (IVAR.EQ.7) THEN DO 701 I=1,NUMNOD VM(1,I)=VNODE(1,I) VM(2,I)=VNODE(2,I) 701 CONTINUE C C 8=CRUSTAL VELOCITY C ELSE IF (IVAR.EQ.8) THEN DO 801 I=1,NUMNOD VC(1,I)=VNODE(1,I) VC(2,I)=VNODE(2,I) 801 CONTINUE C C 9=CRUSTAL DILATATION C ELSE IF (IVAR.EQ.9) THEN DO 950 I=1,NUMEL I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) I4=NODES(4,I) I5=NODES(5,I) I6=NODES(6,I) IF (CHANGE(I1).OR.CHANGE(I2).OR.CHANGE(I3).OR. + CHANGE(I4).OR.CHANGE(I5).OR.CHANGE(I6)) THEN DO 940 M=1,7 DIL=0.0 DO 910 K=1,6 DIL=DIL+CONDNS(NODES(K,I))*PHI(K,M) 910 CONTINUE DILNOW=ESUMM(1,1,M,I)*ESUMM(2,2,M,I)- + ESUMM(1,2,M,I)*ESUMM(2,1,M,I) IF (DILNOW.GT.0.0) THEN DILFAC=SQRT(DIL/DILNOW) ELSE DILFAC=SQRT(-DIL/DILNOW) ENDIF ESUMM(1,1,M,I)=DILFAC*ESUMM(1,1,M,I) ESUMM(1,2,M,I)=DILFAC*ESUMM(1,2,M,I) ESUMM(2,1,M,I)=DILFAC*ESUMM(2,1,M,I) ESUMM(2,2,M,I)=DILFAC*ESUMM(2,2,M,I) 940 CONTINUE ENDIF 950 CONTINUE C C 10=CRUSTAL DILATATION C ELSE IF (IVAR.EQ.10) THEN DO 1050 I=1,NUMEL I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) I4=NODES(4,I) I5=NODES(5,I) I6=NODES(6,I) IF (CHANGE(I1).OR.CHANGE(I2).OR.CHANGE(I3).OR. + CHANGE(I4).OR.CHANGE(I5).OR.CHANGE(I6)) THEN DO 1040 M=1,7 DIL=0.0 DO 1010 K=1,6 DIL=DIL+CONDNS(NODES(K,I))*PHI(K,M) 1010 CONTINUE DILNOW=ESUMC(1,1,M,I)*ESUMC(2,2,M,I)- + ESUMC(1,2,M,I)*ESUMC(2,1,M,I) IF (DILNOW.GT.0.0) THEN DILFAC=SQRT(DIL/DILNOW) ELSE DILFAC=SQRT(-DIL/DILNOW) ENDIF ESUMC(1,1,M,I)=DILFAC*ESUMC(1,1,M,I) ESUMC(1,2,M,I)=DILFAC*ESUMC(1,2,M,I) ESUMC(2,1,M,I)=DILFAC*ESUMC(2,1,M,I) ESUMC(2,2,M,I)=DILFAC*ESUMC(2,2,M,I) 1040 CONTINUE ENDIF 1050 CONTINUE ENDIF RETURN END C C C SUBROUTINE INLAND (INPUT,VNODE,NUMEL,NUMNOD, + NELROW,XNODC,YNODC, + OUTPUT,FROMW) C C COMPUTE DISTANCE INLAND FROM COLUMN-1 EDGE OF CRUSTAL GRID C NOTE: MODIFIED TO COMPUTE AT NODES, INSTEAD OF AT INT. POINTS C DIMENSION FROMW(NUMNOD),VNODE(2,NUMNOD),XNODC(NUMNOD), + YNODC(NUMNOD) NELCOL=NUMEL/(2*NELROW) NCOLN=2*NELCOL+1 NLL=NUMNOD-NCOLN+1 DO 90 I=1,NUMNOD X=VNODE(1,I) Y=VNODE(2,I) FROMW(I)= + HOWFAR(X,Y,NELROW,NCOLN,NLL,XNODC,YNODC) 90 CONTINUE RETURN END C C C REAL FUNCTION HOWFAR (X,Y, + NELROW,NCOLN,NLL,XNODC,YNODC) C C COMPUTES ORTHOGONAL DISTANCE FROM FREE EDGE OF CRUSTAL GRID. C DIMENSION XNODC(NLL),YNODC(NLL) C D2M=9.99E39 DO 10 N=1,NELROW NM=NCOLN*(2*N-1)+1 R2=(X-XNODC(NM))**2+(Y-YNODC(NM))**2 IF (R2.LT.D2M) THEN D2M=R2 NMID=NM ENDIF 10 CONTINUE NTOP=NMID-NCOLN NBOT=NMID+NCOLN R2TOP=(X-XNODC(NTOP))**2+(Y-YNODC(NTOP))**2 R2BOT=(X-XNODC(NBOT))**2+(Y-YNODC(NBOT))**2 X1=X Y1=Y IF (R2TOP.LE.R2BOT) THEN X2=XNODC(NTOP) Y2=YNODC(NTOP) X3=XNODC(NMID) Y3=YNODC(NMID) ELSE X2=XNODC(NMID) Y2=YNODC(NMID) X3=XNODC(NBOT) Y3=YNODC(NBOT) ENDIF BASE=((X2-X3)**2+(Y2-Y3)**2)**0.5 AREA=0.5*(X1*Y2-X2*Y1 + +X2*Y3-X3*Y2 + +X3*Y1-X1*Y3) HOWFAR=MAX(0.,2.*AREA/BASE) RETURN END C C C SUBROUTINE NUMBER (INPUT,LENGTH,STRING, + MODIFY,IPOINT, + OUTPUT,X) C C READ CHARACTER 'STRING' OF LENGTH 'LENGTH', C BEGINNING AT POSITION 'IPOINT' AND RETURNS FLOATING-POINT C RESULT 'X' IF SOME COMBINATION OF 0-9, ., +, AND/OR - IS C FOUND. (ELSE, RETURNS 0.) C AT RETURN, IPOINT POINTS TO THE NEXT (NON-NUMERIC) CHARACTER. C CHARACTER*80 STRING CHARACTER*1 Z LOGICAL BEGUN,POINT IF (LENGTH.GT.80) THEN WRITE(6,1) LENGTH 1 FORMAT(' VARIABLE STRING IS SUBPROGRAM NUMBER ', + 'NEEDS TO BE LONGER:',I5) STOP ENDIF X=0. BEGUN=.FALSE. POINT=.FALSE. SIGN=+1. FACTOR=1. IPOINT=MIN(IPOINT,LENGTH) JS=IPOINT DO 100 J=IPOINT,LENGTH Z=STRING(J:J) IF (ICHAR(Z).GE.240.AND.ICHAR(Z).LE.249) THEN BEGUN=.TRUE. JS=JS+1 N=ICHAR(Z)-240 IF (POINT) THEN FACTOR=FACTOR*0.10 X=X+SIGN*N*FACTOR ELSE X=X*10.+SIGN*N ENDIF ELSE IF (Z.EQ.'+') THEN IF (BEGUN) THEN GO TO 101 ELSE JS=JS+1 SIGN= +1. ENDIF ELSE IF (Z.EQ.'-') THEN IF (BEGUN) THEN GO TO 101 ELSE JS=JS+1 SIGN= -1. ENDIF ELSE IF (Z.EQ.'.') THEN IF (POINT) THEN GO TO 101 ELSE JS=JS+1 POINT=.TRUE. ENDIF ELSE IF ((.NOT.BEGUN).AND.(Z.EQ.' ')) THEN JS=JS+1 ELSE GO TO 101 ENDIF 100 CONTINUE 101 IPOINT=JS RETURN END C C C SUBROUTINE FIXDOT (INPUT,INODE,NODES,NUMNOD,NUMEL,VNODE, + MODIFY,OUTSCA,OUTVEC) C C CORRECTS TEMPORARY ARRAYS OUTSCA = DETERMINANT OF JACOBIAN C AND OUTVEC = POSITIONS OF INTEGRATION POINTS C FOLLOWING THE MOVEMENT OF A NODE. C THE NODE THAT MOVED IS NUMBER INODE, AND ALL NODE POSITIONS ARE C STORED IN VNODE. C DOUBLE PRECISION PHI,POINTS COMMON /PHITAB/ PHI COMMON /L1L2L3/ POINTS DIMENSION B(4),C(4),DN(6,2), + NODES(6,0:NUMEL), + OUTSCA(7,NUMEL), + OUTVEC(2,7,NUMEL), + PHI(6,7),POINTS(5,7), + VNODE(2,NUMNOD), + X(6),Y(6) C DO 100 I=1,NUMEL I1=NODES(1,I) I2=NODES(2,I) I3=NODES(3,I) I4=NODES(4,I) I5=NODES(5,I) I6=NODES(6,I) IF ((I1.EQ.INODE).OR.(I2.EQ.INODE).OR.(I3.EQ.INODE).OR. + (I4.EQ.INODE).OR.(I5.EQ.INODE).OR.(I6.EQ.INODE)) THEN DO 10 K=1,6 X(K)=VNODE(1,NODES(K,I)) Y(K)=VNODE(2,NODES(K,I)) 10 CONTINUE AREA= 0.5*(X(1)*Y(2)-X(2)*Y(1) + +X(2)*Y(3)-X(3)*Y(2) + +X(3)*Y(1)-X(1)*Y(3)) AI2=1./(2.*AREA) B(1)=Y(2)-Y(3) B(2)=Y(3)-Y(1) B(3)=Y(1)-Y(2) B(4)=B(1) C(1)=X(3)-X(2) C(2)=X(1)-X(3) C(3)=X(2)-X(1) C(4)=C(1) DO 90 M=1,7 OUTVEC(1,M,I)=0. OUTVEC(2,M,I)=0. DO 20 K=1,6 OUTVEC(1,M,I)=OUTVEC(1,M,I)+X(K)*PHI(K,M) OUTVEC(2,M,I)=OUTVEC(2,M,I)+Y(K)*PHI(K,M) 20 CONTINUE DO 30 J=1,3 DN(J,1)=AI2*B(J)*(4.*POINTS(J,M)-1.) DN(J+3,1)=AI2*4.*(B(J)*POINTS(J+1,M) + +B(J+1)*POINTS(J,M)) DN(J,2)=AI2*C(J)*(4.*POINTS(J,M)-1.) DN(J+3,2)=AI2*4.*(C(J)*POINTS(J+1,M) + +C(J+1)*POINTS(J,M)) 30 CONTINUE AJ11=0. AJ12=0. AJ21=0. AJ22=0. DO 40 J=1,6 AJ11=AJ11+DN(J,1)*X(J) AJ12=AJ12+DN(J,1)*Y(J) AJ21=AJ21+DN(J,2)*X(J) AJ22=AJ22+DN(J,2)*Y(J) 40 CONTINUE DETJAC=AJ11*AJ22-AJ12*AJ21 OUTSCA(M,I)=DETJAC 90 CONTINUE ENDIF 100 CONTINUE RETURN END /* //SYSLIN DD DISP=(MOD,PASS),DSN=&&LOADSET, // UNIT=VIO,SPACE=(3040,(40,40),,,ROUND), // DCB=(BLKSIZE=3040,BUFNO=1) //SYSTERM DD SYSOUT=* //SYSPRINT DD SYSOUT=* //LKED EXEC PGM=IEWL,REGION=3500K,COND=(4,LT,FORT),PARM='MAP,LIST' //SYSLIB DD DISP=(SHR,PASS),DSN=APP1.FORTVS.LIBRARY (FORTRAN) // DD DISP=(SHR,PASS),DSN=APP1.ESSLV (LINEAR EQUATION SOLVER) // DD DISP=(SHR,PASS),DSN=APP1.GDDM4.LOAD (GRAPHICS) //SYSLIN DD DSN=&&LOADSET,DISP=(OLD,DELETE) // DD DDNAME=SYSIN //SYSLMOD DD DISP=(NEW,CATLG),UNIT=DATA, CREATE LOAD MODULE // SPACE=(TRK,(12,4,1)),DSN=EFF9GPB.FIXERMOD(FIX) //SYSIN DD * ENTRY MAIN INCLUDE SYSLIB(FSINN) INCLUDE SYSLIB(ADMLSYS1) /* //SYSPRINT DD SYSOUT=* //SYSUT1 DD UNIT=VIO,SPACE=(TRK,(5,5),,,ROUND) //