C DATA SET COLOCATE AT LEVEL 012 AS OF 09/29/80 C ******************************************************************* 00001*11 C * * 00002*11 C * * 00003*11 C * * 00004*11 C * PROGRAM 'C O L O C A T E' * 00005*11 C * ========================= * 00006*11 C * * 00007*11 C * * 00008*11 C * * 00009*11 C * THE COMPUTER PROGRAN PACKAGE 'COLOCATE' IS DESIGNED TO * 00010*11 C * PREDICT FOR GRAVITY ANOMALY AT THE CENTRE OF A ONE DEGREE * 00011*11 C * SQUARE BLOCK USING THE METHOD OF LEAST-SQUARES COLLOCATION. * 00012*11 C * * 00013*11 C * FOR THE MATHEMATICAL MODELS, SEE: KASSIM (1980). "AN EVA- * 00014*11 C * LUATION OF THREE TECHNIQUES FOR THE PREDICTION OF GRAVITY ANO- * 00015*11 C * MALIES IN CANADA". M.SC. THESIS, DEPARTMENT OF SURVEYING * 00016*11 C * ENGINEERING, U.N.B., FREDERICTON. * 00017*11 C * * 00018*11 C * THE MAIN PROGRAM 'COLOCATE' READS THE LATITUDE AND * 00019*11 C * LONGITUDE OF THE CENTRE OF A ONE DEGREE SQUARE BLOCK, SEEKS * 00020*11 C * THE APPROPRIATE DATA FILE AND STORES ALL AVAILABLE DATA FOR THE* 00021*11 C * BLOCK , COMPUTING AT THE SAME TIME THE X,Y COORDINATES OF EACH * 00022*11 C * DATA POINT WITH THE PREDICTION POINT AS ORIGIN. THE CLOSEST * 00023*11 C * DATA TO THE PREDICTION POINT ARE LATER SELECTED FOR THE PREDI- * 00024*11 C * CTION. PROVISION IS MADE FOR THE DATA TO BE SCREENED FOR * 00025*11 C * CLUSTERS OF POINTS WHICH ARE REPLACED BY A WEIGHTED MEAN. * 00026*11 C * * 00027*11 C * THE PROGRAM IS DESIGNED TO PREDICT FOR ONE DEGREE SQUARE * 00028*11 C * BLOCK AT A TIME. HENCE, THE NUMBER OF BLOCKS THAT CAN BE PRO- * 00029*11 C * CESSED IS UNLIMITED SO LOG AS THERE IS ENOUGH TIME * 00030*11 C * * 00031*11 C * PROVISION IS MADE FOR THE RESULTS OF A SERIES OF * 00032*11 C * PREDICTIONS TO BE TESTED STATISTICALLY. THE TESTS INCLUDE: * 00033*11 C * CHI-SQUARE GOODNESS OF FIT TEST, MEASURES OF SKEWNESS AND * 00034*11 C * KURTOSIS, AND TESTS ON MEANS AND VARIANCES. THE HISTOGRAM OF * 00035*11 C * STANDARDIZED OBSERVATIONS IS ALSO PLOTTED. * 00036*11 C * * 00037*11 C * THE COVARIANCE FUNCTIONS USED IN THIS PROGRAM WERE * 00038*11 C * COMPUTED BY SCHWARZ AND LACHAPELLE (1979), "LOCAL CHARACTERIS- * 00039*11 C * TICS OF THE GRAVITY ANOMALY COVARIANCE FUNCTION". BULLETIN * 00040*11 C * GEODESIQUE 54 (1980) PP. 21 - 36.; * 00041*11 C * * 00042*11 C * "EMPIRICAL DETERMINATION OF THE GRAVITY ANOMALY COVARIANCE * 00043*11 C * FUNCTION IN MOUNTAINOUS AREAS". TO BE PUBLISHED IN SEPTEMBER * 00044*11 C * 1980 ISSUE OF THE CANADIAN SURVEYOR. * 00045*11 C * * 00046*11 C * * 00047*11 C * * 00048*11 C * REQUIRED INPUT DATA: * 00049*11 C * -------------------- * 00050*11 C * * 00051*11 C * * 00052*11 C * 1. NAME - DESCRIBING THE TERRAIN TYPE - READ AS A CHARACTER * 00053*12 C * STRING * 00054*12 C * FORMAT(25A1) * 00055*11 C * * 00056*11 C * 2. PARAMETERS SPECIFYING: * 00057*11 C * - WHETHER (OR NOT) MEAN MEAN ANOMALY OF DATA SHOULD BE * 00058*11 C * COMPUTED; * 00059*11 C * - WHETHER (OR NOT) TO PLOT DATA TO SHOW THEIR DISTRIBUTION * 00060*11 C * AROUND THE RPEDICTION POINT; * 00061*11 C * - WHETHER (OR NOT) TO USE BOUGUER ANOMALY * 00062*11 C * - NUMBER OF POINTS (+1) DEFINING POLYGON USED TO SELECT DATA* 00063*11 C * THROUGH SUBPROGRAM 'INPOLY'; * 00064*11 C * - WHETHER (OR NOT) TO PRINT INTERMEDIATE RESULTS; * 00065*11 C * - WHETHER (OR NOT) TO SCREEN DATA; * 00066*11 C * - WHETHER (OR NOT) TO USE REGRESSED DATA IN THE MOUNTAINOUS * 00067*11 C * TERRAIN. * 00068*11 C * * 00069*11 C * FORMAT(7I2) * 00070*11 C * * 00071*11 C * 3. PARAMETERS FOR SUBROUTINE 'COVAX' : * 00072*11 C * THE VARIOUS VALUES OF ESSENTIAL PARAMETERS DERIVED FROM * 00073*11 C * CANADIAN GRAVITY DATA BY SCHWARZ AND LACHAPELLE (1979) * 00074*11 C * ARE: * 00075*11 C * * 00076*11 C * GLOBAL: * 00077*11 C * 0.999617D00; 425.28D00; 24; 0; T; 2; 2; F; 2; 0.0D00; F * 00078*11 C * FLAT TERRAIN: * 00079*11 C * 0.99895D00; 181.19D00; 24; 0; T; 2; 2; F; 2; 0.0D00; F * 00080*11 C * MOUNTAINOUS TERRAIN (REGRESSED DATA) : * 00081*11 C * 0.999351D00; 131.00D00; 24; 0; T; 2; 2; F; 2; 0.0D00; F * 00082*11 C * * 00083*11 C * FORMAT(2D14.7,2I5,L2,I5,D14.7,L2) * 00084*11 C * * 00085*11 C * 4. PARAMETERS SPECIFYING: * 00086*11 C * - EXTENT OF BLOCK FOR PREDICTION (1 DEGREE, SAY) * 00087*11 C * - RADIAL INCREMENT FOR THE SELECTION OF CLOSEST DATA * 00088*11 C * - SCALE FOR X,Y COORDINATES * 00089*11 C * * 00090*11 C * FORMAT(3D9.2) * 00091*11 C * * 00092*11 C * 5. REGRESSION PARAMETERS: * 00093*11 C * - SLOPE OF REGRESSION LINE * 00094*11 C * - INTERCEPT OF LINE * 00095*11 C * - MEAN OF HEIGHT FREE GRAVITY ANOMALIES * 00096*11 C * * 00097*11 C * FORMAT(3D14.7) * 00098*11 C * * 00099*11 C * 6. VARIABLES FOR CHI-SQUARE TEST * 00100*11 C * - NUMBER OF CLASSES/SEGMENTS * 00101*11 C * - NUMBER OF POPULATION PARAMETERS KNOWN * 00102*11 C * - CHI-SQUARE VALUE FROM TABLES * 00103*11 C * * 00104*11 C * FORMAT(2I2,D10.3) * 00105*11 C * * 00106*11 C * 7. FACTOR FOR WITHIN CONTEXT TEST OF OUTLIERS * 00107*11 C * * 00108*11 C * FORMAT(D12.3) * 00109*11 C * * 00110*11 C * 8. SCALE FACTOR FOR VARIANCE OF PREDICTED ANOMALY TO CHECK * 00111*11 C * OPTIMISM OR PESSIMISM * 00112*11 C * * 00113*11 C * FORMAT(D12.5) * 00114*11 C * * 00115*11 C * 9. LIMIT OF SEPARATION BETWEEN PAIRS OF DATA FOR SCREENING * 00116*11 C * PURPOSE * 00117*11 C * * 00118*11 C * FORMAT(D10.5) * 00119*11 C * * 00120*11 C * 10. * 00121*11 C * - IDENTIFICATION NUMBER OF PREDICTION POINT * 00122*11 C * - POSITION COORDINATES OF:- LATITUDE (+VE N), LONGITUDE * 00123*11 C * (+VE W) IN DECIMAL DEGREE * 00124*11 C * - OBSERVED FREE-AIR ANOMALY (MGAL) * 00125*11 C * - OBSERVED BOUGUER ANOMALY (MGAL) * 00126*11 C * - HEIGHT (METRES) * 00127*11 C * - VARIANCE OF FREE-AIR ANOMALY (MGAL**2) * 00128*11 C * - VARIANCE OF BOUGUER ANOMALY (MGAL**2) * 00129*11 C * - VARIANCE OF HEIGHT (METRE**2) * 00130*11 C * * 00131*11 C * FORMAT(I5,5F10.5,2F6.2,F8.2) * 00132*11 C * * 00133*11 C * 11. * 00134*11 C * GRAVITY DATA FILES SPECIFIED BY THEIR APPROPRIATE JCL CARDS* 00135*11 C * * 00136*11 C * THE DATA USED IN THIS PROGRAM ARE RESIDING ON SEVEN DATA SETS, * 00137*12 C * THE DETAILS OF WHICH WOULD BE FOUND IN NASSAR AND VANICEK, * 00138*12 C * (1975). * 00139*12 C * "LEVELLING AND GRAVITY". LECTURE NOTES NO. 33, DEPARTMENT OF * 00140*11 C * SURVEYING ENGINEERING, U.N.B., FREDERICTON. * 00141*11 C * * 00142*11 C * * 00143*11 C * OUTPUT: * 00144*11 C * ------- * 00145*11 C * * 00146*11 C * 1. * 00147*11 C * - IDENTIFICATION NUMBER OF PREDICTION POINT * 00148*11 C * - POSITION COORDINATES - LATITUDE (+VE N), LONGITUDE (+VE * 00149*11 C * - W) IN DECIMAL DEGREE * 00150*11 C * - OBSERVED ANOMALY (MGAL) * 00151*11 C * - PREDICTED ANOMALY (MGAL) * 00152*11 C * - STANDARD DEVIATION OF OBSERVED ANOMALY (MGAL) * 00153*11 C * - STANDARD DEVIATION OF PREDICTED ANOMALY (MGAL) * 00154*11 C * - DIFFERENCE BETWEEN OBSERVED AND PREDICTED ANOMALY (MGAL) * 00155*11 C * - RATIO OF DIFFERENCE TO QUADRATIC SUM OF STANDARD DEVIA- * 00156*11 C * TIONS OF OBSERVED AND PREDICTED ANOMALIES * 00157*11 C * - NUMBER OF POINT GRAVITY DATA USED FOR THE PREDICTION. * 00158*11 C * * 00159*11 C * 2. AT THE END OF A SERIES OF PREDICTIONS: * 00160*11 C * - MEAN OF DIFFERENCES .... MD * 00161*11 C * - STANDARD DEVIATION OF THE MEAN .... S * 00162*11 C * * 00163*11 C * 3. RESULTS OF: * 00164*11 C * - TEST OF OUTLIERS (PERHAPS WITH NEW MD & S) * 00165*11 C * - TEST ON MEAN DIFFERENCE * 00166*11 C * - CHI-SQUARE GOODNESS OF FIT TEST (PRECEDED BY HISTOGRAM) * 00167*11 C * - TESTS FOR SKEWNESS AND KURTOSIS * 00168*11 C * - TEST FOR OPTIMISM OR PESSIMISM * 00169*11 C * * 00170*11 C * 4. PLOT OF DATA TO CHECK THEIR DISTRIBUTION (IF REQUIRED) * 00171*11 C * * 00172*11 C * 5. INTERMEDIATE RESULTS FROM SUBROUTINES 'GRAV' (IF REQUIRED) * 00173*11 C * * 00174*11 C * A COMMON BLOCK CONTAINING THE ESSENTIAL PARAMETERS OF * 00175*11 C * COVARIANCE FUNCTION ((INPUT FOR SUBROUTINE COVAX) IS REQUIRED * 00176*11 C * BETWEEN THE MAIN PROGRAM AND SUBROUTINE 'GRAV'. IN ORDER THAT * 00177*11 C * THE SUBROUTINE COVAX MAY BE USED IN THIS PROGRAM, IT HAS TO BE * 00178*11 C * LINKED WITH THE LIBRARY PACKAGE IN WHICH 'COVAX' IS STORED. * 00179*11 C * DETAILS ABOUT THE NECESSARY JCL CARDS ARE OBTAINABLE FROM * 00180*11 C * THE DEPARTMENT'S COMPUTER LIBRARY. * 00181*11 C * * 00182*11 C ******************************************************************* 00183*11 C 00184*11 C 00185 C 00186 IMPLICIT REAL*8(A-H,O-Z),LOGICAL (L) 00187 COMMON /CMCOV/CI(12),CR(51),SIGMA0(300),SIGMA(300), 00188 *KI(25),N1,LOCAL 00189 DIMENSION DLT(50),DLN(50),X(50),VDG(50),CSX(50),CXX(2550), 00190**3 *CX(50,50),CX2(50,50),CSXA(1,50),CSXB(1,50),CX1(50,50),CSX1(50,1) 00191**2 DIMENSION XI(9),YI(9),ALATT(1500),ALOGG(1500),DGA(1500),HTA(1500),00192 *DGAC(1500),DIS(1500),TX(5),TY(5),ALT(1500),ALG(1500),DG(1500), 00193**3 *DGC(1500),TDS(1500),ALTT(1500),ALGG(1500),DGG(1500),DGGC(1500), 00194**3 *TDSS(1500),DG3(300),GD(50),PWT(50),XFY(1500),YAL(1500) 00195**3 *,DG4(300),SEG(10),SNORM(300),Z(300) 00196**9 DATA GM,RE/3.98D14,6371.0D3/, 00197 *D0,D1,D2/0.0D0,1.0D0,2.0D0/,PI/3.1415926535D0/ 00198 *,ANAD,BNAD/6378.2064D0,6356.5838D0/ 00199 REAL AIR(8) 00200**2 INTEGER HVEC(20),MINE(25),IQ(4),NB(11),ICHI(10),JNC(50) 00201**9 C 00202 C... CHECK MEAN ANOMALY ALONG HIGHWAY: 00203 C READ PARAMETER ICHECK. IF EQUAL 1,CHECK MEAN ANOMALIES. 00204 C IF EQUAT 0,NO CHECK NECESSARY: 00205 C 00206 READ(5,31) MINE 00207**2 31 FORMAT(25A1) 00208**2 WRITE(6,32) (MINE(I),I=1,25) 00209**2 32 FORMAT(/,50X,25A1//) 00210**2 READ(5,500) ICHECK,IPLOT,IBG,NN1,IPRINT,ISCRN,IRGR 00211**4 500 FORMAT(7I2) 00212**4 WRITE(6,503) ICHECK,IPLOT,IBG,IPRINT,ISCRN,IRGR 00213**4 503 FORMAT(//,10X,'MEAN ANOMALY CHECK=',I5/,24X,'PLOT=',I5,/,10X,'CONV00214 *ERT TO BOUGUER=',I5/,13X,'PRINT DATA USED=',I5/,17X,'SCREEN DATA='00215**4 *,I5//,1X,'USE HEIGHT INDEPENDENT DATA=',I5/,5X,18('=')//) 00216**4 C READ THE PARAMETERS FOR SUBROUTINE COVV: 00217 C 00218 103 READ(5,21) S,A,KI(3),KI(4),LOCAL,N,KT,LSUM,N2,HMAX,LAST1 00219 21 FORMAT(2D14.7,2I5,L2,2I5,L2,I5,D14.7,L2) 00220 READ(5,22) EXTENT,EXTEND,SCALXY 00221 22 FORMAT(3D9.2) 00222 CONVRT=110.4656D 00 00223**2 EXTEND=EXTEND*CONVRT 00224**2 EXTEND=EXTEND*0.5 00225**3 IF(N.LT.2) N2=2 00226 IF(N2.GT.2001) N2=2000 00227 N2=N2+1 00228 KI(5)=KT 00229 WRITE(6,23) S,A,KI(3),KI(4),N,KT 00230 23 FORMAT('PARAMETERS SPECIFYING THE MODEL DEGREE VARIANCES:',/, 00231 *' S,A= ',2D14.7,/,' K0-K3,N,KT= -2 -1 ',4I5/) 00232 READ(5,28) B0,AH0,HFMNAN 00233**4 28 FORMAT(3D14.7) 00234**4 WRITE(6,29) B0,AH0,HFMNAN 00235**4 29 FORMAT(/,1X,'REGRESSION PARAMETERS:'/,2X,'B0,AH0=',2D14.7/, 00236**4 *2X,'MEAN H/F ANOMALY=',D14.7/) 00237**4 C... READ VARIABLES FOR CHI-SQUARE TEST: 00238**9 READ(5,1000) NKI,NKI1,CHI 00239**9 1000 FORMAT(2I2,D10.3) 00240**9 IF(NKI.EQ.2) GO TO 1001 00241**9 SEG(1)=0.431D0 00242**9 SEG(2)=-0.431D0 00243**9 GO TO 1002 00244**9 1001 SEG(1)=0.0D0 00245**9 1002 CONTINUE 00246**9 C... FACTOR FOR REJECTION OF OUTLIERS: 00247**9 READ(5,9) REJ 00248**9 9 FORMAT(D12.3) 00249**9 C... SCALING FACTOR FOR COMPUTED VARIANCE OF DIFFERENCE: 00250**9 READ(5,1500) FACTOR 00251**9 1500 FORMAT(D12.5) 00252**9 C 00253**9 READ(5,20) DLIMIT 00254**8 20 FORMAT(D10.5) 00255**8 IF(ICHECK.EQ.1) GO TO 510 00256**4 FEETMS=0.3048D0 00257 HF=0.1119D0 00258 TESTV=1.0D-16 00259**3 AREA=EXTENT/2.D0 00260 RHODEG=PI/180.0 00261 ECC1SQ=(ANAD+BNAD)*(ANAD-BNAD)/(ANAD*ANAD) 00262 C 00263 501 CONTINUE 00264 C SET-UP UPPER&LOWER LIMITS IN LATITUDE & LONGITUDE: 00265 C FOR DATA SETS: 00266 C 00267 PHIL21=41.5D 0 00268 PHIU21=63.5D 0 00269 ALGL21=52.5D 0 00270 ALGU21=73.5D 0 00271 PHIL22=43.5D 0 00272 PHIU22=63.5D 0 00273 ALGL22=73.5D 0 00274 ALGU22=91.5D 0 00275 PHIL23=48.5D 0 00276 PHIU23=63.5D 0 00277 ALGL23=91.5D 0 00278 ALGU23=118.5D 0 00279 PHIL24=48.5D 0 00280 PHIU24=63.5D 0 00281 ALGL24=118.5D 0 00282 ALGU24=140.5D 0 00283 PHIL25=63.5D 0 00284 PHIU25=70.5D 0 00285 ALGL25=61.5D 0 00286 ALGU25=142.5D 0 00287 PHIL16=70.5D 0 00288 PHIU16=80.5D 0 00289 ALGL16=66.5D 0 00290 ALGU16=142.5D 0 00291 PHIL17=80.5D 0 00292 PHIU17=82.5D 0 00293 ALGL17=51.5D 0 00294 ALGU17=78.5D 0 00295 C 00296 IF(ICHECK.EQ.1) GO TO 502 00297 DO 24 I=1,300 00298 24 SIGMA0(I)=0.0D0 00299 RB2=RE*RE*S 00300 LMODEL=N.EQ.0 00301 CI(8)=A*RB2*1.0D-10 00302 CI(10)=S 00303 IF(N.NE.0) GO TO 25 00304 LOCAL=.TRUE. 00305 N=2 00306 25 N1=N+1 00307 IF(.NOT.LOCAL) READ(5,26) (SIGMA0(I),I=1,N1) 00308 26 FORMAT(9F8.5) 00309 COKO=1.0D0 00310 DO 121 I=1,N1 00311 121 SIGMA0(I)=SIGMA0(I)*COKO 00312 IF(.NOT.LOCAL) WRITE(6,27) (SIGMA0(I),I=1,N1) 00313 27 FORMAT(' EMPIRICAL ANOMALY DEGREE VARIANCES IN UNITS OF MGAL**2:',00314 */,25(12F6.2/)) 00315 N2=2001 00316 CALL COVAX 00317 KI(6)=3 00318 KI(7)=3 00319 CALL COVBX 00320 CR(2)=0.0D0 00321 CR(3)=0.0D0 00322 CR(10)=GM/(RE+CR(2))**2 00323 CR(11)=GM/(RE+CR(3))**2 00324 ALAX=0.0 00325 ALAN=0.0 00326 ALOX=0.0 00327 ALN=0.0 00328 IUNIT1=40 00329 C 00330 C PRINT-UUT SOME HEADINGS FOR THE JOB: 00331 C 00332 WRITE(6,8) 00333 8 FORMAT(//,35X,'GRAVITY ANOMALY PREDICTION AT BENCH MARKS'/,35X, 00334 *41('=')//,10X,'(1)',5X,'(2)',7X,'(3)',8X,'(4)',7X,'(5)',8X,'(6)', 00335 *7X,'(7)',8X,'(8)',7X,'(9)',/,10X,'BM-#',2X,'LATITUDE',1X,'LONGITUD00336 *E',2X,'OBSERVED',2X,'PREDICTED',2X,'SIG(OBS)',1X,'SIG(PRED)',2X,'(00337 *4)-(5)',3X,'(8)/((7)**.5)') 00338 IIV=0 00339 IJ=0 00340 ITML=0 00341**8 XTRM1=120.0D0 00342**8 XTRM2=129.5D0 00343**8 999 CONTINUE 00344 IF(ICHECK.EQ.1) GO TO 502 00345 C READ-IN CENTRE POIT OF AREA UNDER INVESTIATION: 00346 READ(5,30,END=2000) ID,(AIR(I),I=1,8) 00347**2 30 FORMAT(I5,5F10.5,2F6.2,F8.2) 00348**2 XLAT=AIR(1) 00349**2 XLOG=AIR(2) 00350**2 DGK=AIR(3) 00351**2 DGBK=AIR(4) 00352**2 HT1=AIR(5) 00353**2 WFAINV=AIR(6) 00354**2 WFBINV=AIR(7) 00355**2 VVX=AIR(8) 00356**2 367 FORMAT(//) 00357 IF(IPLOT.EQ.0) GO TO 87 00358 C... GENERATE ELEMENTS FOR SUBROUTINE GC-PLOT: 00359 XLATR=XLAT*RHODEG 00360 DEN0=1.0D0-ECC1SQ*DSIN(XLATR)**2 00361 RCN=ANAD/DSQRT(DEN0) 00362 RCM=RCN*(1.0D0-ECC1SQ)/DEN0 00363 ROM=DSQRT(RCN*RCM) 00364 FYFACT=ROM*RHODEG 00365 ALFACT=DCOS(XLATR)*FYFACT 00366 SCALE=0.25 00367 87 CONTINUE 00368 IEXT=0 00369 88 CONTINUE 00370 IEXT=IEXT+1 00371 AREA=EXTENT/2. 00372 ALATMX =XLAT+AREA 00373 ALATMN =XLAT-AREA 00374 ALONMX =XLOG+AREA 00375 ALONMN=XLOG-AREA 00376 C 00377 C CHECK GIVEN LOCATION AGAINST LIMITS AND CHOSA THE RELEVANT FILE 00378 C FOR INTERPOLATION: 00379 C 00380 IF((XLOG.GE.ALGL21).AND.(XLOG.LE.ALGU21)) GO TO 1 00381 IF((XLOG.GE.ALGL22).AND.(XLOG.LE.ALGU22)) GO TO 2 00382 IF((XLOG.GE.ALGL23).AND.(XLOG.LE.ALGU23)) GO TO 3 00383 IF((XLOG.GE.ALGL24).AND.(XLOG.LE.ALGU24)) GO TO 4 00384 C 00385 C PRINT-OUT MESSAGE: FOR BENCH MARKS OUTSIDE RANGE OF OBSERVED GRAVITY: 00386 IF(IEXT.GT.1) GO TO 90 00387 WRITE(6,10) 00388 10 FORMAT(///,10X,'NO G/ANOMALY VALUE FOUND:'/,10X,'EXTENT IS DOUBLED00389 1FOR THE PREDICTION.'//) 00390 EXTENT=2*EXTENT 00391 GO TO 88 00392 90 WRITE(6,102) 00393 102 FORMAT(/,32X,25('*')//) 00394 GO TO 9999 00395 1 IF((XLAT.GE.PHIL21).AND.(XLAT.LE.PHIU21)) GO TO 11 00396 IF((XLAT.GE.PHIL25).AND.(XLAT.LE.PHIU25)) GO TO 15 00397 IF((XLAT.GE.PHIL16).AND.(XLAT.LE.PHIU16)) GO TO 6 00398 IF((XLAT.GE.PHIL17).AND.(XLAT.LE.PHIU17)) GO TO 7 00399 WRITE(6,10) 00400 GO TO 9999 00401 2 IF((XLAT.GE.PHIL22).AND.(XLAT.LE.PHIU22)) GO TO 12 00402 IF((XLAT.GE.PHIL25).AND.(XLAT.LE.PHIU25)) GO TO 15 00403 IF((XLAT.GE.PHIL16).AND.(XLAT.LE.PHIU16)) GO TO 6 00404 IF((XLAT.GE.PHIL17).AND.(XLAT.LE.PHIU17)) GO TO 7 00405 WRITE(6,10) 00406 GO TO 9999 00407 3 IF((XLAT.GE.PHIL23).AND.(XLAT.LE.PHIU23)) GO TO 13 00408 IF((XLAT.GE.PHIL25).AND.(XLAT.LE.PHIU25)) GO TO 15 00409 IF((XLAT.GE.PHIL16).AND.(XLAT.LE.PHIU16)) GO TO 6 00410 WRITE(6,10) 00411 GO TO 9999 00412 4 IF((XLAT.GE.PHIL24).AND.(XLAT.LE.PHIU24)) GO TO 14 00413 IF((XLAT.GE.PHIL25).AND.(XLAT.LE.PHIU25)) GO TO 15 00414 IF((XLAT.GE.PHIL16).AND.(XLAT.LE.PHIU16)) GO TO 6 00415 WRITE(6,10) 00416 GO TO 9999 00417 C 00418 502 CONTINUE 00419 C 00420 C 00421 C SPECIFY DATA FILE NUMBERS: 00422 C 00423 11 IUNIT=21 00424 GO TO 5 00425 12 IUNIT=22 00426 GO TO 5 00427 13 IUNIT=23 00428 GO TO 5 00429 14 IUNIT=24 00430 GO TO 5 00431 15 IUNIT=25 00432 GO TO 5 00433 6 IUNIT=16 00434 GO TO 5 00435 7 IUNIT=17 00436 5 CONTINUE 00437 C 00438 IF(ICHECK.EQ.1) GO TO 510 00439 600 CONTINUE 00440 C 00441 C... REMOVE SYSTEMATIC EFFECT FROM ANOMALY: 00442 C 00443 SDG=-12.74337 00444**3 IF(IUNIT.EQ.24) SDG=-7.17452 00445 C 00446 DO 51 J=1,4 00447 IQ(J)=0 00448 51 CONTINUE 00449 C 00450 RADIUS=EXTENT/2.D0 00451 RR=RADIUS/DSQRT(2.D0) 00452 PLOTFY=RADIUS*FYFACT 00453**3 PLOTAL=RADIUS*ALFACT 00454**3 XFYMIN=-PLOTFY 00455**3 YALMIN=-PLOTAL 00456**3 XFYMAX=+PLOTFY 00457**3 YALMAX=+PLOTAL 00458**3 C... 00459 XI(1)=XLAT 00460 XI(2)=XLAT+RR 00461 XI(3)=XLAT+RADIUS 00462 XI(4)=XI(2) 00463 XI(5)=XI(1) 00464 XI(6)=XLAT-RR 00465 XI(7)=XLAT-RADIUS 00466 XI(8)=XI(6) 00467 XI(9)=XI(1) 00468 C... 00469 YI(1)=XLOG-RADIUS 00470 YI(2)=XLOG-RR 00471 YI(3)=XLOG 00472 YI(4)=XLOG+RR 00473 YI(5)=XLOG+RADIUS 00474 YI(6)=YI(4) 00475 YI(7)=YI(3) 00476 YI(8)=YI(2) 00477 YI(9)=YI(1) 00478 C... 00479 REWIND IUNIT 00480 J=0 00481 50 READ(IUNIT,57,END=300) (AIR(I),I=1,8) 00482**2 57 FORMAT(5F12.5,2F6.2,F8.2) 00483**2 IF((AIR(6).EQ.0.).OR.(AIR(7).EQ.0.).OR.(AIR(8).EQ.0.)) GO TO 50 00484**2 ALAT=AIR(1) 00485**2 ALON=AIR(2) 00486**2 FA=AIR(3) 00487**2 BA=AIR(4) 00488**2 AHT=AIR(5) 00489**2 VFA=AIR(6) 00490**2 VBA=AIR(7) 00491**2 VAHT=AIR(8) 00492**2 IF(IRGR.EQ.1) GO TO 33 00493**4 GO TO 34 00494**4 33 CONTINUE 00495**4 FA=FA-B0*AHT 00496*10 VFA=VFA-VAHT*B0**2 00497**4 34 CONTINUE 00498**4 C 00499 C SELECT DATA WITHIN AREA DEFINED BY LATMAX/MIN & LONMAX/MIN: 00500 C 00501 IF((ALAT.EQ.XLAT).AND.(ALON.EQ.XLOG)) GO TO 50 00502**2 C 00503 CALL INPOLY(XI,YI,NN1,ALAT,ALON,IN) 00504 C 00505 IF(IN.EQ.1) GO TO 53 00506 GO TO 50 00507 53 CONTINUE 00508 J=J+1 00509 ALATT(J)=ALAT 00510 ALOGG(J)=ALON 00511 DGA(J)=FA 00512**2 DGAC(J)=VFA 00513**2 DIFF1=ALAT-XLAT 00514**2 DIFF2=XLOG-ALON 00515**2 DS1=DIFF1*CONVRT 00516**2 DS2=DIFF2*CONVRT*DCOS(XLATR) 00517**2 DS1SQ=DS1*DS1 00518**2 DS2SQ=DS2*DS2 00519**2 DIS(J)=DSQRT(DS1SQ+DS2SQ) 00520**2 IF(IPRINT.GT.1) GO TO 240 00521**3 GO TO 241 00522**3 240 CONTINUE 00523**3 XFY(J)=DS1 00524**3 YAL(J)=DS2 00525**3 WRITE(6,150) J,ALATT(J),ALOGG(J),DGA(J),DGAC(J),XFY(J),YAL(J),DIS(00526**4 *J) 00527**4 241 CONTINUE 00528**3 GO TO 50 00529 300 CONTINUE 00530 152 FORMAT(100X,2I10) 00531 153 FORMAT(100X,I10) 00532 154 FORMAT(100X,F10.2,I10) 00533 IF((IPLOT.EQ.0).OR.(IPRINT.LT.2)) GO TO 155 00534**4 CALL GCPLOT(YAL,XFY,J,YALMIN,XFYMIN,YALMAX,XFYMAX) 00535**3 CMFLSC=EXTENT*60./15. 00536**3 WRITE(6,554) CMFLSC 00537**3 155 CONTINUE 00538**3 WRITE(6,153) J 00539 IF(ISCRN.EQ.0) GO TO 236 00540**4 C 00541**3 C... SCREEN DATA COLLECTED; WHERE CLUSTER OF POINTS EXISTS, REPLACE 00542**3 C WITH THE STRAIGHT MEAN: 00543**3 C 00544**3 IQR=1 00545**3 KL=0 00546**3 NII=5 00547**3 221 CONTINUE 00548**3 C... 00549**3 TX(1)=XLAT 00550**3 TX(2)=TX(1) 00551**3 TX(3)=TX(1)+RADIUS 00552**3 IF((IQR.EQ.3).OR.(IQR.EQ.4)) TX(3)=TX(1)-RADIUS 00553**3 TX(4)=TX(3) 00554**3 TX(5)=TX(1) 00555**3 C... 00556**3 TY(1)=XLOG 00557**3 TY(2)=TY(1)+RADIUS 00558**3 IF((IQR.EQ.2).OR.(IQR.EQ.3)) TY(2)=TY(1)-RADIUS 00559**3 TY(3)=TY(2) 00560**3 TY(4)=TY(1) 00561**3 TY(5)=TY(1) 00562**3 C... 00563**3 IT=0 00564**3 DO 222 I=1,J 00565**3 IF(ALATT(I).EQ.0.) GO TO 222 00566**3 AAL=ALATT(I) 00567**3 AAG=ALOGG(I) 00568**3 CALL INPOLY(TX,TY,NII,AAL,AAG,IN) 00569**3 IF(IN.EQ.1) GO TO 223 00570**3 GO TO 222 00571**3 223 CONTINUE 00572**3 IT=IT+1 00573**3 ALT(IT)=ALATT(I) 00574**3 ALG(IT)=ALOGG(I) 00575**3 DG(IT)=DGA(I) 00576**3 DGC(IT)=DGAC(I) 00577**3 TDS(IT)=DIS(I) 00578**3 ALATT(I)=0. 00579**3 222 CONTINUE 00580**3 IF(IT.EQ.0) GO TO 229 00581**3 ITT=1 00582**3 ICT=1 00583**3 K=1 00584**3 IF(IT.EQ.1) GO TO 224 00585**3 KP1=K+1 00586**3 GD(ICT)=DG(K) 00587**3 IF(DGC(K).LE.TESTV) GO TO 231 00588**3 GO TO 232 00589**3 231 PWT(ICT)=1.D 16 00590**3 232 PWT(ICT)=1.D 0/DGC(K) 00591**3 DSM1=ALT(K)*PWT(ICT) 00592**3 DSM2=ALG(K)*PWT(ICT) 00593**3 DSM3=DG(K)*PWT(ICT) 00594**3 DSM4=PWT(ICT) 00595**3 DSM5=DGC(K) 00596**3 228 CONTINUE 00597**3 DO 224 I=KP1,IT 00598**3 IF(ALT(K).EQ.0.) GO TO 227 00599**3 IF(ALT(I).EQ.0.)GO TO 224 00600**3 DT1=(ALT(K)-ALT(I))*CONVRT 00601**3 DT2=(ALG(I)-ALG(K))*CONVRT*DCOS(XLATR) 00602**3 DTEST=DSQRT(DT1**2+DT2**2) 00603**3 IF(DTEST.LE.DLIMIT) GO TO 225 00604**3 IF((ITT.GT.1).AND.(I.EQ.IT)) GO TO 226 00605**3 GO TO 224 00606**3 225 CONTINUE 00607**3 ICT=ICT+1 00608**3 IF(DGC(I).LE.TESTV) GO TO 233 00609**3 GO TO 234 00610**3 233 PWT(ICT)=1.D 16 00611**3 234 PWT(ICT)=1.D 0/DGC(I) 00612**3 GD(ICT)=DG(I) 00613**3 C 00614**3 DSM1=DSM1+ALT(I)*PWT(ICT) 00615**3 DSM2=DSM2+ALG(I)*PWT(ICT) 00616**3 DSM3=DSM3+DG(I)*PWT(ICT) 00617**3 DSM4=DSM4+PWT(ICT) 00618**3 DSM5=DSM5+DGC(I) 00619**3 C 00620**3 ALT(I)=0. 00621**3 ITT=ITT+1 00622**3 IF(I.EQ.IT) GO TO 226 00623**3 GO TO 224 00624**3 C... WEIGHTED MEAN: 00625**3 226 CONTINUE 00626**3 KL=KL+1 00627**3 ALTT(KL)=DSM1/DSM4 00628**3 ALGG(KL)=DSM2/DSM4 00629**3 DGG(KL)=DSM3/DSM4 00630**3 GDSM=0. 00631**3 DO 230 IM=1,ICT 00632**3 230 GDSM=GDSM+PWT(IM)*(GD(IM)-DGG(KL))**2 00633**4 COUNT=DFLOAT(ICT-1) 00634**3 DGGC(KL)=GDSM/(DSM4*COUNT) 00635**3 DT3=(ALTT(KL)-XLAT)*CONVRT 00636**3 DT4=(XLOG-ALGG(KL))*CONVRT*DCOS(XLATR) 00637**3 TDSS(KL)=DSQRT(DT3**2+DT4**2) 00638**3 GO TO 227 00639**3 224 CONTINUE 00640**3 KL=KL+1 00641**3 ALTT(KL)=ALT(K) 00642**3 ALGG(KL)=ALG(K) 00643**3 DGG(KL)=DG(K) 00644**3 DGGC(KL)=DGC(K) 00645**3 TDSS(KL)=TDS(K) 00646**3 IF(KP1.EQ.IT) GO TO 235 00647**3 227 CONTINUE 00648**3 K=K+1 00649**3 KP1=K+1 00650**3 ITT=1 00651**3 ICT=1 00652**3 GD(ICT)=DG(K) 00653**3 PWT(ICT)=1.D0/DGC(K) 00654**3 DSM1=ALT(K)*PWT(ICT) 00655**3 DSM2=ALG(K)*PWT(ICT) 00656**3 DSM3=DG(K)*PWT(ICT) 00657**3 DSM4=PWT(ICT) 00658**3 DSM5=DGC(K) 00659**3 IF(K.LT.IT) GO TO 228 00660**3 235 CONTINUE 00661**3 IF(ALT(IT).EQ.0.) GO TO 229 00662**3 KL=KL+1 00663**3 ALTT(KL)=ALT(IT) 00664**3 ALGG(KL)=ALG(IT) 00665**3 DGG(KL)=DG(IT) 00666**3 DGGC(KL)=DGC(IT) 00667**3 TDSS(KL)=TDS(IT) 00668**3 229 CONTINUE 00669**3 IQR=IQR+1 00670**3 IF(IQR.LE.4) GO TO 221 00671**3 J=KL 00672**3 C 00673**3 C... NOW SELECT DATA FOR PREDICTION: 00674**3 C 00675**3 WRITE(6,153) J 00676**3 236 CONTINUE 00677**4 TESTER=5.0D0/60.0D0 00678**2 TESTER=TESTER*0.5 00679**3 TESTER=TESTER*CONVRT 00680**2 NDG=0 00681 54 CONTINUE 00682 IF(ISCRN.EQ.0) GO TO 237 00683**4 DO 55 I=1,J 00684 ALATT(I)=ALTT(I) 00685**3 ALOGG(I)=ALGG(I) 00686**3 DGA(I)=DGG(I) 00687**3 DGAC(I)=DGGC(I) 00688**3 DIS(I)=TDSS(I) 00689**3 55 CONTINUE 00690**4 237 CONTINUE 00691**4 DO 58 I=1,J 00692**4 IF(ALATT(I).EQ.0.0D0) GO TO 58 00693**4 IF(DIS(I).GT.TESTER) GO TO 58 00694**4 NDG=NDG+1 00695 C... 00696 DLT(NDG)=ALATT(I)*RHODEG 00697 DLN(NDG)=ALOGG(I)*RHODEG 00698 IF(IRGR.EQ.1) GOTO 302 00699**4 X(NDG)=DGA(I)-SDG 00700**4 GO TO 303 00701**4 302 CONTINUE 00702**4 X(NDG)=DGA(I) 00703**4 303 CONTINUE 00704**4 VDG(NDG)=DGAC(I) 00705**4 C... 00706 IF((ALATT(I).GE.XLAT).AND.(ALOGG(I).LT.XLOG)) IQ(1)=IQ(1)+1 00707**2 IF((ALATT(I).GE.XLAT).AND.(ALOGG(I).GE.XLOG)) IQ(2)=IQ(2)+1 00708**2 IF((ALATT(I).LT.XLAT).AND.(ALOGG(I).GE.XLOG)) IQ(3)=IQ(3)+1 00709**2 IF((ALATT(I).LT.XLAT).AND.(ALOGG(I).LT.XLOG)) IQ(4)=IQ(4)+1 00710 C... 00711 IF(IPLOT.EQ.0) GO TO 56 00712 XFY(NDG)=(ALATT(I)-XLAT)*FYFACT*SCALXY 00713 YAL(NDG)=(XLOG-ALOGG(I))*ALFACT*SCALXY 00714 56 CONTINUE 00715 IF(IPRINT.GT.0) WRITE(6,150) NDG,ALATT(I),ALOGG(I),DGA(I),DGAC(I) 00716**3 *,XFY(NDG),YAL(NDG),DIS(I) 00717**2 150 FORMAT(9X,I5,8F10.5) 00718**2 ALATT(I)=0.0D0 00719 IF(NDG.EQ.50) GO TO 42 00720**9 58 CONTINUE 00721**4 IF(NDG.LT.10) GO TO 41 00722**9 GO TO 42 00723 41 CONTINUE 00724 IF(NDG.EQ.J) GO TO 42 00725 TESTER=TESTER+EXTEND 00726 GO TO 237 00727**4 42 CONTINUE 00728 WRITE(6,154) TESTER,NDG 00729 ITEST=IQ(1)*IQ(2)*IQ(3)*IQ(4) 00730 IF(ITEST.EQ.0) GO TO 990 00731 301 CONTINUE 00732 NBRR=(NDG*(NDG+1))/2 00733 XLAT=XLAT*RHODEG 00734 XLOG=XLOG*RHODEG 00735 IF(IPLOT.EQ.0) GO TO 43 00736 WRITE(6,367) 00737 EXTD=EXTENT 00738 C***********************************************************************00739 C 00740 C CALL THE SUBROUTINE 'GCPLOT' TO PLOT A SKETCH SHOWING THE APPROXIMATE 00741 C DISTRIBUTION OF THE USED DATA POINTS WITHIN THE ONE-AREA-SQUARED 00742 C SPECIFIED AROUND THE POINT OF INTEREST, FOR THE INTERPOLATION PURPOSES00743 C 00744 C 00745 CALL GCPLOT(YAL,XFY,NDG,YALMIN,XFYMIN,YALMAX,XFYMAX) 00746 C 00747 CMFLSC=( EXTD*60.)/15. 00748 WRITE(6,554) CMFLSC 00749 554 FORMAT(5X,'SCALE OF THE ABOVE PLOT IS: 1.0 CENTIMETER REPRESENTS',00750 1F5.2,2X,'MINUTES (IN BOTH LAT.& LONG.)'/) 00751 43 CONTINUE 00752 C 00753 C CALL SUBROUTINE GRAV TO PREDICT G/ANOMALY AT B/MARKS: 00754 C 00755 CALL CPUTIM(ITI,JJR) 00756 CALL GRAV(NDG,VDG,DLT,DLN,X,CXX,CSX,CSXA,CSXB,XLAT,XLOG, 00757 *DGP,VAR,NBRR,CX,CX2,GG,CX1,CSX1,IPRINT) 00758**3 C 00759 CALL CPUTIM(ITI,JJR) 00760 C 00761**8 C... CPUTIME SUMMATION: 00762**8 C 00763**8 ITML=ITML+ITI 00764**8 XLAT=XLAT/RHODEG 00765 XLOG=XLOG/RHODEG 00766 IF(IRGR.EQ.1) GO TO 304 00767**4 DGP=DGP+SDG 00768**4 GO TO 305 00769**4 304 CONTINUE 00770**4 IF(IPRINT.GT.0) WRITE(6,402) ID,XLAT,XLOG,DGK,DGP,VAR 00771**4 DGP=DGP+B0*HT1 00772*10 VAR=VAR+VVX*B0**2 00773**4 305 CONTINUE 00774**4 IIV=IIV+1 00775 IF(IBG.EQ.1) GO TO 403 00776 DG3(IIV)=DGK-DGP 00777 SVAR=DSQRT(VAR) 00778 SVK=DSQRT(WFAINV) 00779 FACTSQ=FACTOR*FACTOR 00780**9 SS=DSQRT(FACTSQ*VAR+WFAINV) 00781**9 SNORM(IIV)=DG3(IIV)/SS 00782**9 TEST1=DABS(SNORM(IIV)) 00783**9 RATIO=SVAR/SVK 00784**4 IF(TEST1.LT.1.) IJ=IJ+1 00785 WRITE(6,402) ID,XLAT,XLOG,DGK,DGP,SVK,SVAR,DG3(IIV),TEST1,NDG 00786 *,RATIO 00787**4 IF(IPRINT.GT.0) WRITE(6,367) 00788**4 GO TO 404 00789 403 CONTINUE 00790 DGP=DGP-HF*HT1 00791 DG3(IIV)=DGBK-DGP 00792 SVK=DSQRT(WFBINV) 00793 SVAR=DSQRT(VAR+VVX*HF**2) 00794 VAR=SVAR*SVAR 00795 SS=DSQRT(WFBINV+VAR) 00796 TEST1=DABS(DG3(IIV))/SS 00797 IF(TEST1.LT.1.) IJ=IJ+1 00798 WRITE(6,402) ID,XLAT,XLOG,DGBK,DGP,SVK,SVAR,DG3(IIV),TEST1,NDG 00799 402 FORMAT( 9X,I5,8F10.5,I5,F6.1) 00800**4 404 CONTINUE 00801 ALAX=ALATMX 00802 ALAN=ALATMN 00803 ALOX=ALONMX 00804 ALN=ALONMN 00805 IUNIT1=IUNIT 00806 GO TO 9999 00807 C***********************************************************************00808 510 CONTINUE 00809 IUNIT=24 00810**4 WRITE(6,516) IUNIT 00811**4 516 FORMAT(//,25X,'***FILE#',I4,'***'//) 00812**4 ALATSM=0. 00813**4 ALONSM=0. 00814**4 FASM=0. 00815**4 IKOUNT=0 00816**4 REWIND IUNIT 00817**4 511 READ(IUNIT,512,END=513) (AIR(I),I=1,8) 00818**4 512 FORMAT(2F10.5,2F11.5,F8.2,3F9.4) 00819**4 AALAT=AIR(1) 00820**4 AALON=AIR(2) 00821**4 FA=AIR(3) 00822**4 AHT=AIR(5) 00823**4 FA=FA-B0*(AHT-AH0) 00824**4 C... 00825**4 IKOUNT=IKOUNT+1 00826**4 ALATSM=ALATSM+AALAT 00827**4 ALONSM=ALONSM+AALON 00828**4 FASM=FASM+FA 00829**4 GO TO 511 00830**4 513 CONTINUE 00831**4 ACOUNT=DFLOAT(IKOUNT) 00832**4 CLT=ALATSM/ACOUNT 00833**4 CLN=ALONSM/ACOUNT 00834**4 CFA=FASM/ACOUNT 00835**4 WRITE(6,514) CLT,CLN,CFA 00836**4 514 FORMAT(//,25X,'MEAN LATITUDE=',F12.5/,25X,'MEAN LONGITUDE=',F12.5/00837**4 *,25X,'MEAN HEIGHT INDEPENDENT ANOMALY=',F12.5/,30X,18('=')//) 00838**4 WRITE(6,515) IKOUNT 00839**4 515 FORMAT(25X,'TOTAL NUMBER OF ANOMALIES USED=',I10//) 00840**4 GO TO 99 00841**4 2000 CONTINUE 00842 SUMM=0.0D0 00843 SUMM1=0.0D0 00844 VSM=0.D0 00845**8 IPCT=0 00846**8 DO 543 IB=1,IIV 00847 SUMM=SUMM+DG3(IB) 00848 VSM=VSM+DG3(IB)**2 00849**8 543 CONTINUE 00850 AIV=DFLOAT(IIV) 00851**9 VSM1=VSM 00852**9 VSM=VSM/AIV 00853**9 SUMM2=SUMM*SUMM/AIV 00854**9 DMEAN=SUMM/AIV 00855**9 RMEAN=(VSM1-SUMM2)/(AIV-1.) 00856**9 RMSQ=DSQRT(RMEAN) 00857 RMS2=DSQRT(VSM) 00858**8 WRITE(6,544) DMEAN 00859 544 FORMAT(//,47X,'MEAN OF DIFFERENCES=',F10.5) 00860 WRITE(6,545) SUMM1,RMEAN,RMSQ,RMS2 00861**8 545 FORMAT(//,47X,'SUM SQUARE ERROR=',F12.5/,46X,'MEAN SQUARE ERROR=',00862 *F12.5/,41X,'ROOT MEAN SQUARE ERROR=',F12.5//,40X,'RMS ERROR WRPT P00863**8 *OP-MEAN=',F12.5) 00864**8 DO 546 I=1,IIV 00865 TEST4=DG3(I)/RMS2 00866**8 DG4(I)=TEST4 00867**8 IF(TEST4.LT.1) IPCT=IPCT+1 00868**8 546 CONTINUE 00869 TEST3=(IJ/DFLOAT(IIV))*100. 00870 TEST5=(IPCT/DFLOAT(IIV))*100. 00871**8 WRITE(6,548) IIV,IJ,TEST3,TEST5 00872**8 548 FORMAT(//,41X,'TOTAL # OF TEST POINTS=',I5/,52X,'PASS POINTS=',I5/00873**8 *,53X,'PERCENTAGE=',F8.2//,42X,'PERCENT WRPT POP-MEAN=',F8.2) 00874**8 C 00875**3 C... TOTAL CPUTIME IN SECONDS: 00876**8 TIME=DFLOAT(ITML)/1.0D 04 00877**8 WRITE(6,549) TIME 00878**8 549 FORMAT(//,33X,'TOTAL CPU TIME FOR PREDICTIONS=',F12.5//,33X,43('*'00879**8 *)) 00880**8 CALL REJECT(DG3,IIV,RMSQ,REJ,JNC,NLCN,50) 00881**9 IF(NLCN.EQ.0) GO TO 808 00882**9 WRITE(6,807) 00883**9 DO 805 I=1,NLCN 00884**9 IRJ=JNC(I) 00885**9 WRITE(6,806) IRJ,DG3(IRJ) 00886**9 SUMM=SUMM-DG3(IRJ) 00887**9 VSM1=VSM1-DG3(IRJ)**2 00888*10 805 CONTINUE 00889**9 IIV=IIV-NLCN 00890**9 AIV=DFLOAT(IIV) 00891**9 SUMM2=SUMM*SUMM/AIV 00892**9 DMEAN=SUMM/AIV 00893**9 RMEAN=(VSM1-SUMM2)/(AIV-1.) 00894**9 RMSQ=DSQRT(RMEAN) 00895**9 WRITE(6,809) 00896**9 WRITE(6,810) DMEAN,RMSQ 00897**9 808 CONTINUE 00898**9 806 FORMAT(30X,'OBSERVATION #',I4,2X,'(VALUE:',F7.2,') IS REJECTED') 00899**9 807 FORMAT(//,30X,'REJECTED OBSERVATIONS'/,30X,21('=')//) 00900**9 809 FORMAT(//,40X,'NEW MEAN AND RMS VALUES:'/,40X,23('=')//) 00901**9 810 FORMAT(40X,'MEAN DIFFERENCE=',F10.5/,40X,'STANDARD DEVIATION=', 00902**9 *F10.5//) 00903**9 C 00904**8 C... HISTOGRAM PLOT ROUTINE: 00905**3 C 00906**3 C 00907**9 C... STUDENTS-T TEST ON MEAN DIFFERENCE: 00908**9 CALL TONMN(DMEAN,RMSQ,IIV,1) 00909**9 DO 601 I=1,IIV 00910**9 Z(I)=DG3(I)/RMSQ 00911**9 601 CONTINUE 00912**9 CALL HSTPLT(Z,IIV,NCLASS,HVEC,NB) 00913*12 WRITE(6,129) 00914**9 C 00915**9 C... CHI-SQUARE GOODNESS OF FIT TEST: 00916**9 C 00917**9 CALL CHISQ(Z,IIV,SEG,ICHI,NKI,CHIVAL,CHI,NKI1) 00918**9 C 00919**9 IAN=IIV 00920**9 AM1=0.D0 00921**9 AM2=0.D0 00922**9 AM3=0.D0 00923**9 AM4=0.D0 00924**9 DO 1007 I=1,IAN 00925**9 DELMN=DG3(I)-DMEAN 00926**9 DMN2=DELMN*DELMN 00927**9 DMN3=DMN2*DELMN 00928**9 DMN4=DMN3*DELMN 00929**9 AM1=AM1+DABS(DELMN) 00930**9 AM2=AM2+DMN2 00931**9 AM3=AM3+DMN3 00932**9 AM4=AM4+DMN4 00933**9 1007 CONTINUE 00934**9 RIAM=DFLOAT(IAN) 00935**9 AM0=DSQRT(AM2) 00936**9 AM2=AM2/RIAM 00937**9 AM3=AM3/RIAM 00938**9 AM4=AM4/RIAM 00939**9 ALPHA3=AM3/((DSQRT(AM2))**3) 00940**9 ALPHA4=AM4/(AM2*AM2) 00941**9 A4=AM1/(RIAM*AM0) 00942**9 WRITE(6,1008) ALPHA3,ALPHA4,A4 00943**9 1008 FORMAT(///,40X,'MEASURE OF SKEWNESS=',F12.3/,40X,'MEASURE OF KURTO00944**9 *SIS=',F12.3,5X,F12.3//) 00945**9 J0=IIV/10.+1 00946**9 J1=1 00947**9 J2=J1+9 00948**9 DO 1003 I=1,J0 00949**9 WRITE(6,1004) (SNORM(J),J=J1,J2) 00950**9 J1=J2+1 00951**9 J2=J1+9 00952**9 1003 CONTINUE 00953**9 1004 FORMAT(10X,10F10.3) 00954**9 CALL HSTPLT(SNORM,IIV,NCLASS,HVEC,NB) 00955*12 WRITE(6,137) 00956**9 CALL AMNSDR(SNORM,IIV,AMEAN,SDEV) 00957**9 WRITE(6,138) AMEAN,SDEV 00958**9 CALL CHISQ(SNORM,IIV,SEG,ICHI,NKI,CHIVAL,CHI,NKI1) 00959**9 C 00960**9 C... TEST ON MEAN: 00961**9 CALL TONMN(AMEAN,SDEV,IIV,0) 00962**9 C 00963**9 C... TEST ON VARIANCE: 00964**9 C 00965**9 CHI0=1.D0 00966**9 CHI1=83.3D0 00967**9 CHI2=40.48D0 00968**9 SDEVSQ=SDEV*SDEV 00969**9 ANVAR=IIV*SDEVSQ 00970**9 ANVAR1=ANVAR/CHI1 00971**9 ANVAR2=ANVAR/CHI2 00972**9 WRITE(6,124) 00973**9 WRITE(6,125) ANVAR1,CHI0,ANVAR2 00974**9 IF((CHI0.GT.ANVAR1).AND.(CHI0.LT.ANVAR2)) GO TO 1005 00975**9 WRITE(6,126) 00976**9 GO TO 1006 00977**9 1005 WRITE(6,127) 00978**9 1006 CONTINUE 00979**9 124 FORMAT(//,40X,'TEST ON VARIANCE'/,40X,14('=')//) 00980**9 125 FORMAT(40X,F12.5,2X,'<',F12.5,2X,'<',F12.5) 00981**9 126 FORMAT(//,40X,'** TEST FAILS **'//) 00982**9 127 FORMAT(//,40X,'** TEST PASSES **'//) 00983**9 129 FORMAT(//,40X,'HISTOGRAM PLOT OF DIFFERENCES'/,40X,29('=')//) 00984**9 130 FORMAT(//,30X,F5.2,4X,F5.2,6X,I2,12X,I2,12X,F5.1,8X,F5.2) 00985**9 131 FORMAT(//,39X,F5.2,6X,I2,12X,I2,12X,F5.1,8X,F5.2) 00986**9 135 FORMAT(//,91X,F5.2) 00987**9 136 FORMAT(//,37X,'SEGMENT',3X,'OBSERVED',6X,'EXPECTED',8X,'(O-E)',6X,00988**9 *'(O-E)**2/E') 00989**9 137 FORMAT(//,40X,'HISTOGRAM PLOT OF STANDARDIZED DIFFERENCES'/,40X,4200990**9 *('=')//) 00991**9 138 FORMAT(40X,'MEAN VALUE=',F10.2/,40X,'STD. DEV=',F10.2//) 00992**9 GO TO 99 00993 990 WRITE(6,991) NDG,(IQ(J),J=1,4),ID 00994 991 FORMAT(//,10X,'DATA DISTRIBUTION AROUND INTERPOLATION POINT IS NOT00995 * UNIFORM:'/,10X,'TOTAL # OF DATA POINTS =',I5,':',5X,4(I3,5X),': 00996 *BM-#',I5/) 00997 JJ=0 00998 DO 992 J=1,4 00999 IF(IQ(J).EQ.0) JJ=JJ+1 01000 992 CONTINUE 01001 IF(JJ.EQ.1) GO TO 301 01002 9999 GO TO 999 01003 99 STOP 01004 END 01005 C 01006 C***********************************************************************01007 SUBROUTINE INPOLY(XPOLY,YPOLY,NP1,XI,YI,IN) 01008 C 01009 C THIS SUBROUTINE RETURNS A VALUE OF IN=1 IF POINT (XI,YI) IS 01010 C WITHIN THE POLYGON SPECIFIED BY XPOLY AND YPOLY. A VALUE OF 01011 C IN=0 IS RETURNED IF THE POINT IS NOT WITHIN THE POLYGON. 01012 C 01013 C NP1 = NUMBER OF POLYGON POINTS + 1 (THE POLYGON MUST BE CLOSED , 01014 C I.E. THE FIRST AND LAST POINT MUST BE THE SAME IN THE DATA 01015 C VECTORS XPOLY AND YPOLY) 01016 C 01017 C********************************************************************* 01018*11 C 01019*11 IMPLICIT REAL*8(A-H,O-Z) 01020 DIMENSION XPOLY(NP1),YPOLY(NP1) 01021 IN=0 01022 N6 = NP1 - 1 01023**3 DO 625 I=1,N6 01024**3 IF(YI.LE.YPOLY(I))GO TO 623 01025 622 IF(YI.GT.YPOLY(I+1)) GO TO 625 01026 624 IF((XI-XPOLY(I) -(YI-YPOLY(I))*(XPOLY(I+1)-XPOLY(I))/(YPOLY(I+1) 01027 & -YPOLY(I))).LT.0.0) IN=1-IN 01028 GO TO 625 01029 623 IF(YI.GT.YPOLY(I+1)) GO TO 624 01030 625 CONTINUE 01031 RETURN 01032 END 01033 C ******************************************************************* 01034*11 C * * 01035*11 C * * 01036*11 C * SUBROUTINE 'GRAV' IS DESIGNED TO PREDICT GRAVITY ANOMALY * 01037*11 C * BY THE METHOD OF LEAST-SQUARES COLLOCATION FOR DETAILS ABOUT * 01038*11 C * THE MATHEMATICAL MODEL SEE REFERENCES IN MAIN PROGRAM 'COLO- * 01039*11 C * CATE'. * 01040*11 C * * 01041*11 C * * 01042*11 C * INPUT DATA: * 01043*11 C * ----------- * 01044*11 C * * 01045*11 C * X, VDG - OBSERVATIONS AND THEIR VARIANCES * 01046*12 C * NDG - NUMBER OF OBSERVATIONS * 01047*12 C * DLT, DLN - POSITION COORDINATES OF EACH OBSERVATION * 01048*12 C * XLAT,XLOG - POSITION COORDINATES OF PREDICTION POINT * 01049*12 C * NBRR - NUMBER OF ELEMENTS IN UPPER TRIANGLE AND DIAGONAL * 01050*12 C * OF COVARIANCE MATRIX (CX) OF OBSERVATIONS * 01051*12 C * IPRINT - PARAMETER SPECIFYING WHETHER (OR NOT) TO PRINT * 01052*12 C * INTERMEDIATE RESULTS * 01053*12 C * * 01054*11 C * * 01055*11 C * OUTPUT: * 01056*11 C * ------- * 01057*11 C * * 01058*11 C * * 01059*11 C * DGP - VALUE OF PREDICTED ANOMALY * 01060*12 C * VAR - VARIANCE OF PREDICTED ANOMALY * 01061*12 C * ELEMENTS OF COVARIANCE MATRICES, VECTORS (IF REQUIRED) : * 01062*12 C * * 01063*12 C * CXX - VECTOR OF ALL ELEMENTS OF UPPER TRIANGLE AND * 01064*12 C * DIAGONAL OF COVARIANCE MATRIX (CX) OF MEASURE- * 01065*12 C * MENTS TAKING THE ROWS SUCCESSIVELY; SEE BELOW: * 01066*12 C * * 01067*12 C * 1 2 3 * 01068*12 C * 4 5 ---> ( 1 2 3 4 5 6 ) * 01069*12 C * 6 * 01070*12 C * CX,CX1,CX2 - COVARIANCE MATRICES OF MEASUREMENTS * 01071*12 C * CX1,CX2 WERE USED TO COPY CX. THEY COULD * 01072*12 C * BE DISCARDED FROM THE SUBROUTINE AND MAIN * 01073*12 C * PROGRAM * 01074*12 C * CSX,CSXA - SIGNAL COVARIANCE MATRIX AND ITS TRANSPOSE * 01075*12 C * ( CSX1 ) CSX1 WAS USED TO COPY CSX. * 01076*12 C * CSXB - SOLUTION MATRIX ( CSXB = CSX * CX(INVERSE) ) * 01077*12 C * GG - AUTO-COVARIANCE OF MEASUREMENTS. * 01078*12 C * * 01079*11 C ******************************************************************* 01080*11 C 01081*11 C 01082 C SUBROUTINE GRAV 01083 C =============== 01084 01085 SUBROUTINE GRAV(NDG,VDG,DLT,DLN,X,CXX,CSX,CSXA,CSXB, 01086 *XLAT,XLOG,DGP,VAR,NBRR,CX,CX2,GG,CX1,CSX1,IPRINT) 01087**3 C...PREPARED BY G. LACHAPELLE, TU GRAZ, OCTOBER 75. 01088 C...MODIFIED IN FEBRUARY 75. 01089 IMPLICIT REAL*8(A-H,O-Z) 01090 COMMON /CMCOV/CI(12),CR(51),SIGMA0(300),SIGMA(300), 01091 *KI(25),N1,LOCAL 01092 DIMENSION DLT(50),DLN(50),X(50),CXX(2550),CSX(50),VDG(50), 01093**3 *CSXA(1,50),CSX1(50,1),CSXB(1,50),CX(50,50),CX1(50,50),CX2(50,50), 01094**2 * NIK(50),NEK(50) 01095**2 C 01096 10 FORMAT(I4,16F9.2) 01097**4 15 FORMAT(//,4X,4F12.5) 01098**3 16 FORMAT(I4,16F9.5) 01099**4 20 FORMAT(//) 01100**3 30 FORMAT(//,4X,'CXX-MATRIX'/,4X,10('=')/) 01101**3 40 FORMAT(//,4X,'CXX-INVERSE'/,4X,11('=')/) 01102**3 50 FORMAT(//,4X,'DETERMINANT=',F7.4,'D',I4) 01103**4 60 FORMAT(//,4X,'CSX-TRANSPOSE'/,4X,13('=')/) 01104**4 70 FORMAT(//,4X,'CSXB=CSXT*CXX(INVERSE)'/,4X,22('-')/) 01105**4 CR(4)=DSIN(DLT(1)) 01106 CR(6)=DCOS(DLT(1)) 01107 CR(5)=CR(4) 01108 CR(7)=CR(6) 01109 CR(8)=0.0D0 01110 CR(9)=1.0D0 01111 CR(1)=CR(4)*CR(5)+CR(6)*CR(7)*CR(9) 01112 CALL COVCX(GG) 01113 NDG1=NDG-1 01114 IV=0 01115 DO 120 II=1,NDG1 01116 IV=IV+1 01117 CXX(IV)=GG+VDG(II) 01118 K=II+1 01119 DO 120 J=K,NDG 01120 IV=IV+1 01121 CR(4)=DSIN(DLT(II)) 01122 CR(6)=DCOS(DLT(II)) 01123 CR(5)=DSIN(DLT(J)) 01124 CR(7)=DCOS(DLT(J)) 01125 CR(8)=DSIN(DLN(II)-DLN(J)) 01126 CR(9)=DCOS(DLN(II)-DLN(J)) 01127 IF((DABS(DLT(II)-DLT(J)).LT.1.0D-5).AND.(DABS(DLN(II)-DLN(J)).LT. 01128 *1.0D-5)) GO TO 115 01129 CR(1)=CR(4)*CR(5)+CR(6)*CR(7)*CR(9) 01130 GO TO 120 01131 115 CR(1)=1.0D0 01132 120 CALL COVCX(CXX(IV)) 01133 IV=IV+1 01134 CXX(IV)=GG+VDG(NDG) 01135 C...VECTOR CSX %CSXA # CSX< 01136 CR(4)=DSIN(XLAT) 01137 CR(6)=DCOS(XLAT) 01138 DO 150 II=1,NDG 01139 CR(5)=DSIN(DLT(II)) 01140 CR(7)=DCOS(DLT(II)) 01141 CR(8)=DSIN(DLN(II)-XLOG) 01142 CR(9)=DCOS(DLN(II)-XLOG) 01143 IF((DABS(XLAT-DLT(II)).LT.1.0D-5).AND.(DABS(XLOG-DLN(II)) 01144 *.LT.1.0D-5)) GO TO 140 01145 CR(1)=CR(4)*CR(5)+CR(6)*CR(7)*CR(9) 01146 GO TO 145 01147 140 CR(1)=1.0D0 01148 145 CALL COVCX(CSX(II)) 01149 CSX1(II,1)=CSX(II) 01150 150 CSXA(1,II)=CSX(II) 01151 N=1 01152 DO 162 J=1,NDG 01153 CX(1,J)=CXX(N) 01154 N=N+1 01155 162 CONTINUE 01156 DO 157 I=2,NDG 01157 DO 157 K=I,NDG 01158 CX(I,K)=CXX(N) 01159 N=N+1 01160 157 CONTINUE 01161 DO 158 I=2,NDG 01162 I1=I-1 01163 DO 158 K=1,I1 01164 158 CX(I,K)=CX(K,I) 01165 IF(IPRINT.LT.2) GO TO 161 01166**3 WRITE(6,30) 01167**3 DO 159 I=1,NDG 01168**3 WRITE(6,10) I,(CX(I,J),J=1,NDG) 01169**3 159 CONTINUE 01170**3 WRITE(6,20) 01171**3 161 CONTINUE 01172**3 C... SCALING OF CXX(CX) AGAINST UNDERFLOW OR OVERFLOW DURING INVERSION 01173 C... DIVIDE EACH ELEMENT BY S=GG=588.13 GGAL**2. 01174 SS1=GG 01175**4 DO 163 I=1,NDG 01176 DO 163 J=1,NDG 01177 CX(I,J)=CX(I,J)/SS1 01178 163 CONTINUE 01179 C...SOLUTION 01180 CALL SPIN(CX,NDG,50,DET,IDEXP) 01181**3 C 01182**3 IF(IPRINT.EQ.2) WRITE(6,50) DET,IDEXP 01183**3 C... REFILL THE MATRIX CX AFTER INVERSION: 01184 DO 164 I=2,NDG 01185 I1=I-1 01186 DO 164 J=1,I1 01187 CX(I,J)=CX(J,I) 01188 164 CONTINUE 01189 C... RESCALE THE INVERTED CX MATRIX: 01190 DO 165 I=1,NDG 01191 DO 165 J=1,NDG 01192 CX(I,J)=CX(I,J)/SS1 01193 165 CONTINUE 01194 C... COMPUTE CSXB=CSXA*CXX(INV): 01195 DO 170 I=1,NDG 01196 CSXB(1,I)=0.0D0 01197 170 CONTINUE 01198 DO 171 I=1,NDG 01199 DO 171 J=1,NDG 01200 CSXB(1,I)=CSXB(1,I)+CSXA(1,J)*CX(J,I) 01201 171 CONTINUE 01202 168 FORMAT(10X,I5,F9.2) 01203 DGP=0.0 01204 SOM=0.0 01205 DO 200 II=1,NDG 01206 DGP=DGP+CSXB(1,II)*X(II) 01207 SOM=SOM+CSXB(1,II)*CSX1(II,1) 01208 200 CONTINUE 01209 VAR=GG-SOM 01210 IF(IPRINT.LT.2) GO TO 167 01211**3 WRITE(6,40) 01212**3 DO 166 I=1,NDG 01213**3 166 WRITE(6,16) I,(CX(I,J),J=1,NDG) 01214**4 WRITE(6,60) 01215**3 I=1 01216**3 WRITE(6,10) I,(CSXA(1,J),J=1,NDG) 01217**3 WRITE(6,70) 01218**3 WRITE(6,10) I,(CSXB(1,J),J=1,NDG) 01219**3 WRITE(6,15) DGP,GG,SOM,VAR 01220**3 WRITE(6,20) 01221**3 167 CONTINUE 01222**3 RETURN 01223 END 01224 C***********************************************************************01225 C* 01226 C* SUBROUTINE 'G C P L O T' :- ( BY : M.M. NASSAR, JANUARY 1975.) 01227 C* ------------------------ 01228 C* 01229 C* THE PURPOSE OF THIS SUBROUTINE IS TO PLOT01230 C* TWO PERPENDICULAR AXES: X (HORIZONTAL) AND Y (VERTICAL), INTERSECTING01231 C* AT THE CENTRE OF A SQUARE 6 BY 6 INCHES; AND THEN PLOT 'N' NUMBER OF 01232 C* POINTS DEFINED BY THEIR X & Y COORDINATES RELATIVE TO THE CENTRE OF 01233 C* THE SQUARE AND GIVEN IN THE USER'S UNITS. 01234 C* 01235 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01236 C* 01237 C* NOTE THAT THE TRANSFORMATION OF COORDINATES FROM THE USER UNITS TO TH01238 C* PLOTTING UNITS IS DONE AUTOMATICALLY VIA THE PLOTTING ROUTINES CALLED01239 C* BY THIS SUBROUTINE. FOR MORE DETAILS CONCERNING THE PLOTTING ROUTI01240 C* AVAILABLE AT THE U.N.B. COMPUTING CENTRE, SEE: ''COMPUTER PLOTTING'',01241 C* A USER GUIDE MANUAL BY: U.G. GUJAR, CL # 278, COMPUTING CENTRE, U.N.B01242 C* FREDERICTON, N.B., OCTOBER, 1972. 01243 C* 01244 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01245 C* 01246 C* USAGE :- CALL GCPLOT(X,Y,N,XMIN,YMIN,XMAX,YMAX) 01247 C* 01248 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01249 C* 01250 C* REQUIRED INPUT DATA :- (VIA THE CALLING ROUTINE) 01251 C* ------------------- 01252 C* 01253 C* X : VECTOR OF X COORDINATES OF THE N-POINTS, IN USER'S UNITS. 01254 C* Y : VECTOR OF Y COORDINATES OF THE N-POINTS, IN USER'S UNITS. 01255 C* N : TOTAL NUMBER OF POINTS TO BE PLOTTED. 01256 C* XMIN : MINIMUM X COORDINATE IN THE USER'S UNITS. 01257 C* YMIN : MINIMUM Y COORDINATE IN THE USER'S UNITS. 01258 C* XMAX : MAXIMUM X COORDINATE IN THE USER'S UNITS. 01259 C* YMAX : MAXIMUM Y COORDINATE IN THE USER'S UNITS. 01260 C* 01261 C* THE OUTPUT WILL BE THE PLOTTED SQUARE, THE TWO PERPENDICULAR AXES, AN01262 C* THE N-POINTS. FINALLY, THIS SUBROUTINE PRINTS-OUT THE TOTAL NUMBER '01263 C* OF PLOTTED POINTS NO LIMITATIONS ON THE NUMBER OF POINTS. 01264 C* 01265 C***********************************************************************01266 C 01267 SUBROUTINE GCPLOT(X1,Y1,N,XMIN,YMIN,XMAX,YMAX) 01268 IMPLICIT REAL*8(A-H,O-Z) 01269 DIMENSION X1(1500),Y1(1500) 01270**3 XLNGTH=6. 01271 YLNGTH=6. 01272 CALL AREA(XLNGTH,YLNGTH) 01273 CALL SETPLT(XMIN,YMIN,XMAX,YMAX) 01274 CALL PRNTCH('+') 01275 RECXL=2.0D 0*XMIN 01276 CALL RECT(XMIN,YMIN,RECXL) 01277 CALL AXES(XMIN,YMIN,XMAX,YMAX,0.,0.,-1.,-1.) 01278 CALL PRNTCH('.') 01279 DO 30 K=1,N 01280 CALL NOWPLT(0,X1(K),Y1(K)) 01281 30 CONTINUE 01282 CALL ENDPLT 01283 WRITE(6,19) N 01284 19 FORMAT(/,5X,'TOTAL NUMBER OF PLOTTED POINTS SURROUNDING THE CENTRA01285 1L-ORIGIN IS =',I6,2X,'POINTS') 01286 RETURN 01287 END 01288 C 01289 C 01290 C***********************************************************************01291 C 01292 SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) 01293**3 C 01294**3 C SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC 01295**3 C POSITIVE-DEFINITE MATRICES. THE MATRIX INVERTED IS THE UPPER 01296**3 C N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM 01297**3 C IN THE CALLING ROUTINE. 01298**3 C 01299**3 C INPUT: 01300**3 C Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE 01301**3 C MATRIX TO BE INVERTED. 01302**3 C 01303**3 C N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER) 01304**3 C OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST 01305**3 C NOT BE LARGER THAN MM) . 01306**3 C MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. 01307**3 C 01308**3 C 01309**3 C OUTPUT: 01310**3 C 01311**3 C Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF 01312**3 C THE INPUT UPPER LEFT N BY N PORTION. 01313**3 C 01314**3 C DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE 01315**3 C INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE 01316**3 C IDEXP BELOW. 01317**3 C 01318**3 C IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED 01319**3 C UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN 01320**3 C TWO PARTS CORRESPONDING TO 01321**3 C DETERMINANT = DET * 10 ** IDEXP . 01322**3 C THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE 01323**3 C COMPUTATION OF THE DETERMINANT. TO PRINT THE DETERM- 01324**3 C INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS; 01325**3 C (FOR EXAMPLE) 01326**3 C PRINT 10,DET,IDEXP 01327**3 C 10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4) 01328**3 C 01329**3 C R.R. STEEVES 01330**3 C SEPT., 1979 01331**3 C 01332**3 C 01333**3 C********************************************************************* 01334*11 C 01335*11 REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART 01336**3 DET=0.D0 01337**3 DO 4 J=1,N 01338**3 DO 4 I=1,J 01339**3 IF(I.EQ.1) GO TO 2 01340**3 M=I-1 01341**3 SUM=0.0D0 01342**3 DO 1 K=1,M 01343**3 1 SUM=SUM+Q(K,I)*Q(K,J) 01344**3 Q(I,J)=Q(I,J)-SUM 01345**3 2 IF(I.EQ.J) GO TO 3 01346**3 Q(I,J)=Q(I,J)/Q(I,I) 01347**3 GO TO 4 01348**3 3 CONTINUE 01349**3 DET=DET+DLOG10(Q(I,I)) 01350**3 Q(I,I)=DSQRT(Q(I,I)) 01351**3 4 CONTINUE 01352**3 IDEXP=DET 01353**3 RPART=DET-IDEXP 01354**3 APART=DABS(RPART) 01355**3 IF(APART.LT.1.D-20) DET=1.D0 01356**3 IF(APART.LT.1.D-20) GO TO 10 01357**3 DET=10**RPART 01358**3 10 CONTINUE 01359**3 DO 7 J=1,N 01360**3 DO 7 I=1,J 01361**3 IF(I.LT.J) GO TO 5 01362**3 Q(J,J)=1.0D0/Q(J,J) 01363**3 GO TO 7 01364**3 5 SUM=0.0D0 01365**3 M=J-1 01366**3 DO 6 K=I,M 01367**3 6 SUM=SUM-Q(I,K)*Q(K,J) 01368**3 Q(I,J)=SUM/Q(J,J) 01369**3 7 CONTINUE 01370**3 DO 9 J=1,N 01371**3 DO 9 I=1,J 01372**3 SUM=0.0D0 01373**3 DO 8 K=J,N 01374**3 8 SUM=SUM+Q(I,K)*Q(J,K) 01375**3 Q(I,J)=SUM 01376**3 IF(I.EQ.J) GO TO 9 01377**3 Q(J,I)=SUM 01378**3 9 CONTINUE 01379**3 RETURN 01380**3 END 01381**3 C***********************************************************************01382**9 C ******************************************************************* 01383*11 C * * 01384*11 C * SUBROUTINE 'CHISQ' IS DESIGNED TO CARRY OUT CHI-SQUARE * 01385*11 C * GOODNESS OF FIT TEST ON A DATA SERIES. * 01386*11 C * * 01387*11 C * * 01388*11 C * INPUT DATA: * 01389*11 C * ----------- * 01390*11 C * * 01391*11 C * DG4 - ARRAY OF STANDARDIZED OBSERVATIONS * 01392*12 C * NV - NUMBER OF OBSERVATIONS * 01393*12 C * NCHI - NUMBER OF CLASSES * 01394*12 C * SEG - DELIMITERS OF CLASSES * 01395*12 C * CHI - CHI-SQUARE VALUE FROM TABLES * 01396*12 C * NCHI - NUMBER OF KNOWN POPULATION PARAMETERS * 01397*12 C * * 01398*11 C * * 01399*11 C * OUTPUT: * 01400*11 C * ------- * 01401*11 C * * 01402*11 C * CHISQD - COMPUTED CHI-SQUARE VALUE * 01403*12 C * NDOF - DEGREE OF FREEDOM * 01404*12 C * ICHI - COUNTER FOR OBS. IN EACH CLASS * 01405*12 C * * 01406*11 C ******************************************************************* 01407*11 C 01408*11 SUBROUTINE CHISQ(DG4,NV,SEG,ICHI,NCHI,CHISQD,CHI,NCHI1) 01409**9 IMPLICIT REAL*8(A-H,O-Z) 01410**9 DIMENSION DG4(NV),SEG(NCHI) 01411**9 INTEGER ICHI(NCHI),NV,ICOT,NSEG,NSEG1 01412**9 NSEG=NCHI-1 01413**9 DO 10 K=1,NCHI 01414**9 10 ICHI(K)=0 01415**9 J=1 01416**9 ICOT=0 01417**9 DO 20 I=1,NV 01418**9 IF(DG4(I).GE.SEG(J)) ICHI(J)=ICHI(J)+1 01419**9 20 CONTINUE 01420**9 ICOT=ICOT+ICHI(J) 01421**9 J=J+1 01422**9 JM1=J-1 01423**9 25 IF(J.GT.NSEG) GO TO 35 01424**9 DO 30 I=1,NV 01425**9 IF((DG4(I).LT.SEG(JM1)).AND.(DG4(I).GE.SEG(J))) ICHI(J)=ICHI(J)+1 01426**9 30 CONTINUE 01427**9 ICOT=ICOT+ICHI(J) 01428**9 IF(ICOT.EQ.NV) GO TO 100 01429**9 JM1=J 01430**9 J=J+1 01431**9 GO TO 25 01432**9 35 CONTINUE 01433**9 ICHI(J)=NV-ICOT 01434**9 100 EXPD=NV/DFLOAT(NCHI) 01435**9 CHISQD=0.D0 01436**9 DO 102 K=1,NCHI 01437**9 IF(ICHI(K).EQ.0) GO TO 101 01438**9 CHISQD=CHISQD+((DFLOAT(ICHI(K))-EXPD)**2)/EXPD 01439**9 GO TO 102 01440**9 101 CHISQD=CHISQD+EXPD 01441**9 102 CONTINUE 01442**9 C 01443**9 IF(NCHI1.EQ.0) NCHI2=3 01444**9 IF(NCHI1.EQ.1) NCHI2=2 01445**9 IF(NCHI1.EQ.2) NCHI2=1 01446**9 NDOF=NCHI-NCHI2 01447**9 WRITE(6,128) NDOF,CHISQD,CHI 01448**9 128 FORMAT(///,50X,'CHI-SQUARE GOODNESS OF FIT TEST:'//,50X,'DEGREE OF01449**9 * FREEDOM=',I6/,42X,'COMPUTED CHI-SQUARE VALUE=',F12.5/,39X,'CHI-SQ01450**9 *UARE VALUE FROM TABLES=',F12.5) 01451**9 RETURN 01452**9 END 01453**9 C ******************************************************************* 01454*11 C * * 01455*11 C * SUBROUTINE 'AMNSDR' IS DESIGNED TO COMPUTE THE MEAN * 01456*11 C * AND STANDARD DEVIATION OF A SERIES OF OBSERVATIONS. * 01457*11 C * * 01458*11 C * * 01459*11 C * INPUT: * 01460*11 C * * 01461*11 C * RVAR - ARRAY OF OBSERVATIONS * 01462*12 C * NRVAR - NUMBER OF OBSERVATIONS * 01463*12 C * * 01464*11 C * * 01465*11 C * OUTPUT: * 01466*11 C * * 01467*11 C * * 01468*11 C * AMN - MEAN OF OBSERVATIONS * 01469*12 C * SDR - STANDARD DEVIATION OF THE MEAN * 01470*12 C * * 01471*11 C ******************************************************************* 01472*11 C 01473*11 SUBROUTINE AMNSDR(RVAR,NRVAR,AMN,SDR) 01474**9 C 01475**9 IMPLICIT REAL*8(A-H,O-Z) 01476**9 DIMENSION RVAR(NRVAR) 01477**9 RVSUM=0.D0 01478**9 RVSQSM=0.D0 01479**9 DO 10 I=1,NRVAR 01480**9 RVSUM=RVSUM+RVAR(I) 01481**9 RVSQSM=RVSQSM+RVAR(I)**2 01482**9 10 CONTINUE 01483**9 ANRV=DFLOAT(NRVAR) 01484**9 AMN=RVSUM/ANRV 01485**9 SDRSQ=(RVSQSM*ANRV-RVSUM**2)/(ANRV*(ANRV-1.)) 01486**9 SDR=DSQRT(SDRSQ) 01487**9 RETURN 01488**9 END 01489**9 C ******************************************************************* 01490*11 C * * 01491*11 C * SUBROUTINE 'HSTPLT' IS DESIGNED TO PLOT THE HISTOGRAM * 01492*11 C * OF A SERIES OF OBSERVATIONS. IT ASSEMLES THE OBSERVATIONS * 01493*11 C * INTO THEIR RESPECTIVE CLASSES AND CALLS SUBROUTINE 'PLOT' * 01494*11 C * TO PLOT THE HISTOGRAM. * 01495*11 C * * 01496*11 C * * 01497*11 C * INPUT DATA: DDD - ARRAY OF OBSERVATIONS * 01498*12 C * NDDD- NUMBER OF OBSERVATIONS. * 01499*12 C * * 01500*11 C * OUTPUT: - HISTOGRAM PLOT. * 01501*11 C * * 01502*11 C ******************************************************************* 01503*11 C 01504*11 SUBROUTINE HSTPLT(DDD,NDDD,NCLASS,HVEC,NB) 01505*12 IMPLICIT REAL*8(A-H,O-Z) 01506**9 INTEGER HVEC 01507**9 DIMENSION DDD(NDDD),HVEC(20),NB(11) 01508**9 NCLASS=20 01509**9 ACLASS=DFLOAT(NCLASS) 01510**9 DO 120 J=1,NCLASS 01511**9 120 HVEC(J)=0 01512**9 20 CALL FAMNMX(DDD,NDDD,GMN,GMX) 01513**9 GM1=DABS(GMN)+0.5 01514**9 GN=DABS(GMX)+0.5 01515**9 GR=GM1/GN 01516**9 IF(GR.GE.1) GO TO 559 01517**9 XX=-GN 01518**9 AXX=2.*GN/ACLASS 01519**9 GO TO 560 01520**9 559 CONTINUE 01521**9 XX=-GM1 01522**9 AXX=2.*GM1/ACLASS 01523**9 560 CONTINUE 01524**9 XX1=XX+AXX 01525**9 NB(1)=XX 01526**9 NC=NCLASS/2+1 01527**9 DO 561 I=2,NC 01528**9 IM1=I-1 01529**9 NB(I)=XX+IM1*AXX*2. 01530**9 561 CONTINUE 01531**9 K=1 01532**9 121 CONTINUE 01533**9 DO 122 J=1,NDDD 01534**9 IF((DDD(J).LT.XX).OR.(DDD(J).GT.XX1)) GO TO 122 01535**9 HVEC(K)=HVEC(K)+1 01536**9 122 CONTINUE 01537**9 K=K+1 01538**9 XX=XX1 01539**9 XX1=XX+AXX 01540**9 IF(K.GT.NCLASS) GO TO 123 01541**9 GO TO 121 01542**9 123 CONTINUE 01543**9 CALL PLOT(NCLASS,HVEC,NB) 01544**9 RETURN 01545**9 END 01546**9 C ******************************************************************* 01547*11 C * SUBROUTINE 'TONMN' IS DESIGNED TO PERFORM THE STUDENT'S T * 01548*11 C * TEST ON THE MEAN OF A SERIES OF OBSERVATIONS. * 01549*11 C * * 01550*11 C * * 01551*11 C * INPUT: RMN - MEAN VALUE * 01552*12 C * ST - STANDARD DEVIATION * 01553*12 C * NST - NUMBER OF OBSERVATIONS * 01554*12 C * NDEX - PARAMETER SPECIFYING WHETHER(OR NOT) * 01555*12 C * TO USE FACTOR FOR STANDARD NORMAL CURVE * 01556*12 C * * 01557*11 C * OUTPUT: - MESSAGE INDICATING WHETHER TEST * 01558*11 C * PASSED OR FAILED. * 01559*11 C * * 01560*11 C ******************************************************************* 01561*11 C 01562*11 SUBROUTINE TONMN(RMN,ST,NST,NDEX) 01563**9 IMPLICIT REAL*8(A-H,O-Z) 01564**9 TDFALF=1.98 01565**9 ANRM=1.96D0 01566**9 AST=DFLOAT(NST) 01567**9 AMULT=ST/DSQRT(AST) 01568**9 IF(NDEX.EQ.0) PROD=AMULT*ANRM 01569**9 PROD=AMULT*TDFALF 01570**9 EX1=RMN-PROD 01571**9 EX2=RMN+PROD 01572**9 IX=NST-1 01573**9 WRITE(6,124) IX,EX1,EX2 01574**9 124 FORMAT(///,50X,'STUDENTS-T TEST ON MEAN DIFFERENCE:'//,50X,'DEGREE01575**9 * OF FREEDOM=',I6/,51X,'CONFIDENCE LEVEL= 95%'//,50X,F10.5,' < 0 01576**9 *< ',F10.5) 01577**9 IF((EX1.LT.0.).AND.(EX2.GT.0.)) GO TO 125 01578**9 WRITE(6,126) 01579**9 GO TO 99 01580**9 125 WRITE(6,127) 01581**9 126 FORMAT(//,50X,'TEST FAILS'/,50X,10('*')//) 01582**9 127 FORMAT(//,50X,'==========> TEST PASSES <=========='//) 01583**9 99 RETURN 01584**9 END 01585**9 C ******************************************************************* 01586*11 C * * 01587*11 C * SUBROUTINE 'REJECT' IS DESIGNED TO PERFORM THE TEST OF * 01588*11 C * OUTLIERS ON A SERIES OF OBSERVATIONS. * 01589*11 C * INPUT: * 01590*11 C * * 01591*11 C * OBS - ARRAY OF OBSERVATIONS * 01592*12 C * RMS - STANDARD DEVIATION * 01593*12 C * IOBS - NUMBER OF OBSERVATIONS * 01594*12 C * NLOC1 - MAXIMUM NUMBER OF OBSERVATIONS DECLARED * 01595*12 C * IN THE MAIN PROGRAM * 01596*12 C * RJCT - FACTOR (MULTIPLIER) * 01597*12 C * * 01598*11 C * OUTPUT: * 01599*11 C * * 01600*11 C * ILOC - LOCATIONS OF THE REJECTED OBSERVATIONS IN THE * 01601*12 C * SERIES * 01602*12 C * NLOC - NUMBER OF REJECTED OBSERVATIONS * 01603*12 C * * 01604*11 C ******************************************************************* 01605*11 C 01606*11 SUBROUTINE REJECT(OBS,IOBS,RMS,RJCT,ILOC,NLOC,NLOC1) 01607**9 IMPLICIT REAL*8(A-H,O-Z) 01608**9 DIMENSION OBS(IOBS),ILOC(NLOC1) 01609**9 J=0 01610**9 RJCT=RJCT*RMS 01611**9 Y1=-RJCT 01612**9 Y2=RJCT 01613**9 DO 10 I=1,IOBS 01614**9 IF((OBS(I).GT.Y1).AND.(OBS(I).LT.Y2)) GO TO 10 01615**9 J=J+1 01616**9 ILOC(J)=I 01617**9 10 CONTINUE 01618**9 NLOC=J 01619**9 RETURN 01620**9 END 01621**9 SUBROUTINE SORT(VCLS,NOR,NRES) 01622**2 C***********************************************************************01623**2 C* 01624**2 C* SORT REARRANGES THE ELEMENTS OF VCLS IN ORDER OF INCREASING MAGNITUD01625**2 C* FOR USE IN SEPARATING THESE ELEMENTS INTO CLASS INTERVALS FOR THE CH01626**2 C* SQUARE GOODNESS OF FIT TEST AND HISTOGRAM PLOT. 01627**2 C* 01628**2 C* 01629**2 C* INPUT: 01630**2 C* NRES- NUMBER OF RESIDUALS IN VCLS 01631**2 C* 01632**2 C* 01633**2 C* WRITTEN BY: 01634**2 C* LAURIE PACH, JULY, 1978 01635**2 C* 01636**2 C***********************************************************************01637**2 IMPLICIT REAL*8(A-H,O-Z) 01638**2 DIMENSION VCLS(NOR) 01639**2 NOUT=NRES*(NRES+1)/2 01640**2 NINS=NRES-1 01641**2 DO 12 J=1,NOUT 01642**2 IFLG=0 01643**2 DO 11 I=1,NINS 01644**2 IF(VCLS(I).GT.VCLS(I+1))GOTO10 01645**2 GOTO11 01646**2 10 IFLG=1 01647**2 TEMP=VCLS(I) 01648**2 VCLS(I)=VCLS(I+1) 01649**2 VCLS(I+1)=TEMP 01650**2 11 CONTINUE 01651**2 IF(IFLG.EQ.0)GOTO13 01652**2 12 CONTINUE 01653**2 13 RETURN 01654**2 END 01655**2 SUBROUTINE ILPRNT(NVAL,MAX,PLOTV,RVEC,WINT,N,KK) 01656**2 C***********************************************************************01657**2 C* 01658**2 C* LPRNT CONTROLS THE LINE SPACING FOR THE NORMAL-HISTOGRAM PLOT PRINTI01659**2 C* 01660**2 C* WRITTEN BY: 01661**2 C* LAURIE PACH, JULY, 1978 01662**2 C* 01663**2 C***********************************************************************01664**2 INTEGER WINT,PLOTV,SV,RVEC 01665**2 DIMENSION PLOTV(110),SV(22),RVEC(WINT) 01666**2 DATA SV/' ',' ',' ',' ','R','E','L','A','T','I','V','E',' ','F', 01667**2 @ 'R','E','Q','U','E','N','C','Y'/ 01668**2 B=.3 01669**2 C=.2 01670**2 D=.1 01671**2 DO 2 I=1,100 01672**2 DO 3 N=1,WINT 01673**2 IF(RVEC(N).EQ.MAX)RETURN 01674**2 3 CONTINUE 01675**2 N=N-1 01676**2 IF(NVAL.EQ.MAX)RETURN 01677**2 PRINT101,(PLOTV(L),L=1,110) 01678**2 IF(MAX.EQ.25)PRINT113,B 01679**2 IF(MAX.EQ.17)PRINT113,C 01680**2 IF(MAX.EQ.9)PRINT113,D 01681**2 MAX=MAX-1 01682**2 IF(MAX.GT.32.OR.KK.EQ.23)GOTO2 01683**2 PRINT114,SV(KK) 01684**2 KK=KK+1 01685**2 2 CONTINUE 01686**2 101 FORMAT(' ',6X,110A1) 01687**2 113 FORMAT('+',3X,F3.1,'-') 01688**2 114 FORMAT('+',1X,A1) 01689**2 RETURN 01690**2 END 01691**2 SUBROUTINE PBLANK(PLOTV) 01692**2 C***********************************************************************01693**2 C* 01694**2 C* PBLANK CLEARS (SETS ELEMENTS TO BLANKS) VECTOR PLOTV, WHICH IS 01695**2 C* USED IN SUBROUTINE PLOT. 01696**2 C* 01697**2 C* 01698**2 C* WRITTEN BY: 01699**2 C* LAURIE PACH, JULY, 1978 01700**2 C* 01701**2 C***********************************************************************01702**2 INTEGER PLOTV 01703**2 DIMENSION PLOTV(110) 01704**2 DATA BLNK/' '/ 01705**2 DO 1 I=1,110 01706**2 PLOTV(I)=BLNK 01707**2 1 CONTINUE 01708**2 RETURN 01709**2 END 01710**2 SUBROUTINE PLOT(WINT,HVEC,NB) 01711**4 C***********************************************************************01712**2 C* 01713**2 C* PLOT PLOTS THE STANDARD NORMAL CURVE OVERLAYED WITH THE HISTOGRAM OF01714**2 C* STANDARD RESIDUALS. 01715**2 C* 01716**2 C* 01717**2 C* INPUT: 01718**2 C* WINT- NUMBER OF HISTOGRAM INTERVALS 01719**2 C* HVEC- VECTOR CONTAINING THE NUMBER OF RESIDUALS IN EACH HISTOG01720**2 C* INTERVAL 01721**2 C* 01722**2 C* 01723**2 C* WRITTEN BY: 01724**2 C* LAURIE PACH, JULY, 1978 01725**2 C* 01726**2 C***********************************************************************01727**2 INTEGER WINT,NVEC,RVEC,HVEC,STAR,HLINE,VLINE,PLOTV,PLOTL,PLOTH 01728**2 DIMENSION PLOTV(110),NVEC(53),PLOTL(110),PLOTH(110),RVEC(20), 01729**2 @ HVEC(20),NB(11) 01730**4 DATA STAR,HLINE,VLINE/'.','-','|'/ 01731**2 DATA NVEC/22*0,4*1,2,2,3,4,4,5,6,8,9,10,12,14,16,17,19,21,23,25, 01732**2 @27,28,29,31,31,32,32,32,31/ 01733**2 MAX=50 01734**2 A=.4 01735**2 K=2 01736**2 I=52 01737**2 NUMR=0 01738**2 KK=1 01739**2 NVAL=32 01740**2 NFLG=0 01741**2 IF(WINT.EQ.20)INT=5 01742**2 IF(WINT.EQ.10.OR.WINT.EQ.9)INT=10 01743**2 IF(WINT.EQ.2)INT=50 01744**2 IF(WINT.EQ.3)INT=30 01745**2 IF(WINT.EQ.5)INT=20 01746**2 IF(WINT.EQ.4)INT=25 01747**2 CALL PBLANK(PLOTV) 01748**2 CALL PBLANK(PLOTH) 01749**2 CALL PBLANK(PLOTL) 01750**2 DO 2 JJ=1,WINT 01751**2 NUMR=NUMR+HVEC(JJ) 01752**2 2 CONTINUE 01753**2 DO 33 JJ=1,WINT 01754**2 RVEC(JJ)=(80*HVEC(JJ)/NUMR)+.5 01755**2 33 CONTINUE 01756**2 PRINT103 01757**2 DO 29 JJ=1,50 01758**2 DO 28 N=1,WINT 01759**2 IF(RVEC(N).GE.MAX)GOTO19 01760**2 28 CONTINUE 01761**2 PRINT104 01762**2 IF(N.EQ.(WINT+1))N=WINT 01763**2 IF(MAX.EQ.32)GOTO21 01764**2 MAX=MAX-1 01765**2 29 CONTINUE 01766**2 19 RVEC(N)=0 01767**2 MM=(INT*(N-1))+1 01768**2 IF(WINT.EQ.9.OR.WINT.EQ.3)MM=(INT*(N-1))+6 01769**2 PLOTL(MM)=VLINE 01770**2 III=INT-1 01771**2 DO 12 JJ=1,III 01772**2 PLOTH(JJ+MM)=HLINE 01773**2 12 CONTINUE 01774**2 PLOTL(MM+INT)=VLINE 01775**2 DO 32 N=1,WINT 01776**2 IF(RVEC(N).EQ.MAX.OR.RVEC(N).GT.50)GOTO19 01777**2 32 CONTINUE 01778**2 IF(NFLG.EQ.1)GOTO25 01779**2 PRINT101,(PLOTH(L),L=1,110) 01780**2 CALL PBLANK(PLOTH) 01781**2 CALL ILPRNT(NVAL,MAX,PLOTL,RVEC,WINT,N,KK) 01782**2 IF(RVEC(N).EQ.MAX.AND.MAX.NE.32)GOTO19 01783**2 21 DO 31 L=2,100 01784**2 PLOTV(I)=STAR 01785**2 IF(NVEC(I).NE.NVEC(I-1))GOTO4 01786**2 K=K+1 01787**2 I=I-1 01788**2 31 CONTINUE 01789**2 4 I=I-1 01790**2 NFLG=1 01791**2 IF(RVEC(N).EQ.MAX)GOTO19 01792**2 25 DO 26 L=1,110 01793**2 IF(PLOTH(L).EQ.HLINE)GOTO36 01794**2 26 CONTINUE 01795**2 GOTO38 01796**2 36 DO 39 L=1,110 01797**2 IF(PLOTH(L).EQ.HLINE)PLOTV(L)=PLOTH(L) 01798**2 39 CONTINUE 01799**2 38 CALL PBLANK(PLOTH) 01800**2 11 PRINT102,(PLOTV(L),L=1,110) 01801**2 IF(MAX.EQ.32)PRINT113,A 01802**2 IF(I.EQ.1)GOTO20 01803**2 NVAL=NVEC(I) 01804**2 IF(NVAL.EQ.0)GOTO20 01805**2 CALL PBLANK(PLOTV) 01806**2 DO 37 L=1,110 01807**2 IF(PLOTL(L).EQ.VLINE)PLOTV(L)=PLOTL(L) 01808**2 37 CONTINUE 01809**2 CALL ILPRNT(NVAL,MAX,PLOTL,RVEC,WINT,N,KK) 01810**2 CALL PBLANK(PLOTV) 01811**2 DO 3 L=2,100 01812**2 PLOTV(I)=STAR 01813**2 PLOTV(I+K)=STAR 01814**2 K=K+2 01815**2 IF(I.EQ.1)GOTO11 01816**2 IF(NVEC(I).NE.NVEC(I-1))GOTO4 01817**2 I=I-1 01818**2 3 CONTINUE 01819**2 GOTO4 01820**2 20 PRINT111 01821**2 PRINT112,(NB(L),L=1,11) 01822**4 IF(WINT.EQ.20)PRINT107,(HVEC(L),L=1,20) 01823**2 IF(WINT.EQ.10)PRINT106,(HVEC(L),L=1,10) 01824**2 IF(WINT.EQ.4)PRINT108,(HVEC(L),L=1,4) 01825**2 IF(WINT.EQ.2)PRINT109,(HVEC(L),L=1,2) 01826**2 IF(WINT.EQ.9)PRINT116,(HVEC(L),L=1,9) 01827**2 IF(WINT.EQ.5)PRINT114,(HVEC(L),L=1,5) 01828**2 IF(WINT.EQ.3)PRINT115,(HVEC(L),L=1,3) 01829**2 101 FORMAT(' ',6X,110A1) 01830**2 102 FORMAT('+',6X,110A1) 01831**2 103 FORMAT('1') 01832**2 104 FORMAT(' ') 01833**2 106 FORMAT(' ',8X,10(I4,6X),/) 01834**2 107 FORMAT(' ',7X,20(I4,1X),/) 01835**2 108 FORMAT(' ',15X,I4,3(21X,I4),/) 01836**2 109 FORMAT(' ',34X,I4,36X,I4,/) 01837**2 111 FORMAT(' ',6X,20('|----'),'|') 01838**2 112 FORMAT(' ',5X,10(I3,7X),I3/) 01839**4 113 FORMAT('+',3X,F3.1,'-') 01840**2 114 FORMAT(' ',13X,I4,4(16X,I4),/) 01841**2 115 FORMAT(' ',23X,3(I4,26X),/) 01842**2 116 FORMAT(' ',14X,9(I4,6X),/) 01843**2 RETURN 01844**2 END 01845**2 C ******************************************************************* 01846*11 C * * 01847*11 C * SUBROUTINE 'FAMNMX' IS DESIGNED TO DETERMINE THE * 01848*11 C * MAXIMUM AND MINIMUM VALUES OF A DATA SERIES. * 01849*11 C * * 01850*11 C * * 01851*11 C * INPUT: * 01852*11 C * * 01853*11 C * AIN - ARRAY OF DATA * 01854*12 C * NAIN - NUMBER OF DATA IN ARRAY * 01855*12 C * * 01856*11 C * OUTPUT: * 01857*11 C * * 01858*11 C * BOUT,AOUT - MAXIMUM AND MINIMUN OF DATA SERIES * 01859*12 C * * 01860*11 C ******************************************************************* 01861*11 C 01862*11 SUBROUTINE FAMNMX(AIN,NAIN,AOUT,BOUT) 01863**4 IMPLICIT REAL*8(A-H,O-Z) 01864**4 DIMENSION AIN(NAIN) 01865**4 AOUT=1000. 01866**4 BOUT=-1000. 01867**4 DO 20 I=1,NAIN 01868**4 IF(AIN(I).LT.AOUT) AOUT=AIN(I) 01869**4 IF(AIN(I).GT.BOUT) BOUT=AIN(I) 01870**4 20 CONTINUE 01871**4 RETURN 01872**4 END 01873**4