C P R O G R A M 4D-VCM C C PURPOSE:COMPUTES A LEAST SQUARES UPLIFT SURFACE FIT TO RELEVELLED C SEGMENTS AND/OR TIDE GAUGE UPLIFTS SCATTERED IN SPACE AND TIME. C THE SURFACE IS OF NX DEGREE IN X,NY DEGREE IN Y AND NPOT DE- C GREE IN TIME,AND MAY CONTAIN PRESCRIBED EPISODES OF ACCELER- C ATED UPLIFT (SUBSIDENCE) THE MOVEMENTS CAN BE EITHER RELATIVE C TO THE CENTROID OF SEGMENTS OR ABSOLUTE,REFERED TO TIDE GAUGES C .THE SURFACE IS COMPUTED FOR ANY NUMBER OF PRESCRIBED DATES. C C C DESCRIPTION OF INPUT PARAMETERS C C SPECIFIED IN PROGRAM : C IR= READ CODE OF COMPUTER C IW= PRINT CODE OF COMPUTER C STD1ST= ASSUMED STANDARD DEVIATION OF 1KM OF 1-ST ORDER LEVEL. C STEPFI= SPACING OF PREDICTION POINTS IN PHI C STEPLA= SPACING OF PREDICTION POINTS IN LAMBDA C C INPUT DATA C INDEX2= STATISTICAL FILTERING OF INSIGNIFICANT FOURIER COEF- C FICIENTS.PROVIDES THE LEVEL OF SIGNIFICANCE IN MULTIPL- C ES OF STANDARD DEVIATION. C INDEX3= FLAG INDICATING WHETHER INPUT HEIGHT DIFFERENCES ARE IN C FEET (1) OR IN METRES (0). C INDEX4= FLAG INDICATING WHETHER WEIGHTS (1) OR WEIGHT FACTORS C (0) ARE INPUT C NX= DEGREE OF SURFACE IN X C NY= DEGREE OF SURFACE IN Y C NPOT= DEGREE OF SURFACE IN TIME. MUST BE GE.1 C YFS= STIPULATED DATE OF ZERO UPLIFT (YEAR OF FLAT SURFACE) C TE(1,I),TE(2,I)= PRESPECIFIED DATES OF BEGINNING AND END OF UP C TO NEPI=12 EPISODES IN MOVEMENTS.THESE ARE C OPTIONAL AND MAY BE LEFT OUT. C TPRD(I)= UP TO NTPRD=50 DATES FOR WHICH PREDICTION IS SOUGHT. C THESE ARE OPTIONAL.DEFAULT VALUE IS YFS + 100. C ISEPLT = FLAG INDICATING IF PLOT OF SEGMENTS IS DESIRED.IF C ISEPLT=1 SEGMENTS GET PLOTTED C IMAP= FLAG INDICATING IF PLOT OF UPLIFT MAPS IS DESIRED.IF IMAP C = 1,MAPS GET PLOTTED. C IPUNCH = FLAG INDICATING IF PUNCHED OUTPUT OF PREDICTIONS IS C DESIRED. IF IPUNCH = 1, PREDICTIONS ARE PUNCHED. C C THE FOLLOWING INFORMATION IS OPTIONAL.IF NOT SUPPLIED C THE PROGRAM PICKS DEFAULT VALUES C THE DEFAULT VALUES FOR XINCH AND YINCH (USED FOR C SEGMENT PLOT) ARE THE SAME AS THE WIDTH AND HEIGHT OF C THE UPLIFT CONTOUR MAPS. C C XINCH= WIDTH OF PLOTS IN INCHES. C YINCH= HEIGHT OF PLOTS IN INCHES. C XINCH AND YINCH WILL BE CALCULATED RELATIVE TO THE C AREA INCLUDED IN THE PLOTS I.E. MAXIMUM - MINIMUM PHI, C MAXIMUM - MINIMUM LAMBDA. C CONINT= CONTOUR INTERVAL FOR PLOTS IN CENTIMETRES PER CENTURY. C DEFAULT VALUE = 10 CM/CENT. C VCMMIN= LOWEST CUT-OFF CONTOUR FOR PLOTTING IN CM/CENT.DEFAULT C VALUE = -200 CM/CENT. C VCMMAX= LARGEST CUT-OFF CONTOUR FOR PLOTTING IN CM/CENT.DEFAULT C VALUE = 200 CM/CENT. C STDMIN= LOWEST CUT-OFF CONTOUR FOR PLOTTING STANDARD DEVIATION C SURFACE IN CM/CENT. DEFAULT VALUE = 0 CM/CENT. C STDMAX= LARGEST CUT-OFF CONTOUR FOR PLOTTIG STANDARD DEVIATION C SURFACE IN CM/CENT.DEFAULT VALUE = 200 CM/CENT. C C THE FOLLOWING INFORMATION PERTAINS TO A POLYGON SPECIFYING THE C AREA FROM WHICH DATA AND TIDE GAUGES IS DESIRED. C N = THE NUMBER OF POINTS IN THE POLYGON. THIS MUST BE A C CLOSED POLYGON SO THE FIRST AND LAST SETS OF COORDINATES C MUST BE THE SAME. C FI = LATITUDE OF POLYGON POINTS IN DECIMAL DEGREES. C YL = LONGITUDE OF POLYGON POINTS IN DECIMAL DEGREES. C C THE FOLLOWING INFORMATION PERTAINS TO INFORMATION FOR C NTID TIDE GAUGES PICKED FROM WITHIN THE SPECIFIED C POLYGON. THE LATITUDE AND LONGITUDE OF THE TIDE GAUGES C ARE READ FROM LOGICAL FILE 11, ALONG WITH THE UPLIFT C AND STANDARD DEVIATION FOR THE YEAR 2000. THIS DATE C (TTGU(I)), MUST BE SPECIFIED IN THE INPUT DECK. C PREDICTED UPLIFTS FROM PREVIOUS RUNS MAY ALSO BE USED C AS TIDE GAUGES. THESE ARE STORED ON LOGICAL FILE 12 C IN THE SAME FORMAT AS THE REAL TIDE GAUGES. A FLAG C (NPT) INDICATES WHETHER THESE PREDICTED TIDE GAUGES C ARE DESIRED OR NOT. C C PREDICTIONS ARE STORED ON LOGICAL FILE 13 FOR C A NARROW STRIP STRETCHING FROM B.C. TO ONTARIO. C THESE PREDICTIONS ARE USED AS TIDE GAUGES FOR ALL RUNS. C C GFD(I)= LATITUDE OF TIDE GAUGE - DEGREES C GFM(I)= LATITUDE OF TIDE GAUGE - MINUTES C GLD(I)= LONGITUDE OF TIDE GAUGE - DEGREES C GLM(I)= LONGITUDE OF TIDE GAUGE - MINUTES C TTGU(J)= UP TO ND=18 DATES FOR WHICH TIDE GAUGE UPLIFT IS INPUT C UTG(I,J)= UPLIFT OF I-TH TIDE GAUGE FOR J-TH DATE.MUST BE C TAKEN RELATIVE TO YFS,I.E. UTG(I) FOR TTGU(J)=YFS C MUST EQUAL TO ZERO. C SIGUTG(I,J)= ITS STANDARD DEVIATION C C C THE FOLLOWING INFORMATION PERTAINS TO THE LEVELLING DATA C ALSO PICKED FROM WITHIN THE SPECIFIED POLYGON. DATA C FOR THESE RELEVELLED SEGMENTS IS READ FROM LOGICAL FILE C 10. C C T(1,I)= DATE OF FIRST LEVELLING C DHT1(I)= OBSERVED HEIGHT DIFFERENCE BETWEEN A AND B (HEIGHT B C MINUS HEIGHT A) AT T(1,I) C T(2,I)= DATE OF SECOND LEVELLING C DHT2(I)= OBSERVED HEIGHT DIFFERENCE AT T(2,I) C AFD(I)= LATITUDE OF A - DEGREES C AFM(I)= LATITUDE OF A - MINUTES C ALD(I)= LONGITUDE OF A - DEGREES C ALM(I)= LONGITUDE OF A - MINUTES C BFD(I)= LATITUDE OF B - DEGREES C BFM(I)= LATITUDE OF B - MINUTES C BLD(I)= LONGITUDE OF B - DEGREES C BLM(I)= LONGITUDE OF B - MINUTES C W(I)= ACTUAL WEIGHT OF THE HEIGHT DIFFERENCE DIFFERENCE (DHT2- C DHT1) OR WEIGHT FACTOR (2 FOR TWO FIRST ORDERS, 5 FOR C FIRST AND SECOND ORDERS, 8 FOR TWO SECOND ORDERS, 17 FOR C FIRST AND THIRD ORDERS, 20 FOR SECOND AND THIRD ORDERS, C 32 FOR TWO THIRD ORDERS). C NO(I)= SEQUENCE NUMBER OF SEGMENT C C C IFLAGT = FLAG INDICATING WHETHER PREDICTIONS FROM THIS C RUN SHOULD BE ADDED TO THE PREDICTION FILE(12). C C C C DESCRIPTION OF FILES USED BY THE PROGRAM C C FILE 10 = THE INPUT LEVELLING DATA C FILE 11 = THE INPUT TIDE GAUGES C FILE 12 = THE OUTPUT PREDICTIONS . ALSO USED FOR INPUT C TO LATER RUNS IF NPT = 1. C FILE 13 = A NARROW BAND OF PREDICTIONS STRETCHING ACROSS C THE COUNTRY. AS OF NOW, THIS FILE IS ALWAYS USED AS C INPUT. C FILE 23 = THE PLOT FILE ON WHICH THE PLOT IS STORED. C C C DATA DECK ASSEMBLY INSTRUCTIONS C C CARD NUMBER COLUMNS VARIABLE FORMAT UNITS C C 1 C 1 - 2 INDEX2 I2 C 3 - 4 INDEX3 I2 C 5 - 6 INDEX4 I2 C 2 C 1 - 2 NX I2 C 3 - 4 NY I2 C 5 - 6 NPOT I2 C 7 - 16 YFS F10.3 YEARS C C THE FOLLOWING NEPI CARDS ARE OPTIONAL AND MAY BE LEFT OUT C C 3,...,2+NEPI C 1 - 10 TE(1,I) F10.3 YEARS C 11 - 20 TE(2,I) F10.3 YEARS C C 3+NEPI BLANK C C C THE FOLLOWING NTPRD CARDS ARE OPTIONAL AND MAY BE LEFT OUT C C C 4+NEPI,...,3+NEPI+NTPRD C 1- 10 TPRD(I) F10.3 YEARS C C 4+NEPI+NTPRD BLANK C C 5+NEPI+NTPRD C 1 - 2 ISEPLT I2 C 3 - 4 IMAP I2 C 5 - 6 IPUNCH I2 C C THE FOLLOWING CARD IS OPTIONAL.MAY BE REPLACED BY A BLANK CARD C 6+NEPI+NTPRD C C 1 - 8 XINCH F8.2 INCHES C 9 - 16 YINCH F8.2 INCHES C 17 - 24 CONINT F8.2 CM/CENT C 25 - 32 VCMMIN F8.2 CM/CENT C 33 - 40 VCMMAX F8.2 CM/CENT C 41 - 48 STDMIN F8.2 CM/CENT C 49 - 56 STDMAX F8.2 CM/CENT C C C 7+NEPI+NTPRD C 1-5 N I5 C FOR N<=8 , 1 CARD OF THE FOLLOWING FORMAT IS REQUIRED, FOR N>8, C EXTRA CARDS ARE NEEDED. COORDINATES OF 8 POINTS WILL FIT ON EACH C CARD. C 8+NEPI+NTPRD,...,8+NEPI+NTPRD+N/16 C 1-5 FI(1) F5.1 DECIMAL DEG. C 6-10 YL(1) F5.1 DECIMAL DEG. C . . . C . . . C . . . C 71-75 FI(8) F5.1 DECIMAL DEG. C 75-80 YL(8) F5.1 DECIMAL DEG. C C 9+NEPI+NTPRD+N/16 C 1-5 NPT I5 C C 10+NEPI+NTPRD+N/16,...,10+NEPI+NTPRD+N/16+ND C 1 - 10 TTGU(J) F10.3 YEARS C C 11+NEPI+NTPRD+N/16+ND BLANK C C 12+NEPI+NTPRD+N/16+ND C 1-2 IFLAGT I2 C C C EMERGENCY STOPS C C 0 TIME ASSIGNED TO TIDE GAUGE UPLIFT COINCIDES WITH YFS C 1 NO TIMES ASSIGNED TO TIDE GAUGE UPLIFTS C 2 NON-POSITIVE DEGREES OF FREEDOM (IF L.LE.(NX+NY+NX*NY)*(NPOT+NEPI) C AND NTID*ND. LE. (NX+1)*(NY+1)*(NPOT+NEPI)) C 3 ALL COEFFICIENTS DISCARDED.TOO HIGH A LEVEL OF PROBABILITY C REQUESTED C 4 IF ANY OF THE MAXIMUM DIMENSIONS (MX...) IS EXCEDED C 5 TOO MANY SEGMENTS REJECTED ON PRESCRIBED REJECTION LEVEL. C C REMARKS C C 0 IF MAXIMUM DIMENSIONS SPECIFIED IN THE PROGRAM ARE FOUND IN- C SUFFICIENT,CHANGE DIMENSION STATEMENTS AND MAXIMUM DIMENSI- C ON OF VARIABLE STATEMENTS.NO OTHER CHANGES TO THE PROGRAM C ARE REQUIRED. C 1 UPLIFT SURFACE IS PREDICTED FOR POINTS ON REGULAR GRID C STEPLA*STEPFI IN THE MINIMAL RECTANGLE COVERING ALL THE INPUT C RELEVELLED SEGMENTS AND TIDE GAUGES. C 2 LOCAL (X,Y) SYSTEM HAS ORIGIN AT CENTROID OF REVELLED C SEGMENTS AND TIDE GAUGES AND SCALE = 2/((XMAX-XMIN)+(YMAX-YMIN)) C 3 IF THE LENGTH OF INPUT SEGMENT IS ZERO,IT IS SET TO 0.1 KM C 4 IF NO TIDE GAUGES GIVEN,THE UPLIFT IS REFERED TO THE CENTROID C HELD AT ZERO IN TIME C 6 OUTPUT IS SELFEXPLANATORY C 7 PROGRAM USES SUBROUTINES : C PHI,PREMAT,PREDIC,ORTHO,MATRIT,SEGMNT,DRIVER C 8 THE PROGRAM WAS WRITTEN BY PETR VANICEK IN JANUARY 1977 ON C THE BASIS OF A 3D PROGRAM DEVELOPPED AT THE UNIVERSITY OF C NEW BRUNSWICK,FREDERICTON,N.B.,CANADA C IMPLICIT REAL*8(A-H,O-Z) C REAL XTRA,YTRA REAL VCMMIN,VCMMAX,STDMIN,STDMAX REAL XINCH,YINCH,CONINT REAL CTRMAT,STDMTR REAL UPLIFT,STADEV REAL WRK C C FOLLOWING DIMENSIONS PERTAIN TO MAXIMUM NUMBER OF SEGMENTS MXLL C DIMENSION X(2,2500),Y(2,2500),T(2,2500),D(2500),W(2500),V(2500), & DIST(2500),XYS(2500),NO(2500) DIMENSION XTRA(2,2500),YTRA(2,2500) C C FOLLOWING DIMENSINS PERTAIN TO MAXIMUM NUMBER OF PREDICTED POINTS C MXNPRD C DIMENSION XP(3600),YP(3600),PU(3600),SIGP(3600) DIMENSION WRK(3600),IDUM1(1) C C FOLLOWING DIMENSIONS PERTAIN TO MAXIMUM DIMENSION OF MATRIX OF PRE- C DICTED POINTS. MAXDIM = 60. NOTE THAT MAXDIM**2 = MXNPRD. C DIMENSION UPLIFT(60),STADEV(60),CTRMAT(60,60),STDMTR(60,60) C C FOLLOWING DIMENSIONS PERTAIN TO MAXIMUM NUMBER OF DATES FOR WHICH C THE UPLIFT SURFACE IS REQUIRED C DIMENSION TPRD(2) C C FOLLOWING DIMENSIONS PERTAIN TO MAIMUM NUMBER OF COEFFICIENTS REQUI- C RED MXM.MUST BE LARGER OR EQUAL TO (NX+1)*(NY+1)*(NPOT+NEPI) C DIMENSION C(30),FOUR(30),VAFO(30),XXXX(30),SC2(30),STDP(30) DIMENSION SIGMAC(30,30),ALPHA(30,30) C C FOLLOWING DIMENSIONS PERTAIN TO MAXIMUM NUMBER OF TIDE GAUGES C MXNTID AND DATES OF TIDE GAUGE UPLIFT MXND. DIMENSION TTGU(1),UTG(3000,1),SIGUTG(3000,1),XTID(3000),YTID(3000) DIMENSION NUMTID(3000) C C C FOLLOWING DIMENSIONS PERTAIN TO MAXIMUM NUMBER OF EPISODES MXNEPI C DIMENSION TE(2,1),DDD(1),FUT(1) C C C FOLLOWING DIMENSIONS PERTAIN TO MAXIMUM NUMBER OF POINTS IN C THE POLYGON SPECIFYING LIMITS OF AREA OF CONCERN. DIMENSION FI(100),YL(100),DATA(13) COMMON/BLOCK1/M,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND COMMON/CTRLIM/ VCMMIN,VCMMAX,STDMIN,STDMAX COMMON/XYINCH/ XINCH,YINCH COMMON/IDPLOT/ ID COMMON/EXTREM/ YMIN,YMAX,XMIN,XMAX COMMON/IBOX/ IBOX COMMON/MAXDIM/ MAXDIM COMMON/LIMITS/ NEPI,NN,NJ,NT COMMON/INVLIM/ VNJ,VNT COMMON / XYTID/ NDNTID,XTID,YTID,ND,NTID,L COMMON X,Y,T,TE,TTGU FC(A1,A2)=A1+A2/60 C C C MAXIMUM DIMENSIONS OF VARIABLES C MXNT=2 MXNEPI=1 MXM = 30 MAXDIM=60 MXNPRD=3600 MXLL=2500 MXNTID =3000 MXND=1 MXTPRD = 2 C C SPECIFY READ AND WRITE CODES OF COMPUTER C IR=5 IRD = 5 IW=6 C C READ FLAGS C READ(IR,1) INDEX2,INDEX3,INDEX4 1 FORMAT(3I2) C C READ MAXIMUM DEGREE OF SURFACE IN X,Y,TIME AND THE INITIAL YEAR C FOR WHICH ZERO UPLIFT IS STIPULATED (YEAR OF FLAT SURFACE). C READ(IR,2) NX,NY,NPOT,YFS 2 FORMAT(3I2,F10.3) YFSC=1D-2*YFS C C READ STARTING AND FINISHING DATES (IN YEARS) OF UP TO 12 EPISODIC C MOVEMENTS.IF NO EPISODES SUSPECTED INSERT A BLANK CARD. C DO 3 I=1,1000 IF(I.GT.MXNEPI) GO TO 200 READ(IR,4) TE(1,I),TE(2,I) 4 FORMAT(2F10.3) IF(TE(1,I).EQ.0D0) GO TO 5 TE(1,I)=1D-2*TE(1,I)-YFSC TE(2,I)=1D-2*TE(2,I)-YFSC GO TO 3 200 WRITE (IW,201) MXNEPI,I 201 FORMAT(//,5X,'AUGMENT DIMENSIONS MXNEPI (',I3,') TO AT LEAST', 1 I4,' AND TRY AGAIN') GO TO 999 3 CONTINUE 5 NEPI=I-1 C C DETERMINE NUMBER OF SELECTED FUNCTIONS. C NT=NPOT+NEPI NJ=NT*(NY+1) NN=NJ+2*NT VNJ=1D0/NJ VNT=1D0/NT IF(NT.LE.MXNT) GO TO 202 WRITE(IW,203) MXNT,NT 203 FORMAT(//,5X,'AUGMENT DIMENSIONS MXNT (',I3,') TO AT LEAST',I4,' A &ND TRY AGAIN') GOTO 999 202 CONTINUE M=(NX*NY+NX+NY)*NT + NT IF(M.LE.MXM) GO TO 204 WRITE(IW,205) MXM,NT 205 FORMAT(//,5X,'AUGMENT DIMENSIONS MXM (',I4,' ) TO AT LEAST',I5,' A &ND TRY AGAIN') GO TO 999 204 CONTINUE C C C READ UP TO 50 DESIRED DATES (IN YEARS) FOR WHICH PREDICTION IS C WANTED.TO BE FOLLOWED BY A BLANK CARD. C DO 7 I=1,1000 IF(I.GT.MXTPRD) GO TO 206 READ(IR,8) TPRD(I) 8 FORMAT(F10.3) IF(TPRD(I).EQ.0D0) GO TO 9 TPRD(I)=1D-2*TPRD(I)-YFSC GO TO 7 206 WRITE(IW,207) MXTPRD,I 207 FORMAT(//,5X,'AUGMENT DIMENSIONS MXTPRD (',I4,' ) TO AT LEAST', 1 I3,'AND TRY AGAIN') GO TO 999 7 CONTINUE 9 NTPRD=I-1 IF(NTPRD.EQ.0) TPRD(1) = 1D0 IF(NTPRD.EQ.0) NTPRD=1 C C READ SPECIFICATIONS FOR PLOTTING. C IF NO PLOTTING IS DESIRED, INSERT A BLANK CARD FOR ISEPLT, C IMAP, IPUNCH. C READ(IR,220) ISEPLT,IMAP,IPUNCH 220 FORMAT(3I2) IF(ISEPLT.EQ.1.OR.IMAP.EQ.1) CALL PLOTS(IDUM1,IDUM2) READ(IR,221) XINCH,YINCH,CONINT,VCMMIN,VCMMAX,STDMIN,STDMAX 221 FORMAT(7F8.2) IF(XINCH.LE.0.0) XINCH=10.0 IF(YINCH.LE.0.0) YINCH=10.0 IF(CONINT.LE.0.0) CONINT=10.0 IF(VCMMIN.GE.0.0) VCMMIN=-200.0 IF(VCMMAX.LE.0.0) VCMMAX=200.0 IF(STDMIN.GE.0.0) STDMIN=0.0 IF(STDMAX.LE.0.0) STDMAX=200.0 C C READ PHI AND LAMBDA OF N POINTS OF POLYGON FROM WHICH DATA IS C DESIRED C READ(5,701)N,(FI(I),YL(I),I=1,N) 701 FORMAT(I5/(16F5.1)) C C BLOCK THAT READS TIDE GAUGE INFORMATION FOLLOWS. C C WHEN TIDE GAUGES ARE PRINTED, THE NUMBER 0 BESIDE THEM INDICATES C PREDICTIONS FROM PREVIOUS RUNS, THE NUMBER -1 INDICATES C PREDICTIONS FROM THE NARROW STRIP, AND ANY OTHER NUMBER C INDICATES A REAL TIDE GAUGE, WITH THE NUMBER REPRESENTING C ITS POSITION ON THE FILE OF TIDE GAUGES. C C READ PHI , LAMBDA, UPLIFT FOR YEAR 2000,AND STD. DEV. FROM C LOGICAL FILE 11. IF TIDE GAUGE IS WITHIN POLYGON SPECIFIED, USE IT. C N1=N-1 NTID=0 DO 10 I=1,1000 READ(11,11,END=12)TFD,TFM,TLD,TLM,TUTG,TSIUTG 11 FORMAT(4F5.1,2F10.3) IF(NTID.GT.MXNTID) GO TO 208 FA=FC(TFD,TFM) YYA=FC(TLD,TLM) IN1=0 CALL INPOLY(FI,YL,N1,FA,YYA,IN1) IF(IN1.EQ.0) GO TO 10 NTID=NTID+1 NUMTID(NTID)=I XTID(NTID)=TLD*6D1+TLM YTID(NTID)=TFD*6D1+TFM UTG(NTID,1)=TUTG SIGUTG(NTID,1)=TSIUTG 10 CONTINUE 12 REWIND 13 C C READ PREDICTIONS FROM FILE 13 (NARROW BAND) AND C IF WITHIN POLYGON, USE AS TIDE GAUGES C DO 401 I=1,50000 READ(13,11,END=794)PFD,PFM,PLD,PLM,PUTG,PSIUTG IF(NTID.GT.MXNTID)GO TO 208 FB=FC(PFD,PFM) YYB=FC(PLD,PLM) IN1=0 CALL INPOLY(FI,YL,N1,FB,YYB,IN1) IF(IN1.EQ.0)GO TO 401 NTID=NTID+1 NUMTID(NTID)=-1 XTID(NTID)=PLD*6D1+PLM YTID(NTID)=PFD*6D1+PFM UTG(NTID,1)=PUTG SIGUTG(NTID,1)=PSIUTG 401 CONTINUE 794 READ(5,310)NPT 310 FORMAT(I5) IF(NPT.EQ.0)GO TO 350 C C IF NPT=1, READ PREDICTIONS FROM PREVIOUS RUNS(LOGICAL FILE 12) C AND USE AS TIDE GAUGES IF WITHIN POLYGON SPECIFIED. C REWIND 12 DO 312 I=1,50000 READ(12,11,END=350)PFD,PFM,PLD,PLM,PUTG,PSIUTG IF(NTID.GT.MXNTID) GO TO 208 FB=FC(PFD,PFM) YYB=FC(PLD,PLM) IN1=0 CALL INPOLY(FI,YL,N1,FB,YYB,IN1) IF(IN1.EQ.0) GO TO 312 NTID=NTID+1 NUMTID(NTID)=0 XTID(NTID)=PLD*6D1+PLM YTID(NTID)=PFD*6D1+PFM UTG(NTID,1)=PUTG SIGUTG(NTID,1)=PSIUTG 312 CONTINUE 208 WRITE(IW,209) MXNTID,I 209 FORMAT(//,5X,'AUGMENT DIMENSIONS MXNTID (',I4,')TO AT LEAST', 1 I4,'AND TRY AGAIN') GO TO 999 C C IF NO TIDE GAUGE USED SKIP THE REST OF THIS BLOCK C 350 IF(NTID.EQ.0)GO TO 22 WRITE(6,790)NTID 790 FORMAT(' ','NTID=',I5) C C READ UP TO 18 DATES (IN YEARS) FOR WHICH UPLIFT DERIVED FROM C TIDE GAUGES IS AVAILABLE. C DO 13 I=1,1000 READ(IR,14) TEMP 14 FORMAT(F10.3) IF(TEMP.EQ.0D0) GO TO 15 IF(I.GT.MXND) GO TO 210 TTGU(I) = 1D-2*TEMP-YFSC GO TO 13 210 WRITE(IW,211) MXND,I 211 FORMAT(//,5X,'AUGMENT DIMENSIONS MXND (',I3,') TO AT LEAST',I3, 1 ' AND TRY AGAIN') GO TO 999 13 CONTINUE 15 ND=I-1 IF(ND.EQ.1.AND.TTGU(1).EQ.0D0) WRITE(IW,170) 170 FORMAT(///,10X,'GIVE ANOTHER TIDE GAUGE UPLIFT FOR A DATE DIFFEREN &T FROM YFS') IF(ND.EQ.1.AND.TTGU(1).EQ.0D0) GO TO 999 NDNTID=ND*NTID C C CHECK IF TIDE GAUGE INFORMATION MEANINGFUL. C IF(NTID.NE.0.AND.ND.EQ.0) GO TO 16 GO TO 17 16 WRITE(IW,18) 18 FORMAT(1H,57HNO TIMES ASSIGNED TO TIDE GAUGE DERIVED UPLIFTS.TRY A &GAIN) GO TO 999 17 CONTINUE C 22 CONTINUE C C READ TIDE GAUGE DERIVED UPLIFTS (IN CM) AND THEIR STANDARD DEVIATION C D (IN CM).CONVERT DATES FROM YEARS TO CENTURIES. C DATA READ FROM LOGICAL FILE 10. INCLUDE ONLY IF WITHIN POLYGON. C L=0 702 READ(10,24,END=25)DATA,NUM 24 FORMAT(F6.1,F9.4,F7.1,F9.4,8F5.1,F3.0,I7) FA=FC(DATA(5),DATA(6)) YYA=FC(DATA(7),DATA(8)) FB=FC(DATA(9),DATA(10)) YYB=FC(DATA(11),DATA(12)) IN1=0 CALL INPOLY(FI,YL,N1,FA,YYA,IN1) IN2=0 CALL INPOLY(FI,YL,N1,FB,YYB,IN2) IF(IN1*IN2.EQ.0)GO TO 702 L=L+1 IF(L.GT.MXLL) GO TO 212 IF(DABS(D(L)).GT.0.5)WRITE(IW,51)DATA(2),DATA(4),NUM 51 FORMAT(//,10X,'THE DIFFERENCE IN HEIGHT DIFFERENCE EXCEEDS 0.5.', 1//,10X,'THE VALUES OF DHT1, DHT2, AND NO(I) ARE:',2X,F9.4,3X,F9.4, 23X,I7) IF(DABS(DATA(5)-DATA(9)).GT.5.0.OR.DABS(DATA(7)-DATA(11)).GT. & 5.0)WRITE(IW,52)DATA(5),DATA(9),DATA(7),DATA(11),NUM 52 FORMAT(//,10X,'THERE IS A 5.0 DEGREE DIFFERENCE IN LATITUDE OR L 1ONGITUDE. THE DATA WAS: ',//,10X,'AFD= ',F5.1,3X,'BFD= ',F5.1,3X, 2'ALD= ',F5.1,3X,'BLD= ',F5.1,3X,'NO(I)= ',I7) NO(L)=NUM W(L)=DATA(13) X(1,L)=DATA(7)*6D1+DATA(8) X(2,L)=DATA(11)*6D1+DATA(12) Y(1,L)=DATA(5)*6D1+DATA(6) Y(2,L)=DATA(9)*6D1+DATA(10) T(1,L)=DATA(1)*1D-2-YFSC T(2,L)=DATA(3)*1D-2-YFSC D(L)=(DATA(4)-DATA(2))*1D2 IF(INDEX3.EQ.1)D(L)=D(L)*0.3048D0 GO TO 702 212 WRITE(IW,213)MXLL,L 213 FORMAT(//,5X,'AUGMENT DIMENSIONS MXLL (',I5,') TO AT LEAST ', 1 I5,' AND TRY AGAIN') GO TO 999 25 CONTINUE WRITE(6,791)L 791 FORMAT(' ','L=',I5) C C THIS WAS THE LAST READ STATEMENT IN THE PROGRAM. C C DETERMINATION OF ORIGIN AND SCALE OF (X,Y) SYSTEM FOLLOWS. C XMAX=0D0 YMAX=0D0 XMIN=1D20 YMIN=1D20 XMEAN=0D0 YMEAN=0D0 IF(L.EQ.0) GO TO 55 DO 26 I=1,L DO 26 J=1,2 IF(X(J,I).GT.XMAX) XMAX=X(J,I) IF(Y(J,I).GT.YMAX) YMAX=Y(J,I) IF(X(J,I).LT.XMIN) XMIN=X(J,I) IF(Y(J,I).LT.YMIN) YMIN=Y(J,I) XMEAN=XMEAN+X(J,I) YMEAN=YMEAN+Y(J,I) 26 CONTINUE 55 IF(NTID.EQ.0) GO TO 57 DO 56 I=1,NTID IF(XTID(I).GT.XMAX) XMAX=XTID(I) IF(YTID(I).GT.YMAX) YMAX = YTID(I) IF(XTID(I).LT.XMIN) XMIN = XTID(I) IF(YTID(I).LT.YMIN) YMIN = YTID(I) XMEAN = XMEAN + XTID(I) YMEAN =YMEAN+YTID(I) 56 CONTINUE 57 CONTINUE C C XMEAN = XMEAN/(2*L+NTID) YMEAN = YMEAN/(2*L+NTID) Q=YMEAN*2.908882D-4 COF0=DCOS(Q) IF(L.EQ.0.AND.NTID.EQ.1) SCALE=1.0D0 IF(L.NE.0.OR.NTID.GT.1) SCALE=2.0D0/((XMAX-XMIN)+(YMAX-YMIN)*COF0) C C PLOTTING OF SEGMENTS TAKES PLACE HERE IF ISEPLT=1. DY=.33334 DX=DY/COF0 STEPFI=DY*6D1 STEPLA=DX*6D1 NDLA=(XMAX-XMIN)/STEPLA+2 NDFI=(YMAX-YMIN)/STEPFI+2 NDFIM=NDFI-1 NDLAM=NDLA-1 SY=ALOGAR(FLOAT(NDFIM)*DY/10.) XINCH=FLOAT(NDLAM)*DX/SY YINCH=FLOAT(NDFIM)*DY/SY IF(ISEPLT.EQ.1) CALL SEGMNT(X,Y,XTRA,YTRA,L,NDFI,NDLA,DY,DX) C COMPUTATION OF X,Y OF TIDE GAUGES FOLLOWS. C IF(NTID.EQ.0) GO TO 28 58 DO 29 I=1,NTID XTID(I)=(XMEAN-XTID(I))*COF0*SCALE 29 YTID(I)=(YTID(I)-YMEAN)*SCALE IF(NTID.NE.1.OR.XTID(1).NE.0D0.OR.YTID(1).NE.0D0) GO TO 28 XMEAN = XMEAN+ 0.1/(COF0*SCALE) YMEAN = YMEAN+0.1/SCALE GO TO 58 C C COMPUTATION OF X,Y OF SEGMENT ENDS FOLLOWS. C 28 CONTINUE IF(L.EQ.0) GO TO 62 DO 27 I=1,L DO 27 J=1,2 X(J,I)=(XMEAN-X(J,I))*COF0*SCALE Y(J,I)=(Y(J,I)-YMEAN)*SCALE 27 CONTINUE C C C RECAPITULATION OF INPUT PARAMETERS C 62 WRITE(IW,100) 100 FORMAT(1H1) WRITE(IW,101) 101 FORMAT(35X,'* * I N P U T P A R A M E T E R S * *',////) IF(INDEX3.EQ.1) WRITE(IW,102) IF (INDEX3.EQ.0) WRITE(IW,103) 102 FORMAT(10X,'LEVELLED HEIGHT DIFFERENCES GIVEN IN FEET',//) 103 FORMAT(10X,'LEVELLED HEIGHT DIFFERENCES GIVEN IN METRES',//) IF(INDEX4.EQ.1) WRITE(IW,104) IF(INDEX4.EQ.0) WRITE(IW,105) 104 FORMAT(10X,'FINAL WEIGHTS OF SEGMENTS SUPPLIED',///) 105 FORMAT(10X,'ONLY WEIGHT FACTORS SUPPLIED',///) WRITE(IW,106) NX,NY,NPOT 106 FORMAT(10X,'DESIRED DEGREE OF UPLIFT SURFACE : IN X =',I2,' IN Y &=',I2,' IN TIME = ',I2,///) WRITE(IW,158) YFS 158 FORMAT(10X,'DATE FOR WHICH UPLIFT EQUAL TO ZERO IS STIPULATED :', &F10.3,///) IF(NEPI.EQ.0) GO TO 107 WRITE(IW,108) NEPI 108 FORMAT(10X,I2,' ENFORCED EPISODES OF MOVEMENT WITH FOLLOWING BEGI &NNINGS AND ENDS',//) DO 109 I=1,NEPI BEG=1D2*TE(1,I)+YFS EN=1D2*TE(2,I)+YFS 109 WRITE(IW,110) I,BEG,EN 110 FORMAT(25X,I2,10X,F10.3,5X,F10.3) 107 CONTINUE WRITE(IW,111) M 111 FORMAT(///,10X,'TOTAL NUMBER OF COEFFICIENTS SOUGHT =',I3,///) WRITE(IW,112) 112 FORMAT(10X,'UPLIFT SURFACE WANTED FOR FOLLOWING DATES',//) DO 113 I=1,NTPRD Q=1D2*TPRD(I)+YFS 113 WRITE(IW,114) I,Q 114 FORMAT(25X,I2,10X,F10.3) WRITE(IW,115) 115 FORMAT(///,10X,'ORIGIN OF LOCAL (X,Y) COORDINATE SYSTEM',//) IPHOD=YMEAN/6D1 LAMOD=XMEAN/6D1 PHOM=YMEAN-DFLOAT(IPHOD)*6D1 AMBOM=XMEAN-DFLOAT(LAMOD)*6D1 WRITE(IW,116) IPHOD,PHOM,LAMOD,AMBOM 116 FORMAT(10X,'PHI = ',I3,F7.3,',LAMBDA = ',I4,F7.3,///) WRITE(IW,146) XMAX,XMIN,YMAX,YMIN,SCALE 146 FORMAT(10X,'XMAX= ',F16.8,3X,'XMIN= ',F16.8,//,10X,'YMAX= ', 1F16.8,3X,'YMIN= ',F16.8,//,10X,'SCALE= ',F16.8,//) IF(NTID.NE.0) GO TO 117 C C IF NO TIDE GAUGES INPUTTED, THEN 1 MUST BE SIMULATED WITH THE # OF C DATES EQUAL TO NT C XTID(1) = 0.0D0 YTID(1) = 0.0D0 NTID = 1 ND = NT DO 147 I =1, ND UTG(1,I) = 0.0D0 SIGUTG(1,I) = 1.0D0 147 TTGU(I) = 0.1 * I NDNTID=ND*NTID WRITE(IW,118) 118 FORMAT(10X,'NO TIDE GAUGE INFORMATION SUPPLIED. UPLIFT SURFACES CO &MPUTED RELATIVE TO ORIGIN.',///) C C END OF RECAPITULATION OF INPUT PARAMETERS. C C TEST FOR DEGREES OF FREEDOM. C 117 LNTID = L + NDNTID IF(M.LT.LNTID) GO TO 53 WRITE(IW,31) 31 FORMAT(10X,'NON-POSITIVE DEGREES OF FREEDOM.TRY AGAIN.') GO TO 999 C C FOLLOWING IS THE ASSUMED STANDARD DEVIATION OF 1 KM LINE IN 1-ST C ORDER LEVELLING (IN CM).MAY BE REPLACED IF BETTER VALUE IS KNOWN. C 53 STD1ST = 3.79D-1 Q=1.853249D0/SCALE IF(L.EQ.0) GO TO 148 DO 32 I=1,L DIST(I)=Q*DSQRT((X(2,I)-X(1,I))**2+(Y(2,I)-Y(1,I))**2) IF (DIST(I).LT.1D-1) DIST(I)=1D-1 32 CONTINUE C C FORM WEIGHT VECTOR IF ONLY WEIGHT FACTORS (2,5,8,17,20,32) WERE C INPUT INSTEAD OF THE FULL WEIGHTS,I.E. INDEX4 = 0. C IF(INDEX4.NE.0.OR.L.EQ.0) GO TO 148 DO 33 I=1,L IF(W(I).EQ.0.0D0) W(I)=1D10 IF(W(I).NE.1D10) W(I)=1D0/(DIST(I)*STD1ST**2*W(I)) 33 CONTINUE 148 K = L + 1 DO 149 I = 1,NTID DO 149 J =1 ,ND D(K) = UTG(I,J) IF(SIGUTG(I,J).EQ.0.0D0) W(K)=1D10 IF(W(K).EQ.1D10) GO TO 149 W(K) = 1.0D0 /(SIGUTG(I,J)**2) 149 K = K + 1 C C SOLVE FOR THE UNKNOWN COEFFICIENTS OF THE MODEL. C CALL ORTHO(LNTID,1D0,SIGMAF,VFAC,NFC,INDEX2,V,SIGMAC,D,W,C,FOUR, 1 VAFO,IW,ALPHA,SC2,STDP) DO 50 I=1,M DO 50 J=1,M Q=DABS(SIGMAC(I,J)) 50 IF(Q.LT.1D-10) SIGMAC(I,J)=0D0 C C CHECK IF ALL COEFFICIENTS ARE DISCARDED ON THE PRESCRIBED LEVEL C OF PROBABILITY. C IF(NFC.GT.0) GO TO 34 WRITE(IW,35) 35 FORMAT(10X,'ALL COEFFICIENTS DISCARDED.TRY AGAIN WITH LOWER LEVEL &OF PROBABILITY') GO TO 999 34 CONTINUE C C RECAPITULATION OF INPUT LEVELLING AND TIDE GAUGE DATA AND ITS RESIDUALS C WRITE(IW,100) WRITE (IW,143) 143 FORMAT(35X,'* * TIDE GAUGE DATA * *') KK = L + 1 DO 145 I=1,NTID PHIM = YMEAN + YTID(I) /SCALE AMBM = XMEAN - XTID(I) / (COF0 * SCALE) IPHID = PHIM / 6D1 LAMBD = AMBM / 6D1 PHIM = PHIM - DFLOAT(IPHID)*6D1 AMBM = AMBM - DFLOAT(LAMBD)*6D1 WRITE(IW,122)NUMTID(I),IPHID, PHIM,LAMBD,AMBM 122 FORMAT(/// 10X,I2,10X,'PHI = ',I3,F6.2,' LAMBDA = ',I4,F6.2//) WRITE(IW,123) 123 FORMAT(//,10X,'YEAR',8X,'UPLIFT',3X,'ST.DEV.',2X,'RESIDUAL'//) DO 144 J=1,ND TTT = 1D2 * TTGU(J) + YFS WRITE(IW,125) TTT,UTG(I,J),SIGUTG(I,J),V(KK) 125 FORMAT(6X,F8.2,4X,3F10.2) 144 KK = KK + 1 145 CONTINUE SCALIN=1D0/SCALE COSCIN=SCALIN/COF0 IF(L.EQ.0) WRITE(IW,159) 159 FORMAT(//35X,'* * NO LEVELLING DATA SUPPLIED * *'/) IF(L.EQ.0) GO TO 63 WRITE(IW,127) 127 FORMAT(///35X,'* * L E V E L L I N G D A T A * *',/////) WRITE(IW,128) 128 FORMAT(2X,'SEGMENT NO',4X,'LOCATION',9X,'LENGTH',4X,'AZIMUTH',5X, &'T1',6X,'T2',6X,'WEIGHT',6X,'DDH',5X,'RESIDUAL',5X,'LINK',7X,'SEQ. & NO.',//) WRITE(IW,131) 131 FORMAT(12X,'PHI',7X,'LAMBDA',6X,'KM',6X,'DEGREES',4X,'YEAR',4X,'YE &AR',18X,'CM',8X,'CM',//) DO 129 I=1,L QQ=DABS(V(I)) LINK=0 IP1=I+1 DELX=DABS(X(2,I)-X(1,IP1)) DELY=DABS(Y(2,I)-Y(1,IP1)) IF(DELX.LT.1D-6.AND.DELY.LT.1D-6) LINK=1 PHIM=YMEAN+Y(1,I)*SCALIN AMBM=XMEAN-X(1,I)*COSCIN IPHID=PHIM*1.666666D-2 LAMBD=AMBM*1.666666D-2 PHIM=PHIM-DFLOAT(IPHID)*6D1 AMBM=AMBM-DFLOAT(LAMBD)*6D1 DELX=X(2,I)-X(1,I) DELY=Y(2,I)-Y(1,I) AZIM=DATAN2(DELX,DELY) LAZIM=AZIM*57.29578D0 T1=1D2*T(1,I)+YFS T2=1D2*T(2,I)+YFS 129 WRITE(IW,130) I,IPHID,PHIM,LAMBD,AMBM,DIST(I),LAZIM,T1,T2,W(I),D(I &),V(I),LINK,NO(I) 130 FORMAT(2X,I4,4X,I3,F6.2,2X,I3,F6.2,4X,F5.1,5X,I4,4X,F6.1,2X,F6.1, &4X,F7.4,4X,F7.2,4X,F7.2,8X,I1,7X,I7) C C PRINT OUT RESULTS OF THE SURFACE FITTING. C 63 WRITE(IW,100) WRITE(IW,132) 132 FORMAT(35X,'* * S U R F A C E P A R A M E T E R S * *',////) WRITE(IW,133) INDEX2 133 FORMAT(10X,'REQUESTED LEVEL OF PROBABILITY FOR FILTERING OUT INSIG &NIFICANT FOURIER COEFFICIENTS IS',I2,'-TIMES ST.DEVIATION',//) MDIS=M-NFC WRITE(IW,134) MDIS,M 134 FORMAT(10X,I3,' OUT OF ',I3,' COEFFICIENTS WERE FILTERED OUT ON TH &E PRESCRIBED PROBABILITY LEVEL',//) SIGMAF=DSQRT(SIGMAF) WRITE(IW,135) SIGMAF 135 FORMAT(10X,'A POSTERIORI VARIANCE FACTOR BEFORE FILTERING =',F10.3 &,//) VFAC=DSQRT(VFAC) WRITE(IW,136) VFAC 136 FORMAT(10X,'A POSTERIORI VARIANCE FACTOR AFTER FILTERING =',F10.3, &//) WRITE(IW,137) 137 FORMAT(10X,'COEFFICIENTS OF ORTHOGONALISED POLYNOMIAL AFTER FILTER &ING',//,16X,'FOURIER COEFFICIENTS',8X,'STANDARD DEVIATIONS',//) DO 138 I=1,M STD=DSQRT(VAFO(I)) IF(DABS(FOUR(I)).LT.1.0D-6)FOUR(I)=0.0D0 IF(STD.LT.1.0D-6) STD = 0.0D0 138 WRITE(IW,139) I,FOUR(I),STD 139 FORMAT(18X,I3,2X,D12.5,13X,D12.5) WRITE(IW,140) 140 FORMAT(///,10X,'COEFFICIENTS OF ORIGINAL POLYNOMIAL AFTER FILTERIN &G',//) NXP1=NX+1 NYP1=NY+1 IF(NEPI.EQ.0) GO TO 151 DO 150 K=1,NEPI WRITE(IW,152) K 152 FORMAT(//,20X,'COEFFICIENTS FOR F(T) = ',I2,'-TH EPISODE',//) WRITE(IW,153) 153 FORMAT(10X,'COLUMN NO. INDICATES POWER OF X + 1, ROW NO. INDICAT &ES POWER OF Y + 1',//) DO 154 I=1,NYP1 DO154 J=1,NXP1 IWORK = K+J*NJ+I*NT-NN+NT IF(DABS(C(IWORK)).LT.1.0D-6) C(IWORK)=0.0D0 154 ALPHA(I,J)=C(IWORK) CALL MATRIT(ALPHA,NYP1,NXP1,0,IW) 150 CONTINUE 151 CONTINUE DO 155 KK=1,NPOT WRITE(IW,156) KK 156 FORMAT(//,20X,'COEFFICIENTS FOR F(T) = T**',I2,//) WRITE(IW,153) K=KK+NEPI DO 157 I=1,NYP1 DO 157 J=1,NXP1 IWORK = K+J*NJ+I*NT-NN+NT IF(DABS(C(IWORK)).LT.1.0D-6) C(IWORK)=0.0D0 157 ALPHA(I,J)=C(IWORK) CALL MATRIT(ALPHA,NYP1,NXP1,0,IW) 155 CONTINUE WRITE(IW,142) 142 FORMAT(//,10X,'COVARIANCE MATRIX OF COEFFICIENTS',//) WRITE(IW,141) 141 FORMAT(5X,'POSITION FUNCTION OF COEFFICIENTS IS L+J*NT*(NY+1)+(I-1 &)*NT, L=1,NT; J=0,NX; I=0,NY',//) CALL MATRIT(SIGMAC,M,M,1,IW) C C DETERMINATION OF X,Y OF NPRD PREDICTION POINTS. C MINLA=XMIN*1.666666D-2 MINFI=YMIN*1.666666D-2 IF(NDLA.LE.MAXDIM) GO TO 214 WRITE(IW,215) MAXDIM,NDLA 215 FORMAT(//,5X,'AUGMENT DIMENSIONS MAXDIM (',I3,') TO AT LEAST',I3,' &AND TRY AGAIN') GO TO 999 214 CONTINUE IF(NDFI.LE.MAXDIM) GO TO 216 WRITE(IW,217) MAXDIM,NDFI 217 FORMAT(//,5X,'AUGMENT DIMENSIONS MAXDIM (',I3,') TO AT LEAST',I3,' &AND TRY AGAIN') GO TO 999 216 CONTINUE STELAD=STEPLA/6.D1 STEFID=STEPFI/6.D1 DIFXD=XMIN/6.D1-MINLA DIFYD=YMIN/6.D1-MINFI IDX=DIFXD/STELAD IDY=DIFYD/STEFID ALAMIN=MINLA+IDX*STELAD PHIMIN=MINFI+IDY*STEFID XMIN=(XMEAN-ALAMIN*6.D1)*COF0*SCALE YMIN=(PHIMIN*6.D1-YMEAN)*SCALE STEPX=-STEPLA*COF0*SCALE STEPY=STEPFI*SCALE NPRD=0 DO 39 I=1,NDFI DO 39 J=1,NDLA NPRD=NPRD+1 XP(NPRD)=XMIN+(J-1)*STEPX YP(NPRD)=YMIN+(I-1)*STEPY 39 CONTINUE IF(NPRD.LE.MXNPRD) GO TO 218 WRITE(IW,219) MXNPRD,NPRD 219 FORMAT(//,5X,'AUGMENT DIMENSIONS MXNPRD (',I4,') TO AT LEAST',I4,' & AND TRY AGAIN') GO TO 999 218 CONTINUE WRITE(IW,100) C C COMPUTATION OF PREDICTED UPLIFT SURFACES FOR ALL REQUIRED DATES. C WRITE(IW,43) 43 FORMAT(//,35X,'* * P R E D I C T E D U P L I F T S * *',///) ID=9 DO 40 IT=1,NTPRD TT=TPRD(IT)*1D2+YFS CALL PREDIC(NPRD,NX,NY,NPOT,NEPI,NT,DDD,TE,XP,YP,TPRD(IT),C, 1 SIGMAC,PU,SIGP,XXXX,XYS,FUT) C C PRINT OUT PREDICTED UPLIFT FOR DATE TPRD(IT). C WRITE(IW,42) TT 42 FORMAT(10X,'UPLIFT FOR DATE = ',F8.3,//) WRITE(IW,46) 46 FORMAT(3X,'GRID NUMBERS',8X,'PHI(DEG.)',6X,'LAMBDA(DEG.)',6X,'UP &LIFT(CM)',3X,'ST.DEV.(CM)',//) IPRD=0 C REWIND 3 REWIND 8 REWIND 12 C C IF IFLAGT=1, ADD PREDICTIONS TO PREDICTION FILE(12). ONLY THOSE C WITHIN SPECIFIED POLYGON ARE PRINTED OR ADDED TO FILE. C READ(5,313)IFLAGT 313 FORMAT(I2) C DO 41 I=1,NDFI DO 47 J=1,NDLA IPRD=IPRD+1 ALAD=ALAMIN+(J-1)*STELAD PHID=PHIMIN+(I-1)*STEFID IN1=0 CALL INPOLY(FI,YL,N1,PHID,ALAD,IN1) IF(IN1.EQ.0) GO TO 777 WRITE(IW,44) I,J,PHID,ALAD,PU(IPRD),SIGP(IPRD) IF(IPUNCH.EQ.1)WRITE(7,44) I,J,PHID,ALAD,PU(IPRD),SIGP(IPRD) 44 FORMAT(5X,I3,3X,I3,10X,F6.3,10X,F7.3,10X,F9.2,5X,F7.2) IF(IFLAGT.EQ.0)GO TO 777 IALAD=ALAD TALAM=(ALAD-IALAD)*6D1 TALAD=IALAD IPHID=PHID TPHIM=(PHID-IPHID)*6D1 TPHID=IPHID WRITE(12,715)TPHID,TPHIM,TALAD,TALAM,PU(IPRD),SIGP(IPRD) 715 FORMAT(4F5.1,2F10.3) 777 K=IPRD UPLIFT(J)=PU(K) STADEV(J)=SIGP(K) 47 CONTINUE C WRITE(3) UPLIFT WRITE(8) STADEV C 41 CONTINUE C C PLOTTING OF UPLIFT SURFACES FOR TPRD(I) TAKES PLACE HERE IF IMAP = 1 C IF(IMAP.EQ.1) CALL DRIVER(WRK,NDFI,NDLA,UPLIFT,CTRMAT,STADEV, & STDMTR,CONINT,STEPFI,STEPLA) C ID=ID+2 40 CONTINUE WRITE(IW,100) WRITE(IW,45) 45 FORMAT(35X,'* * E N D O F C O M P U T A T I O N * *') 999 CONTINUE IF(ISEPLT.EQ.1.OR.IMAP.EQ.1) CALL PLOT(0.0,0.0,999) STOP END SUBROUTINE INPOLY(XPOLY,YPOLY,N1,XI,YI,IN) C C THIS SUBROUTINE RETURNS A VALUE OF IN=1 IF POINT (XI,YI) IS C WITHIN THE POLYGON SPECIFIED BY XPOLY AND YPOLY. A VALUE OF C IN=0 IS RETURNED IF THE POINT IS NOT WITHIN THE POLYGON. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION XPOLY(N1),YPOLY(N1) IN=0 DO 625 I=1,N1 IF(YI.LE.YPOLY(I))GO TO 623 622 IF(YI.GT.YPOLY(I+1)) GO TO 625 624 IF((XI-XPOLY(I) -(YI-YPOLY(I))*(XPOLY(I+1)-XPOLY(I))/(YPOLY(I+1) & -YPOLY(I))).LT.0.0) IN=1-IN GO TO 625 623 IF(YI.GT.YPOLY(I+1)) GO TO 624 625 CONTINUE RETURN END FUNCTION PHI(I,J) C C C C THIS FUNCTION COMPUTES ELEMENTS OF DESIGN MATRIX FOR A GIVEN PAIR C OF SUBSCRIPTS I,J. THE FIRST SUBSCRIPT REFERS TO THE NUMBER OF C RELEVELLED SEGMENT,OR DATE OF TIDE GAUGE UPLIFT.THE SECOND GIVES C THE SEQUENCE NUMBER OF BASE FUNCTION (I.E. J=IT+IX*NJ+IY*NT FOR C LEVELLING SEGMENTS). THE TWO MATRICES ARE CONCATENATED WITH THAT C PERTAINING TO THE LEVELLING SEGMENTS FIRST WITH THE TIDE GAUGE C MATRIX AFTERWARDS. C IMPLICIT REAL*8(A-H,O-Z) COMMON/LIMITS/ NEPI,NN,NJ,NT COMMON/INVLIM/ VNJ,VNT DIMENSION FUT(2) COMMON /IXYT/ IX,IY,IT COMMON /IKJCJJ/ II,K,JC,JJ COMMON /XYTID/ NDNTID,XTID,YTID,ND,NTID,L COMMON X(2,2500),Y(2,2500),T(2,2500),TE(2,1),TTGU(1) DIMENSION XTID(3000),YTID(3000) C JJ = J - NJ IF(JJ.LE.0) GO TO 5 IX=(JJ-1)*VNJ+1D-6 JJ=JJ-IX*NJ 1 IY = (JJ-1)*VNT + 1D-6 IT=JJ-IY*NT IX=IX+1 C C IF I EXCEEDS L(NUMBER OF SEGMENTS) THEN THE MATRIX PERTAINING TO C THE TIDE GAUGES MUST BE OUTPUT(I.E.ONE ELEMENT OF IT) C IF(I.GT.L) GO TO 8 IF(J.LE.NT) GO TO 4 IF((IX.EQ.0).AND.(IY.EQ.0).AND.(IT.NE.0)) GO TO 4 IF(IY.EQ.0) B=X(2,I)**IX-X(1,I)**IX IF(IX.EQ.0) B=Y(2,I)**IY-Y(1,I)**IY IF(IX.NE.0.AND.IY.NE.0)B=X(2,I)**IX* Y(2,I)**IY 1 - X(1,I)**IX*Y(1,I)**IY IF(IT.GT.NEPI) GO TO 3 DO 2 JJ=1,2 IF(T(JJ,I).LE.TE(1,IT)) FUT(JJ)=0D0 IF(T(JJ,I).GE.TE(2,IT)) FUT(JJ)=1D0 2 IF(T(JJ,I).GT.TE(1,IT).AND.T(JJ,I).LT.TE(2,IT)) FUT(JJ)=(T(JJ,I)- & TE(1,IT))/(TE(2,IT)-TE(1,IT)) PHI=B*(FUT(2)-FUT(1)) RETURN 3 JJ = IT-NEPI PHI=B*(T(2,I)**JJ-T(1,I)**JJ) RETURN C C FOR THE FIRST NT(NPOT + NEPI) COLUMNS OF THE LEVELLING SEGMENT C MATRIX THE VALUE RETURNED IS 0.0 . C 4 PHI = 0D0 RETURN 5 IX = -1 JJ = J GO TO 1 8 N = I - L C C THIS CODE FINDS THE APPROPRIATE TIDE GAUGE AND MATCHES IT WITH C THE CORRECT DATE. (THE DATE CHANGES FASTER THAN THE TIDE GAUGE) C DO 10 K=1,NTID LL = K*ND 10 IF(N.LE.LL) GO TO 20 WRITE(6,500) 500 FORMAT(//9X,'*** ERROR *** I IS OUT OF RANGE OF LNTID IN FUNC. P 1HI') 20 II = N -((K-1)*ND) IF(JJ.GT.NEPI) GO TO 30 IF(TTGU(II).LT.TE(1,JJ)) PHI = 0D0 IF(TTGU(II).GT.TE(2,JJ)) PHI = 1D0 IF(TTGU(II).GE.TE(1,JJ).AND.TTGU(II).LE.TE(2,JJ)) 1 PHI = (TTGU(II) -TE(1,JJ)) / (TE(2,JJ)- TE(1,JJ)) IF(XTID(K).EQ.0.0D0.AND. IX.EQ.0) XTIX = 1.0D0 IF(YTID(K).EQ.0.0D0. AND.IY.EQ.0) YTIY = 1.0D0 IF(XTID(K).NE.0.0D0.OR. IX.NE.0) XTIX = XTID(K)**IX IF(YTID(K).NE.0.0D0.OR.IY.NE.0) YTIY = YTID(K)**IY PHI = PHI*XTIX*YTIY RETURN 30 JJ = JJ - NEPI IF(XTID(K).EQ.0.0D0.AND. IX.EQ.0) XTIX = 1.0D0 IF(YTID(K).EQ.0.0D0. AND.IY.EQ.0) YTIY = 1.0D0 IF(XTID(K).NE.0.0D0.OR. IX.NE.0) XTIX = XTID(K)**IX IF(YTID(K).NE.0.0D0.OR.IY.NE.0) YTIY = YTID(K)**IY PHI = TTGU(II)**JJ*XTIX*YTIY RETURN END SUBROUTINE PREDIC(IPRD,NX,NY,NPOT,NEPI,NT,DDD,TE,X,Y,T,C,SIGC,PU, 1 SIGP,XYPR,XYS,FUT) C C C C THIS SUBROUTINE COMPUTES PREDICTED UPLIFT PU AND ITS STANDARD C DEVIATIN SIGP FOR THE SET OF DESIRED POINTS (X,Y) AT DATE T. IT USES C SUBROUTINE PREMAT. C C INPUT PARAMATERS : C IPRD= NUMBER OF POINTS FOR WHICH PREDICTION IS WANTED C NX= HIGHEST DEGREE OF SURFACE IN X C NY= HIGHEST DEGREE OF SURFACE IN Y C NPOT= HIGHEST DEGREE OF SURFACE IN TIME C NEPI= NUMBER OF EPISODES ENFORCED C NT= NPOT + NEPI C TE= DATES OF BEGINNING AND END OF EPISODES IN CENTURIES C X,Y= COORDINATES OF POINTS FOR WHICH PREDICTIOS ARE WANTED C T= DATE FOR WHICH PREDICTION IS WANTED IN CENTURIES C C= COEFFICIENTS OF SURFACE C SIGC= THEIR COVARIANCE MATRIX C P1= PREDICTED UPLIFT C SIGP= ITS STANDARD DEVIATIONS C DDD,XYPR,XYS,FUT= DUMMY VARIABLES PASSED ON FOR REASONS OF DIMENSION C ING C IMPLICIT REAL*8(A-H,O-Z) COMMON/BLOCK1/MM,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND C DIMENSION DDD(MXNEPI),TE(2,MXNEPI),X(MXNPRD),Y(MXNPRD),C(MXM), & SIGC(MXM,MXM),PU(MXNPRD), & SIGP(MXNPRD),XYPR(MXM),XYS(MXLL),FUT(MXNEPI) C C COMPUTE PREDICTIONS AND THEIR STANDARD DEVIATIONS C REWIND 2 CALL PREMAT (IPRD,NX,NY,NPOT,NEPI,TE,X,Y,T,XYPR,2,DDD) REWIND 2 DO 6 K=1,IPRD READ (2) XYPR PU(K)=0D0 SIGP(K)=0D0 DO 8 J=1,MM PU(K)=PU(K)+XYPR(J)*C(J) XYS(J)=0D0 DO 7 I=1,MM XYS(J)=XYS(J)+XYPR(I)*SIGC(I,J) 7 CONTINUE SIGP(K)=SIGP(K)+XYS(J)*XYPR(J) 8 CONTINUE SIGP(K)=DSQRT(SIGP(K)) 6 CONTINUE RETURN END C SUBROUTINE MATRIT(A,NR,NC,ISYM,IW) C C THIS SUBROUTINE PRINTS DOUBLE PRECISION MATRICES,7 COLUMNS AND NR C ROWS.IN CASE NC>7,THE NEXT COLUMNS ARE PRINTED BELOW. C C INPUT PARAMETERS: C C A= MATRIX TO BE PRINTED C NR= ITS ROW DIMENSION C NC= ITS COLUMN DIMENSION C ISYM= FLAG IDENTIFYING A SYMMETRIC MATRIX,0 FOR NON-SYMMETRIC, C 1 FOR SYMMETRIC C IW= WRITE CODE OF COMPUTER. C IMPLICIT REAL*8(A-H,O-Z) COMMON/BLOCK1/M,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND C DIMENSION A(MXM,MXM) C INIT=1 NEND=0 1 NEND=NEND+7 IF(NC.LT.NEND) NEND=NC WRITE(IW,100) (I,I=INIT,NEND) IF(ISYM.EQ.1) GO TO 2 DO 3 I=1,NR WRITE(IW,101) I,(A(I,J),J=INIT,NEND) 3 CONTINUE 4 IF(NC.EQ.NEND) RETURN INIT=INIT+7 WRITE(IW,102) GO TO 1 2 DO 5 I=INIT,NR IF(I.GT.NEND) IN=NEND IF(I.LE.NEND) IN=I 5 WRITE(IW,101) I,(A(I,J),J=INIT,IN) GO TO 4 100 FORMAT(15X,I3,6(12X,I3)) 101 FORMAT(7X,I3,7(3X,D12.5)) 102 FORMAT(////) END C SUBROUTINE ORTHO(N,SIGMA, SIGMAF,VFC,NPC,INDEX,V,SUMD,F,W,D,C, &SUMC,IW, ALPHA,SC2,STDP) C 00017390 C THIS SUBROUTINE ORTHOGONALIZES THE MATRIX PHI USING THE GRAM-SCHMIDT 00007390 C METHOD, COMPUTES THE FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX00007400 C DERIVES THE COEFFICIENTS OF PHI,COMPUTES THE VARIANCES OF THE FOURIER 00007410 C COEFFICIENTS AND THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 00007420 C INPUTS : 00007430 C 1. PHI(OPTIONAL - COULD BE FUNCTION SUBPROGRAM INSTEAD) - AN N BY M 00007440 C CONTAINING THE BASE FUNCTIONS EVALUATED FOR EACH OBSERVATION 00007450 C 2. N - THE NUMBER OF OBSERVATIONS 00007460 C 3. M - THE NUMBER OF BASE FUNCTIONS (EQUAL OR GREATER THAN 2) 00007470 C 4. W - A VECTOR OF LENGTH N CONTAINING THE COMPUTED WEIGHT FUNCTIONS00007480 C 5. F - FUNCTIONAL VALUES 00007490 C 6. SIGMA - THE A PRIORI VARIANCE FACTOR 00007500 C 9. INDEX - PERMITS OPTIONAL TEST FOR STATISTICAL SIGNIFICANCE 00007530 C OF FOURIER COEFFICIENTS.... 00007540 C IF 0,STATISTICAL TEST FOR FOURIER COEFFICIENTS ABANDONED 00007550 C IF 1,TESTS AGAINST ONE TIME ITS STANDARD DEVIATION 00007560 C IF 2,TESTS AGAINST TWICE ITS ST.DEVIATION 00007570 C IF 3,TESTS AGAINST THREE TIMES ITS ST.DEVIATION 00007580 C 10. IW - WRITE CODE OF THE COMPUTER 00007590 C OUTPUTS : 00007600 C 1. ALPHA - AN MRD BY M MATRIX CONTAINING THE ALPHA'S USED IN COMPUTI00007610 C THE ORTHOGONALIZED MATRIX AND IN COMPUTING THE COEFFICIENTS OF PH00007620 C 2. C - THE M FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX 00007630 C 3. D - THE M COEFFICIENTS OF THE INPUT MATRIX PHI 00007640 C 4. SUMC - THE VARIANCES OF THE FOURIER COEFFICIENTS 00007650 C 5. SUMD - THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 00007660 C 6. SC2 - THE SQUARES OF THE NORMS OF THE ORTHOGONALIZED MATRIX 00007670 C 7. SIGMAF - THE FOURIER POLYNOMIAL A POSTEORI VARIANCE FACTOR 00007680 C 8. V - THE N RESIDUALS 00007690 C 9. VFC - THE ORIGINAL POLYNOMIAL A POSTEORI VARIANCE FACTOR 00007700 C 10. NPC - NUMBER OF THE COEFFICIENTS OF THE ORIGINAL POLYNOMIAL 00007710 C AFTER THE STATISTICAL TEST IS PERFORMED 00007720 C 11. STDP - VECTOR AGAINST WHICH THE ABSOLUTE VALUES OF FOURIER 00007730 C COEFFICIENTS ARE TESTED 00007740 IMPLICIT REAL*8(A-H,O-Z) 00007750 COMMON/BLOCK1/M,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND C DIMENSION V(MXLL),SUMD(MXM,MXM),F(MXLL),W(MXLL),D(MXM),C(MXM), & SUMC(MXM), ALPHA(MXM,MXM),SC2(MXM),STDP(MXM) C C TEST FOR NEGATIVE DEGREES OF FREEDOM 00007800 IF (N.LT.M) GO TO 100 00007810 IF(M.EQ.1) GO TO 50 K=1 00007820 ALPHA(M,M)=1.D0 00007830 C DETERMINE THE ALPHA'S FOR COMPUTATION OF ORTHOGONALIZED MATRIX 00007840 10 DO 3 J=K,M 00007850 IF(J.NE.K) GO TO 6 00007860 ALPHA(K,K)=1.D0 00007870 GO TO 3 00007880 6 SC1=0.D0 00007890 SC2(K)=0.D0 00007900 SC3=0.D0 00007910 DO 2 I=1,N 00007920 P=PHI(I,K) IF(K.EQ.1) GO TO 4 00007950 K1=K-1 00007960 DO 5 J1=1,K1 00007970 5 P=P+ALPHA(J1,K)*PHI(I,J1) 4 SC1=SC1+W(I)*PHI(I,J)*P SC3=SC3+F(I)*W(I)*P 00008000 2 SC2(K)=SC2(K)+W(I)*P**2 00008010 ALPHA(J,K)=-SC1/SC2(K) 00008020 ALPHA(K,J)=ALPHA(J,K) 00008030 3 CONTINUE 00008040 C DETERMINE THE FOURIER COEFFICIENTS FOR THE ORTHOGONALIZED MATRIX 00008050 C(K)=SC3/SC2(K) 00008060 K=K+1 00008070 IF(M.EQ.2) GO TO 34 00008080 IF(K.LT.3) GO TO 10 00008090 C DETERMINE THE ALPHA'S USED IN COMPUTING THE COEFFICIENTS OF PHI 00008100 JK=K-1 00008110 9 JL=K 00008120 JK=JK-1 00008130 JJ=K-JK-1 00008140 DO 8 LM=1,JJ 00008150 JL=JL-1 00008160 8 ALPHA(JK,K)=ALPHA(JK,K)+ALPHA(JK,JL)*ALPHA(K,JL) 00008170 IF(JK.NE.1) GO TO 9 00008180 IF(K.LT.M) GO TO 10 00008190 C DETERMINE THE LAST FOURIER COEFFICIENT 00008200 34 SC2(K)=0.D0 00008210 SC3=0.D0 00008220 DO 7 I=1,N 00008230 P=PHI(I,K) K1=K-1 00008260 DO 1 J=1,K1 00008270 1 P=P+ALPHA(J,K)*PHI(I,J) SC2(K)=SC2(K)+W(I)*P**2 00008290 7 SC3=SC3+F(I)*W(I)*P 00008300 C(K)=SC3/SC2(K) 00008310 C DETERMINE THE COEFFICIENTS OF PHI 00008320 IDEKT=1 00008330 ICOUNT=0 00008340 1000 CONTINUE 00008350 DO 13 I=1,M 00008360 D(I)=C(I) 00008370 IF(I.EQ.M) GO TO 13 00008380 II=I+1 00008390 DO 14 J=II,M 00008400 14 D(I)=D(I)+ALPHA(I,J)*C(J) 00008410 13 CONTINUE 00008420 C COMPUTE THE VARIANCE OF THE FOURIER COEFFICIENTS AND THE VARIANCE-COVA00008430 C MATRIX OF THE COEFFICIENTS 00008440 DO 15 I=1,M 00008450 DO 15 J=1,M 00008460 15 SUMD(I,J)=0.D0 00008470 SC4=0.D0 00008480 DO 22 I=1,N 00008490 PN=0.D0 00008510 DO 21 J=1,M 00008520 21 PN=PN+D(J)*PHI(I,J) V(I)=F(I)-PN 00008540 V2=V(I)**2 00008550 22 SC4=SC4+V2*W(I) 00008560 SIGMAF=SC4/(N-M+ICOUNT)*SIGMA 00008570 VFC=SIGMAF 00008580 IF(IDEKT.EQ.2) VFC=SC4/(N-NPC)*SIGMA 00008590 IF(INDEX.EQ.0) NPC=M 00008600 DO 28 I=1,M 00008610 SUMC(I)=SIGMAF/SC2(I) 00008620 IF(IDEKT.EQ.1) GO TO 28 00008630 IF(C(I).EQ.0D0) SUMC(I)=0D0 00008640 28 CONTINUE 00008650 CALL ERRSET(208,256) DO 23 I=1,M 00008660 DO 23 J=1,I 00008670 DO 23 K=J,I 00008680 23 SUMD(J,K)=SUMD(J,K)+ALPHA(J,I)*ALPHA(K,I)*SUMC(I) 00008690 CALL ERRSET(208,1) DO 24 I=1,M 00008700 IT=I+1 00008710 IF(IT.GT.M) GO TO 30 00008720 DO 24 J=IT,M 00008730 24 SUMD(J,I)=SUMD(I,J) 00008740 C OPTIONAL CHECK FOR STATISTICALLY SIGNIFICANT FOURIER COEFFICIENTS 00008750 30 CONTINUE 00008760 IF(INDEX.EQ.0) GO TO 40 00008770 IF(IDEKT.EQ.2) GO TO 40 00008780 PINDEX=DFLOAT(INDEX) 00008790 DO 31 I=1,M 00008800 STDP(I)=PINDEX*DSQRT(SUMC(I)) 00008810 IF(DABS(C(I)).LT.STDP(I)) GO TO 32 00008820 GO TO 31 00008830 32 C(I)=0D0 00008840 ICOUNT=ICOUNT+1 00008850 SUMC(I)=0D0 00008860 31 CONTINUE 00008870 NPC=0 00008880 DO 33 I=1,M 00008890 IF(C(I).NE.0D0) NPC=I 00008900 33 CONTINUE 00008910 IDEKT=2 00008920 GO TO 1000 00008930 50 ALPHA(1,1) = 1.0D0 D(1) = 0.0D0 D(2) = 0.0D0 DO 55 I=1,N D(1) = PHI(I,M) *W(I) *F(I) + D(1) 55 D(2)=PHI(I,M) * W(I) * PHI(I,M) + D(2) D(1) = D(1) /D(2) C(1) = D(1) SC2(1) = 0.0D0 DO 60 I=1,N V(I) = PHI(I,M) * D(1) -F(I) 60 SC2(1) = SC2(1) + V(I)*W(I)*V(I) SUMC(1) = SC2(1) / ((N-1)*D(2)) SUMD(1,1) = SUMC(1) SIGMAF = SC2(1) /(N-1) VFC = SIGMAF NPC = 1 STDP(1) = SUMC(1) **0.5 40 RETURN 00008940 100 WRITE(IW,102) 00008950 102 FORMAT('0','*ERROR* NEGATIVE DEGREES OF FREEDOM') 00008960 RETURN 00008970 END 00008980 C 00020050 SUBROUTINE DRIVER (WRK,NROWS,NCOLS,UPLIFT,CTRMAT,STADEV,STDMTR, & CONINT,STEPFI,STEPLA) C C THIS SUBROUTINE SETS UP THE STAGE FOR CONTOURING THE UPLIFT AND STAN C DARD DEVIATION SURFACES GIVEN THE DATA UPLIFT AND STADEV IN MATRIX C FORM.IT USES SUBROUTINES SETLIM,CONBLK AND CONMAP. C C INPUT PARAMETERS: C C NROWS= NUMBER OF ROWS IN PLOTTED MATRICES C NCOLS= NUMBER OF COLUMNS IN PLOTTED MATRICES C UPLIFT= MATRIX OF UPLIFT TO BE PLOTTED.PASSED ON IN FILE 5 C STADEV= MATRIX OF STANDARD DEVIATION TO BE PLOTTED.PASSED ON IN FILE C 6 C CONINT= COTOUR INTERVAL IN CM/CENT C STEPFI= SPACING OF DATA POINTS IN LATITUDE IN MINS OF ARC C STEPLA= SPACING OF DATA POINTS IN LONGITUDE IN MINS OF ARC C WRK,CTRMAT,STDMTR= DUMMY VARIABLES PASSED ON FOR REASONS OF DIMENSI C ONING C DIMENSION WRK(MXNPRD),UPLIFT(MAXDIM),CTRMAT(MAXDIM,MAXDIM), &STADEV(MAXDIM),STDMTR(MAXDIM,MAXDIM) COMMON/MAXDIM/ MAXDIM COMMON/EXTREM/ YMIN,YMAX,XMIN,XMAX COMMON/IDPLOT/ ID COMMON/XYINCH/ XINCH,YINCH COMMON/CTRLIM/ VCMMIN,VCMMAX,STDMIN,STDMAX COMMON/BLOCK1/ MM,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND COMMON/MEASUR/WIDTH,HEIGHT DOUBLE PRECISION STEPFI,STEPLA DOUBLE PRECISION YMIN,YMAX,XMIN,XMAX REWIND 3 REWIND 8 DO 1 I=1,NROWS READ(3) UPLIFT READ(8) STADEV DO 1 J=1,NCOLS K=NCOLS-J+1 CTRMAT(I,J)=UPLIFT(K) STDMTR(I,J)=STADEV(K) 1 CONTINUE ENDFILE 3 ENDFILE 8 DY=STEPFI/60 DX=STEPLA/60 CALL SETLIM(VCMMIN,VCMMAX) CALL CONBLK (0,30,0,0) CALL CONMAP (CTRMAT,MAXDIM,NROWS,NCOLS,1,1,CONINT,5,1,DY,DX,WRK, & MXNPRD) CALL PLOT(WIDTH/1000.+2.,0.0,-3) ID=ID+1 CALL SETLIM (STDMIN,STDMAX) CALL CONBLK (0,30,0,0) CALL CONMAP (STDMTR,MAXDIM,NROWS,NCOLS,1,1,CONINT,5,1,DY,DX,WRK, & MXNPRD) CALL PLOT(WIDTH/1000.+2.,0.0,-3) RETURN END C***********************************************************************00020040 C 00020050 C 00020060 SUBROUTINE PLTRTN(IC,IX,IY) 00020070 C 00020080 C 00020090 C THIS SUBROUTINE REPLACES THE U.N.B. PLOTTING PACKAGES' VERSION OF00020100 C PLTRTN WHICH WOULD NORMALLY STORE A PLOT OF THE CONTOUR MAP 00020110 C PRODUCED BY CONMAP. IT SETS UP THE DATA IN THE CAM FORMAT WHICH IS 00020120 C AS FOLLOWS: 00020130 C 00020140 C COL. FORMAT VARIABLE ITEM 00020150 C 00020160 C 1-6 I6 LSNUM LINE SEGMENT NUMBER 00020170 C 7-26 E20.8 LATR LATITUDE IN RADIANS 00020180 C 27-46 E20.8 LONGR LONGITUDE IN RADIANS 00020190 C 47-50 4X BLANK 00020200 C 51-52 I2 LATD LATITUDE-DEGREE PART (ABSOLUTE VALUE) 00020210 C 53-54 I2 LATM LATITUDE-MINUTE PART (ABSOLUTE VALUE) 00020220 C 55-56 I2 LATS LATITUDE-SECOND PART (ABSOLUTE VALUE) 00020230 C 57 A1 ALAT N OR S FOR NORTH OR SOUTH 00020240 C 58-60 I3 LONGD LONGITUDE-DEGREE PART (ABSOLUTE VALUE) 00020250 C 61-62 I2 LONGM LONGITUDE-MINUTE PART (ABSOLUTE VALUE) 00020260 C 63-64 I2 LONGS LONGITUDE-SECOND PART (ABSOLUTE VALUE) 00020270 C 65 A1 ALONG E OR W FOR EAST OR WEST 00020280 C 66-71 6X BLANK 00020290 C 72-80 I9 SNUM RECORD SEQUENCE NUMBER 00020300 C 00020310 C 00020320 C PASSED VARIABLES: 00020330 C 00020340 C IC GIVES THE POSITION OF THE PLOTTING PENCIL 00020350 C IX THE X COORDINATE OF THE POINT 00020360 C IY THE Y COORDINATE OF THE POINT 00020370 C 00020380 C 00020390 C VARIABLES: 00020400 C 00020410 C ALAMDA THE DECIMAL DEGREE FORM OF THE LONGITUDE 00020420 C ALAT,ALONG N OR S AND E OR W RESPECTIVELY 00020430 C APHI THE DECIMAL DEGREE FORM OF THE LATITUDE 00020440 C CEILFD MAXIMUM PHI VALUE TO BE COVERED BY CONTOURING 00020450 C CEILLD MAXIMUM LAMBDA VALUE TO BE COVERED BY CONTOURING 00020460 C FLORFD MINIMUM PHI VALUE TO BE COVERED BY CONTOURING 00020470 C FLORLD MINIMUM LAMBDA VALUE TO BE COVERED BY CONTOURING 00020480 C IBOX USED TO INSURE THAT THE 5 VECTORS FORMING A BOX AROUND 00020490 C THE CONTOUR MAP DO NOT GET INTO THE CAM DATA BANK 00020500 C ID IDENTIFIES THE PLOT FILE TO BE WRITTEN 00020510 C HEIGHT DESIRED HEIGHT OF PLOT IN 1000THS OF INCHES 00020520 C LATD LATITUDE DEGREES(INTEGER) 00020530 C LATM LATITUDE MINUTES(INTEGER) 00020540 C LATR DECIMAL RADIAN FORM OF THE LATITUDE 00020550 C LATS LATITUDE SECONDS(INTEGER) 00020560 C LSNUM LINE SEGMENT NUMBER 00020570 C LONGD LONGITUDE DEGREES (INTEGER) 00020580 C LONGM LONGITUDE MINUTES (INTEGER) 00020590 C LONGR DECIMAL RADIAN FORM OF THE LONGITUDE 00020600 C LONGS LONGITUDE SECONDS (INTEGER) 00020610 C SNUM SEQUENCE NUMBER OF DATA GOING INTO CAM DATA BANK 00020620 C WIDTH DESIRED WIDTH OF PLOT IN 1000THS OF INCHES 00020630 C 00020640 C 00020650 C TEMPORARY VARIABLES: 00020660 C 00020670 C REM USED IN FINDING INTEGER VALUES MINUTES AND SECONDS 00020680 C SIXREM OF LATITUDE AND LONGITUDE 00020690 C 00020700 C ENTRY POINTS: 00020710 C 00020720 C ENDPLT 00020730 C LCTRTN ENTRY POINTS AVAILABLE IN THE U.N.B. PLOTTING PACKAGE 00020740 C PLOTID 00020750 C SFTRTN 00020760 00020770 00020780 00020790 CC 00020800 C***********************************************************************00020810 INTEGER SNUM 00020820 DOUBLE PRECISION LONGR,LATR 00020830 COMMON /IDPLOT/ ID 00020840 COMMON /EXTREM/ FLORFD,CEILFD,FLORLD,CEILLD 00020850 COMMON/MEASUR/ WIDTH,HEIGHT 00020860 COMMON /NUMBRS/ LSNUM,SNUM 00020870 COMMON /IBOX/ IBOX 00020880 DATA ALAT,ALONG /'N','W'/ 00020890 C IF A POINT HAS AN IC VALUE OF < 0, THEN IT IS THE BEGINNING OF A 00020900 C NEW LINE AND SO THE LINE SEQUENCE NUMBER IS INCREASED. 00020910 IF(IC.LE.0) LSNUM = LSNUM+1 00020920 IF(IBOX.GT.0) GO TO 200 00020930 C TO CONVERT X & Y TO DECIMAL LONG. AND LAT. 00020940 ALAMDA = -FLORLD-(WIDTH-IX)*(CEILLD-FLORLD)/WIDTH 00020950 APHI = FLORFD + IY*(CEILFD-FLORFD)/HEIGHT 00020960 C TO CONVERT X TO INTEGER DEGREES,INTEGER MINUTES,& INTEGER SECONDS 00020970 C OF LONGITUDE. 00020980 LONGD = INT(ALAMDA) 00020990 REM = ALAMDA-FLOAT(LONGD) 00021000 SIXREM = REM * 60. 00021010 LONGM = INT (SIXREM) 00021020 REM = SIXREM -FLOAT(LONGM) 00021030 LONGS = INT(REM*60.) 00021040 LONGR = ALAMDA * 3.14159265/180. 00021050 C TOCONVERT Y TO INTEGER DEGREES, INTEGER MINUTES, AND INTEGER 00021060 C SECONDS OF LATITUDE. 00021070 LATD = INT(APHI) 00021080 REM = APHI-FLOAT(LATD) 00021090 SIXREM = REM * 60. 00021100 LATM = INT(SIXREM) 00021110 REM = SIXREM-FLOAT(LATM) 00021120 LATS = INT(REM*60.) 00021130 LATR = APHI*3.14159265/180. 00021140 C TO KEEP TRACK OF THE SEQUENCE NUMBER OF EACH SEGMENT 00021150 SNUM = SNUM + 1 00021160 LONGD = IABS(LONGD) 00021170 LONGM = IABS(LONGM) 00021180 LONGS = IABS(LONGS) 00021190 IF(ID.EQ.1) GO TO 20 00021200 IF (ID.EQ.3) GO TO 30 00021210 C THE FIRST WRITES THE CONTOUR DATA FILE FOR VCM 00021220 C THE SECOND WRITES THE LINE SEGMENT DATA FILE 00021230 C THE LAST WRITES THE CONTOUR DATA FILE FOR THE STANDARD DEVIATION 00021240 WRITE ( 4, 100) LSNUM,LATR,LONGR,LATD,LATM,LATS,ALAT,LONGD, 00021250 & LONGM,LONGS,ALONG,SNUM 00021260 RETURN 00021270 20 WRITE(11,100) LSNUM,LATR,LONGR,LATD,LATM,LATS,ALAT,LONGD,LONGM, 00021280 C LONGS,ALONG,SNUM 00021290 RETURN 00021300 30 WRITE (13,100) LSNUM,LATR,LONGR,LATD,LATM,LATS,ALAT,LONGD, 00021310 & LONGM,LONGS,ALONG,SNUM 00021320 RETURN 00021330 100 FORMAT(I6,E20.8,E20.8,4X,I2,I2,I2,A1,I3,I2,I2,A1,6X,I9) 00021340 C 00021350 C THE PURPOSE OF THE VARIABLE IBOX IS TO PREVENT THE OUTLINING BOX 00021360 C OF THE CONTOUR MAP (5 PAIR OF COORDINATES) FROM BEING PLOTTED. 00021370 C 00021380 200 IBOX = IBOX + 1 00021390 IF(IBOX.GT.5) IBOX = 0 00021400 RETURN 00021410 ENTRY LCTRTN(I,J) 00021420 RETURN 00021430 ENTRY PLOTID(IARRAY) 00021440 RETURN 00021450 ENTRY ENDPLT 00021460 RETURN 00021470 ENTRY SFTRTN(K,L) 00021480 RETURN 00021490 END 00021500 SUBROUTINE SEGMNT(X,Y,XTRA,YTRA,ILIM,NDFI,NDLA,DY,DX) C C THIS SUBROUTINE IS DESIGNED TO PLOT THE INPUT SEGMENTS.IT USES C SUBROUTINES PLTRTN,SYMBOL,PLOT. C C INPUT PARAMETERS: C C X(1,I)= VECTOR OF LONGITUDES OF BEGINNING BENCHMARKS IN MINS OF ARC C X(2,I)= VECTOR OF LONGITUDES OF END BENCHMARKS IN MINS OF ARC C Y(1,I)= VECTOR OF LATITUDES OF BEGINNING BENCHMARKS IN MINS OF ARC C Y(2,I)= VECTOR OF LATITUDES OF END BENCHMARKS IN MINS OF ARC C ILIM= NUMBER OF SEGMENTS TO BE PLOTTED C XTRA,YTRA= DUMMY VARIABLES PASSED ON FOR REASONS OF DIMENSIONING C COMMON/NUMBRS/ LSNUM,NUM COMMON /IDPLOT/ ID 00019500 COMMON /MEASUR/ WIDTH,HEIGHT 00019510 COMMON /EXTREM/ YMIN,YMAX,XMIN,XMAX 00019520 COMMON/XYINCH/ XINCH,YINCH 00019550 COMMON/BLOCK1/ MM,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND REAL XINCH,YINCH 00019540 DOUBLE PRECISION X,Y,YMIN,YMAX,XMIN,XMAX DIMENSION X(2,MXLL),Y(2,MXLL) DIMENSION XTRA(2,MXLL),YTRA(2,MXLL) WRITE(6,101) XINCH,YINCH 101 FORMAT(//,5X,'XINCH= ',F10.4,2X,'YINCH= ',F10.4) WRITE(6,102) XMAX, YMAX,XMIN,YMIN 102 FORMAT(//,5X,'XMAX= ',F11.4,2X,'YMAX= ',F11.4,2X,'XMIN= ',F11.4, 1 2X,'YMIN= ',F11.4) WIDTH= XINCH*1000. 00019630 HEIGHT = YINCH*1000. 00019640 LSNUM = 100000 00019650 NUM=0 ID=4 DO 6 I=1,ILIM 00019880 XTRA(1,I)=((X(1,I)-XMAX)/(XMIN-XMAX))*XINCH XTRA(2,I)=((X(2,I)-XMAX)/(XMIN-XMAX))*XINCH YTRA(1,I)=((Y(1,I)-YMIN)/(YMAX-YMIN))*YINCH YTRA(2,I)=((Y(2,I)-YMIN)/(YMAX-YMIN))*YINCH 6 CONTINUE 00019930 DO 8 J = 1,ILIM 00019950 CALL SYMBOL (XTRA(1,J),YTRA(1,J),.1,1,0.,-1) CALL PLOT (XTRA(2,J),YTRA(2,J),2) CALL SYMBOL(XTRA(2,J),YTRA(2,J),.1,1,0.,-1) 8 CONTINUE 00019990 CALL PLOT(WIDTH/1000.+2.,0.0,-3) RETURN 00020000 END 00020030 *PROCESS OPT=1 SUBROUTINE PREMAT(IPRD,NX,NY,NPOT,NEPI,TE,X,Y,T,XYPR,IFILE,D) C C THIS SUBROUTINE COMPUTES PREDICTION MATRIX XYPR FOR PREDICTIONS PU: C PU = XYPR*C FOR GIVEN SET OF POINTS (X,Y) AND A GIVEN DATE T. XYPR C IS STORED AS IPRD VECTORS OF LENGTH (MM = (NX*NY+NX+NY)*(NPOT+NEPI). C C INPUT PARAMETERS : C IPRD= NUMBER OF POINTS FOR WHICH PREDICTIONS ARE WANTED C NX= HIGHEST DEGREE OF SURFACE IN X C NY= HIGHEST DEGREE OF SURFACE IN Y C NPOT= HIGHEST DEGREE OF SURFACE IN TIME C NEPI= NUMBER OF EPISODES ENFORCED C TE= DATES OF BEGINNINGS AND ENDS OF ENFORCED EPISODES IN CENTURIES C X,Y= COORDINATES OF POINTS FOR WHICH PREDICTION IS WANTED C T= DATE FOR WHICH PREDICTION IS WANTED IN CENTURIES C XYPR= PREDICTION MATRIX RETURNED BY THIS SUBROUTINE C IFILE= NUMBER OF FILE IN WHICH XYPR IS PASSED C D= DUMMY VARIABLE PASSED ON FOR REASONS OF DIMENSIONING C IMPLICIT REAL*8(A-H,O-Z) COMMON/BLOCK1/MM,MXNT,MXNEPI,MXM,MXNPRD,MXLL,MXNTID,MXND C DIMENSION TE(2,MXNEPI),XYPR(MXM),X(MXNPRD),Y(MXNPRD),D(MXNEPI) NXP1=NX+1 NYP1=NY+1 IF(NEPI.EQ.0) GO TO 6 DO 1 J=1,NEPI 1 D(J)=1D0/(TE(2,J)-TE(1,J)) 6 CONTINUE DO 2 I=1,IPRD M=0 DO 8 L=1,NXP1 DO 8 K=1,NYP1 IF((L+K).EQ.2) XY= 1.0D0 IF((L+K).EQ.2) GO TO 7 IF(L.EQ.1) XY=Y(I)**(K-1) IF(L.EQ.1) GO TO 7 IF(K.EQ.1) XY=X(I)**(L-1) IF(K.EQ.1) GO TO 7 XY =X(I)**(L-1)*Y(I)**(K-1) 7 IF(NEPI.EQ.0) GO TO 5 DO 3 J=1,NEPI FUT=D(J)*(T-TE(1,J)) IF(T.LE.TE(1,J)) FUT=0D0 IF(T.GE.TE(2,J)) FUT=1D0 M=M+1 3 XYPR(M)=XY*FUT 5 DO 4 J=1,NPOT M=M+1 XYPR(M)=XY*T**J 4 CONTINUE 8 CONTINUE WRITE(IFILE) XYPR 2 CONTINUE RETURN END C CONMAP -- PLOTTER CONTOUR MAP OF GRIDDED DATA 01 0010 C 01 0020 C SUBROUTINE CONMAP (DATA,NDIM,NY,MX,NYS,MXS,CONINT, NDARK,NLABEL, 01 0030 C 1 DY,DX, WRK,LWRK) C 01 0050 C THIS PROGRAM PRODUCES A CALCOMP CONTOUR MAP FROM A TWO-DIMENSIONAL01 0060 C MATRIX OF FUNCTION VALUES. 01 0070 C 01 0080 C DATA IS (REAL) ARRAY OF VALUES TO BE CONTOURED. COLUMNS 01 0090 C LENGTH NY CONTAIN DATA FOR INCREASING Y-VALUES, ROWS OF LENGTH 01 0100 C MX CONTAIN DATA FOR INCREASING X-VALUES. NDIM IS THE DIMENSIONED 01 0110 C LENGTH OF A DATA COLUMN (NDIM.GE.NY). THUS THE FIRST ROW OF THE 01 0120 C MATRIX FORMS THE BOTTOM EDGE OF THE MAP, STARTING FROM THE LEFT 01 0130 C (AT THE CURRENT PLOTTER ORIGIN). 01 0140 C 01 0150 C THE DATA GRID MAY BE REFINED USING TWO-DIMENSIONAL FIFTH-ORDER 01 0160 C SPLINE INTERPOLATION. NYS & MXS ARE THE NUMBER OF SUB-INTERVALS 01 0170 C DESIRED FOR EACH GRID RECTANGE IN Y & X DIRECTIONS, RESPECTIVELY, 01 0180 C NO INTERPOLATION IS DONE IF NYS=1=MXS. THE ACTUAL CONTOURING IS 01 0190 C DONE USING LINEAR INTERPOLATION IN THE (REFINED) GRID. 01 0200 C 01 0210 C CONINT IS CONTOUR INTERVAL, IN SAME UNITS AS THE DATA VALUES. 01 0230 C IF CONINT IS ZERO, AN APPROPRIATE INTERVAL IS SELECTED BY CONMAP. 01 0240 C IF CONINT IS NEGATIVE, ITS MAGNITUDE IS IGNORED AND A COMMON BLOCK01 0250 C MUST BE PROVIDED CONTAINING A LIST OF SPECIFIC CONTOUR LEVELS (IN 01 0260 C ANY ORDER): 01 0270 C COMMON /CONLEV/ NLEV,CLEVS(NLEV) 01 0280 C IN ADDITION, THE RANGE OF CONTOURING MAY BE LIMITED BY CODING: 01 0290 C CALL SETLIM (CMIN,CMAX) 01 0300 C BEFORE THE CALL TO CONMAP, WHERE CMIN & CMAX ARE THE LOWEST AND 01 0310 C HIGHEST VALUES OF THE DATA RANGE TO BE CONTOURED. 01 0320 C 01 0330 C EVERY NDARKTH CONTOUR IS MADE DARKER UNLESS NDARK=0. 01 0340 C EVERY NLABELTH CONTOUR IS LABELLED WITH ITS VALUE UNLESS NLABEL=0.01 0350 C 01 0360 C DY & DX ARE RELATIVE SPACINGS OF SUPPLIED GRID; MAP WILL BE SCALED01 0370 C TO ABOUT TEN INCHES HEIGHT. IF THE SPACINGS ARE NEGATIVE, THEIR 01 0380 C MAGNITUDES ARE USED AS ACTUAL GRID DIMENSIONS IN INCHES, WITH 01 0390 C NO SCALING. 01 0400 C 01 0410 C WRK IS WORKSPACE OF DIMENSION LWRK. 01 0420 C NOTE THAT MAP WILL BE SEGMENTED IF NECESSARY FOR CONTOURING 01 0430 C WITHIN GIVEN WORKSPACE. IF YOU WISH TO ALLOW NI*MI SEGMENTS, 01 0440 C WORKSPACE SIZE IS DETERMINED GENEROUSLY BY: 01 0450 C NYW = NYS*(NY+4)/NI IF NI*MI=1, USE: 01 0460 C MXW = MXS*(MX+4)/MI NYW=NYS*(NY-1)+1 01 0470 C NPTCON=3*(NYW+MXW) MXW=MXS*(MX-1)+1 01 0480 C LWRK = NYW*MXW + 2*NPTCON 01 0490 C THE VALUES PRINTED BY THE PROGRAM MAY BE USED IN SUBSEQUENT RUNS. 01 0500 C ANY OR ALL OF THESE VALUES MAY BE SPECIFIED BEFORE CALLING CONMAP,01 0510 C USING: 01 0520 C CALL SPECWK (NYW,MXW,NPTCON) 01 0530 C PUT ZEROS FOR VALUES TO BE DETERMINED AUTOMATICALLY: VALUES OF 01 0540 C NYW & MXW, IF DESIRED, SHOULD BE SPECIFIED AS A PAIR (GIVING JUST 01 0550 C ONE HAS NO EFFECT). 01 0560 C 01 0570 C IF THE DATA GRID IS NOT TO BE REGINED BEFORE CONTOURING, THE DATA 01 0580 C MATRIX MAY BE CONTOURED IN PLACE IN A SINGLE SEGMENT (WITH A 01 0590 C LARGE SAVING IN MEMORY) BY EQUIVALENCING DATA & WRK, AND USING: 01 0600 C CALL SPECWK (NDIM,MX, 3*(NDIM+MX)) 01 0610 C WITH LWRK SET ACCORDING TO THE ABOVE FORMULA. 01 0620 C 01 0630 C 01 0640 C ADDITIONAL DISPLAY OPTIONS 01 0650 C 01 0660 C BLANKING: 01 0670 C 01 0680 C DESIGNATED AREAS OF THE MAP MAY BE LEFT BLANK BY SUTABLE DATA 01 0690 C PREPARATION, AND INSERTION OF: 01 0700 C CALL CONBLK (1,0,0,0) 01 0710 C TO PREPARE THE DATA, FOR EACH GRID POINT CODE 01 0720 C CALL SAVBLK (DATA(I,J),IND) 01 0730 C WHERE IND=0 FOR POINTS IN AREAS TO BE BLANKED, 01 0740 C =1 FOR POINTS IN CONTOURED AREAS. 01 0750 C THE SAVBLK ROUTINE ADJUSTS LOW-ORDER DATA BITS FOR USE AS FLAGS. 01 0760 C THE PERIMETER OF ANY AREA TO BE BLANKED SHOULD CONTAIN AT LEAST 01 0770 C ONE ROW OF REASONABLE DATA VALUES, TO ENSURE UNDISTORTED 01 0780 C INTERPOLATION IN NEARBY AREAS OF THE MAP. 01 0790 C 01 0800 C SUPPRESSION OF CLOSELY SPACED CONTOURS: 01 0810 C 01 0820 C IF A MIXTURE OF LIGHT AND DARK CONTOURS IS BEGIN PLOTTED (NDARKS1)01 0830 C THE LIGHT CONTOURS MAY BE SUPPRESSED IN AREAS OF HIGH CONTOUR 01 0840 C DENSITY. BEFORE CALLING CONMAP, 01 0850 C CALL CONBLK (IBLANK,NCINCH,0,0) 01 0860 C WHERE NCINCH IS THE NUMBER OF CONTOURS PER INCH DENSITY ABOVE 01 0870 C WHICH LIGHT CONTOURS ARE TO BE SUPPRESSED (E.G.,NCINCH=50). 01 0880 C IBLANK IS ZERO UNLESS THE BLANKING FEATURE IS BEING USED. 01 0890 C 01 0900 C PARTIAL MAP: 01 0910 C 01 0920 C THE GIVEN ARRAY WILL BE CONSIDERED AS A BLOCK OF A LARGER ARRAY 01 0930 C (WHICH MAY NOT FIT IN MEMORY) IF THE LAST TWO ARGUMENTS OF 01 0940 C THE CALL TO CONBLK ARE NON-ZERO. BEFORE CALLING CONMAP, 01 0950 C CALL CONBLK (IBLANK,NCINCH, LOWUPR,LEFTRT) 01 0960 C IBLANK IS ZERO UNLESS THE BLANKING FEATURE IS BEING USED. 01 0970 C NCINCH IS ZERO UNLESS THE SUPPRESSION FEATURE IS BEING USED. 01 0980 C LOWUPR & LEFTRT DESCRIBE THE POSITION OF THE PRESENT DATA BLOCK 01 0990 C WITHIN THE LARGER DATA ARRAY; VALUES ARE: 01 1000 C LOWUPR=+2 ON UPPER EDGE BUT NOT ON LOWER EDGE 01 1010 C +1 NOT ON EITHER UPPER OR LOWER EDGE 01 1020 C -1 ON LOWER EDGE BUT NOT ON UPPER EDGE 01 1030 C -2 ON BOTH EDGES (BLOCK SPANS MAP FROM TOP TO BOTTOM) 01 1040 C LEFTRT=+2 ON RIGHT SIDE BUT NOT LEFT 01 1050 C +1 NEITHER SIDE 01 1060 C -1 ON LEFT SIDE BUT NOT RIGHT 01 1070 C -2 ON BOTH SIDES 01 1080 C NOTE THAT AN ADDITIONAL ROW OR COLUMN OF GRID VALUES MUST BE 01 1090 C INCLUDED IN THE DATA ARRAY ON ANY SIDE WHICH IS NOT AN EDGE OF THE01 1100 C LARGER MAP, IN ORDER TO PROVIDE DATA OVERLAP FOR CONTINUITY OF 01 1110 C CONTOURS WHEN THE MAP BLOCKS ARE JOINED. 01 1120 C 01 1130 C 01 1140 C SAMPLE MINIMAL CALL: 01 1150 C 01 1160 C REAL DATA(8,8), WRK(650) 01 1170 C CALL CONMAP (DATA,8,6,7, 3,4, 0.,0,1, O.,4., WRK,650) 01 1180 C 01 1190 C THIS CALL CONTOURS A SAMLL 6 X 7 MAP IN A DATA MATRIX WITH COLUMNS01 1200 C OF LENGTH 8. EACH GRID IS DIVIDED INTO 3X4=12 SUBGRIDS. 01 1210 C A CONTOUR INTERVAL WILL BE SELECTED AUTOMATICALLY, AND CONTOURS 01 1220 C WILL BE LABELLED BUT NOT BOLDFACE. THE GRIDS HAVE SIDES IN THE 01 1230 C RATIO OF 3. TO 4. (NOTE THIS GIVES SQUARE SUB-GRIDS). THE 01 1240 C WORKSPACE IS LARGE ENOUGH TO AVOID SEGMENTATION. 01 1250 C 01 1260 C SUBROUTINES WITHIN THE PACKAGE ARE: ISOMAP,ISOVAL,ROW,ALOGAR 01 1270 C CONBLK,SPECWK,SAVBLK,SETLIM 01 1280 C COMMON BLOCKS INCLUDED ARE: CONLEV,CONSCL 01 1290 C PACKAGE SIZE IS ABOUT 18K BYTES. 01 1300 C 01 1310 C OTHER ROUTINES CALLED ARE: 01 1320 C CINSPL IN SAME LIBRARY AS CONMAP 01 1330 C DIPOL ENTRY TO CINSPL 01 1340 C SQRT,ATAN,IBCOM# STANDARD FORTRAN ROUTINES 01 1360 C PLOT,NUMBER CALCOMP ROUTINES 01 1350 C SQRT,ATAN,IBCOM# STANDARD FORTRAN ROUTINES 01 1360 C 01 1370 C AT U. OF T., CONMAP PACKAGE IS IN LIBRARYD DSN=USER.QQGP.GEOLIB 01 1380 C 01 1390 C CODED BY YVES LAMONTAGNE, 1972 01 1400 C WRITE-UP BY T. CLEE, APR. 1973 01 1410 C 01 1420 C 01 1430 SUBROUTINE CONMAP (DATA,NDIM,NY,MX,NYS,MXS,CONINT,NDARK,NLABEL, 01 1440 1 DY,DX,WRK,LWRK) 01 1450 LOGICAL NCRITE 01 1460 REAL DATA(NDIM,NDIM), WRK(LWRK) 01 1470 EXTERNAL CINSPL 01 1480 COMMON /CONSCL/ CFILL(6),IDENC 01 1490 DATA IM1,IM2,IB,KMAX,NW,MW,IDEN/-2,-2,5*0/ 01 1500 IDENC=IABS(IDEN) 01 1510 NS=MAX0(NYS,1) 01 1520 MS=MAX0(MXS,1) 01 1530 C 01 1540 C WORKSPACE OPTIMIZATION 01 1550 IF (NW*MW.NE.0) GO TO 3 01 1560 NA=NS*(NY-1)+1 01 1570 MA=MS*(MX-1)+1 01 1580 KKW=5*(NA+MA) 01 1590 IF (KMAX.NE.0) KKW=KMAX 01 1600 IF (LWRK.LT.NA*MA+KKW) GO TO 4 01 1610 NW=NA 01 1620 MW=MA 01 1630 GO TO 3 01 1640 4 LW=SQRT(FLOAT(LWRK)) 01 1650 KTOP=KMAX 01 1660 IF (KMAX.EQ.0) KTOP=5*LW 01 1670 KKW=LWRK-2*KTOP 01 1680 IF (NW*MW.EQ.0) GO TO 2 01 1690 IF (NW*MW.LE.KKW) GO TO 3 01 1700 2 LW=SQRT(FLOAT(KKW)) 01 1710 CALL DIPOL(IPOL) 01 1720 NC=(LW-1)/NS-IPOL 01 1730 NI=(NY+NC-2)/NC 01 1740 NC=(NY+NI-2)/NI 01 1750 MMC=(KKW/(NS*(NC+IPOL)+1)-1)/MS-IPOL 01 1760 MMI=(MX+MMC-2)/MMC 01 1770 MC=(LW-1)/MS-IPOL 01 1780 MI=(MX+MC-2)/MC 01 1790 MC=(MX+MI-2)/MI 01 1800 NNC=(KKW/(MS*(MC+IPOL)+1)-1)/NS-IPOL 01 1810 NNI=(NY+NNC-2)/NNC 01 1820 NCRITE=NI*MMI.LE.MI*NNI 01 1830 IF (NCRITE) MC=MMC 01 1840 IF (.NOT.NCRITE) NC=NNC 01 1850 NW=NS*(NC+IPOL)+1 01 1860 MW=MS*(MC+IPOL)+1 01 1870 3 IF (KMAX .EQ. 0) KMAX=3*(NW+MW) 01 1880 KMAX=MIN0(KMAX,(LWRK-NW*MW)/2) 01 1890 IX=NW*MW+1 01 1900 C 01 1910 CALL ISOMAP (DATA,NDIM,NY,MX, WRK,NW,MW, CINSPL,KMAX,NS,MS, 01 1920 1 WRK(IX), WRK(IX+KMAX), WRK(IX),WRK,CONINT, NDARK,NLABEL, 01 1930 2 DX,DY, IM1,IM2, IB,DATA) 01 1940 KMAX=0 01 1950 IB=0 01 1960 NW=0 01 1970 MW=0 01 1980 1 RETURN 01 1990 C 01 2000 ENTRY SPECWK (NYW,MXW,NPTCON) 01 2010 NW=NYW 01 2020 MW=MXW 01 2030 KMAX=NPTCON 01 2040 GO TO 1 01 2050 C 01 2060 ENTRY CONBLK (IBLANK,NCINCH,LOWUPR,LEFTRT) 01 2070 IDEN=NCINCH 01 2080 IM1=LOWUPR 01 2090 IM2=LEFTRT 01 2100 IF (IM1.EQ.0) IM1=-2 01 2110 IF (IM2.EQ.0) IM2=-2 01 2120 IB=IBLANK 01 2130 GO TO 1 01 2140 C 01 2150 ENTRY SAVBLK(IVAL,IND) 01 2160 DATA MAGIC/Z08100008/ 01 2170 C NOTE MAGIC HAS LAST TWO BITS ZERO. AND IS NEARLY ZERO (WITHOUT 01 2180 C UNDERFLOW) WHEN TREATED AS FLOATING POINT. 01 2190 IF (IVAL.EQ.0) IVAL=MAGIC 01 2200 IVAL=((IVAL-3)/4)*4-IND 01 2210 GO TO 1 01 2220 END 01 2230 FUNCTION ALOGAR(XX) 01 2240 DIMENSION A(10) 01 2250 DATA A/1.25,1.5,2.,2.5,3.,4.,5.,6.,8.,10./ 01 2260 F=1. 01 2270 X=XX 01 2280 IF(X.EQ.0.)GO TO 7 01 2290 GO TO 2 01 2300 1 X=X*10. 01 2310 F=F/10. 01 2320 2 IF(X.LE.1.)GO TO 1 01 2330 GO TO 4 01 2340 3 X=X/10. 01 2350 F=F*10. 01 2360 4 IF(X.GT.10.)GO TO 3 01 2370 DO 5 I=1,10 01 2380 IF(X.LE.A(I))GO TO 6 01 2390 5 CONTINUE 01 2400 6 ALOGAR=A(I)*F 01 2410 RETURN 01 2420 7 ALOGAR=0. 01 2430 RETURN 01 2440 END 01 2450 SUBROUTINE ISOMAP(V,ND1,NG,MG,W,ND,MD,INTERP,KMAX,NS,MS,X,Y,U,IFF,01 2460 1ENT,NR,NV,EX,EY,IM1,IM2,IB,IDATA) 01 2470 DIMENSION V(ND1,ND1),W(ND,MD),IFF(ND,MD),X(KMAX),Y(KMAX),U(MG), 01 2480 1 IDATA(ND1,ND1) 01 2490 COMMON /CONSCL/ DX,DY,TX,TY, ISEG,NSEG 01 2500 COMMON /IBOX/ IBOX 01 2501 DATA MAGIC/Z08100008/ 01 2510 INTEGER SNUM 01 2511 COMMON /NUMBRS/ LSNUM,SNUM 01 2512 COMMON/MEASUR/ WIDTH,HEIGHT 01 2513 LSNUM = 100000 01 2514 SNUM = 0 01 2515 IBOX = 1 01 2516 CALL DIPOL(IPOL) 01 2520 NP1=NG 01 2530 NP2=MG 01 2540 IF(IM1.EQ.-2) GO TO 41 01 2550 IF(IM1.GT.0)NP1=NP1-IPOL/2 01 2560 IF(IABS(IM1).NE.2)NP1=NP1-IPOL/2 01 2570 41 IF(IM2.EQ.-2) GO TO 42 01 2580 IF(IM2.GT.0)NP2=NP2-IPOL/2 01 2590 IF(IABS(IM2).NE.2)NP2=NP2-IPOL/2 01 2600 42 NG1=NP1-1 01 2610 MG1=NP2-1 01 2620 NC=(ND-1)/NS-IPOL 01 2630 NI=(NG1+NC-1)/NC 01 2640 IF (NS*(NG-1).LT.ND) NI=1 01 2650 NC=(NG1+NI-1)/NI 01 2660 NW=NS*((NG-2)/NI+IPOL+1)+1 01 2670 IF (NI.EQ.1) NW=NS*(NG-1)+1 01 2680 NND=NS*NC+1 01 2690 MC=(MD-1)/MS-IPOL 01 2700 MI=(MG1+MC-1)/MC 01 2710 IF (MS*(MG-1).LT.MD) MI=1 01 2720 MC=(MG1+MI-1)/MI 01 2730 MW=MS*((MG-2)/MI+IPOL+1)+1 01 2740 IF (MI.EQ.1) MW=MS*(MG-1)+1 01 2750 NSEG=NI*MI 01 2760 MMD=NW*MW+2*KMAX 01 2770 WRITE (6,10) NSEG,NW,MW,KMAX,MMD,ND,MD 01 2780 10 FORMAT ('0ISOMAP: THIS MAP WILL BE CONTOURED IN', I3, 01 2790 1 ' SEGMENTS, USING WORKSPACE SIZE ',I4,' *',I4,' + 2 *',I4,' =', 01 2800 1 I6/ 10X, 'ALLOCATED RECTANGULAR WORKSPACE DIMENSIONS ARE ', 01 2810 1 I4, ' BY ', I4 ) 01 2820 SY=-1. 01 2830 IF (EY.GT.0.) SY=ALOGAR(FLOAT(NG1)*EY/10.) 01 2840 12 DX=SY*FLOAT(MS)/EX 01 2850 DY=SY*FLOAT(NS)/EY 01 2860 IF(IM1.NE.-2.OR.IM2.NE.-2)GO TO 43 01 2870 X(1)=0. 01 2880 X(2)=FLOAT(MG1)*EX/SY 01 2890 X(3)=FLOAT(MG1)*EX/SY 01 2900 X(4)=0. 01 2910 X(5)=0. 01 2920 Y(1)=0. 01 2930 Y(2)=0. 01 2940 Y(3)=FLOAT(NG1)*EY/SY 01 2950 Y(4)=FLOAT(NG1)*EY/SY 01 2960 Y(5)=0. 01 2970 WRITE (6,13) X(3),Y(3) 01 2980 WIDTH = X(3) * 1000. 01 2981 HEIGHT = Y(3) * 1000. 01 2982 13 FORMAT ('0',9X,'MAP SIZE: WIDTH=',F6.2,' IN., HEIGHT=',F6.2, 01 2990 1 ' IN.') 01 3000 IPEN=3 01 3010 DO 11 IN=1,5 01 3020 CALL PLOT (X(IN),Y(IN),IPEN) 01 3030 11 IPEN=2 01 3040 43 MMD=MS*MC+1 01 3050 ISEG=0 01 3060 DO 1 IN=1,NI 01 3070 IP=0 01 3080 IF(IN.NE.1.AND.IM1.LT.0)IP=IPOL/2 01 3090 NI1=(IN-1)*NC-IP 01 3100 NIS=NI1+1 01 3110 NI2=IN*NC+IPOL+1-IP 01 3120 ITN=NS 01 3130 IFN=1 01 3140 IF(IN.NE.1.OR.IM1.GT.0)GO TO 4 01 3150 IFN=-1 01 3160 ITN=0 01 3170 4 NN=NND 01 3180 IF(IN.NE.NI)GO TO 5 01 3190 NI2=NG 01 3200 NN=NS*(NP1-NC*(IN-1)-1)+1 01 3210 IF(IABS(IM1).NE.2)GO TO 5 01 3220 IFN=2*IFN 01 3230 5 NJ=NI2-NI1 01 3240 NDJ=(NJ-1)*NS+1 01 3250 TY=(IN-1)*NC*NS 01 3260 NNI=1+ITN*IPOL/2 01 3270 DO 1 IM=1,MI 01 3280 ISEG=ISEG+1 01 3290 TX=(IM-1)*MC*MS 01 3300 JP=0 01 3310 IF(IM.NE.1.AND.IM2.LT.0)JP=IPOL/2 01 3320 MI1=(IM-1)*MC-JP 01 3330 ITM=MS 01 3340 IFM=1 01 3350 MIS=MI1+1 01 3360 MI2=IM*MC+IPOL+1-JP 01 3370 IF(IM.NE.1.OR.IM2.GT.0)GO TO 6 01 3380 IFM=-1 01 3390 ITM=0 01 3400 6 MM=MMD 01 3410 IF(IM.NE.MI)GO TO 7 01 3420 MI2=MG 01 3430 MM=MS*(NP2-MC*(IM-1)-1)+1 01 3440 IF(IABS(IM2).NE.2)GO TO 7 01 3450 IFM=2*IFM 01 3460 7 MJ=MI2-MI1 01 3470 MDJ=(MJ-1)*MS+1 01 3480 CALL INTERP(V(NIS,MIS),W,NJ,MJ,NDJ,MDJ,ND1,MD1,ND,MD,U,IPOL, 01 3490 1IFN,IFM) 01 3500 MMI=1+ITM*IPOL/2 01 3510 DO 48 J=1,MDJ 01 3520 DO 48 I=1,NDJ 01 3530 IF (IFF(I,J).EQ.0) IFF(I,J)=MAGIC 01 3540 48 IFF(I,J)=(IFF(I,J)/4)*4 - 2 01 3550 IF(IB.NE.1)GO TO 46 01 3560 J1=MIS+ITM*IPOL/MS/2 01 3570 J2=J1+(MM-1)/MS 01 3580 I1=NIS+IPOL*ITN/NS/2 01 3590 I2=I1+(NN-1)/NS 01 3600 NT=2*NS-1 01 3610 MT=2*MS-1 01 3620 K2=MMI-2*MS 01 3630 DO 45 J=J1,J2 01 3640 K2=K2+MS 01 3650 K1=NNI-2*NS 01 3660 DO 45 I=I1,I2 01 3670 K1=K1+NS 01 3680 IF((IDATA(I,J)/4)*4.NE.IDATA(I,J))GO TO 45 01 3690 DO 47 M=1,MT JJ=K2+M 01 3710 IF(JJ.LE.0)JJ=1 01 3720 DO 47 N=1,NT 01 3730 II=K1+N 01 3740 IF(II.LE.0)II=1 01 3750 47 IFF(II,JJ)=((IFF(II,JJ)-3)/4)*4 01 3760 45 CONTINUE 01 3770 46 CONTINUE 01 3780 1 CALL ISOVAL(W(NNI,MMI),NN,MM,KMAX,X,Y,ND,MD,IFF(NNI,MMI),ENT,NR, 01 3790 1NV) 01 3800 RETURN 01 3810 END 01 3820 SUBROUTINE ROW(W,IG1,INN,IMM,NG,MG,IFF,C,JN,JM,JJ,KMAX,ND,MD,X,Y, 01 5070 1 IR,IS,ENT) 01 5080 DIMENSION W(ND,MD),IFF(ND,MD),X(KMAX),Y(KMAX) 01 5090 COMMON /CONSCL/ XI,YI,TX,TY,ISEG,NSEG,NCINCH 01 5100 DATA XP,YP/0.,0./ 01 5110 CALL ERRSET(208,256,10,1,0,0) 01 5111 DEN1=(FLOAT(NCINCH)*ENT)**2 01 5120 H=0.07 01 5130 DEN=ENT*ENT/H/H/3.5 01 5140 RAT=XI/YI 01 5150 IBL=0 01 5160 PI=3.141593 01 5170 NG1=NG-1 01 5180 MG1=MG-1 01 5190 KN=INN*INN 01 5200 KM=IMM*IMM 01 5210 N=JN*KN 01 5220 M=JM*KM 01 5230 DO 1 J=1,IG1 01 5240 IM=IMM 01 5250 IN=INN 01 5260 N=N+KM 01 5270 M=M+KN 01 5280 IST=0 01 5290 ICOND=1 01 5300 N1=N+(1-IM)/2 01 5310 N2=N+(1+IM)/2 01 5320 NN=N2 01 5330 M1=M+(1-IN)/2 01 5340 M2=M+(1+IN)/2 01 5350 MM=M2 01 5360 K=0 01 5370 IF((IFF(N,M)/2)*2.NE.IFF(N,M))GO TO 1 01 5380 IF (ABS(W(N1,M1)-W(N2,M2)).LT.1.E-30) GO TO 1 01 5390 D=(C-W(N2,M2))/(W(N1,M1)-W(N2,M2)) 01 5400 IF(D.GT.1.0.OR.D.LE.0.0)GO TO 1 01 5410 IF(D.NE.1.)GO TO 19 01 5420 IF((C-W(N1,M1))/(W(N2,M2)-W(N1,M1)).LT.0.)GO TO 1 01 5430 19 ICO=(INN+2*IMM)*(MM+(NN-1)*MG) 01 5440 NNN=N 01 5450 MMM=M 01 5460 2 IF((IABS(IFF(NNN,MMM))/4)*4.GE.IABS(IFF(NNN,MMM))-1)GO TO 110 01 5470 K=K+1 01 5480 X(K)=FLOAT(M2-1)+D*FLOAT(M1-M2) 01 5490 Y(K)=FLOAT(N2-1)+D*FLOAT(N1-N2) 01 5500 IF(K.GE.KMAX)GO TO 8 01 5510 GO TO 111 01 5520 110 IF(K.EQ.0)GO TO 111 01 5530 IBL=1 01 5540 IF(IST.EQ.ICOND)IBL=0 01 5550 GO TO 8 01 5560 9 IBL=0 01 5570 K=0 01 5580 111 IF(IST.EQ.ICOND)GO TO 8 01 5590 III=0 01 5600 M1=MM 01 5610 N1=NN 01 5620 M2=M1+IM 01 5630 N2=N1+IN 01 5640 I=1 01 5650 GO TO 7 01 5660 4 M1=M2 01 5670 N1=N2 01 5680 M2=M1-IN 01 5690 N2=N1-IM 01 5700 I=2 01 5710 GO TO 7 01 5720 5 M1=M2 01 5730 N1=N2 01 5740 M2=M1-IM 01 5750 N2=N1-IN 01 5760 I=3 01 5770 7 IF (ABS(W(N1,M1)-W(N2,M2)).LT.1.E-30) GO TO (4,5,8),I 01 5780 D=(C-W(N2,M2))/(W(N1,M1)-W(N2,M2)) 01 5790 IF(D.GT.1.0.OR.D.LT.0.0)GO TO (4,5,8),I 01 5800 IF(D.NE.1.)GO TO 17 01 5810 IF((C-W(N1,M1))/(W(N2,M2)-W(N1,M1)).LT.0.)GO TO (4,5,8),I 01 5820 17 IMS=IM 01 5830 INS=IN 01 5840 IM=(N1-N2)*(IMS*IMS-INS*INS) 01 5850 IN=(M1-M2)*(INS*INS-IMS*IMS) 01 5860 MMM=(M1+M2)/2 01 5870 NNN=(N1+N2)/2 01 5880 IF((IFF(NNN,MMM)/2)*2.EQ.IFF(NNN,MMM))IFF(NNN,MMM)=IFF(NNN,MMM) 01 5890 1+ISIGN(1,IFF(NNN,MMM)) 01 5900 NN=(N1+N2+IM)/2 01 5910 MM=(M1+M2+IN)/2 01 5920 GO TO (11,12),JJ 01 5930 11 IST=(IM*MM+IN*NN) 01 5940 ICOND=IN*((1+IN)*NG+1-IN)/2+IM*((1+IM)*MG+1-IM)/2 01 5950 GO TO 2 01 5960 12 IST=(IN+2*IM)*(MM+(NN-1)*MG) 01 5970 ICOND=ICO 01 5980 GO TO 2 01 5990 8 IF(K.EQ.0)GO TO 1 01 6000 NCON=IR*3 01 6010 KKK=K 01 6020 IBACK=0 01 6030 IF(IS.EQ.0)GO TO 61 01 6040 IIS=IS 01 6050 IF(IS.GT.0)GO TO 31 01 6060 IIS=-IS 01 6070 31 CC=ABS(C) 01 6080 NDIG=1 01 6090 33 IF(CC.LT.1.)GO TO 32 01 6100 CC=CC/10. 01 6110 NDIG=NDIG+1 01 6120 GO TO 33 01 6130 32 NDEC=IIS-1 01 6140 IF(NDIG.EQ.1)NDIG=NDIG+1 01 6150 IF(NDEC.EQ.0)NDEC=-1 01 6160 IF(C.LT.0.)NDIG=NDIG+1 01 6170 NSP=NDIG+NDEC+1 01 6180 GX=H*FLOAT(NSP)*AMIN1(XI,YI)+1. 01 6190 NSPG=GX 01 6200 IF(K.LT.NSPG*3+5)GO TO 61 01 6210 NPEF=1.1*GX 01 6220 PN=NPEF 01 6230 MK=K 01 6240 COURS=1.E30 01 6250 T1=(Y(2)-Y(1))/(X(2)-X(1)+1.E-30)*RAT 01 6260 T2=(Y(3)-Y(2))/(X(3)-X(2)+1.E-30)*RAT 01 6270 T3=(Y(4)-Y(3))/(X(4)-X(3)+1.E-30)*RAT 01 6280 T4=(Y(5)-Y(4))/(X(5)-X(4)+1.E-30)*RAT 01 6290 C1=ABS((T2-T1)/(1.+T2*T1+1.E-30)) 01 6300 C2=ABS((T3-T2)/(1.+T3*T2+1.E-30)) 01 6310 C3=ABS((T4-T3)/(1.+T4*T3+1.E-30)) 01 6320 COUR=C1+C2+C3 01 6330 K1=K-2 01 6340 DO 73 IK=6,K1 01 6350 T5=(Y(IK)-Y(IK-1))/(X(IK)-X(IK-1)+1.E-30)*RAT 01 6360 C4=ABS((T5-T4)/(1.+T5*T4+1.E-30)) 01 6370 COUR=COUR+C4 01 6380 IF(COUR.GE.COURS)GO TO 72 01 6390 MY=Y(IK)+2. 01 6400 MX=X(IK)+2. 01 6410 IF(MX.GT.MG.OR.MY.GT.NG)GO TO 72 01 6420 IF(((XI*(W(MY,MX-1)-W(MY,MX)))**2+(YI*(W(MY-1,MX)-W(MY,MX)))**2) 01 6430 1.GT.DEN)GO TO 72 01 6440 MK=IK 01 6450 COURS=COUR 01 6460 72 COUR=COUR-C1 01 6470 C1=C2 01 6480 C2=C3 01 6490 C3=C4 01 6500 73 T4=T5 01 6510 IF(MK.EQ.K)GO TO 61 01 6520 MK=MIN0(MK,K-NPEF/2-1) 01 6530 K2=MK+NPEF/2 01 6540 K1=K2-NSPG+1 01 6550 IF(K1.LE.0)K1=1 01 6560 IK=NSPG 01 6570 52 IF(SQRT((X(K1)-X(K2))*(X(K1)-X(K2))/XI/XI+(Y(K1)-Y(K2))*(Y(K1)- 01 6580 1Y(K2))/YI/YI).GT. H*FLOAT(NSP))GO TO 50 01 6590 IF(IK.GT.K/3)GO TO 61 01 6600 IK=IK+1 01 6610 K1=K1-1 01 6620 IF(JJ.EQ. 1)GO TO 47 01 6630 GO TO 49 01 6640 47 IF(K1.LE.0)GO TO 61 01 6650 IF(K1.GT.K)GO TO 61 01 6660 49 IF(K1.LE.0)K1=K1+K-1 01 6670 IF(K1.GT.K)K1=K1+1-K 01 6680 GO TO 52 01 6690 50 KL=K1 01 6700 KR=K2 01 6710 IF(X(K1).GT.X(K2))GO TO 48 01 6720 GO TO 51 01 6730 48 KL=K2 01 6740 KR=K1 01 6750 51 RX=X(KR)-X(KL) 01 6760 IF(RX.EQ.0.)GO TO 61 01 6770 RY=Y(KR)-Y(KL) 01 6780 RR=SQRT(RX*RX+RY*RY) 01 6790 CS=RX/RR 01 6800 SN=RY/RR 01 6810 TH=ATAN(RY/RX*RAT) 01 6820 THETA=TH*180./PI 01 6830 X0=(X(KL)+TX)/XI+.7*H*CS+.5*H*SN 01 6840 Y0=(Y(KL)+TY)/YI+.7*H*SN-.5*H*CS 01 6850 G=C*(1.+1.E-6) 01 6860 CALL NUMBER(X0,Y0, H ,G,THETA,NDEC) 01 6870 IF(K1.GT.K2)GO TO 64 01 6880 KKK=K1 01 6890 IBACK=1 01 6900 KC=K-K2+1 01 6910 GO TO 61 01 6920 64 DO 65 IK=K2,K1 01 6930 IIK=IK-K2+1 01 6940 X(IIK)=X(IK) 01 6950 65 Y(IIK)=Y(IK) 01 6960 KKK=K1-K2+1 01 6970 GO TO 61 01 6980 63 DO 60 IK=1,KC 01 6990 IIK=IK+K2-1 01 7000 X(IK)=X(IIK) 01 7010 60 Y(IK)=Y(IIK) 01 7020 KKK=KC 01 7030 IBACK=0 01 7040 61 NI=0 01 7050 C PLOT THE CONTOUR 01 7060 62 IPEN=3 01 7070 D1=(XP-X(1))**2 + (YP-Y(1))**2 01 7080 D2=(XP-X(KKK))**2 + (YP-Y(KKK))**2 01 7090 IF (D2.LT.D1) GO TO 67 01 7100 ISN=+1 01 7110 IBAS=0 01 7120 GO TO 68 01 7130 67 ISN=-1 01 7140 IBAS=KKK+1 01 7150 68 DO 66 IPT=1,KKK 01 7160 L=IBAS+IPT*ISN 01 7170 IF (IR.GE.0) GO TO 69 01 7180 MX=X(L)+2. 01 7190 MY=Y(L)+2. 01 7200 IF (MX.GT.MG .OR. MY.GT.NG) GO TO 69 01 7210 IF (((W(MY,MX)-W(MY,MX-1))*XI)**2+((W(MY,MX)-W(MY-1,MX))*YI)**2 01 7220 $ .LT. DEN1) GO TO 69 01 7230 IPEN=3 01 7240 GO TO 66 01 7250 69 CALL PLOT ((X(L)+TX)/XI,(Y(L)+TY)/YI, IPEN) 01 7260 IPEN=2 01 7270 66 CONTINUE 01 7280 XP=X(L) 01 7290 YP=Y(L) 01 7300 NI=NI+1 01 7310 IF(NI.LT.NCON)GO TO 62 01 7320 IF(IBACK.EQ.1)GO TO 63 01 7330 IF(IBL.EQ.1)GO TO 9 01 7340 IBL=0 01 7350 1 CONTINUE 01 7360 RETURN 01 7370 END 01 7380 SUBROUTINE CINSPL(U,W,NG,MG,II,JJ,NDIM,MDIM,ND,MD,V,IPOL,IN,IM) 01 7390 DIMENSION U(NDIM,MDIM),W(ND,MD),V(MDIM) 01 7400 IPOL=2 01 7410 IF(MG.LE.3.OR.NG.LE.3)GO TO 14 01 7420 NS=(II-1)/(NG-1) 01 7430 MS=(JJ-1)/(MG-1) 01 7440 MS1=MS-1 01 7450 NG1=NG-1 01 7460 MG1=MG-1 01 7470 IONE=2 01 7480 NNG=NG1 01 7490 IF(IN.LT.0)IONE=1 01 7500 IF(IABS(IN).GE.2) NNG=NG 01 7510 IL=(IONE-1)*NS 01 7520 DO 1 N=IONE,NNG 01 7530 IA=-2 01 7540 IF(N.EQ.1)IA=-1 01 7550 IF(N.EQ.NG1)IA=-3 01 7560 IF(N.EQ.NG)IA=-4 01 7570 NIA=N+IA 01 7580 DO 6 JS=1,MG 01 7590 6 V(JS)=U(N,JS) 01 7600 NSS=NS 01 7610 IF(N.EQ.NNG)NSS=1 01 7620 DO 1 L=1,NSS 01 7630 IL=IL+1 01 7640 DI=FLOAT(L)/FLOAT(NS)-FLOAT(2+IA) 01 7650 JONE=2 01 7660 MMG=MG1 01 7670 IF(IM.LT.0)JONE=1 01 7680 IF(IABS(IM).GE.2) MMG=MG 01 7690 K=(JONE-1)*MS 01 7700 DO 21 M=JONE,MMG 01 7710 K=K+1 01 7720 JA=-2 01 7730 IF(M.EQ.1)JA=-1 01 7740 IF(M.EQ.MG1)JA=-3 01 7750 W(IL,K)=V(M) 01 7760 IF (M.EQ.MMG .OR. MS1.EQ.0) GO TO 21 01 7770 MJA=M+JA 01 7780 V1=V(MJA+1) 01 7790 V2=V(MJA+2) 01 7800 V3=V(MJA+3) 01 7810 V4=V(MJA+4) 01 7820 AQ=V2 01 7830 BQ=(V3-V1)/2. 01 7840 CQ=(11.*V1-27.*V2+21.*V3-5.*V4)/12. 01 7850 DQ=(V4-3.*V3+3.*V2-V1)/6. 01 7860 EQ=(5.*V4-15.*V3+15.*V2-5.*V1)/12. 01 7870 DO 11 KK=1,MS1 01 7880 K=K+1 01 7890 DJ=FLOAT(KK)/FLOAT(MS)-FLOAT(2+JA) 01 7900 11 W(IL,K)=AQ+DJ*(BQ+DJ*(CQ+DJ*(DQ+DJ*(EQ-DJ*DQ)))) 01 7910 21 CONTINUE 01 7920 IF(L.EQ.NSS)GO TO 1 01 7930 DO 12 JS=1,MG 01 7940 V1=U(NIA+1,JS) 01 7950 V2=U(NIA+2,JS) 01 7960 V3=U(NIA+3,JS) 01 7970 V4=U(NIA+4,JS) 01 7980 AQ=V2 01 7990 BQ=(V3-V1)/2. 01 8000 CQ=(11.*V1-27.*V2+21.*V3-5.*V4)/12. 01 8010 DQ=(V4-3.*V3+3.*V2-V1)/6. 01 8020 EQ=(5.*V4-15.*V3+15.*V2-5.*V1)/12. 01 8030 12 V(JS)=AQ+DI*(BQ+DI*(CQ+DI*(DQ+DI*(EQ-DI*DQ)))) 01 8040 1 CONTINUE 01 8050 ENTRY DIPOL(IPOL) 01 8060 IPOL=2 01 8070 RETURN 01 8080 14 WRITE (6,100) NG,MG 01 8090 100 FORMAT ('0WORKSPACE TOO SMALL FOR CINSPL. NG=',I2,' MG=',I2) 01 8100 STOP 1 01 8110 END 01 8120 *PROCESS OPT=1 SUBROUTINE ISOVAL(W,NG,MG,KMAX,X,Y,ND,MD,IFF,CONINT,NR,NV) 01 3830 COMMON /CONLEV/ NLEV,CC(2) 01 3840 COMMON /CONSCL/ FILL(4),ISEG,NSEG,NCINCH 01 3850 LOGICAL LIM/.FALSE./ 01 3860 DATA ENT,BMN,BMX/0., 1.E60,-1.E60/, MAGIC/Z08100008/ 01 3870 DIMENSION W(ND,MD),X(KMAX),Y(KMAX),IFF(ND,MD) 01 3880 ITEST=0 01 3890 ISUP=0 01 3900 ISUB=0 01 3910 AM=W(1,1) 01 3920 AN=W(1,1) 01 3930 DO 1 N=1,NG 01 3940 DO 1 M=1,MG 01 3950 IF ((IFF(N,M)/4)*4.EQ.IFF(N,M)) GO TO 1 01 3960 AM=AMAX1(AM,W(N,M)) 01 3970 AN=AMIN1(AN,W(N,M)) 01 3980 1 CONTINUE 01 3990 BMN=AMIN1(BMN,AN) 01 4000 BMX=AMAX1(BMX,AM) 01 4010 IF (LIM) AM=AMIN1(AM,CMAX) 01 4020 IF (LIM) AN=AMAX1 (AN,CMIN) 01 4030 IF (CONINT.NE.0.) GO TO 30 01 4040 IF (ENT.NE.0) GO TO 31 01 4050 ENT=ALOGAR((AM-AN)/20.) 01 4060 WRITE (6,32) ENT 01 4070 32 FORMAT ('0SELECTED CONTOUR INTERVAL IS', G12.3) 01 4080 GO TO 31 01 4090 30 IF (ENT.EQ.0.) WRITE (6,33) CONINT 01 4100 33 FORMAT ('0PROVIDED CONTOUR INTERVAL IS', G12.3) 01 4110 ENT=CONINT 01 4120 31 IF (ENT.LE.0.) GO TO 13 01 4130 EENT=ENT 01 4140 IC=IFIX(AN/ENT)-1 01 4150 IF(AN.GT.0.)IC=IC+1 01 4160 ICR=IFIX(AM/ENT)-IC 01 4170 IF (ICR.GT.1000) GO TO 21 01 4180 20 NG1=NG-1 01 4190 MG1=MG-1 01 4200 DO 2 L=1,ICR 01 4210 IC=IC+1 01 4220 IR=0 01 4230 IF(NR.EQ.0)GO TO 7 01 4240 IF(IC.EQ.NR*(IC/NR))IR=1 01 4250 7 IS=0 01 4260 IF(NV.EQ.0)GO TO 8 01 4270 IF(IC.EQ.NV*(IC/NV))IS=1 01 4280 8 IF (ENT.GT.0.) C=FLOAT(IC)*ENT 01 4290 IF(ENT.LE.0.)C=CC(IC) 01 4300 IF (C.GT.AM) ISUP=1 01 4310 IF (C.LT.AN) ISUB=1 01 4320 IF (C.GT.AM .OR. C.LT.AN) GO TO 2 01 4330 ITEST=1 01 4340 IF (NCINCH.GT.1 .AND. NR.GT.1 .AND. IR.EQ.0) IR=-1 01 4350 IA=-4 01 4360 IF(C.EQ.0.)IA=MAGIC 01 4370 IF(IS.EQ.0)GO TO 11 01 4380 CE=ABS(C)*1.00001 01 4390 NDEC=1 01 4400 IF (ENT.GT.0.) GO TO 10 01 4410 EENT=(AM-AN)/FLOAT(ICR) 01 4420 IF (L.GT.1 .OR. L.LT.ICR) EENT=AMIN1(EENT,ABS(CC(IC+1)-CC(IC-1)) 01 4430 1 /2.) 01 4440 10 ICRIT=CE 01 4450 CRIT=FLOAT(ICRIT)+1.E-2 01 4460 IF(CRIT.GE.CE)GO TO 9 01 4470 CE=CE*10. 01 4480 NDEC=NDEC+1 01 4490 GO TO 10 01 4500 9 IS=NDEC*NV/IABS(NV) 01 4510 11 DO 4 N=1,NG 01 4520 DO 4 M=1,MG 01 4530 IFF(N,M)=(IFF(N,M)/2)*2-4 01 4540 4 IF(W(N,M).EQ.C)IFF(N,M)=IFF(N,M)+IA 01 4550 DO 3 N=1,NG,NG1 01 4560 INN=1 01 4570 IF(N.EQ.NG)INN=-1 01 4580 3 CALL ROW(W,MG1,INN,0,NG,MG,IFF,C,N,0, 1,KMAX,ND,MD,X,Y,IR,IS,EENT)01 4590 DO 5 M=1,MG,MG1 01 4600 IMM=1 01 4610 IF(M.EQ.MG)IMM=-1 01 4620 5 CALL ROW(W,NG1,0,IMM,NG,MG,IFF,C,0,M, 1,KMAX,ND,MD,X,Y,IR,IS,EENT)01 4630 DO 6 N=2,NG1 01 4640 6 CALL ROW(W,MG1,1,0,NG,MG,IFF,C,N,0, 2,KMAX,ND,MD,X,Y,IR,IS,EENT) 01 4650 2 CONTINUE 01 4660 IF (ITEST.EQ.0) GO TO 23 01 4670 25 IF (ISEG.LT.NSEG) RETURN 01 4680 WRITE (6,36) BMN,BMX 01 4690 36 FORMAT ('0DATA RANGE FOR THIS MAP WAS',G12.3,' TO', G12.3) 01 4700 IF (LIM) WRITE (6,37) CMIN,CMAX 01 4710 37 FORMAT ('0SPECIFIED CONTOURING RANGE WAS',G12.3,' TO',G12.3) 01 4720 LIM=.FALSE. 01 4730 ENT=0. 01 4740 BMN=1.E60 01 4750 BMX=-1.E60 01 4760 RETURN 01 4770 13 ICR=MAX0(NLEV,1) 01 4780 ICR=MIN0(ICR,1000) 01 4800 IC=0 01 4810 IF (ISEG.GT.1) GO TO 20 01 4820 WRITE (6,34) NLEV 01 4830 34 FORMAT ('0NUMBER OF PROVIDED CONTOURS IS:',I5) 01 4840 WRITE (6,35) ICR,CC(1),CC(ICR) 01 4850 35 FORMAT ('01ST AND',I4,'TH CONTOUR LEVELS ARE:',2G12.3) 01 4860 GO TO 20 01 4870 C 01 4880 21 WRITE (6,22) ISEG 01 4890 22 FORMAT ('0ISOVAL: DATA RANGE IN SEGMENT',I3,' CONTAINS OVER A THOU01 4900 $SAND CONTOURS; NONE PLOTTED') 01 4910 GO TO 25 01 4920 23 IF (ISUP+ISUB.EQ.2) GO TO 19 01 4930 WRITE (6,24) ISEG 01 4940 24 FORMAT ('0ISOVAL: DATA RANGE IN SEGMENT',I3,' CONTAINS NONE OF THE01 4950 $ SPECIFIED CONTOURS') 01 4960 GO TO 25 01 4970 19 WRITE (6,18) ISEG 01 4980 18 FORMAT ('0ISOVAL: ENTIRE DATA RANGE IN SEGMENT',I3,' FALLS BETWEEN01 4990 $ TWO CONTOUR LEVELS') 01 5000 GO TO 25 01 5010 C 01 5020 ENTRY SETLIM (CMIN,CMAX) 01 5030 LIM=.TRUE. 01 5040 RETURN 01 5050 END 01 5060