PROGRAM BAYTAP C C BAYESIAN TIDAL ANALYSIS PROCEDURE BY GROUPING MODEL (BAYTAP-G) C C DESIGNED AND CODED BY C C MAKIO ISHIGURO(*) AND YOSHIAKI TAMURA(**) C C * THE INSTITUTE OF STATISTICAL MATHEMATICS, C 4-6-7 MINAMI-AZABU, MINATO-KU, TOKYO, 106, JAPAN. C TEL. 03-3446-1501 FAX 03-3446-1695 C E-MAIL ishiguro@ism.ac.jp C ** NATIONAL ASTRONOMICAL OBSERVATORY, MIZUSAWA, C 2-12 HOSHIGAOKA-CHO, MIZUSAWA-SHI, IWATE-KEN, 023, JAPAN. C TEL. 0197-22-7131 FAX 0197-22-3410 C E-MAIL tamura@miz.nao.ac.jp C C THIS PROGRAM REALIZES A DECOMPOSITION OF TIME SERIES Y C INTO THE FORM C Y(I) = TDC(I)+T(I)+S(I)+R(I)+I(I) C WHERE TDC=TIDAL COMPONENT T=SMOOTH TREND S=STEP C R=RESPONSE TO THE ASSOCIATED PHENOMENA C I=IRREGULAR C C THIS PROGRAM IS DESIGNED TO BE APPLICABLE TO DATA OF THE LENGTH C 90000, THAT IS HOURLY SAMPLED DATA OF MORE THAN TEN YEARS. C C THE PROCEDURE IS BASED ON A BAYESIAN MODEL AND ITS C PERFORMANCE IS CONTROLLED BY THE SELECTION OF THE PARAMETERS OF C THE PRIOR DISTRIBUTION. THE CONSTRUCTION OF THE BASIC MODEL IS C DISCUSSED IN THE FOLLOWING PAPERS: C C ABIC AND STATISTICS: C AKAIKE, H. (1979) LIKELIHOOD AND THE BAYES PROCEDURE. A PAPER C PRESENTED AT THE INTERNATIONAL MEETING ON BAYESIAN STATISTICS, C MAY 26 - JUNE 2, 1979, VALENCIA, SPAIN. TO BE PUBLISHED IN C TRABAJOS DE ESTADISTICA. C AKAIKE, H. AND ISHIGURO, M. (1979) TREND ESTIMATION WITH C MISSING OBSERVATIONS. A PAPER PRESENTED AT THE INTERNATIONAL C CONFERENCE IN STATISTICS IN TOKYO, NOVEMBER 28 - 30, 1979, C TOKYO. C ISHIGURO, M. AND AKAIKE, H. (1980) A BAYESIAN APPROACH TO C THE TRADING-DAY ADJUSTMENT OF MONTHLY DATA. A PAPER C PRESENTED AT THE INTERNATIONAL TIME SERIES MEETING, C AUGUST 14 - 15, 1980, HOUSTON, U.S.A. C ISHIGURO, M. (1981) COMPUTATIONALLY EFFICIENT IMPLIMENTATION C OF A BAYESIAN SEASONAL ADJUSTMENT PROCEDURE, RES. MEMO. 204, C INSTITUTE OF STATISTICAL MATHEMATICS, TOKYO C TIDAL ANALYSIS AND PROGRAM: C ISHIGURO, M., AKAIKE, H., OOE, M. AND NAKAI, S. (1981) C A BAYESIAN APPROACH TO THE ANALYSIS OF EARTH TIDES. C PROC. 9TH INT. SYMPOS. ON EARTH TIDES, PP283-292, C HELD IN AUG. 17-22, 1981, IN NEW YORK. C ISHIGURO, M. (1981) A BAYESIAN APPROACH TO THE ANALYSIS OF C THE DATA OF CRUSTAL MOVEMENTS, JOURNAL OF THE GEODETIC C SOCIETY OF JAPAN, 27, PP256-262. C ISHIGURO, M., SATO, T, TAMURA, Y. AND OOE, M. (1984) C TIDAL DATA ANALYSIS --AN INTRODUCTION TO BAYTAP--, C PROC. INSTITUTE OF STATISTICAL MATHEMATICS, 32, PP71-85. C ISHIGURO, M. AND TAMURA, Y. (1985) BAYTAP-G IN TIMSAC 84, C COMPUTER SCIENCE MONOGRAPHS, 22, PP56-117. C TAMURA, Y. (1985) BAYTAP-G USER'S MANUAL. C TAMURA, Y., SATO, T., OOE, M. AND ISHIGURO, M. (1991) C A PROCEDURE FOR TIDAL ANALYSIS WITH A BAYESIAN INFORMATION C CRITERION, GEOPHYSICAL JOURNAL INTERNATIONAL, 104, PP507-516. C RESPONSE METHOD: C MUNK, W. M. AND CARTWRIGHT, D. E. (1966) TIDAL SPECTROSCOPY AND C PREDICTION, PHIL. TRANS. ROY. SOC. LONDON, A, 259, PP533-581. C POTENTIAL: C CARTWRIGHT, D. E. AND TAYLER, R. J. (1971) NEW COMPUTATIONS OF C THE TIDE-GENERATING POTENTIAL, GEOPHYS. J. ROY. ASTRON. SOC., C 23, PP45-74. C CARTWRIGHT, D. E. AND EDDEN, A. C. (1973) CORRECTED TABLES OF C THE TIDAL HARMONICS, GEOPHYS. J. ROY. ASTRON. SOC., 33, C PP253-264. C TAMURA, Y. (1987) A HARMONIC DEVELOPMENT OF THE TIDE-GENERATING C POTENTIAL, MAREES TERRESTRES BULLETIN D'INFORMATION (BIM), 99, C PP6813-6855. C C C THE PRIOR DISTRIBUTION CONTROLS THE SMOOTHNESS C OF THE TREND BY ASSUMING A HOMOGENEOUS GAUSSIAN C DISTRIBUTION OF DIFFERENCE OF COMPONENT. THE C CHOICE OF THE VARIANCE OF THE HOMOGENEOUS GAUSSIAN DISTRIBUTION IS C REALIZED BY MAXIMIZING THE LOG LIKELIHOOD OF THE BAYESIAN MODEL. C FOR THE PURPOSE OF COMPARISON OF MODELS WITH DIFFERENT PARAMETERS C THE CRITERION ABIC IS DEFINED BY C C ABIC = - 2 * LOG( MAXIMUM LIKELIHOOD OF THE MODEL ) C + 2 * ( NUMBER OF ESTIMATED PARAMETERS ) C C SMALLER VALUE OF ABIC MEANS BETTER FIT OF THE MODEL. C THE VALUE OF THE VARIANCE OF THE PRIOR DISTRIBUTION WHICH C MINIMIZES ABIC WITHIN A PRESCRIBED FINITE SET OF POSSIBLE VALUES C IS USED TO PRODUCE THE FINAL OUTPUT FOR A PARTICULAR CHOICE OF C OTHER PARAMETERS. C FOR THE COMPARISON OF THE OVERALL PERFORMANCES OF VARIOUS C MODELS THE AVERAGED ABIC (AVABIC) IS USEFUL. C C THIS PROGRAM REQUIRES THE FOLLOWING PARAMETERS. WITHOUT FURTHER C SPECIFICATION BY THE PROGRAM USER, THEY ARE RESPECTIVELY SET EQUAL C TO THE VALUES GIVEN IN THE PARENTHESES AT THE ENDS OF THEIR C DESCRIPTIONS. C C PARAMETERS: C SPAN: NUMBER OF DATA TO BE PROCESSED AT ONE TIME (0) C SPAN IS AUTOMATICALLY SELECTED WHEN SPAN=0 C SHIFT: NUMBER OF DATA TO BE SHIFTED C SHIFT IS AUTOMATICALLY SELECTED WHEN SPAN=0 C TO DEFINE THE NEW SPAN OF DATA (720) C ORDER: ORDER OF DIFFERENCING OF TREND (2) C IPOTEN: SELECTION OF REFERENCE POTENTIAL (2) C POTENTIAL BY CTE WHEN IPOTEN=1 C POTENTIAL BY TAMURA WHEN IPOTEN=2 C P4FLAG: FLAG FOR 4TH ORDER POTNTIAL (1) C P4FLAG=1 IS VALID WHEN IPOTEN=2 C IAUG: NUMBER OF ASSOCIATED DATASETS (0) C LAGP: MAXIMUM LAG NUMBER (3) C LAGINT: LAG INTERVAL (1) C DMIN: LOWER LIMIT OF HYPERPARAMETER D (0.25) C 4*DMIN IS USED FOR A INITIAL VALUE OF D SEARCH C WEIGHT: CONTROLS THE ADAPTIVITY OF THE INPULSE RESPONSE C MORE DATA-ADAPTIVE WITH SMALLER WEIGHT (1.0) C MAXJMP: MAXIMUM COUNT OF THE STEP WITHIN A SPAN (5) C LPOUT: LP OUTPUT LEVEL (0), LPOUT=1 IS VALID WHEN SPAN>0 C FILOUT: FALG FOR FILE OUTPUT THROUGH IO UNIT NO 16 C FILOUT=1 IS VALID WHEN SPAN>0 C PREPRO: FLAG FOR INTERPOLATION (0) C PREPRO=1, 2 OR 3 IS VALID WHEN LPOUT=1 AND SPAN>0 C C C TO MODIFY THE PARAMETERS FOLLOW THE FOLLOWING EXAMPLE: C IF THE USER WANTS TO SET IAUG=1,WEIGHT=0.5 AND SO ON, C PLACE THE CARDS WHICH CONTAIN THE FOLLOWING C THREE TYPES OF STATEMENTS ON TOP OF THE INPUT DATA( PUNCH ONE C SPACE AT THE FIRST COLUMN OF EACH CARD) : C &PARAM (ON THE FIRST CARD) C IAUG=1, WEIGHT=0.5, AND SO ON, (ON THE SECOND AND LATER CARDS) C &END (ON THE LAST CARD) C NOTE THE COMMA AT THE END OF EACH PARAMETER SPECIFICATION. C C INPUT DATA: C C THE FOLLOWING DATA SET SHOULD BE FED THROUGH THE CARD READER C IN THE FORMATS SHOWN IN THE PARENTHSES. C C C KIND C 1 GRAVITY TIDE C 2 TILT TIDE N-S C 3 TILT TIDE E-W (OR KIND=2, AZ=90) C 4 NS STRAIN C 5 EW STRAIN (OR KIND=4, AZ=90) C 6 NE STRAIN (OR KIND=4, AZ=45) C 7 FOR STRAIN DATA (NONE DIMENSIONAL:POTENTIAL/G/R) C 8 OCEAN TIDE C ORDERING OF PARAMETERS C COS PART OF GROUP 1 C SIN PART OF GROUP 1 C COS PART OF GROUP 2 C SIN PART OF GROUP 2 C AND SO ON. C C NUMERICAL OUTPUTS: C TREND: ESTIMATED TREND C C GRAPHICAL OUTPUTS: C ORIGINAL DATA C TREND + STEP C IRREGULAR C C C STRUCTURE C C BAYTAP :->TJUL C :->PRMSET C :->NAMELS C :->PRMCHK C :->SETGRP C :->TBREAD C :->SETMWV C :->DATAR2 :->FNORMG C : :->TJUL C : :->DATAR :->UPDATE C : C :->GEOGR1 C :->SPHAS1 C :->GEOGR2 C :->SPHAS2 C :->SETCOS C :->TIDE01 C :->TIDE02 C :->RTITLE C :->CALEND C :->SUBSEA :->SETX :->YREADO,YREADT,YREADA C : : :->HUSHLD C : : :->SETD C : : :->INIT C : : :->EXQRCN C : : :->SOLVE C : : C : :->SOLVE C : :->DECODE :->YREADO,YREADT,YREADA C : :->COPY C : C :->FRQRSP :->DSP3 C :->COPY C :->OTLCHK C :->SAUTCO :->SMEADL :->SUMF C : :->CROSCO C : :->CORNOM C : C :->SICP2 C :->SNRASP :->FOUGER C : :->SUBVCP C : :->DSP3 C : C :->GRAPH3 :->UPDATE C : :->KANDO :->TJUL C : C :->KANDO :->TJUL C :->BLOCKD C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C PARAMETER STATEMENTS IN SUBROUTINES SHOULD BE MODIFIED C CONSISTENTLY WITH MAIN PROGRAM. C FOR 1MB MEMORY CC PARAMETER ( MAXH2 = 30000, MAXDC = 2500, MAXOUT = 1000 ) CC PARAMETER ( MAXGRP = 15, MAXGR2 = MAXGRP*2 ) CC PARAMETER ( LBUFF = 500 ) CC PARAMETER ( MAXASC = 1 ) C FOR 2MB MEMORY CC PARAMETER ( MAXH2 = 90000, MAXDC = 10000, MAXOUT = 3000 ) CC PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) CC PARAMETER ( LBUFF = 1250 ) CC PARAMETER ( MAXASC = 3 ) C FOR 3MB MEMORY ( 3 MONTH DATA WITH DRIFT ESTIMATION(SHIFT>0) ) CC PARAMETER ( MAXH2 = 200000, MAXDC = 20000, MAXOUT = 3000 ) CC PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) CC PARAMETER ( LBUFF = 1250 ) CC PARAMETER ( MAXASC = 3 ) C FOR 9MB MEMORY (1 YEAR DATA WITH DRIFT ESTIMATION(SHIFT>0) ) PARAMETER ( MAXH2 = 600000, MAXDC = 30000, MAXOUT = 9000 ) PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( LBUFF = 4500 ) PARAMETER ( MAXASC = 3 ) C FOR 15MB MEMORY (3 YEAR DATA WITH DRIFT ESTIMATION(SHIFT>0) ) C (WITH FILE I/O) CC PARAMETER ( MAXH2 = 1300000, MAXDC = 80000, MAXOUT = 26000 ) CC PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) CC PARAMETER ( LBUFF = 4500 ) CC PARAMETER ( MAXASC = 3 ) C FOR 34MB MEMORY (4 YEAR DATA WITH DRIFT ESTIMATION(SHIFT>0) ) C (WITHOUT FILE I/O) CC PARAMETER ( MAXH2 = 2600000, MAXDC = 110000, MAXOUT = 36000 ) CC PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) CC PARAMETER ( LBUFF = 18000 ) CC PARAMETER ( MAXASC = 3 ) C C PARAMETER ( MAXUNK = 100 ) PARAMETER ( MAXOP2 = MAXOUT + MAXUNK ) PARAMETER ( MAXCLB = 100 ) PARAMETER ( MXLSTA = 50, MXLSTB = 100, MXOTL = 20 ) PARAMETER ( MAXFLE = MAXASC + 2 ) C C INTEGER ORDER, SPAN, SHIFT, FILOUT, SPECTW, PREPRO, S1WAVE INTEGER PAGLEN, P4FLAG, RSPERR INTEGER ZERORM DOUBLE PRECISION IRREG DIMENSION YO(LBUFF,2), YA(LBUFF,2,MAXASC) REAL YT(LBUFF,2,MAXGR2) C DIMENSION H2(MAXH2), DC(MAXDC) DIMENSION TREND(MAXOUT), DATA(MAXOUT), IRREG(MAXOUT), * TIDAL(MAXOP2), RESPNS(MAXOP2), ERR(MAXOP2) DIMENSION AMP(MAXGRP), PHASE(MAXGRP) DIMENSION CENTER(MAXGRP), FCR(MAXGRP), FCI(MAXGRP), CNTR(MAXGRP) DIMENSION FTRN(10), FTRNBK(10), RESPC(MAXUNK) DIMENSION LSMRY(3,MXLSTB),LISTIR(MXOTL),LISTTR(MXOTL),RAT(MAXFLE) DIMENSION JLIST(MXLSTA), QJUMP(MXLSTA), QQJUMP(MXLSTB) DIMENSION AUTCO(101), TAUTCO(101) DIMENSION IG1R(22), IG2R(22), FCRR(11), FCRI(11) DIMENSION COMG(10), COMG0(10), TCOS(10) C CHARACTER MARK(4) CHARACTER*8 CMODEL, CVERSI, SYMBOL CHARACTER*39 OBUF1 CHARACTER*40 CSTNAM, CINSTR(10), CUNIT, CASSOC(MAXASC) CHARACTER*60 CCALCU*60 CHARACTER*72 CHRBUF CHARACTER*80 CTITLE(10) C COMMON H2, DC, TREND, DATA, IRREG, TIDAL, RESPNS, ERR COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /IDATA/ ORDER,LAGDBL,LAGP,NJUMP,SPAN,LAG, * NCOEF,ITH,ITHP1,LAGINT,IFCR COMMON /RDATA/ WEIGHT,DMIN,MAXITR COMMON /LPCNT/ LPOUT COMMON /JSTUT/ TIMSYS COMMON /STNGRP/ IG1(36),IG2(36),NGROUP,IMWAVE(36) COMMON /STNGRC/ SYMBOL(36) COMMON /VERNAM/ CMODEL, CVERSI, CCALCU COMMON /VALIDD/ NUMUSE, MAXJMP COMMON /INFORM/ IFORMT, NETALT, ETALT(MAXCLB), TJET(MAXCLB) COMMON /BUFOTA/ YO, YT, YA, INCORE COMMON /PRSTEP/ STEP00 COMMON /PHASE0/ PHCORR(4,4) COMMON /SAVEWS/ WS000(919), AMPC0(919), DIFC0(919), CON00(919), * AMPS0(919), DIFS0(919) COMMON /DTDMUT/ TDMUT(130) C DATA IG1R / 1, 66, 92,114,117,124,127,138,140,147,169, * 1,144,202,250,257,267,271,289,293,306,346 / DATA IG2R / 65, 91,113,116,123,126,137,139,146,168,205, * 143,201,249,256,266,270,288,292,305,345,450 / C CCC EXTERNAL TJUL C C C ALLOCATION OF CONTROL DATA AND LP OUTPUT OPEN( 5, FILE='input05.dat', STATUS='OLD' ) OPEN( 6, FILE='output06.dat', STATUS='NEW' ) C C ALLOCATION OF OPTIONAL OUTPUT FILE C OPEN( 16, FILE='output16.dat', STATUS='NEW' ) OPEN( 16, FILE='output16.dat' ) C C ALLOCATION OF TIDAL DATA AND ASSOCIATED DATASETS C I/O UNIT NUMBERS ARE DEFINED BY CONTROL DATA C OPEN( 12, FILE='input12.dat', STATUS='OLD' ) C OPEN( 13, FILE='input13.dat', STATUS='OLD' ) C OPEN( 14, FILE='input14.dat', STATUS='OLD' ) C OPEN( 15, FILE='input15.dat', STATUS='OLD' ) OPEN( 12, FILE='input12.dat' ) OPEN( 13, FILE='input13.dat' ) OPEN( 14, FILE='input14.dat' ) OPEN( 15, FILE='input15.dat' ) C C ALLOCATION OF WORKING AREA OPEN( 30, STATUS='SCRATCH', FORM='UNFORMATTED' ) OPEN( 31, STATUS='SCRATCH', FORM='UNFORMATTED' ) OPEN( 32, STATUS='SCRATCH', FORM='UNFORMATTED' ) CC OPEN( 33, STATUS='SCRATCH', FORM='UNFORMATTED' ) C C C PRESET CONSTANTS C DDMIN = 1.0D0 NFILE = 29 CALL BLOCKD C C C HEADDER PRINTOUT C NF0 = NFILE + 1 NF1 = NFILE + MAXFLE WRITE(6, 100) CMODEL WRITE(6, 200) CVERSI WRITE(6, 300) MAXH2, MAXDC, MAXOUT, MAXGRP, MAXASC, LBUFF, MAXUNK WRITE(6, 400) NF0,NF1 100 FORMAT(1H /' ***** ',A8,' *****') 200 FORMAT(1H ,' VERSION ',A8) 300 FORMAT(1H ,' WITH CONSTANTS'/ * 1H ,' MAXH2 =',I9/ * 1H ,' MAXDC =',I9/ * 1H ,' MAXOUT =',I9/ * 1H ,' MAXGRP =',I9/ * 1H ,' MAXASC =',I9/ * 1H ,' LBUFF =',I9/ * 1H ,' MAXUNK =',I9) 400 FORMAT(1H ,' TEMPORAL FILE =',I3,' TO',I3) C C C PRESET PARAMETERS C CALL PRMSET( ITH, ORDER, SPAN, SHIFT, RLIM, LPOUT, FILOUT, * SPECTW, PREPRO, WEIGHT, IAUG, LAGP, LAGINT, MAXJMP, * KIND, AZ, DMIN, TIMSYS, IFCR, QVAL, MAXITR, * ICETGR, S1WAVE, PAGLEN, IPOTEN, P4FLAG, LOVENM, * RSPERR, SIGMA1, SIGMA2, ZERORM ) C C PARAMETER MODIFICATION AND DATA INPUT C CALL NAMELS( ITH, ORDER, SPAN, SHIFT, RLIM, LPOUT, FILOUT, * SPECTW, PREPRO, WEIGHT, IAUG, LAGP, LAGINT, MAXJMP, * KIND, AZ, DMIN, TIMSYS, IFCR, QVAL, MAXITR, * ICETGR, S1WAVE, PAGLEN, IPOTEN, P4FLAG, LOVENM, * RSPERR, SIGMA1, SIGMA2, ZERORM ) C C PARAMETER CHECK C CALL PRMCHK( ITH, ORDER, SPAN, SHIFT, RLIM, LPOUT, FILOUT, * SPECTW, PREPRO, WEIGHT, IAUG, LAGP, LAGINT, MAXJMP, * KIND, AZ, DMIN, TIMSYS, IFCR, QVAL, MAXITR, * ICETGR, S1WAVE, PAGLEN, IPOTEN, P4FLAG, LOVENM, * RSPERR, SIGMA1, SIGMA2, ZERORM ) C IF( IPOTEN.EQ.1 ) THEN IWND = 205 IWNS = 360 IWNT = 377 ELSE IF( IPOTEN.EQ.2 ) THEN IWND = 450 IWNS = 827 IWNT = 909 ELSE WRITE(6,405) 405 FORMAT(1H /' ERROR ILLEGAL IPOTEN') STOP END IF C C READ(5,410) JA,JB,JC,HR,NDATA,DELT 410 FORMAT(3I5,F5.0,I5,F10.0) INCORE = 0 IF( NDATA.LE.2*LBUFF ) INCORE = 1 JA = MOD( JA, 1900 ) NGROUP=0 NCOSGR=0 IG1(1)=1 DO 430 I = 1, MIN0( MAXGRP+1, 36 ) C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) DO 412 L = 1, 72 412 CHRBUF(L:L) = ' ' READ(5,415,ERR=414) (CHRBUF(L:L),L=1,72) 414 IF( CHRBUF(1:3) .EQ. 'COS' ) GO TO 436 READ(CHRBUF,420) IG2I, SYMBOL(I) 415 FORMAT(72A1) 420 FORMAT(I5,2X,A8) IF( IG2I.LE.6 .OR. IG2I.GT.IWNT ) GO TO 440 IG2(I) = IG2I IG1(I+1) = IG2I+1 NGROUP = NGROUP+1 IF( IG2I.GT.IWNT .OR. IG2I.LE.IG1(I) ) THEN WRITE(6,425) I, IG2I 425 FORMAT(1H /' INVALID GROUPING I =',I3, * 6X,'WAVE NO.=', I8) STOP END IF 430 CONTINUE WRITE(6,435) 435 FORMAT(1H /' ERROR TOO MANY GROUPS') STOP C 436 IF( CHRBUF(1:5) .EQ. 'COS ' ) THEN NCOSGR = NCOSGR + 1 IF( NCOSGR.GT.10 ) WRITE(6,435) IF( NCOSGR.GT.10 ) STOP READ(CHRBUF,437) COMG(NCOSGR), IYY,IMM,IDD,HHH, COMG0(NCOSGR) 437 FORMAT(5X,D15.8,3I5,F5.1,D15.8) TCOS(NCOSGR) = TJUL(IYY,IMM,IDD,HHH) END IF C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) DO 452 L = 1, 72 452 CHRBUF(L:L) = ' ' READ(5,415,ERR=454) (CHRBUF(L:L),L=1,72) 454 IF( CHRBUF(1:3) .EQ. 'COS' ) GO TO 436 READ(CHRBUF,438) IG2I 438 FORMAT(I5) C 440 CONTINUE IF( NGROUP.EQ.0 ) THEN OBSLEN = DBLE(SPAN)*DELT/24.D0 IF( SPAN .EQ.0 ) OBSLEN = DBLE( MIN0(744,NDATA) )*DELT/24.D0 IF( SHIFT.EQ.0 ) OBSLEN = DBLE( NDATA )*DELT/24.D0 CALL SETGRP( OBSLEN, IG2I, ICETGR, IPOTEN ) END IF CALL TBREAD( IPOTEN ) CALL SETMWV( IPOTEN ) C IF( ITH.EQ.0 ) THEN NGROUP = 0 IFCR = 0 END IF C C NGROUP : NUMBER OF TIDAL GROUPS C NCOSGR : NUMBER OF COSINE WAVES C NNNGRP : TOTAL NNNGRP = NGROUP + NCOSGR IF( NNNGRP .GT. MAXGRP ) WRITE(6,435) IF( NNNGRP .GT. MAXGRP ) STOP C C CALL DATAR2(RAT,H2,0,NDATA,RLIM,JA,JB,JC,HR,DELT,NFILE, * SLAT,SLONG,SH,GE, KIND,AZ, ZERORM) C C WRITE(6,700) SPAN, ITH, SHIFT, IPOTEN, LPOUT, P4FLAG, * FILOUT, S1WAVE, MAXJMP, ICETGR, ORDER, IFCR, * WEIGHT, QVAL, IAUG, RLIM, LAGP, MAXITR, * LAGINT, PREPRO, DMIN, PAGLEN, TIMSYS, SPECTW WRITE(6,710) KIND, LOVENM IF( KIND.EQ.3 ) AZ = 90.D0 IF( KIND.EQ.5 ) AZ = 90.D0 IF( KIND.EQ.6 ) AZ = 45.D0 IF( KIND.GE.2 .AND. KIND.LE.6 ) WRITE(6,720) AZ WRITE(6,725) SIGMA1, SIGMA2, ZERORM WRITE(6,730) JA, JB, JC, HR, NDATA, DELT 700 FORMAT(1H // *' ********** PARAMETER SETTING **********'/ *' SPAN =',I9, 10X, 'ITH =',I9/ *' SHIFT =',I9, 10X, 'IPOTEN =',I9/ *' LPOUT =',I9, 10X, 'P4FLAG =',I9/ *' FILOUT =',I9, 10X, 'S1WAVE =',I9/ *' MAXJMP =',I9, 10X, 'ICETGR =',I9/ *' ORDER =',I9, 10X, 'IFCR =',I9/ *' WEIGHT =',F9.3, 10X, 'QVAL =',F9.1/ *' IAUG =',I9, 10X, 'RLIM =',F9.1/ *' LAGP =',I9, 10X, 'MAXITR =',I9/ *' LAGINT =',I9, 10X, 'PREPRO =',I9/ *' DMIN =',F9.3, 10X, 'PAGLEN =',I9/ *' TIMSYS =',F9.1, 10X, 'SPECTW =',I9 ) 710 FORMAT(1H /' THEORETICAL VALUE'/1H ,' KIND =',I7, * 12X, 'LOVENM =',I7) 720 FORMAT(1H ,' AZ =',F7.1) 725 FORMAT(1H /' DATA REJECTION OPTIONS'/ * 1H ,' SIGMA1 =',1PD9.2, 8X, 'SIGMA2 =',1PD9.2 / * 1H ,' ZERORM =',I5 ) 730 FORMAT(1H //' DATE OF THE FIRST DATA'/ * ' YEAR MONTH DAY HOUR DATA COUNT', * ' SAMPLING' / I5, 2I10, F11.2, I10, 4X, F8.2 ) C C DEG = DELT I61 = 60*SPECTW + 1 IF(I61.LT. 61) I61 = 61 IF(I61.GT.361) I61 = 361 C IF( S1WAVE.EQ.1 ) THEN IF( IPOTEN.EQ.1 ) THEN AVC(125) = 15.00000000D0 AMC(124) = 0.D0 CCC AMC(125) = -0.00416D0 AMC(126) = 0.D0 IDN(125,6) = 0 ELSE IF( IPOTEN.EQ.2 ) THEN AVC(269) = 15.00000000D0 AMC(267) = 0.D0 AMC(268) = 0.D0 CCC AMC(269) = -0.004145D0 AMC(270) = 0.D0 IDN(269,6) = 0 END IF END IF C C A C RESONANCE = 1 + -------------------------------- C OMG - ( OMG0 + I * OMG0 / 2Q ) C CCC QVAL = 1150.D0 OMEGA0 = 15.0737D0 Q1 = 0.5D0/QVAL IF( IFCR.EQ.1 ) THEN DO 850 N = 1, 11 IF( IPOTEN.EQ.1 ) THEN NG1R = IG1R(N) NG2R = IG2R(N) ELSE IF( IPOTEN.EQ.2 ) THEN NG1R = IG1R(N+11) NG2R = IG2R(N+11) END IF SUM1 = 0.D0 SUM = 0.D0 DO 840 I = NG1R, NG2R TEM = DABS( AMC(I) ) SUM = SUM + TEM SUM1 = SUM1 + AVC(I)*TEM 840 CONTINUE CENTF = SUM1/SUM DEN = CENTF/OMEGA0 DNU = DEN - 1.D0 DEN = DNU**2 + Q1**2 FCRR(N) = DNU/DEN FCRI(N) = Q1/DEN 850 CONTINUE END IF C C DO 880 N = 1, NGROUP NG1 = IG1(N) NG2 = IG2(N) SUM1 = 0.D0 SUM = 0.D0 DO 870 I = NG1, NG2 TEM = DABS( AMC(I) ) SUM = SUM + TEM SUM1 = SUM1 + AVC(I)*TEM 870 CONTINUE CENTER(N) = SUM1/SUM DEN = CENTER(N)/OMEGA0 DNU = DEN - 1.D0 DEN = DNU**2 + Q1**2 FCR(N) = DNU/DEN FCI(N) = Q1/DEN IF( CENTER(N).GT.20.D0 ) THEN FCR(N) = 0.D0 FCI(N) = 0.D0 END IF 880 CONTINUE DO 882 N = 1, NCOSGR CENTER(N+NGROUP) = COMG(N) FCR(N+NGROUP) = 0.D0 FCI(N+NGROUP) = 0.D0 882 CONTINUE C DO 888 I = 1, NGROUP-1 IF( CENTER(I+1) - CENTER(I) .GT. 8.D0 ) THEN CNTR(I) = 0.05D0 ELSE CNTR(I) = 1.D0 END IF 888 CONTINUE DO 890 I = MAX0(NGROUP,1), NNNGRP 890 CNTR(I) = 0.001D0 CC890 CNTR(I) = 0.D0 C C IF( SHIFT .EQ. 0 ) SPAN = MIN0( 744, NDATA ) IF( NDATA .LT. SPAN ) SPAN = NDATA IF( SPAN .EQ. 0 ) THEN IF( NDATA .LE. 960 ) THEN SPAN = NDATA SHIFT = NDATA ELSE SHIFT = 720 IX = (NDATA - 24)/SHIFT JX = NDATA - 24 - IX*SHIFT IF( JX .LE. 12 ) THEN SPAN = SHIFT + JX + 24 ELSE KX = JX/IX /6*6 SHIFT = SHIFT + KX SPAN = SHIFT + JX - KX*IX + 24 END IF END IF END IF WRITE(6,1391) SPAN, SHIFT 1391 FORMAT(1H /' *** SPAN =',I5,' SHIFT =',I5,' ***') C C C WORK AREA CHECK C CCC IF( ITH.NE.0 ) ITH = ( NNNOUP + IFCR )*2 ITH = ( NNNGRP + IFCR )*2 C IORDP1 = ORDER + 1 LAGDBL = 1 LAG = 0 ITHP1 = ITH + 1 ITHAUG = ITH + IAUG IAUGP2 = IAUG + 2 LAGPP1 = LAGP + 1 NCOEF = LAGDBL*ITH + LAGPP1*IAUG N2MAX = NCOEF + IORDP1 + MAXJMP LIMIT = N2MAX + SPAN ITEM1 = IORDP1*LIMIT ITEM2 = N2MAX*LIMIT WRITE(6,1500) WRITE(6,1510) SPAN, MAXOUT WRITE(6,1520) NDATA, MAXH2-11 WRITE(6,1530) IAUGP2, MAXASC+2 WRITE(6,1540) N2MAX, MAXUNK WRITE(6,1550) LIMIT, MAXOP2 WRITE(6,1560) ITEM1, MAXDC WRITE(6,1570) ITEM2, MAXH2 WRITE(6,1580) ILL = 0 IF(SPAN .GT. MAXOUT) ILL = 1 IF(NDATA .GT. MAXH2-11) ILL = 1 IF(IAUGP2 .GT. MAXASC+2) ILL = 1 IF(N2MAX .GT. MAXUNK) ILL = 1 IF(LIMIT .GT. MAXOP2) ILL = 1 IF(ITEM1 .GT. MAXDC) ILL = 1 IF(ITEM2 .GT. MAXH2) ILL = 1 C EMERGENCY STOP IF( ILL.NE.0 ) THEN WRITE(6,1600) STOP END IF 1500 FORMAT(1H /' ******** WORK AREA CHECK ********'/ * 21X,'COUNT',21X,'LIMIT') 1510 FORMAT(1H ,'SPAN ',I9,8X,'MAXOUT ',I9) 1520 FORMAT(1H ,'DATA LENGTH ',I9,8X,'MAXH2-11 ',I9) 1530 FORMAT(1H ,'IAUG + 2 ',I9,8X,'MAXASC+2 ',I9) 1540 FORMAT(1H ,'N2MAX ',I9,8X,'MAXUNK ',I9) 1550 FORMAT(1H ,'LIMIT ',I9,8X,'MAXOP2 ',I9) 1560 FORMAT(1H ,'(ORDER+1)*LIMIT ',I9,8X,'MAXDC ',I9) 1570 FORMAT(1H ,'N2MAX*LIMIT ',I9,8X,'MAXH2 ',I9) 1580 FORMAT(1H ,5X,'N2MAX = (NO. OF GROUP)*2 + (LAGP+1)*IAUG + ORDER +' * ,' MAXJMP + 1'/ * 1H ,5X,'LIMIT = N2MAX + SPAN') 1600 FORMAT(1H /' ***** INSUFFICENT WORKING AREA *****') C C C DATA INPUT (CONTINUED) AND THEORETICAL VALUE COMPUTATION C 2400 CONTINUE JY=JA JM=JB JD=JC HRT=HR NCALC=NDATA + LAG*2 TEM=HR-LAG*DELT NFILE = NFILE + 1 C IF( ITH.EQ.0 ) GO TO 2490 C TJ0 = TJUL(JA,JB,JC,HR) IF( IPOTEN.EQ.1 ) THEN CALL GEOGR1( KIND, SLAT, SH, GE, AZ, LOVENM ) CALL SPHAS1( TJ0, HR, SLONG ) ELSE IF( IPOTEN.EQ.2 ) THEN TC0 = TJ0 + DBLE(NDATA)*DELT/24.D0*0.5D0 CALL GEOGR2( KIND, SLAT, SH, GE, AZ, LOVENM ) CALL SPHAS2( TJ0, HR, SLONG, TC0 ) END IF IF( NCOSGR.GT.0 ) * CALL SETCOS( NCOSGR, COMG, COMG0, TCOS, TJ0, IPOTEN ) NBLC = (NCALC - 1)/LBUFF + 1 C DO 2480 NB = 1, NBLC JEND = LBUFF IF( NB.EQ.NBLC ) JEND = MOD( NCALC-1, LBUFF ) + 1 JSTART = LBUFF*(NB - 1) DO 2470 N = 1, NNNGRP NG1 = IG1(N) NG2 = IG2(N) DO 2460 J = 1, JEND JCOPY = J T = DELT*(J-1+JSTART) IF( IPOTEN.EQ.1 ) THEN CALL TIDE01( NG1, NG2, T, CT, ST, DELT, JCOPY ) ELSE IF( IPOTEN.EQ.2 ) THEN CALL TIDE02( NG1, NG2, T, CT, ST, DELT, JCOPY, * P4FLAG, TC0 ) END IF H2(J) = CT H2(J+LBUFF) = ST 2460 CONTINUE IF( INCORE.EQ.1 ) THEN DO 2465 J = 1, JEND YT(J+JSTART,1,2*N-1) = H2(J) YT(J+JSTART,1,2*N ) = H2(J+LBUFF) 2465 CONTINUE ELSE WRITE(31) (SNGL(H2(J)),J=1,2*LBUFF) END IF 2470 CONTINUE IF(IFCR .EQ. 1) THEN DO 2410 I = 1, 2*LBUFF 2410 H2(I) = 0.D0 DO 2474 N = 1, 11 IF( IPOTEN.EQ.1 ) THEN NG1R = IG1R(N) NG2R = IG2R(N) ELSE IF( IPOTEN.EQ.2 ) THEN NG1R = IG1R(N+11) NG2R = IG2R(N+11) END IF DO 2472 J = 1, JEND JCOPY = J T = DELT*(J-1+JSTART) IF( IPOTEN.EQ.1 ) THEN CALL TIDE01( NG1R, NG2R, T, CT, ST, DELT, JCOPY ) ELSE IF( IPOTEN.EQ.2 ) THEN CALL TIDE02( NG1R, NG2R, T, CT, ST, DELT, JCOPY, * P4FLAG, TC0 ) END IF H2(J) = H2(J) + (FCRR(N)*CT+FCRI(N)*ST) H2(J+LBUFF) = H2(J+LBUFF) + (FCRR(N)*ST-FCRI(N)*CT) 2472 CONTINUE 2474 CONTINUE IF( INCORE.EQ.1 ) THEN DO 2475 J = 1, JEND YT(J+JSTART,1,2*NNNGRP+1) = H2(J) YT(J+JSTART,1,2*NNNGRP+2) = H2(J+LBUFF) 2475 CONTINUE ELSE WRITE(31) (SNGL(H2(J)),J=1,2*LBUFF) END IF END IF 2480 CONTINUE C IF( INCORE.NE.1 ) REWIND 31 2490 CONTINUE C C DO 2500 K = 1, IAUG CALL DATAR2(RAT,H2,1,NDATA,RLIM,JA,JB,JC,HR,DELT,NFILE, * SLAT2,SL2,SH2,GE2, KIND,AZ, ZERORM) IF( NDATA.LT.SPAN ) THEN WRITE(6,2495) 2495 FORMAT(1H /' LENGTH OF ASSOCIATED DATASET IS TOO SHORT') STOP END IF 2500 CONTINUE C C C C READING TITLS, STATION NAME, INSTRUMENT NAME, UNIT AND C ASSOCIATED DATA NAME C CALL RTITLE( NUMTIT, NUMINS, CTITLE, CSTNAM, * CINSTR, CUNIT, CASSOC, IAUG ) C C C C**** INITIALIZATION **** C RLIM=1.D25 DO 2700 I=1,ORDER 2700 FTRN(I) = 0.D0 DO 2750 I=1,101 2750 AUTCO(I)=0.D0 NSUM=0 N=NDATA IF(SHIFT .NE. 0) N=SPAN AVSD=0.D0 AVABIC=0.D0 COUNT=0.D0 AJUMP=0.D0 LP=0 C C C C ****************** C * MAIN BLOCK * C ****************** C DO 5400 IIII=1,100 JLISTP=0 INTEST=0 IF(IIII.EQ.1) INTEST=ORDER NJUMP=INTEST N2=NCOEF+1+NJUMP ISTART=1+(IIII-1)*SHIFT IEND=ISTART-1+N C EXIT OF THE MAIN BLOCK IF(ISTART-1+N .GT. NDATA) GO TO 5500 CALL CALEND(JA,JB,JC,HR,ISTART,DELT,JY,JM,JD,HRT) CALL CALEND(JY,JM,JD,HRT,N,DELT,JY1,JM1,JD1,HRT1) WRITE(6, 2801) WRITE(6, 3372) WRITE(6, 2900) WRITE(6, 3000) WRITE(6, 3100) JY,JM,JD,HRT,JY1,JM1,JD1,HRT1 WRITE(6, 3000) WRITE(6, 2900) WRITE(6, 3372) 2800 FORMAT(1H /1H ) 2801 FORMAT(1H1) 2900 FORMAT(1H ,'**************************************************', *'**************') 3000 FORMAT(1H ,'* ', *' *') 3100 FORMAT(1H ,'* START =',3I4,F5.1,3X, * '; END =',3I4,F5.1,' *') C C C C**** MAIN PROCEDURE **** C DO 90110 IRTY0 = 1, ORDER 90110 FTRNBK(IRTY0) = FTRN(IRTY0) C CALL SUBSEA(DDD,RESPC,ITHAUG,ABIC,TREND,RESPNS,IRREG,TIDAL, * FTRN,N,ISTART,RLIM,DC,IORDP1,H2,N2,N2MAX, * DATA,INTEST,LIMIT,SHIFT,ERR,SQE,CNTR,SD, * MAXDC,MAXH2) C C C C /* TEST VERSION, INCORE=1, SHIFT>0 ONLY ! */ C IF( ( SIGMA1.GT.1.D25 .AND. SIGMA2.GT.1.D25 ) .OR. * INCORE.NE.1 .OR. SHIFT.EQ.0 ) GO TO 90300 C DO 90200 IRTY9 = 1, 2 C RTYERR = SIGMA1 IF( IRTY9.EQ.2 ) THEN IF( SIGMA2.GE.SIGMA1 ) GO TO 90200 RTYERR = SIGMA2 END IF IF( RTYERR.GT.1.D25 ) GO TO 90200 C IRETRY = 0 DO 90010 IRTY1 = 1, N RTYD1 = DABS(IRREG(IRTY1)) IF( RTYD1.GT.RTYERR .AND. RTYD1.LT.1.D25 ) THEN IRETRY = IRETRY + 1 IF( IRETRY.EQ.1 ) WRITE(6,90020) RTYERR IRTY2 = IRTY1 + ISTART - 1 WRITE(6,90030) IRTY2, YO(IRTY2,1), IRREG(IRTY1) YO(IRTY2,1) = 1.D30 END IF 90010 CONTINUE 90020 FORMAT(1H ///' MAXIMUM ERROR IS LARGER THAN THE LIMIT', F9.3/) 90030 FORMAT(1H ,' POSITION, YOBS, IRREG =', I6, F13.3, F11.3) C C SUBSEA AGAIN C IF( IRETRY.NE.0 ) THEN WRITE(6,90040) IRETRY 90040 FORMAT(1H /' TOTAL', I5, ' DATA', 9X,'ABOVE DATA MIGHT', * ' BE REMOVED FROM INPUT DATASET' /// ) C INTEST = 0 IF( IIII.EQ.1 ) INTEST = ORDER NJUMP = INTEST N2 = NCOEF + 1 + NJUMP DO 90120 IRTY0 = 1, ORDER 90120 FTRN(IRTY0) = FTRNBK(IRTY0) C CALL SUBSEA(DDD,RESPC,ITHAUG,ABIC,TREND,RESPNS,IRREG,TIDAL, * FTRN,N,ISTART,RLIM,DC,IORDP1,H2,N2,N2MAX, * DATA,INTEST,LIMIT,SHIFT,ERR,SQE,CNTR,SD, * MAXDC,MAXH2) END IF C 90200 CONTINUE 90300 CONTINUE C C C AVABIC=AVABIC+ABIC AVSD=AVSD+SD*N COUNT=COUNT+N C SKIP TO THE OUTPUT ROUTINE IF( ITH.EQ.0 ) GO TO 3600 C DO 3300 J = 1, NCOEF 3300 ERR(J) = DSQRT( ERR(J)*SQE ) C WRITE(6,2800) WRITE(6,3400) WRITE(6,3402) ( RESPC(J), J = 1, 2*NNNGRP ) WRITE(6,3404) ( ERR(J), J = 1, 2*NNNGRP ) IF( IFCR.EQ.1 ) THEN WRITE(6,3410) WRITE(6,3412) RESPC(2*NNNGRP+1), RESPC(2*NNNGRP+2) WRITE(6,3414) ERR (2*NNNGRP+1), ERR (2*NNNGRP+2) END IF DO 3305 I = 1, IAUG WRITE(6,3420) I JE = MIN0( LAGP, 5 ) WRITE(6,3421) ( J, J = 0, JE ) FNM = RAT(I+1) JS = (NNNGRP+IFCR)*2 + I*LAGPP1 JE = JS - LAGP WRITE(6,3422) ( RESPC(J)*FNM, J = JS, JE, -1 ) WRITE(6,3424) ( ERR(J)*FNM, J = JS, JE, -1 ) 3305 CONTINUE WRITE(6,3430) SQE C RAD = 3.141592653589793D0/180.D0 C WRITE(6,2801) WRITE(6,2800) C DO 3310 II = 1, NUMTIT 3310 WRITE(6,3435) CTITLE(II) WRITE(6,3440) CSTNAM WRITE(6,3445) SLONG WRITE(6,3450) SLAT WRITE(6,3455) SH WRITE(6,3460) CINSTR(1) DO 3320 II = 2, NUMINS 3320 WRITE(6,3465) CINSTR(II) WRITE(6,3470) CUNIT WRITE(6,3475) JY,JM,JD,HRT, JY1,JM1,JD1,HRT1 WRITE(6,3476) NUMUSE, N WRITE(6,3480) CMODEL, CVERSI IF( NGROUP.GE.1 ) THEN IF( IPOTEN.EQ.1 ) THEN WRITE(6,3485) ELSE IF( IPOTEN.EQ.2 ) THEN IF( P4FLAG.EQ.0 ) THEN WRITE(6,3486) ELSE WRITE(6,3487) END IF END IF ELSE WRITE(6,3488) END IF WRITE(6,3490) CCALCU C WRITE(6,2800) WRITE(6,3372) C WRITE(6,3380) NPP1 = NNNGRP*2 + 1 NPP2 = NPP1 + 1 NP = 0 DO 3360 I = 1, NGROUP NP = NP + 1 SUMR = RESPC(NP) IF(IFCR .EQ. 1) SUMR=SUMR+RESPC(NPP1)*FCR(I)-RESPC(NPP2)*FCI(I) ERRR = ERR(NP)**2 IF(IFCR .EQ. 1) ERRR = ERRR + (ERR(NPP1)*FCR(I))**2 * + (ERR(NPP2)*FCI(I))**2 NP = NP + 1 SUMI = -RESPC(NP) IF(IFCR .EQ. 1) SUMI = SUMI - RESPC(NPP2)*FCR(I) * - RESPC(NPP1)*FCI(I) ERRI = ERR(NP)**2 IF(IFCR .EQ. 1) ERRI = ERRI + (ERR(NPP2)*FCR(I))**2 * + (ERR(NPP1)*FCI(I))**2 AMP(I) = DSQRT(SUMI**2+SUMR**2) EAMP = DSQRT( ERRR*SUMR**2 + ERRI*SUMI**2 )/AMP(I) PHASE(I) = DATAN2(SUMI,SUMR) / RAD EPH = DSQRT( ERRI*SUMR**2 + ERRR*SUMI**2 )/AMP(I)**2 / RAD C AMPLMW = 0.D0 IF(IG2(I).LE.IWND) AMPLMW = DABS( GC(2,1)*AMP(I)*AMC(IMWAVE(I)) ) IF(IG1(I).GT.IWND .AND. IG2(I).LE.IWNS) * AMPLMW = DABS( GC(2,2)*AMP(I)*AMC(IMWAVE(I)) ) IF(IG1(I).GT.IWNS) AMPLMW = DABS( GC(3,3)*AMP(I)*AMC(IMWAVE(I)) ) AMPLER = AMPLMW*EAMP/AMP(I) C WRITE(6,3390) I, IG1(I), IG2(I), SYMBOL(I), * AMP(I), EAMP, PHASE(I), EPH, AMPLMW, AMPLER 3360 IF( IG2(I).EQ.IWND .OR. IG2(I).EQ.IWNS ) WRITE(6,3372) C IF( NCOSGR.GT.0 ) WRITE(6,2800) DO 3370 I = NGROUP+1, NNNGRP NP = NP + 1 SUMR = RESPC(NP) ERRR = ERR(NP)**2 NP = NP + 1 SUMI = -RESPC(NP) ERRI = ERR(NP)**2 AMP(I) = DSQRT( SUMI**2 + SUMR**2 ) EAMP = DSQRT( ERRR*SUMR**2 + ERRI*SUMI**2 )/AMP(I) PHASE(I) = DATAN2( SUMI, SUMR ) / RAD EPH = DSQRT( ERRI*SUMR**2 + ERRR*SUMI**2 )/AMP(I)**2 / RAD WRITE(6,3392) I, CENTER(I), AMP(I), EAMP, PHASE(I), EPH 3370 CONTINUE C WRITE(6,2800) WRITE(6,3510) ABIC WRITE(6,3515) SD WRITE(6,3520) WEIGHT WRITE(6,3525) DDD IF(IFORMT.NE.1) WRITE(6,3530) DELT IF(IFORMT.EQ.1) WRITE(6,3531) DELT IF(SHIFT.GT.0) WRITE(6,3535) SPAN WRITE(6,3540) SHIFT IF(IAUG.EQ.0) THEN WRITE(6,3545) ELSE WRITE(6,3372) DO 3371 II=1,IAUG 3371 WRITE(6,3550) CASSOC(II) WRITE(6,3555) LAGP, LAGINT END IF C 3372 FORMAT(1H ) 3380 FORMAT(1H ,' GROUP',11X,'SYMBOL',7X, * 'FACTOR (RMSE) PHASE (RMSE) AMPLITUDE', * ' (RMSE)' /1H ,49X,'(LOCAL, LAG:NEGATIVE)' ) 3390 FORMAT(1H ,I5,' (',I3,' -',I3,' : ',A8,')', * F9.4,' (',F7.4,' )',F10.2,' (',F6.2,' )', * F10.3,' (',F6.3,' )' ) 3392 FORMAT(1H ,I5,' (',F16.8,' )', * F9.4,' (',F7.4,' )',F10.2,' (',F6.2,' )' ) C 3400 FORMAT(1H ,'TIDAL FACTORS (COS AND SIN PARTS)') 3402 FORMAT(1H ,' FACT. = ',1P2D13.5,2X,2D13.5,2X,2D13.5) 3404 FORMAT(1H ,' ERROR = ',1P2D13.5,2X,2D13.5,2X,2D13.5) 3410 FORMAT(1H /' FCR FACTORS') 3412 FORMAT(1H ,' FACT. =',1P2D14.5) 3414 FORMAT(1H ,' ERROR =',1P2D14.5) 3420 FORMAT(1H /' RESPONSE WEIGHTS',8X,'ASSOCIATED DATASET NO. =', I3 ) 3421 FORMAT(1H ,' LAG =',8X,6(I1,14X) ) 3422 FORMAT(1H ,' RESPC =',1P6D15.5) 3424 FORMAT(1H ,' ERROR =',1P6D15.5) 3430 FORMAT(1H //' SQE =', 1PD15.5) C 3435 FORMAT(1H ,A80) 3440 FORMAT(1H ,44X,'STATION ',4X,A40) 3445 FORMAT(1H ,44X,' LON.= ',4X,F9.3,' E') 3450 FORMAT(1H ,44X,' LAT.= ',4X,F9.3,' N') 3455 FORMAT(1H ,44X,' H = ',4X,F8.2,' M') 3460 FORMAT(1H ,44X,'INSTRUMENT',4X,A40) 3465 FORMAT(1H ,44X,' ',4X,A40) 3470 FORMAT(1H ,44X,'AMPLITUDE ',4X,A40) 3475 FORMAT(1H /45X,'PERIOD 19',I2,2I3,F5.1,' - 19',I2,2I3,F5.1) 3476 FORMAT(1H /45X,'AVAILABLE / BLOCK LEN.',I9,' /',I6) 3480 FORMAT(1H /' ANALYZED BY ',A8,' VERSION ',A8) 3485 FORMAT(1H ,' POTENTIAL : CARTWRIGHT, TAYLER AND EDDEN') 3486 FORMAT(1H ,' POTENTIAL : TAMURA, Y., 1987, BIM, 99. ', * '(WITHOUT P4)' ) 3487 FORMAT(1H ,' POTENTIAL : TAMURA, Y., 1987, BIM, 99. (WITH P4)') 3488 FORMAT(1H ,' POTENTIAL : NONE') 3490 FORMAT(1H ,' CALCULATION : ',A60) C 3510 FORMAT(1H /' MINIMUM ABIC =',F11.2) 3515 FORMAT(1H ,' STANDARD DEVIATION =',F13.4) 3520 FORMAT(1H /' HYPERPARAMETER WEIGHT =',F12.3) 3525 FORMAT(1H ,' " D =',1PD12.4) 3530 FORMAT(1H /' SAMPLING INTERVAL ORIGINAL = 0.5 HOUR' / * 1H ,' RESAMPLING =',F5.1,' HOUR') 3531 FORMAT(1H /' SAMPLING INTERVAL (ORIGINAL) =',F5.1,' HOUR') 3535 FORMAT(1H ,' SPAN =',I5) 3540 FORMAT(1H ,' SHIFT =',I5) 3545 FORMAT(1H /' ASSOCIATED DATA : NONE') 3550 FORMAT(1H ,' ASSOCIATED DATA : ',A40) 3555 FORMAT(1H ,' MAXIMUM LAG =',I3 / * 1H ,' LAG INTERVAL =',I3 ) C C C FOURIER TRANSFORM OF RESPONSE WEIGHT C C IF( IAUG .GE. 1 ) THEN KRESPC = (NNNGRP+IFCR)*2 + 1 CALL FRQRSP( RAT, RESPC(KRESPC), ERR(KRESPC), LAGINT, LAGPP1, * IAUG, DEG, PAGLEN, RSPERR ) END IF C C C IN CASE C THEORETICAL TIDE : NONE ( ITH.EQ.0 ) C ASSOCIATED DATA : EXIST ( IAUG.GE.1 ) C 3600 IF( ITH.EQ.0 .AND. IAUG.GE.1 ) THEN DO 3610 J = 1, NCOEF 3610 ERR(J) = DSQRT( ERR(J)*SQE ) C WRITE(6,2800) DO 3615 I = 1, IAUG WRITE(6,3420) I JE = MIN0( LAGP, 5 ) WRITE(6,3421) ( J, J = 0, JE ) FNM = RAT(I+1) JS = IAUG*LAGPP1 JE = JS - LAGP WRITE(6,3422) ( RESPC(J)*FNM, J = JS, JE, -1 ) WRITE(6,3424) ( ERR(J)*FNM, J = JS, JE, -1 ) 3615 CONTINUE WRITE(6,3430) SQE C CALL FRQRSP( RAT, RESPC(1), ERR(1), LAGINT, LAGPP1, * IAUG, DEG, PAGLEN, RSPERR ) END IF C C C NUMERICAL OUTPUT (IF FILOUT.EQ.1) C C EXIT OF THE MAIN BLOCK FOR SHIFT = 0 IF(SHIFT .EQ. 0) GO TO 5700 DO 3700 I=1,N 3700 TREND(I)=TREND(I)-AJUMP NOUT=N IF(IEND .LT. NDATA .AND. SHIFT .NE. 0) NOUT=SHIFT C IF(FILOUT.EQ.1) THEN IF( IIII.EQ.1 ) WRITE(16,4000) IF( IIII.EQ.1 ) WRITE(16,4100) DO 3800 I=1,NOUT ISER=ISTART-1+I WRITE(16,4200) ISER, IRREG(I), TREND(I), * RESPNS(I), TIDAL(I), DATA(I) 3800 CONTINUE ELSE IF(FILOUT.EQ.2) THEN WRITE(16) NOUT, ISTART, ISTART+NOUT-1 WRITE(16) ( SNGL( IRREG(I) ), I = 1, NOUT ) WRITE(16) ( SNGL( TREND(I) ), I = 1, NOUT ) WRITE(16) ( SNGL( RESPNS(I) ), I = 1, NOUT ) WRITE(16) ( SNGL( TIDAL(I) ), I = 1, NOUT ) WRITE(16) ( SNGL( DATA(I) ), I = 1, NOUT ) END IF C 3900 FORMAT(1H ,'STEP ',1PD13.4,' AT', I5) 4000 FORMAT(1H ,' *(ORIGINAL DATA)=(TIDAL COMP.)+(TREND)+(RESPONSE)' * ,'+(IRREGULAR)+(STEP)' / * 1H ,' FORMAT IS 1H ,I5,5D13.5, 1.D30 AND -1.D30', * ' ARE SPECIAL MARKS') 4100 FORMAT(1H ,4X,'I',4X,'IRREGULAR',4X, * 'TREND',8X,'RESPONSE',4X, * 'TIDAL COMP. ORIGINAL') 4200 FORMAT(1H ,I5,1P5D13.5) C C C C STEP CORRECTION AND GRAPHICAL OUTPUT (IF LPOUT .GE. 1) C ITEM=1+SHIFT-ORDER CALL COPY(FTRN,TREND,ORDER,ITEM) CALL OTLCHK(TREND,NOUT,1,LISTTR,LTR0,H2,MXOTL) I101=101 IF(I101 .GT. NOUT-1) I101=NOUT CALL SAUTCO(H2,TAUTCO,NOUT-1,I101,LPOUT) IF( SPECTW.GE.1 ) THEN WRITE(6,6410) CALL SICP2(TAUTCO,51,NOUT,H2,L,SGME2,OAIC) CALL SNRASP(H2,DC,H2,SGME2,L,0,I61,DEG,PAGLEN) END IF DO 4350 I = 1, I101 AUTCO(I) = ( AUTCO(I)*DBLE(NSUM) + TAUTCO(I)*DBLE(NOUT) ) * / DBLE( NSUM + NOUT ) 4350 CONTINUE NSUM=NSUM+NOUT CALL OTLCHK(IRREG,NOUT,0,LISTIR,LIR0,H2,MXOTL) ITEM=ISTART-1+SHIFT JTEM=NCOEF + NJUMP + 1 JPCHK = 0 DO 4600 I = 1, NOUT IF( DATA(I).LE.-RLIM .AND. JPCHK.NE.1 ) THEN JTEM = JTEM - 1 AJUMP0 = RESPC(JTEM) JPCHK = 1 END IF IF( JPCHK.NE.0 .AND. DABS(DATA(I)).LE.RLIM ) THEN AJUMP = AJUMP + AJUMP0 JPCHK = 0 IJUMP = ISTART + I - 1 IF( LPOUT.GE.2 ) WRITE(6,3900) AJUMP0, IJUMP JLISTP = JLISTP + 1 IF( JLISTP.GT.MXLSTA ) THEN WRITE(6,4410) 4410 FORMAT(1H /' TOO MANY STEPS STOPPED AT MAIN') STOP END IF JLIST(JLISTP) = I - 1 QJUMP(JLISTP) = AJUMP0 END IF TREND(I) = TREND(I) + AJUMP 4600 CONTINUE C IF( JPCHK.EQ.1 .AND. SPAN.GT.NOUT ) THEN IJUMP = ISTART + NOUT - 1 IF( LPOUT.GE.2 ) WRITE(6,3900) AJUMP0, IJUMP JLISTP = JLISTP + 1 IF( JLISTP.GT.MXLSTA ) THEN WRITE(6,4410) STOP END IF JLIST(JLISTP) = NOUT QJUMP(JLISTP) = AJUMP0 END IF C IF( LPOUT.GE.1 ) * CALL GRAPH3( DATA, TREND, IRREG, RESPNS, TIDAL, NOUT, * PREPRO, RLIM, JY, JM, JD, HRT, DELT ) C C C C ASSESSMENT OF THE RESULT OF THE PRESENT SPAN C C EMERGENCY EXIT (BUFFER OVER FLOW) IF(LP .GE. MXLSTB) GO TO 5200 IF(DDD .LT. DDMIN) THEN LP=LP+1 LSMRY(1,LP)=ISTART LSMRY(2,LP)=IEND LSMRY(3,LP)=0 END IF LIR=1 LTR=1 LJP=1 C EMERGENCY EXIT (BUFFER OVER FLOW) 4800 IF(LP .GE. MXLSTB) GO TO 5200 MIR=1000000 MTR=MIR MJP=MIR IF(LIR .LE. LIR0) MIR=LISTIR(LIR) IF(LTR .LE. LTR0) MTR=LISTTR(LTR) IF(LJP .LE. JLISTP) MJP=JLIST(LJP) MIN=900000 INDXIR=0 INDXTR=0 INDXJP=0 IF(MIN .GE. MIR) THEN MIN=MIR INDXIR=1 END IF IF(MIN .GE. MTR) THEN MIN=MTR INDXTR=1 IF(MTR .LT. MIR) INDXIR=0 END IF IF(MIN .GE. MJP) THEN MIN=MJP INDXJP=1 IF(MJP .LT. MIR) INDXIR=0 IF(MJP .LT. MTR) INDXTR=0 END IF IX=INDXIR + 2*INDXTR + 4*INDXJP C EXIT IF(IX .EQ. 0) GO TO 5200 LP=LP+1 LSMRY(1,LP)=0 LSMRY(2,LP)=MIN + ISTART - 1 LSMRY(3,LP)=IX IF(INDXIR .EQ. 1) LIR=LIR+1 IF(INDXTR .EQ. 1) LTR=LTR+1 IF(INDXJP .EQ. 0) GO TO 4800 QQJUMP(LP)=QJUMP(LJP) LJP=LJP+1 GO TO 4800 C C C C PREPERATION FOR THE NEXT SPAN C 5200 CONTINUE IF( JPCHK .EQ. 1 ) AJUMP = AJUMP + AJUMP0 JPCHK = 0 DO 5300 I=1,ORDER 5300 FTRN(I)=FTRN(I)+AJUMP 5400 CONTINUE C C ***************************** C * END OF THE MAIN BLOCK * C ***************************** C C C C SUMMARY OUTPUT C 5500 CONTINUE IF( IIII.GT.2 ) THEN WRITE(6,6410) CALL SICP2(AUTCO,51,NSUM,TAUTCO,L,SGME2,OAIC) CALL SNRASP(TAUTCO,DC,H2,SGME2,L,0,I61,DEG,PAGLEN) END IF C WRITE(6,5800) WRITE(6,5900) IF( LP.LE.0 ) THEN WRITE(6,5910) ELSE DO 5600 I=1,LP IT0=LSMRY(1,I) IT1=LSMRY(2,I) ID=LSMRY(3,I) IF(IT0 .NE. 0) CALL CALEND(JA,JB,JC,HR,IT0,DELT,JY,JM,JD,HRT) CALL CALEND(JA,JB,JC,HR,IT1,DELT,JY1,JM1,JD1,HRT1) IF(IT0 .NE. 0) WRITE(OBUF1,6000) JY,JM,JD,HRT,JY1,JM1,JD1,HRT1 IF(IT0 .EQ. 0) WRITE(OBUF1,6100) JY1,JM1,JD1,HRT1 MARK(1) = ' ' MARK(2) = ' ' MARK(3) = ' ' ID1=ID/2 INDXIR=ID - ID1*2 INDXJP=ID1/2 INDXTR=ID1 - INDXJP*2 IF(ID .EQ. 0) MARK(1) = '*' IF(INDXTR .EQ. 1) MARK(2) = '*' IF(INDXIR .EQ. 1) MARK(3) = '*' IF(INDXJP .EQ. 0) WRITE(6,6200) OBUF1, (MARK(K),K=1,3) IF(INDXJP. EQ. 1) THEN CALL KANDO( COEF, JY1, JM1, JD1, HRT1, 1 ) QNJUMP = - QQJUMP(I) / COEF WRITE(6,6200) OBUF1, (MARK(K),K=1,3), QQJUMP(I), QNJUMP END IF 5600 CONTINUE END IF C C ENTRY IN CASE OF SHIFT = 0 5700 CONTINUE IF( IIII.LE.2 ) STOP AVABIC=AVABIC/COUNT AVSD=AVSD/COUNT AVABIC=AVABIC*NDATA WRITE(6,6400) AVABIC,AVSD STOP C C 5800 FORMAT(1H1/' *** SUMMARY OF THE ANALYSIS ***') 5900 FORMAT(1H / 7X, 'DATE', 17X, 'DATE', 12X, 'ABIC TREND', * ' IRREGULAR STEP -STEP/FACTOR' ) 5910 FORMAT(1H / 7X, 'NONE' // ) 6000 FORMAT(3I4,F6.1,' -',3I4,F6.1) 6100 FORMAT(3I4,F6.1,21X) 6200 FORMAT(1H ,A39,6X,A1,7X,A1,7X,A1,6X,1PD11.4,D18.5) 6400 FORMAT(1H /' *AVABIC =',F13.2,' AVSD =',1PD12.4 ////) 6410 FORMAT(1H1,'SPECTRUM OF TREND') C END SUBROUTINE PRMSET( ITH, ORDER, SPAN, SHIFT, RLIM, LPOUT, FILOUT, * SPECTW, PREPRO, WEIGHT, * IAUG, LAGP, LAGINT, MAXJMP, * KIND, AZ, DMIN, TIMSYS, IFCR, QVAL, MAXITR, * ICETGR, S1WAVE, PAGLEN, IPOTEN, P4FLAG, LOVENM, * RSPERR, SIGMA1, SIGMA2, ZERORM ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER ORDER, SPAN, SHIFT, FILOUT, SPECTW, PREPRO, S1WAVE INTEGER PAGLEN, P4FLAG, RSPERR INTEGER ZERORM C ITH = 1 ORDER = 2 SPAN = 744 SHIFT = 720 RLIM = 0.D0 LPOUT = 0 FILOUT = 0 SPECTW = 1 PREPRO = 0 WEIGHT = 1.D0 IAUG = 0 LAGP = 3 LAGINT = 1 MAXJMP = 5 KIND = 99 AZ = 0.D0 DMIN = 0.25D0 TIMSYS = 0.D0 IFCR = 0 QVAL = 1150.D0 MAXITR = 21 ICETGR = 1 S1WAVE = 0 PAGLEN = 66 IPOTEN = 2 P4FLAG = 1 LOVENM = 2 RSPERR = 0 SIGMA1 = 1.D30 SIGMA2 = 1.D30 ZERORM = 0 C RETURN END SUBROUTINE NAMELS( ITH, ORDER, SPAN, SHIFT, RLIM, LPOUT, FILOUT, * SPECTW, PREPRO, WEIGHT, * IAUG, LAGP, LAGINT, MAXJMP, * KIND, AZ, DMIN, TIMSYS, IFCR, QVAL, MAXITR, * ICETGR, S1WAVE, PAGLEN, IPOTEN, P4FLAG, LOVENM, * RSPERR, SIGMA1, SIGMA2, ZERORM ) C C SIMULATION OF NAMELIST READ ROUTINE C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C INTEGER ORDER, SPAN, SHIFT, FILOUT, SPECTW, PREPRO, S1WAVE INTEGER PAGLEN, P4FLAG, RSPERR INTEGER ZERORM CHARACTER*6 VNAM, PNAME(2) CHARACTER*8 FMTI, FMTD CHARACTER*21 FMT21 CHARACTER*72 BUFF DATA MONITR / 0 / DATA PNAME / '&PARAM', '¶m' / C C ILL = 0 IST = 0 C 10 CONTINUE C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) DO 12 I = 1, 72 12 BUFF(I:I) = ' ' READ(5,16,END=900,ERR=20) (BUFF(I:I),I=1,72) 16 FORMAT(72A1) 20 CONTINUE C C I = 0 40 I = I + 1 IF( I.GT.72 ) GO TO 10 IF( BUFF(I:I).EQ.' ' .OR. BUFF(I:I).EQ.',' ) GO TO 40 C C FIND ONE SYMBOL IS = I IELIM = IS + 6 50 I = I + 1 IF( I.GT.72 ) GO TO 60 IF( I.GT.IELIM ) THEN WRITE(6,52) BUFF ILL = 1 GO TO 10 END IF IF( BUFF(I:I).NE.'=' .AND. BUFF(I:I).NE.' ' ) GO TO 50 60 IE = I - 1 IL = I - IS VNAM = ' ' VNAM(1:IL) = BUFF(IS:IE) C C IF( VNAM.EQ.PNAME(1) .OR. VNAM.EQ.PNAME(2) ) THEN IST = 1 GO TO 40 ELSE IF( VNAM.EQ.'&END ' .OR. VNAM.EQ.'&end ' ) THEN IF( ILL.NE.0 ) STOP IF( IST.EQ.0 ) THEN WRITE(6,910) PNAME(1), BUFF STOP END IF RETURN ELSE IF( IST.EQ.0 ) THEN WRITE(6,910) PNAME(1), BUFF STOP END IF C C 70 I = I + 1 IF( I.GT.72 ) THEN WRITE(6,72) BUFF ILL = 1 GO TO 10 END IF IF( BUFF(I:I).EQ.' ' .OR. BUFF(I:I).EQ.'=' ) GO TO 70 C C FIND ONE VALUE IS = I 80 I = I + 1 IF( I.GT.72 ) THEN WRITE(6,82) BUFF ILL = 1 GO TO 10 END IF IF( BUFF(I:I).NE.',' .AND. BUFF(I:I).NE.' ' ) GO TO 80 IE = I - 1 IL = I - IS C C MONITOR PRINT OUT IF( MONITR.GE.1 ) THEN FMT21 = '(1H ,A6,1X,I2,2X,A72)' WRITE(FMT21(19:20),'(I2)') IL WRITE(6,FMT21) VNAM, IL, BUFF(IS:IE) END IF C WRITE(FMTI,'(A2,I2,A4)') '(I', IL, ') ' WRITE(FMTD,'(A2,I2,A4)') '(F', IL, '.0) ' IF( VNAM.EQ.'ITH ' ) THEN READ(BUFF(IS:IE),FMTI) ITH ELSE IF( VNAM.EQ.'ORDER ' .OR. VNAM.EQ.'order ' ) THEN READ(BUFF(IS:IE),FMTI) ORDER ELSE IF( VNAM.EQ.'SPAN ' .OR. VNAM.EQ.'span ' ) THEN READ(BUFF(IS:IE),FMTI) SPAN ELSE IF( VNAM.EQ.'SHIFT ' .OR. VNAM.EQ.'shift ' ) THEN READ(BUFF(IS:IE),FMTI) SHIFT ELSE IF( VNAM.EQ.'RLIM ' .OR. VNAM.EQ.'rlim ' ) THEN READ(BUFF(IS:IE),FMTD) RLIM ELSE IF( VNAM.EQ.'LPOUT ' .OR. VNAM.EQ.'lpout ' ) THEN READ(BUFF(IS:IE),FMTI) LPOUT ELSE IF( VNAM.EQ.'FILOUT' .OR. VNAM.EQ.'filout' ) THEN READ(BUFF(IS:IE),FMTI) FILOUT ELSE IF( VNAM.EQ.'SPECTW' .OR. VNAM.EQ.'spectw' ) THEN READ(BUFF(IS:IE),FMTI) SPECTW ELSE IF( VNAM.EQ.'PREPRO' .OR. VNAM.EQ.'prepro' ) THEN READ(BUFF(IS:IE),FMTI) PREPRO ELSE IF( VNAM.EQ.'WEIGHT' .OR. VNAM.EQ.'weight' ) THEN READ(BUFF(IS:IE),FMTD) WEIGHT ELSE IF( VNAM.EQ.'IAUG ' .OR. VNAM.EQ.'iaug ' ) THEN READ(BUFF(IS:IE),FMTI) IAUG ELSE IF( VNAM.EQ.'LAGP ' .OR. VNAM.EQ.'lagp ' ) THEN READ(BUFF(IS:IE),FMTI) LAGP ELSE IF( VNAM.EQ.'LAGINT' .OR. VNAM.EQ.'lagint' ) THEN READ(BUFF(IS:IE),FMTI) LAGINT ELSE IF( VNAM.EQ.'MAXJMP' .OR. VNAM.EQ.'maxjmp' ) THEN READ(BUFF(IS:IE),FMTI) MAXJMP ELSE IF( VNAM.EQ.'KIND ' .OR. VNAM.EQ.'kind ' ) THEN READ(BUFF(IS:IE),FMTI) KIND ELSE IF( VNAM.EQ.'AZ ' .OR. VNAM.EQ.'az ' ) THEN READ(BUFF(IS:IE),FMTD) AZ ELSE IF( VNAM.EQ.'DMIN ' .OR. VNAM.EQ.'dmin ' ) THEN READ(BUFF(IS:IE),FMTD) DMIN ELSE IF( VNAM.EQ.'TIMSYS' .OR. VNAM.EQ.'timsys' ) THEN READ(BUFF(IS:IE),FMTD) TIMSYS ELSE IF( VNAM.EQ.'IFCR ' .OR. VNAM.EQ.'ifcr ' ) THEN READ(BUFF(IS:IE),FMTI) IFCR ELSE IF( VNAM.EQ.'QVAL ' .OR. VNAM.EQ.'qval ' ) THEN READ(BUFF(IS:IE),FMTD) QVAL ELSE IF( VNAM.EQ.'MAXITR' .OR. VNAM.EQ.'maxitr' ) THEN READ(BUFF(IS:IE),FMTI) MAXITR ELSE IF( VNAM.EQ.'ICETGR' .OR. VNAM.EQ.'icetgr' ) THEN READ(BUFF(IS:IE),FMTI) ICETGR ELSE IF( VNAM.EQ.'S1WAVE' .OR. VNAM.EQ.'s1wave' ) THEN READ(BUFF(IS:IE),FMTI) S1WAVE ELSE IF( VNAM.EQ.'PAGLEN' .OR. VNAM.EQ.'paglen' ) THEN READ(BUFF(IS:IE),FMTI) PAGLEN ELSE IF( VNAM.EQ.'IPOTEN' .OR. VNAM.EQ.'ipoten' ) THEN READ(BUFF(IS:IE),FMTI) IPOTEN ELSE IF( VNAM.EQ.'P4FLAG' .OR. VNAM.EQ.'p4flag' ) THEN READ(BUFF(IS:IE),FMTI) P4FLAG ELSE IF( VNAM.EQ.'LOVENM' .OR. VNAM.EQ.'lovenm' ) THEN READ(BUFF(IS:IE),FMTI) LOVENM ELSE IF( VNAM.EQ.'RSPERR' .OR. VNAM.EQ.'rsperr' ) THEN READ(BUFF(IS:IE),FMTI) RSPERR ELSE IF( VNAM.EQ.'SIGMA1' .OR. VNAM.EQ.'sigma1' ) THEN READ(BUFF(IS:IE),FMTD) SIGMA1 ELSE IF( VNAM.EQ.'SIGMA2' .OR. VNAM.EQ.'sigma2' ) THEN READ(BUFF(IS:IE),FMTD) SIGMA2 ELSE IF( VNAM.EQ.'ZERORM' .OR. VNAM.EQ.'zerorm' ) THEN READ(BUFF(IS:IE),FMTI) ZERORM ELSE WRITE(6,92) VNAM, BUFF CCC ILL = 1 END IF C GO TO 40 C 900 WRITE(6,920) BUFF STOP C 52 FORMAT(/' ERROR IN NAMELS TOO LONG SYMBOL '/' REC : ',A72) 72 FORMAT(/' ERROR IN NAMELS = NOT FOUND'/' REC : ',A72) 82 FORMAT(/' ERROR IN NAMELS VALUE NOT FOUND'/' REC : ',A72) 92 FORMAT(/' WARNING IN NAMELS SYMBOL ',A6,' NOT FOUND'/ * ' REC : ',A72) 910 FORMAT(/' ERROR IN NAMELS NO ',A6 /' REC : ',A72) 920 FORMAT(/' ERROR IN NAMELS NO &END' /' REC : ',A72) C END SUBROUTINE PRMCHK( ITH, ORDER, SPAN, SHIFT, RLIM, LPOUT, FILOUT, * SPECTW, PREPRO, WEIGHT, * IAUG, LAGP, LAGINT, MAXJMP, * KIND, AZ, DMIN, TIMSYS, IFCR, QVAL, MAXITR, * ICETGR, S1WAVE, PAGLEN, IPOTEN, P4FLAG, LOVENM, * RSPERR, SIGMA1, SIGMA2, ZERORM ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER ORDER, SPAN, SHIFT, FILOUT, SPECTW, PREPRO, S1WAVE INTEGER PAGLEN, P4FLAG, RSPERR INTEGER ZERORM C IF( ITH .LT. 0 ) ITH = 0 IF( ITH .GT. 1 ) ITH = 1 IF( ORDER .LT. 2 ) ORDER = 2 IF( ORDER .GT. 3 ) ORDER = 3 IF( SHIFT .LT. 0 ) SHIFT = 0 IF( SPAN .NE. 0 ) THEN IF( SPAN .LT. 24 ) SPAN = 24 IF( SHIFT .GT. SPAN ) SHIFT = SPAN END IF IF( RLIM .LT. -1.D25 ) RLIM = -1.D25 IF( RLIM .GT. 1.D25 ) RLIM = 1.D25 IF( LPOUT .LT. 0 ) LPOUT = 0 IF( LPOUT .GT. 2 ) LPOUT = 2 IF( FILOUT .LT. 0 ) FILOUT = 0 IF( FILOUT .GT. 2 ) FILOUT = 2 IF( SPECTW .LT. 0 ) SPECTW = 0 IF( SPECTW .GT. 6 ) SPECTW = 6 IF( PREPRO .LT. 0 ) PREPRO = 0 IF( PREPRO .GT. 3 ) PREPRO = 3 IF( SHIFT .EQ. 0 ) THEN LPOUT = 0 FILOUT = 0 PREPRO = 0 END IF IF( WEIGHT .LT. 0.D0 ) WEIGHT = 0.D0 IF( WEIGHT .GT. 1.D2 ) WEIGHT = 1.D2 IF( IAUG .LT. 0 ) IAUG = 0 IF( LAGP .LT. 0 ) LAGP = 0 IF( LAGINT .LT. 1 ) LAGINT = 1 IF( MAXJMP .LT. 0 ) MAXJMP = 0 IF( ITH .EQ. 0 ) THEN IFCR = 0 IPOTEN = 1 END IF IF( KIND .LT. 1 ) KIND = 99 IF( KIND .GT. 8 ) KIND = 99 IF( DABS(AZ).GT.360.D0 ) WRITE(6,100) AZ IF( DMIN .LE. 0.D0 ) DMIN = 1.D0/256.D0 IF( DMIN .GE. 250.D0 ) DMIN = 128.D0 IF( DABS(TIMSYS).GT.24.D0 ) WRITE(6,110) TIMSYS IF( IFCR .LT. 0 ) IFCR = 0 IF( IFCR .GT. 1 ) IFCR = 1 IF( QVAL .LT.1.D1 ) QVAL = 1.D1 IF( QVAL .GT.1.D5 ) QVAL = 1.D5 IF( MAXITR .LT. 1 ) MAXITR = 1 IF( MAXITR .GT.30 ) MAXITR = 30 IF( ICETGR .LT. 1 ) ICETGR = 1 IF( ICETGR .GT. 2 ) ICETGR = 2 IF( S1WAVE .LT. 0 ) S1WAVE = 0 IF( S1WAVE .GT. 1 ) S1WAVE = 1 IF( PAGLEN .LT.66 ) PAGLEN = 60 IF( PAGLEN .GT.66 ) PAGLEN = 66 IF( IPOTEN .LT. 1 ) IPOTEN = 1 IF( IPOTEN .GT. 2 ) IPOTEN = 2 IF( IPOTEN .EQ. 1 ) THEN P4FLAG = 0 ELSE IF( P4FLAG. LT. 0 ) P4FLAG = 0 IF( P4FLAG. GT. 1 ) P4FLAG = 1 END IF IF( LOVENM .LT. 0 ) LOVENM = 0 IF( LOVENM .GT. 2 ) LOVENM = 0 IF( KIND.GE.4 .AND. KIND.LE.6 ) LOVENM = 1 IF( RSPERR .NE. 0 ) RSPERR = 1 C IF( SIGMA1 .LE. 0.0D0 ) SIGMA1 = 1.D30 IF( SIGMA2 .LE. 0.0D0 ) SIGMA2 = 1.D30 IF( ZERORM .LT. 0 ) ZERORM = 0 IF( ZERORM .GT. 1 ) ZERORM = 1 C RETURN C 100 FORMAT(1H /' INVALID AZIMUTH AZ =',F10.1/) 110 FORMAT(1H /' INVALID LOCAL TIME TIMSYS =',F10.1/) C END SUBROUTINE SETGRP( OBSLEN, IGRPLV, ICETGR, IPOTEN ) C C GROUPING OF WAVES (AUTO) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( MTBLEN = 31 ) C CHARACTER*8 SYMBOL, KSYMBL(MAXGRP,8) DIMENSION KG2(MAXGRP,8) C COMMON /STNGRP/ IG1(36),IG2(36),NGROUP,IMWAVE(36) COMMON /STNGRC/ SYMBOL(36) C C IF( IGRPLV.EQ.0 ) THEN J = 1 IF( OBSLEN.GE. 8.D0 ) J = 2 IF( OBSLEN.GE. 16.D0 ) J = 3 IF( OBSLEN.GE. 180.D0 ) J = 4 IF( OBSLEN.GE. 364.D0 ) J = 5 IF( OBSLEN.GE.2000.D0 ) J = 6 ELSE IF( IGRPLV.GE.-2 .AND. IGRPLV.LE.6 ) THEN J = IGRPLV + 2 IF( IGRPLV.LT.0 ) J = IGRPLV + 3 ELSE WRITE(6,200) IGRPLV 200 FORMAT(1H /' ERROR STOPPED AT SETGRP IGRPLV =', I5) STOP END IF C C IF( ICETGR.GE.1 .AND. ICETGR.LE.2 ) THEN IGX = ICETGR ELSE WRITE(6,210) ICETGR 210 FORMAT(1H /' WARNING IN SETGRP ICETGR =',I8) IGX = 1 END IF C C IF( IPOTEN.EQ.1 .AND. IGX.EQ.1 ) THEN OPEN( 95, FILE = '/home/tamura/baytap/GROUP31.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/GROUP31.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/GROUP31.DAT', C * STATUS = 'OLD' ) ELSE IF( IPOTEN.EQ.1 .AND. IGX.EQ.2 ) THEN OPEN( 95, FILE = '/home/tamura/baytap/GROUP32.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/GROUP32.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/GROUP32.DAT', C * STATUS = 'OLD' ) ELSE IF( IPOTEN.EQ.2 .AND. IGX.EQ.1 ) THEN OPEN( 95, FILE = '/home/tamura/baytap/GROUP41.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/GROUP41.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/GROUP41.DAT', C * STATUS = 'OLD' ) ELSE IF( IPOTEN.EQ.2 .AND. IGX.EQ.2 ) THEN OPEN( 95, FILE = '/home/tamura/baytap/GROUP42.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/GROUP42.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/GROUP42.DAT', C * STATUS = 'OLD' ) END IF C C GROUP NOUMBER DO 10 I = 1, MAXGRP READ(95,220) ( KG2(I,K), K = 1, 8 ) 10 CONTINUE C SKIP DO 12 I = MAXGRP+1, MTBLEN 12 READ(95,220) KG2IK KG2IK = KG2IK C SYMBOL NAME DO 20 I = 1, MAXGRP READ(95,230) ( KSYMBL(I,K), K = 1, 8 ) 20 CONTINUE 220 FORMAT(8I9) 230 FORMAT(8(1X,A8)) C CLOSE( 95 ) C C DO 30 I = 1, MAXGRP IF( KG2(I,J).EQ.0 ) GO TO 40 NGROUP = I 30 CONTINUE 40 CONTINUE C C IF( NGROUP.GT.36 ) THEN WRITE(6,240) 240 FORMAT(1H /' ERROR IN SETGRP TOO MANY GROUPS') STOP END IF C C DO 110 I = 1, NGROUP 110 IG2(I) = KG2(I,J) IG1(1) = 1 DO 120 I = 2, NGROUP 120 IG1(I) = KG2(I-1,J) + 1 DO 130 I = 1, NGROUP 130 SYMBOL(I) = KSYMBL(I,J) C RETURN C END SUBROUTINE TBREAD( IPOTEN ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) C IF( IPOTEN.EQ.1 ) THEN C OPEN( 95, FILE = '/home/tamura/baytap/TABLE3.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/TABLE3.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/TABLE3.DAT', C * STATUS = 'OLD' ) C DO 10 I = 1, 377 READ(95,110) IP(I), (IDN(I,J),J=1,6), AVC(I), AMC(I) 10 CONTINUE 110 FORMAT(7X,I1,1X,6I3,8X,F12.8,F12.6) C ELSE IF( IPOTEN.EQ.2 ) THEN C OPEN( 95, FILE = '/home/tamura/baytap/TABLE4.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/TABLE4.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/TABLE4.DAT', C * STATUS = 'OLD' ) C DO 20 I = 1, 909 READ(95,120) IP(I),(IDN(I,J),J=1,8),AVC(I),AMC(I),AMC2(I) 20 CONTINUE 120 FORMAT(7X,I1,1X,8I3,2X,F12.8,F12.6,F11.6) C ELSE WRITE(6,130) 130 FORMAT(1H /' ERROR IPOTEN?') STOP END IF C CLOSE( 95 ) C RETURN END SUBROUTINE SETMWV( IPOTEN ) C C IMWAVE(I) : MAIN WAVE NO. IN EACH GROUP C EXCEPT( MU2 ---> 2N2 ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /STNGRP/ IG1(36),IG2(36),NGROUP,IMWAVE(36) C C IF( IPOTEN.EQ.1 ) THEN IWN1 = 240 IWN2 = 233 IWN3 = 377 ELSE IF( IPOTEN.EQ.2 ) THEN IWN1 = 534 IWN2 = 520 IWN3 = 909 END IF C C DO 30 I = 1, NGROUP JS = IG1(I) JE = IG2(I) NW = 0 AMP = 0.D0 DO 10 J = JS, JE IF( DABS(AMC(J)).LT.AMP ) GO TO 10 AMP = DABS(AMC(J)) NW = J 10 CONTINUE IF( NW.EQ.IWN1 .AND. JS.LE.IWN2 ) NW = IWN2 IF( NW.LT.1 .OR. NW.GT.IWN3 ) THEN WRITE(6,200) I, NW 200 FORMAT(1H /' ERROR AT SETMWV I, NW =', 2I6 ) STOP END IF IMWAVE(I) = NW 30 CONTINUE C RETURN END SUBROUTINE DATAR2( RAT,Y,IY,NDATA,RLIM,LA,LB,LC,HLR,DELT,NFILE, * SLAT,SLONG,SH,G, KIND,AZ, ZERORM ) C C THIS SUBROUTINE READS DATA Y FROM UNIT NO. MT C INPUTS C IY : 0:TIDAL DATA; 1:ASSOSIATED DATA C MT : INPUT DEVICE SPECIFICATION (UNIT NO.) C IFORM : FORMAT SPECIFICATION C : ( 1: USER DEF.; 0: 4 DEG.; -1: 5 DEG.; 2: 7 DEG. ) C OUTPUTS C Y : DATA VECTOR C NDATA : DATA LENGTH C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( LBUFF = 4500 ) PARAMETER ( MAXASC = 3 ) PARAMETER ( MAXCLB = 100 ) C DIMENSION Y(*), RAT(*) DIMENSION ETAL(MAXCLB), TJE(MAXCLB) DIMENSION JACAL(MAXCLB), JBCAL(MAXCLB), JCCAL(MAXCLB) DIMENSION HDCAL(MAXCLB) DIMENSION DELTCK(6), X1(12) DIMENSION YO(LBUFF,2), YA(LBUFF,2,MAXASC) REAL YT(LBUFF,2,MAXGR2) INTEGER ZERORM C CHARACTER*4 XA(12), DUMMY CHARACTER*10 CCKIND CHARACTER*80 TITLE C SAVE DIF, JABCLS C COMMON /LPCNT/ LPOUT COMMON /INFORM/ IFORMT, NETALT, ETALT(MAXCLB), TJET(MAXCLB) COMMON /BUFOTA/ YO, YT, YA, INCORE COMMON /PRSTEP/ STEP00 C DATA IRAT / 0 / DATA DELTCK / 0.5D0, 1.D0, 1.5D0, 2.D0, 3.D0, 6.D0 / DATA EPS12 / 1.D-12 / C CCC EXTERNAL FNORMG, TJUL C **************************** * PARAMETERS READING * **************************** C C I/O UNIT, FORMAT C C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) ISTFLG = 0 READ(5,5000,ERR=5001) MT, IFORM, ISTFLG 5001 IF( MT.LE.0 .OR. MT.EQ.5 .OR. MT.EQ.16 * .OR. (MT.GE.30 .AND. MT.LE.34) * .OR. MT.EQ.95 .OR. MT.GE.100 ) THEN WRITE(6,6000) MT STOP END IF IF(IFORM.EQ.2) THEN READ(5,5000) MTC IF( MTC.LE.0 .OR. MTC.EQ.5 .OR. MTC.EQ.16 * .OR. (MTC.GE.30 .AND. MTC.LE.34) * .OR. MTC.EQ.95 .OR. MTC.GE.100 ) THEN WRITE(6,6000) MTC STOP END IF END IF IF( IY.EQ.0 ) IFORMT = IFORM 5000 FORMAT(3I5) 6000 FORMAT(/' INVALID I/O UNIT NUMBER UNIT =',I5) C C TITLE OF DATASET C C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) DO 5002 I = 1, 80 5002 TITLE(I:I) = ' ' READ(MT,5010,ERR=5008) (TITLE(I:I),I=1,80) 5008 WRITE(6,6010) TITLE 5010 FORMAT(80A1) 6010 FORMAT(1H //' ***** TOP OF ORIGINAL DATA *****' // 1H ,A80) C C LATITUDE, LONGITUDE, ALTTUDE, GRAVITY, KIND C C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) G = 0.D0 CCKIND = ' ' READ(MT,5020,ERR=5018) SLAT, SLONG, SH, G, CCKIND 5018 IF( G.EQ.0.D0 ) G = FNORMG( SLAT, SH ) WRITE(6,6020) SLAT, SLONG, SH, G 5020 FORMAT(4F10.0,A10) 6020 FORMAT(/' LATITUDE LONGITUDE HEIGHT GRAVITY' / * 1H ,F10.4,F11.4,F10.2,F10.4 / ) C ILL = 0 IF( DABS(SLAT) .GT. 90.D0 ) ILL = 1 IF( DABS(SLONG).GE. 360.D0 ) ILL = 1 IF( DABS(SH) .GT.9999.D0 ) ILL = 1 IF( G.LT.977.D0 .OR. G.GE.984.D0 ) ILL = 1 IF( ILL.NE.0 ) THEN WRITE(6,6030) STOP END IF 6030 FORMAT(/' INPUT GEOGRAPHICAL PARAMETERS MAY BE WRONG') C ************************************* * CASE OF USER DEFINED FORMAT * ************************************* C IRAT = IRAT + 1 IF( IFORM.EQ.1 ) THEN CALL DATAR(RAT,IRAT,Y,IY,NDATA,RLIM,LA,LB,LC,HLR,DELT, * NFILE,MT,DIF, ZERORM) IF( IY.EQ.0 ) STEP00 = 0.D0 RETURN END IF C ********************* * OTHER CASES * ********************* C C KIND OF TIDAL DATASET FOR IFORM = 2 C IF( KIND.EQ.99 .AND. IFORM.EQ.2 .AND. IY.EQ.0 ) THEN IF( CCKIND(1:7).EQ.'GRAVITY' ) KIND = 1 IF( CCKIND(1:5).EQ.'TILT ' ) THEN IF( CCKIND(6:7).EQ.'NS' ) THEN KIND = 2 ELSE IF( CCKIND(6:7).EQ.'EW' ) THEN KIND = 3 ELSE KIND = 2 READ(CCKIND(6:10),'(F5.1)') AZ END IF END IF IF( CCKIND(1:6).EQ.'STRAIN' ) KIND = 7 IF( CCKIND(1:6).EQ.'OCEAN ' ) KIND = 8 END IF IF( KIND.LE.0 .OR. KIND.GE.9 ) THEN WRITE(6,6040) KIND STOP END IF 6040 FORMAT(/' INVALID KIND KIND =',I5) C C START TIME, SAMPLING INTERVAL C IF( HLR.LT.0.D0 .OR. HLR.GT.23.5D0 ) THEN WRITE(6,6050) HLR STOP END IF C ILL = 1 DO 10 I = 1, 6 IF( DABS( DELT - DELTCK(I) ) .LT. EPS12 ) ILL = 0 10 CONTINUE IF( ILL.NE.0 ) THEN WRITE(6,6060) DELT STOP END IF C HLRC = DELT * DBLE( IDINT( (HLR+0.001D0) / DELT ) ) IF( DABS(HLRC-HLR) .GT. EPS12 ) WRITE(6,6070) HLRC HLR = HLRC C IYMDA = LA*10000 + LB*100 + LC IYMDA = IYMDA*10 + IDINT( (HLR+0.001D0)/6.D0 ) + 1 ISAMPL= IDINT( (DELT+0.001D0)/0.5D0 ) C 6050 FORMAT(/' ERROR AT DATAR2 INVALID TIME HOUR =',F8.1) 6060 FORMAT(/' ERROR AT DATAR2 INVALID SAMPLING INTERVAL', * ' DELT =', F8.1 ) 6070 FORMAT(/' WARNING BEGINNING TIME IS EXCHANGED', * ' TIME =',F5.1 /) C C CALIBRATION DATA C IF( IFORM.EQ.2 ) THEN READ(MTC,5060) DUMMY READ(MTC,5060) DUMMY DUMMY = DUMMY I = 1 20 READ(MTC,5080,END=40) JA,JB,JC,HD, ETAL(I) TJE(I) = TJUL(JA,JB,JC,HD) JACAL(I) = JA JBCAL(I) = JB JCCAL(I) = JC HDCAL(I) = HD NETAL = I I = I + 1 IF( I.LE.MAXCLB ) GO TO 20 WRITE(6,6090) I STOP ELSE READ(MT,5070) NETAL DO 30 I = 1, NETAL IF( I.GT.MAXCLB ) THEN WRITE(6,6090) I STOP END IF READ(MT,5080) JA,JB,JC,HD, ETAL(I) TJE(I) = TJUL(JA,JB,JC,HD) JACAL(I) = JA JBCAL(I) = JB JCCAL(I) = JC HDCAL(I) = HD 30 CONTINUE END IF C 40 IF( NETAL.LT.2 ) THEN WRITE(6,6100) NETAL STOP ELSE IF( NETAL.LE.6 ) THEN DO 42 I = 1, NETAL 42 WRITE(6,6080) JACAL(I),JBCAL(I),JCCAL(I),HDCAL(I), ETAL(I) ELSE NETALX = NETAL/2 NETALY = NETAL - NETALX*2 DO 44 I = 1, NETALX J = I + NETALX + NETALY WRITE(6,6082) JACAL(I),JBCAL(I),JCCAL(I),HDCAL(I), ETAL(I), * JACAL(J),JBCAL(J),JCCAL(J),HDCAL(J), ETAL(J) 44 CONTINUE IF( NETALY.EQ.1 ) * WRITE(6,6080) JACAL(I),JBCAL(I),JCCAL(I),HDCAL(I), ETAL(I) END IF C NETAL1 = NETAL - 1 C 5060 FORMAT(A4) 5070 FORMAT(I5) 5080 FORMAT(3I5,F5.0,F10.5) 6080 FORMAT(1H ,'Y,M,D,HR =', I6,I4,I4, F6.1,' CALIB.=',1PD13.5) 6082 FORMAT(1H ,'Y,M,D,HR =', I6,I4,I4, F6.1,' CALIB.=',1PD13.5, * 12X,'Y,M,D,HR =', I6,I4,I4,0PF6.1,' CALIB.=',1PD13.5) 6090 FORMAT(/' TOO MANY CARIBRATION CARDS N =',I5/ * ' MODIFY MAXCLB IN PARAMETER STATEMENTS') 6100 FORMAT(/' TOO FEW CARIBRATION CARDS N =',I5) C C PREPARATION C JMPCHK = 0 DT = 0.5D0/24.D0 IEND = 12/ISAMPL YSUM = 0.D0 YMAX = -1.D25 YMIN = 1.D25 NUM = 0 IPLUS = IDINT( ( DMOD( HLR+0.001D0, 6.D0 ) + 0.001D0 ) / DELT ) KLIM = ( NDATA + IPLUS - 1 ) / IEND + 1 NDATAX = 0 KKK = 0 JABC0 = 0 STEP = 0.D0 NETALK = 1 C ********************** * DATA READING * ********************** C 1000 CONTINUE C C DATA READ C IF( IFORM.EQ.0 ) THEN READ(MT,5200,END=1100) (X1(M),M=1,12),NST,STVAL,JA,JB,JC,IHC, * (XA(M),M=1,12), JABC IF( XA(1).EQ.'9999' ) GO TO 1100 ELSE IF( IFORM.EQ.-1 ) THEN READ(MT,5210,END=1100) (X1(M),M=1,12),NST,STVAL,JA,JB,JC,IHC, * (XA(M),M=1,12), JABC IF( X1(1).GE.9999.0 ) GO TO 1100 ELSE IF( IFORM.EQ.2 ) THEN READ(MT,5220,END=1100) JA,JB,JC,IHC,(X1(M),M=1,12),NST,STVAL, * JABC, (XA(M),M=1,12) IF( JABC.GE.9999000 ) GO TO 1100 ELSE WRITE(6,6200) IFORM STOP END IF C 5200 FORMAT( 12F4.1,I3,F5.1,17X,3I2,I1, T1,12A4,25X,I7 ) 5210 FORMAT( 12F5.1,I3,F5.1,5X,3I2,I1, T1,12(1X,A4),13X,I7 ) 5220 FORMAT( 3I2,I1,12(F7.1,1X),I2,F8.1, T1,I7,3X,12(A4,4X) ) 6200 FORMAT(/' ERROR INVALID FORMAT SPECIFICATION IFORM =',I5) C C CHECK OF DATE SEQUENCE IF( JABC0.NE.0 ) THEN IF( IHC.LE.0 ) THEN WRITE(6,6300) MT, JABC0, JABC STOP ELSE IF( IHC.EQ.1 ) THEN IF( JABC0+7.NE.JABC ) THEN IF( JC.NE.1 ) THEN WRITE(6,6300) MT, JABC0, JABC STOP END IF JALX = JABC0/100000 JABC0 = JABC0 - JALX*100000 JBLX = JABC0/1000 JCLX = (JABC0-JBLX*1000)/10 TIMLX = TJUL( JALX, JBLX, JCLX, 0.D0 ) TIMNW = TJUL( JA, JB, JC, 0.D0 ) IF( DABS(TIMNW-TIMLX-1.D0) .GE. 1.D-12 ) THEN WRITE(6,6300) MT, JABC0, JABC STOP END IF END IF ELSE IF( IHC.LE.4 ) THEN IF( JABC0+1.NE.JABC ) THEN WRITE(6,6300) MT, JABC0, JABC STOP END IF ELSE WRITE(6,6300) MT, JABC0, JABC STOP END IF END IF C IF( NST.LT.0 .OR. NST.GT.12 ) THEN WRITE(6,6310) MT, JABC STOP END IF JABC0 = JABC C 6300 FORMAT(/' SEQUENCE OF DATA IS WRONG', * ' I/O UNIT =',I3,' DATE =',I8,' /',I8) 6310 FORMAT(/' POSITION OF STEP IS WRONG', * ' I/O UNIT =',I3,' DATE =',I8) C C SKIP ? C IF(JABC.LT.IYMDA) THEN IF( NST.EQ.0 ) GO TO 1000 IF( STVAL.EQ.0.D0 ) GO TO 1000 HD = 0.5D0 * DBLE( 12*IHC + NST - 13 ) T = TJUL( JA, JB, JC, HD ) DO 50 K = NETALK, NETAL1 IF( T.GE.TJE(K) .AND. T.LT.TJE(K+1) ) GO TO 60 50 CONTINUE WRITE(6,6500) STOP 60 NETALK = K DET = ( ETAL(K+1) - ETAL(K) ) / ( TJE(K+1) - TJE(K) ) SENSE = DET*( T - TJE(K) ) + ETAL(K) STEP = STEP + STVAL*SENSE GO TO 1000 END IF 6500 FORMAT(/' INVALID CALIBRATION CARD STOPPED AT DATAR2') C C END ?, START DATE CHECK, STEP C IF(KKK.GT.KLIM) GO TO 1100 KKK = KKK + 1 IF( KKK.EQ.1 ) THEN IF( IY.EQ.0 ) STEP00 = STEP IF( JABC.GT.IYMDA ) THEN HLR = DBLE( (IHC-1)*6 ) IPLUS = 0 END IF IF( IY.EQ.0 ) JABCLS = JABC IF( JABCLS.NE.JABC ) THEN WRITE(6,6400) STOP END IF JABCLS = JABC LA = JA LB = JB LC = JC END IF 6400 FORMAT(/' ERROR AT DATAR2 INCOMPATIBLE BEGINNING DATE') C C MULTIPLICATION OF SENSITIVITY C HD = DBLE( (IHC-1)*6 ) TJ0 = TJUL(JA,JB,JC,HD) IF( NST.EQ.0 ) THEN MSKIP = ISAMPL ELSE MSKIP = 1 END IF DO 80 M = 1, 12, MSKIP T = TJ0 + DT*(M-1) DO 70 K = NETALK, NETAL1 IF( T.GE.TJE(K) .AND. T.LT.TJE(K+1) ) GO TO 72 70 CONTINUE WRITE(6,6500) STOP 72 NETALK = K DET = ( ETAL(K+1) - ETAL(K) ) / ( TJE(K+1) - TJE(K) ) SENSE = DET*( T - TJE(K) ) + ETAL(K) IF( M.EQ.NST ) STEP = STEP + STVAL*SENSE X1(M) = X1(M)*SENSE + STEP 80 CONTINUE C C POSITION OF STEP C C NST 0 1 2 3 4 5 6 7 8 9 10 11 12 ORIGINAL POSITION C NSTC 0 1 2 3 4 5 6 7 8 9 10 11 12 FOR ISAMPL = 1 C NSTC 0 1 1 3 3 5 5 7 7 9 9 11 11 = 2 C NSTC 0 1 1 1 4 4 4 7 7 7 10 10 10 = 3 C NSTC 0 1 1 1 1 5 5 5 5 9 9 9 9 = 4 C NSTC 0 1 1 1 1 1 1 7 7 7 7 7 7 = 6 C NSTC 0 1 1 1 1 1 1 1 1 1 1 1 1 = 12 C IF( NST.EQ.0 ) THEN NSTC = 0 ELSE NSTC = NST - MOD( NST-1, ISAMPL ) END IF C C RESAMPLING C DO 180 I = 1, 12, ISAMPL X1I = X1(I) C CC IF( XA(I).NE.' ' ) THEN C NORMAL DATA CC NUM = NUM + 1 CC YSUM = YSUM + X1I CC YMAX = DMAX1( YMAX, X1I ) CC YMIN = DMIN1( YMIN, X1I ) CC ELSE C NO OBSERVATION CC X1I = 1.D30 CC END IF C IF( XA(I).EQ.' ' .OR. * ( ZERORM.EQ.1 .AND. DABS(X1I).LT.1.0D-14 ) ) THEN C NO OBSERVATION X1I = 1.D30 ELSE C NORMAL DATA NUM = NUM + 1 YSUM = YSUM + X1I YMAX = DMAX1( YMAX, X1I ) YMIN = DMIN1( YMIN, X1I ) END IF C IF( IY.EQ.0 ) THEN IF( I.EQ.NSTC .AND. JMPCHK.EQ.1 ) THEN C STEP MARK IF( STVAL.EQ.0.D0 .OR. ISTFLG.EQ.1 ) X1I = -1.D30 END IF ELSE IF( I.EQ.NSTC .AND. JMPCHK.EQ.1 .AND. STVAL.EQ.0.D0 ) THEN WRITE(6,6600) MT, JABC 6600 FORMAT(/' WARNING I/O UNIT =', I3, * ' STEP CORRECTION VALUE OF ', * 'ASSOCIATED DATA IS MISSING DATE =', I8 ) X1I = 1.D30 END IF END IF C NDATAX = NDATAX + 1 IF( X1I.LT.1.D20 .AND. NDATAX.GT.IPLUS ) JMPCHK = 1 Y(NDATAX) = X1I 180 CONTINUE C GO TO 1000 C ****************************** * END OF DATA, SAVING * ****************************** C 1100 CONTINUE IF( NUM.LE.0 ) THEN WRITE(6,6700) STOP END IF 6700 FORMAT(/' ERROR AT DATAR2 NO DATA') C C NORMALIZATION C DIF0 = YMAX - YMIN IF( IY.EQ.0 ) DIF = DIF0 RATIO = DIF/DIF0 RAT(IRAT) = RATIO IF( IY.EQ.0 ) THEN WRITE(6,6800) MT, RATIO ELSE YMEAN = YSUM/NUM DO 1110 I = 1, NDATA IF( DABS(Y(I)) .LE. 1.D25 ) Y(I) = (Y(I)-YMEAN)*RATIO 1110 CONTINUE WRITE(6,6810) MT, RATIO, YMEAN END IF 6800 FORMAT(/' TIDAL DATA I/O UNIT =', I3, * ' DATA RANGE RATIO =',1PD12.4) 6810 FORMAT(/' ASSOCIATED DATA I/O UNIT =', I3, * ' DATA RANGE RATIO =',1PD12.4, * ' MEAN VALUE =',D12.4) C C SAVING C NFILE = NFILE + 1 IF(NDATA .GT. NDATAX-IPLUS) NDATA = NDATAX - IPLUS C IF( INCORE.EQ.1 ) THEN IF(IY.EQ.0) THEN CCC DO NOT CHECK SUBSCRIPT RANGE DO 810 I = 1, NDATA 810 YO(I,1) = Y(I+IPLUS) ELSE CCC DO NOT CHECK SUBSCRIPT RANGE DO 820 I = 1, NDATA 820 YA(I,1,NFILE-31) = Y(I+IPLUS) END IF ELSE NBLOCK = (NDATA-1)/LBUFF + 1 IS = IPLUS + 1 IE = LBUFF + IPLUS DO 1130 J = 1, NBLOCK WRITE(NFILE) ( Y(I), I = IS, IE ) IS = IS + LBUFF IE = IE + LBUFF 1130 CONTINUE REWIND NFILE END IF C IF( IY.EQ.0 ) THEN NETALT = NETAL DO 1200 I = 1, NETAL ETALT(I) = ETAL(I) TJET(I) = TJE(I) 1200 CONTINUE END IF C RETURN C END DOUBLE PRECISION FUNCTION FNORMG( ANORTH, HEIGHT ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DATA R1, R2 / 6378137.D0, 6356752.3141D0 / DATA G1, G2 / 978.03267715D0, 983.21863685D0 / C SLATIT = ANORTH * 3.14159 26535 89793 D0 / 180.D0 C FNORMG = ( R1*G1*DCOS(SLATIT)**2 + R2*G2*DSIN(SLATIT)**2 ) * / DSQRT( ( R1*DCOS(SLATIT) )**2 + ( R2*DSIN(SLATIT) )**2 ) * - 0.3086D-3*HEIGHT C RETURN END DOUBLE PRECISION FUNCTION TJUL( JA, JB, JC, HR ) C C JA : LAST TWO DIGITS OF THE YEAR C JB : MONTH C JC : DAY C HR : HOUR C C JULIAN DATE IN UT FROM 1900. JAN. 0.5 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION NM(12) C DATA NM/0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334/ C JAC = MOD( JA, 1900 ) ILL = 0 IF( JAC.LT.1 .OR. JAC.GT.199) ILL = 1 IF( JB.LT.1 .OR. JB .GT.12 ) ILL = 1 IF( JC.LT.1 .OR. JC .GT.31 ) ILL = 1 CCC IF( HR.LT.0.D0 .OR. HR.GT.24.D0 ) ILL = 1 IF( ILL.EQ.0 ) GO TO 10 WRITE(6,200) JA, JB, JC, HR 200 FORMAT(1H /' INVALID DATE Y,M,D,H =', 3I6, F6.1, * ' STOPPED AT TJUL') STOP C 10 IE = (JAC-1)/4 IF( MOD(JAC,4).EQ.0 .AND. JB.GE.3 ) IE = IE+1 JT = 365*JAC + IE + NM(JB) + JC TJ = DBLE(JT) TJUL = TJ + HR/24.D0 - 0.5D0 C RETURN END SUBROUTINE DATAR(RAT,IRAT,Y,IY,NDATA,RLIM,LA,LB,LC,HLR,DELT, * NFILE,MT,DIF, ZERORM) C C THIS SUBROUTINE READS DATA Y FROM MT. C INPUTS: C IY: 0:TIDAL DATA; 1:ASSOCIATED DATA C MT: INPUT DEVICE SPECIFICATION C OUTPUTS: C Y: DATA VECTOR C NDATA: DATA LENGTH C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( LBUFF = 4500 ) PARAMETER ( MAXASC = 3 ) C DIMENSION Y(*), RAT(*) DIMENSION YO(LBUFF,2), YA(LBUFF,2,MAXASC) REAL YT(LBUFF,2,MAXGR2) INTEGER ZERORM C CHARACTER*80 FMT C COMMON /LPCNT/ LPOUT COMMON /BUFOTA/ YO, YT, YA, INCORE C DATA EPS7 /1.D-7/ C YMDH = LA*10000 + LB*100 + LC + HLR/24.D0 READ(MT,502) JA, JB, JC, HR0 502 FORMAT(3I5,F10.0) JA = MOD( JA, 1900 ) YMDH0 = JA*10000 + JB*100 + JC + HR0/24.D0 READ(MT,501) NDATA0, DELT0, RLIM 501 FORMAT(I5,F5.0,D15.5) KAN = IDINT( DELT/DELT0 + EPS7 ) GOSA = DABS( DELT - DELT0*DBLE(KAN) ) IF( GOSA.GT.EPS7 ) THEN WRITE(6,100) DELT, DELT0, KAN STOP END IF 100 FORMAT(1H /' INVALID SAMPLING INTERVAL STOPPED AT DATAR'/ * 1H ,'DELT, DELT0, KAN =', F10.2, F10.2, I8 ) C C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) DO 512 I = 1, 80 512 FMT(I:I) = ' ' READ(MT,'(72A1)',ERR=518) (FMT(I:I),I=1,72) 518 READ(MT,FMT) (Y(I),I=1,NDATA0) IS = 0 IF( YMDH .LE. YMDH0+EPS7 ) THEN LA = JA LB = JB LC = JC HLR = HR0 IF( YMDH .LT. YMDH0-EPS7 ) WRITE(6,102) LA, LB, LC, HLR ELSE 11 IS = IS + 1 CALL UPDATE(JA,JB,JC,HR0,DELT0) YMDH0 = JA*10000 + JB*100 + JC + HR0/24.D0 IF( YMDH .GT. YMDH0+EPS7 ) GO TO 11 IF( YMDH .LT. YMDH0-EPS7 ) THEN LA = JA LB = JB LC = JC HLR = HR0 WRITE(6,102) LA, LB, LC, HLR END IF END IF 102 FORMAT(1H /' START DATE WAS CHANGED Y,M,D,H = ',3I5,F6.1/) IF( KAN*(NDATA-1)+1 .GT. NDATA0-IS ) NDATA = (NDATA0-IS)/KAN C C RESAMPLING JJ = IS + 1 Y(1) = Y(JJ) IF( KAN.EQ.1 .OR. RLIM.LE.0.D0 ) THEN DO 112 J = 2, NDATA JJ = JJ + KAN 112 Y(J) = Y(JJ) ELSE DO 115 J = 2, NDATA YMIN = 0.D0 DO 114 KK = JJ+1, JJ+KAN-1 114 YMIN = DMIN1( YMIN, Y(KK) ) JJ = JJ + KAN IF( YMIN.LT.-RLIM ) THEN Y(J) = YMIN ELSE Y(J) = Y(JJ) END IF 115 CONTINUE END IF C RNUM = 0.D0 YSUM = 0.D0 JCHK = 0 YMAX = -1.D25 YMIN = 1.D25 DO 14 I = 1, NDATA IF(RLIM.GT.0.D0 .AND. Y(I).LT.-RLIM .AND. JCHK.EQ.1) Y(I) =-1.D30 IF(RLIM.GT.0.D0 .AND. Y(I).GT. RLIM) Y(I) = 1.D30 IF( ZERORM.EQ.1 .AND. DABS(Y(I)).LT.1.0D-14 ) Y(I) = 1.D30 IF( DABS(Y(I)) .LE. 1.D25 ) THEN JCHK = 1 YMAX = DMAX1( YMAX, Y(I) ) YMIN = DMIN1( YMIN, Y(I) ) YSUM = YSUM + Y(I) RNUM = RNUM + 1.D0 END IF 14 CONTINUE C DIF0 = YMAX - YMIN IF( IY.EQ.0 ) DIF = DIF0 RATIO = DIF/DIF0 RAT(IRAT) = RATIO IF( IY.NE.0 ) THEN YMEAN = YSUM/RNUM DO 20 I = 1, NDATA IF( DABS(Y(I)) .LT. 1.D25 ) Y(I) = (Y(I)-YMEAN)*RATIO 20 CONTINUE END IF C IF( IY.EQ.0 ) WRITE(6,6072) MT, RATIO IF( IY.NE.0 ) WRITE(6,6074) MT, RATIO, YMEAN 6072 FORMAT(1H /' TIDAL DATA I/O UNIT =', I3, * ' DATA RANGE RATIO =',1PD12.4) 6074 FORMAT(1H /' ASSOCIATED DATA I/O UNIT =', I3, * ' DATA RANGE RATIO =',1PD12.4, * ' MEAN VALUE =',D12.4) C NFILE = NFILE + 1 C IF( INCORE.EQ.1 ) THEN IF(IY.EQ.0) THEN CCC DO NOT CHECK SUBSCRIPT RANGE DO 810 I = 1, NDATA 810 YO(I,1) = Y(I) ELSE CCC DO NOT CHECK SUBSCRIPT RANGE DO 820 I = 1, NDATA 820 YA(I,1,NFILE-31) = Y(I) END IF ELSE NBLOCK = (NDATA-1)/LBUFF + 1 IS = 1 IE = LBUFF DO 900 J = 1, NBLOCK WRITE(NFILE) ( Y(I), I = IS, IE ) IS = IS + LBUFF IE = IE + LBUFF 900 CONTINUE REWIND NFILE END IF C RETURN C END SUBROUTINE UPDATE( JA, JB, JC, HR, DELT ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION MONTH(12) DATA MONTH / 31,28,31,30,31,30,31,31,30,31,30,31 / DATA EPS9 / 1.D-9 / C HR = HR + DELT IF( HR .LT. 24.D0-EPS9 ) RETURN C 1 HR = HR - 24.D0 IF( DABS(HR) .LT. EPS9 ) HR = 0.D0 JC = JC + 1 IF( HR .GE. 24.D0 ) GO TO 1 C 2 M = MONTH(JB) IF( JB .EQ. 2 .AND. MOD(JA,4) .EQ. 0 ) M = M + 1 IF( JC .LE. M ) RETURN C JC = JC - M JB = JB + 1 IF( JB .LE. 12 ) GO TO 2 C JB = 1 JA = JA + 1 IF( JA .GE. 100 ) JA = MOD( JA, 100 ) GO TO 2 C END SUBROUTINE GEOGR1( IC, PH, ATT, GE, AZ, LOVENM ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DOUBLE PRECISION K2, K3 DOUBLE PRECISION L2, L3 COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /PHASE0/ PHCORR(4,4) C DATA ER / 6378160.0D0 /, EC/3.35289D-03/ DATA K2, K3 / 0.299D0, 0.093D0 / DATA H2, H3 / 0.606D0, 0.291D0 / DATA L2, L3 / 0.0840D0, 0.0149D0 / DATA RAD/0.017453292519943D0/ C C UNIT GRAVITY : MICROGAL TILT : MSECARC STRAIN : 10**(-9) C HIGHT : M GE : CM/(SEC*SEC) C FAI = PH IF( FAI.GE. 90.D0 ) FAI = 89.9999D0 IF( FAI.LE.-90.D0 ) FAI = -89.9999D0 IF( FAI.EQ. 0.D0 ) FAI = 0.0001D0 C FG = FAI*RAD FF = (1.D0-EC)**2*DTAN(FG) FC = DATAN(FF) EP = FG-FC CE = DCOS(EP) SE = DSIN(EP) RC = 1.D0-EC*DSIN(FG)**2+ATT/ER DO 10 J = 1, 3 DO 10 I = 1, 3 GC(I,J) = 0.D0 PHCORR(I,J) = 0.D0 10 CONTINUE CF = DCOS(FC) SF = DSIN(FC) C21 = DSIN(2.D0*FC) C22 = CF*CF C31 = CF*(1.D0-5.D0*SF*SF) C32 = SF*CF*CF C33 = C22*CF S21 = DCOS(2.D0*FC) S31 = SF*(15.D0*C22-4.D0) S32 = CF*(1.D0-3.D0*SF*SF) S33 = C32 C GO TO ( 1, 23, 23, 4, 4, 4, 11, 12 ), IC WRITE(6,200) IC 200 FORMAT(1H /' ERROR AT GEOGR1 INVALID KIND =',I6) STOP C C GRAVITY TIDE 1 C2 = 41.291125D0*RC C3 = C2*RC GC(2,1) = -2.D0*C2*(C21*CE+S21*SE) GC(2,2) = -C2*(2.D0*C22*CE-C21*SE) GC(3,1) = -0.72618D0*C3*(3.D0*C31*CE-S31*SE) GC(3,2) = -2.59808D0*C3*(3.D0*C32*CE+S32*SE) GC(3,3) = -3.D0*C3*(C33*CE-S33*SE) C DFACT2 = 1.D0 + H2 - 1.5D0*K2 DFACT3 = 1.D0 + 2.D0/3.D0*H3 - 4.D0/3.D0*K3 C GO TO 8 C C TILT TIDE (FROM NORTH CLOCKWISELY) 23 C2 = 8.516906D03/GE*RC C3 = C2*RC GC21N = -2.D0*C2*(S21*CE-C21*SE)*(-1.D0) GC22N = C2*(C21*CE+2.D0*C22*SE) *(-1.D0) GC31N = 0.72618D0*C3*(S31*CE-C21*SE)*(-1.0) GC32N = 2.59808D0*C3*(CF*(2.D0-3.D0*C22)*CE+3.D0*C32*SE)*(-1.D0) GC33N = 3.D0*C3*(C32*CE+C33*SE)*(-1.D0) GC21E = 2.D0*C2*SF GC22E = 2.D0*C2*CF GC31E = 0.72618D0*C3*C31/CF GC32E = 2.59808D0*C3*C21 GC33E = 3.D0*C3*C22 IF( IC.EQ.3 ) AZ = 90.D0 COSAZ = DCOS( AZ*RAD ) SINAZ = DSIN( AZ*RAD ) GC(2,1) = DSQRT( (GC21N*COSAZ)**2 + (GC21E*SINAZ)**2 ) GC(2,2) = DSQRT( (GC22N*COSAZ)**2 + (GC22E*SINAZ)**2 ) GC(3,1) = DSQRT( (GC31N*COSAZ)**2 + (GC31E*SINAZ)**2 ) GC(3,2) = DSQRT( (GC32N*COSAZ)**2 + (GC32E*SINAZ)**2 ) GC(3,3) = DSQRT( (GC33N*COSAZ)**2 + (GC33E*SINAZ)**2 ) PHCORR(2,1) = DATAN2( GC21E*SINAZ, GC21N*COSAZ ) / RAD PHCORR(2,2) = DATAN2( GC22E*SINAZ, GC22N*COSAZ ) / RAD PHCORR(3,1) = DATAN2( GC31E*SINAZ, GC31N*COSAZ ) / RAD PHCORR(3,2) = DATAN2( GC32E*SINAZ, GC32N*COSAZ ) / RAD PHCORR(3,3) = DATAN2( GC33E*SINAZ, GC33N*COSAZ ) / RAD C DFACT2 = 1.D0 + K2 - H2 DFACT3 = 1.D0 + K3 - H3 C GO TO 8 C C STRAIN TIDE WITH LOVE'S NUMBER H AND L C 4 IF( IC.EQ.5 ) AZ = 90.D0 IF( IC.EQ.6 ) AZ = 45.D0 COSAZ = DCOS(AZ*RAD) SINAZ = DSIN(AZ*RAD) COSAZ2 = COSAZ*COSAZ SINAZ2 = SINAZ*SINAZ COSIAZ = COSAZ*SINAZ C THETA = 90.D0*RAD - FC CTH = DCOS(THETA) STH = DSIN(THETA) C G = 4.124492D+04*RC/GE GC21H = G*2.D0*CTH*STH * H2 GC22H = G*STH*STH * H2 GC31H = G*0.72618D0*STH*(1.D0-5.D0*CTH*CTH)*RC * H3 GC32H = G*2.59808D0*CTH*STH*STH*RC * H3 GC33H = G*STH*STH*STH*RC * H3 C GC21N = -G*8.D0*CTH*STH * L2 GC22N = G*2.D0*(CTH*CTH-STH*STH) * L2 GC31N = G*0.72618D0*STH*(34.D0-45.D0*STH*STH)*RC * L3 GC32N = G*2.59808D0*CTH*(9.D0*CTH*CTH-7.D0)*RC * L3 GC33N = G*STH*(6.D0-9.D0*STH*STH)*RC * L3 C GC21E = -G*4.D0*CTH*STH * L2 GC22E = G*(2.D0*CTH*CTH-4.D0) * L2 GC31E = G*0.72618D0*STH*(14.D0-15.D0*STH*STH)*RC * L3 GC32E = G*2.59808D0*CTH*(3.D0*CTH*CTH-5.D0)*RC * L3 GC33E = -G*STH*(6.D0+3.D0*STH*STH)*RC * L3 C GC21X = G*4.D0*STH * L2 GC22X = -G*4.D0*CTH * L2 GC31X = -G*0.72618D0*20.D0*CTH*STH*RC * L3 GC32X = G*2.59808D0*4.D0*(STH*STH-CTH*CTH)*RC * L3 GC33X = -G*12.D0*CTH*STH*RC * L3 C PIN21 = ( GC21N + GC21H )*COSAZ2 + ( GC21E + GC21H )*SINAZ2 PIN22 = ( GC22N + GC22H )*COSAZ2 + ( GC22E + GC22H )*SINAZ2 PIN31 = ( GC31N + GC31H )*COSAZ2 + ( GC31E + GC31H )*SINAZ2 PIN32 = ( GC32N + GC32H )*COSAZ2 + ( GC32E + GC32H )*SINAZ2 PIN33 = ( GC33N + GC33H )*COSAZ2 + ( GC33E + GC33H )*SINAZ2 C POUT21 = - GC21X*COSIAZ POUT22 = - GC22X*COSIAZ POUT31 = - GC31X*COSIAZ POUT32 = - GC32X*COSIAZ POUT33 = - GC33X*COSIAZ C GC(2,1) = DSQRT( PIN21*PIN21 + POUT21*POUT21 ) GC(2,2) = DSQRT( PIN22*PIN22 + POUT22*POUT22 ) GC(3,1) = DSQRT( PIN31*PIN31 + POUT31*POUT31 ) GC(3,2) = DSQRT( PIN32*PIN32 + POUT32*POUT32 ) GC(3,3) = DSQRT( PIN33*PIN33 + POUT33*POUT33 ) C PHCORR(2,1) = - DATAN2( POUT21, PIN21 ) / RAD PHCORR(2,2) = - DATAN2( POUT22, PIN22 ) / RAD PHCORR(3,1) = - DATAN2( POUT31, PIN31 ) / RAD PHCORR(3,2) = - DATAN2( POUT32, PIN32 ) / RAD PHCORR(3,3) = - DATAN2( POUT33, PIN33 ) / RAD C GO TO 8 C C POTENTIAL/GRAVITY/RADIUS (NONE DIMENSIONAL, FOR STRAIN TIDE) 11 G = 4.124492D+04*RC/GE GC(2,1) = G*C21 GC(2,2) = G*C22 GC(3,1) = G*0.72618D0*C31*RC GC(3,2) = G*2.59808D0*C32*RC GC(3,3) = G*C33*RC C DFACT2 = 2.D0*H2 - 6.D0*L2 DFACT3 = 2.D0*H3 - 12.D0*L3 C GO TO 8 C C OCEAN TIDE 12 G = 26.75824D0*RC GC(2,1) = G*C21 GC(2,2) = G*C22 GC(3,1) = G*0.72618D0*C31*RC GC(3,2) = G*2.59808D0*C32*RC GC(3,3) = G*C33*RC C DFACT2 = 1.D0 + K2 - H2 DFACT3 = 1.D0 + K3 - H3 C C 8 CONTINUE IF( IC.LE.3 .OR. IC.GE.7 ) THEN IF( LOVENM.EQ.1 ) THEN GC(2,1) = GC(2,1)*DFACT2 GC(2,2) = GC(2,2)*DFACT2 GC(3,1) = GC(3,1)*DFACT3 GC(3,2) = GC(3,2)*DFACT3 GC(3,3) = GC(3,3)*DFACT3 WRITE(6,220) DFACT2 WRITE(6,230) DFACT3 ELSE IF( LOVENM.EQ.2 ) THEN DFACTA = DFACT3/DFACT2 GC(3,1) = GC(3,1)*DFACTA GC(3,2) = GC(3,2)*DFACTA WRITE(6,250) DFACTA END IF END IF 220 FORMAT(//' DELTA FACTOR FOR P21, P22 :', F9.5 ) 230 FORMAT( ' DELTA FACTOR FOR P31, P32, P33 :', F9.5 ) 250 FORMAT(//' NORMALIZATION FACTOR FOR P31, P32 :', F9.5 ) C RETURN C END SUBROUTINE SPHAS1( TJ, HD, ALON ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION T(6) COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /PHASE0/ PHCORR(4,4) COMMON /JSTUT/ TIMSYS DATA DETMUT /51.D0/ C TD = TJ + TIMSYS/24.D0 + DETMUT/86400.D0 FS = TD/36525.0D0 FS2 = FS*FS FS3 = FS2*FS T(2) = 270.434162D0 + 481267.883142D0*FS - 0.0011333D0*FS2 & + 0.000001889D0*FS3 T(3) = 279.696681D0 + 36000.768925D0*FS + 0.0003027D0*FS2 T(4) = 334.329568D0 + 4069.034034D0*FS - 0.0103249D0*FS2 & - 0.0000125D0*FS3 T(5) = 100.843202D0 + 1934.142008D0*FS - 0.002078 D0*FS2 & - 0.000002 D0*FS3 T(6) = 281.220868D0 + 1.719175D0*FS + 0.0004527D0*FS2 & + 0.0000033D0*FS3 C ALON:LONGITUDE (EAST) T(1) = T(3) - T(2) + (HD+TIMSYS)*15.D0 + ALON DO 5 I = 1, 377 R = 0.D0 DO 1 J = 1, 6 F = IDN(I,J) R = R + T(J)*F 1 CONTINUE IS = IDN(I,1) NP = IP(I) R = R + PHCORR( NP, IS ) IF( MOD(NP+IS,2).EQ.1 ) R = R - 90.D0 K = R/360.D0 U = K PHC(I) = R - 360.D0*U 5 CONTINUE C RETURN END SUBROUTINE GEOGR2( KIND, PH, HI, GRAV, AZ, LOVENM ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DOUBLE PRECISION K2, K3, K4 DOUBLE PRECISION L2, L3, L4 COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /PHASE0/ PHCORR(4,4) C C DATA PI / 3.14159 26535 89793 D0/ DATA RE / 6378137.D2 / DATA FLAT / 3.35292D-3 / DATA GE / 3.98600448D20 / DATA FMU / 0.012300034D0 / DATA SINPI / 3422.448D0 / DATA K2, K3, K4 / 0.299D0, 0.093D0, 0.042D0 / DATA H2, H3, H4 / 0.606D0, 0.291D0, 0.178D0 / CCC DATA L2, L3, L4 / 0.0840D0, 0.0149D0, 0.0050D0 / DATA L2, L3, L4 / 0.0840D0, 0.0149D0, 0.0100D0 / C C DEG = 180.D0/PI RAD = PI/180.D0 C C CV = 3.D0 / 4.D0 * FMU * GE * ( SINPI/3600.D0*RAD )**3 / RE**2 C PHI = PH IF( PHI.GE. 90.D0 ) PHI = 89.9999D0 IF( PHI.LE.-90.D0 ) PHI = -89.9999D0 IF( PHI.EQ. 0.D0 ) PHI = 0.0001D0 PHI = PHI * RAD C ALTI = HI * 100.D0 C C COSPHI = DCOS( PHI ) SINPHI = DSIN( PHI ) EP = ( FLAT + FLAT**2/2.D0 ) * DSIN( 2.D0*PHI ) * - FLAT**2/2.D0 * DSIN( 4.D0*PHI ) E2 = 2.D0*FLAT - FLAT**2 ROE = ( 1.D0 / DSQRT( 1.D0-E2*SINPHI**2 ) + ALTI/RE ) * * COSPHI ROP = ( ( 1.D0-E2 ) / DSQRT( 1.D0-E2*SINPHI**2 ) + ALTI/RE ) * * SINPHI ROE2 = ROE**2 ROP2 = ROP**2 RR = RE * DSQRT( ROE2 + ROP2 ) C CCC DO 10 J = 0, 4 DO 10 J = 1, 4 DO 10 I = 1, 4 GC(I,J) = 0.D0 PHCORR(I,J) = 0.D0 10 CONTINUE C IF( KIND.EQ.1 ) THEN C C GRAVITY IN MICRO GAL CV = -1.D6*CV CCC GC(2,0) = CV * ( ROE*COSPHI - 2.D0*ROP*SINPHI ) GC(2,1) = 2.D0 * CV * ( ROP*COSPHI + ROE*SINPHI ) GC(2,2) = 2.D0 * CV * ROE * COSPHI CCC GC(3,0) = 3.D0*DSQRT(5.D0)/2.D0 * CV * ( 2.D0*ROE*ROP CCC * * COSPHI + (ROE2-2.D0*ROP2)*SINPHI ) GC(3,1) = 3.D0*DSQRT(15.D0)/16.D0 * CV * ( (3.D0*ROE2 * -4.D0*ROP2)*COSPHI - 8.D0*ROE*ROP*SINPHI ) GC(3,2) = 3.D0*DSQRT(3.D0)/2.D0 * CV * ROE * * ( 2.D0*ROP*COSPHI + ROE*SINPHI ) GC(3,3) = 3.D0 * CV * ROE2 * COSPHI CCC GC(4,0) = 0.5D0 * CV CCC * * ( 3.D0*ROE*( ROE2-4.D0*ROP2 ) *COSPHI CCC * + 4.D0*ROP*( 2.D0*ROP2-3.D0*ROE2 )*SINPHI ) GC(4,1) = 448.D0 / ( 3.D0+DSQRT(393.D0)) * / DSQRT( 390.D0+2.D0*DSQRT(393.D0) ) * CV * * ( ROP*( 9.D0*ROE2-4.D0*ROP2 )*COSPHI * + 3.D0*ROE*( ROE2-4.D0*ROP2 )*SINPHI ) GC(4,2) = 28.D0/9.D0 * CV * * ( ROE*( ROE2-3.D0*ROP2 )*COSPHI * - 3.D0*ROE2*ROP *SINPHI ) GC(4,3) = 16.D0*DSQRT(3.D0)/9.D0 * CV * * ( 3.D0*ROE2*ROP*COSPHI + ROE**3*SINPHI ) CCC GC(4,4) = 4.D0 * CV * ROE**3 * COSPHI C DFACT2 = 1.D0 + H2 - 1.50D0*K2 DFACT3 = 1.D0 + 2.D0/3.D0*H3 - 4.D0/3.D0*K3 DFACT4 = 1.D0 + 0.5D0*H4 - 1.25D0*K4 C ELSE IF( KIND.EQ.2 .OR. KIND.EQ.3 ) THEN C C TILT IN MILLI SECOND OF ARC FACT = DEG * 3600.D3 / GRAV CV = CV*FACT CCC XX20 = -CV * ( 2.D0*ROP*COSPHI + ROE*SINPHI ) XX21 = 2.D0 * CV * ( ROE*COSPHI - ROP*SINPHI ) XX22 = -2.D0 * CV * ROE * SINPHI CCC XX30 = 3.D0*DSQRT(5.D0)/2.D0 * CV * ( (ROE2 CCC * -2.D0*ROP2)*COSPHI - 2.D0*ROE*ROP*SINPHI ) XX31 = 3.D0*DSQRT(15.D0)/16.D0 * CV * ( -8.D0*ROE*ROP * *COSPHI + ( 4.D0*ROP2-3.D0*ROE2 )*SINPHI ) XX32 = 3.D0*DSQRT(3.D0)/2.D0 * CV * *( ROE2*COSPHI - 2.D0*ROE*ROP*SINPHI ) XX33 = -3.D0 * CV * ROE2 * SINPHI CCC XX40 = 0.5D0 * CV CCC * * ( 4.D0*ROP*( 2.D0*ROP2-3.D0*ROE2 )*COSPHI CCC * + 3.D0*ROE*( 4.D0*ROP2-ROE2 ) *SINPHI ) XX41 = 448.D0 / ( 3.D0+DSQRT(393.D0)) * / DSQRT( 390.D0+2.D0*DSQRT(393.D0) ) * CV * * ( 3.D0*ROE*( ROE2-4.D0*ROP2 )*COSPHI * + ROP*( 4.D0*ROP2-9.D0*ROE2 )*SINPHI ) XX42 = 28.D0/9.D0 * CV * * ( -3.D0*ROE2*ROP *COSPHI * + ROE*(3.D0*ROP2-ROE2)*SINPHI ) XX43 = 16.D0*DSQRT(3.D0)/9.D0 * CV * * ( ROE**3*COSPHI - 3.D0*ROE2*ROP*SINPHI ) CCC XX44 = -4.D0 * CV * ROE**3 * SINPHI C CCC YY20 = 0.D0 YY21 = 2.D0 * CV * ROP YY22 = 2.D0 * CV * ROE CCC YY30 = 0.D0 YY31 = 3.D0*DSQRT(15.D0)/16.D0 * CV * * ( ROE2 - 4.D0*ROP2 ) YY32 = 3.D0*DSQRT(3.D0) * CV * ROE * ROP YY33 = 3.D0 * CV * ROE2 CCC YY40 = 0.D0 YY41 = 448.D0 / ( 3.D0+DSQRT(393.D0) ) * / DSQRT( 390.D0+2.D0*DSQRT(393.D0) ) * CV * * ROP*(3.D0*ROE2-4.D0*ROP2) YY42 = 14.D0/9.D0 * CV * * ROE*(ROE2-6.D0*ROP2) YY43 = 16.D0*DSQRT(3.D0)/3.D0 * CV * ROE2 * ROP CCC YY44 = 4.D0 * CV * ROE**3 C IF( KIND.EQ.3 ) AZ = 90.D0 COSAZ = DCOS( AZ*RAD ) SINAZ = DSIN( AZ*RAD ) C CCC GC(2,0) = DSQRT( (XX20*COSAZ)**2 + (YY20*SINAZ)**2 ) GC(2,1) = DSQRT( (XX21*COSAZ)**2 + (YY21*SINAZ)**2 ) GC(2,2) = DSQRT( (XX22*COSAZ)**2 + (YY22*SINAZ)**2 ) CCC GC(3,0) = DSQRT( (XX30*COSAZ)**2 + (YY30*SINAZ)**2 ) GC(3,1) = DSQRT( (XX31*COSAZ)**2 + (YY31*SINAZ)**2 ) GC(3,2) = DSQRT( (XX32*COSAZ)**2 + (YY32*SINAZ)**2 ) GC(3,3) = DSQRT( (XX33*COSAZ)**2 + (YY33*SINAZ)**2 ) CCC GC(4,0) = DSQRT( (XX40*COSAZ)**2 + (YY40*SINAZ)**2 ) GC(4,1) = DSQRT( (XX41*COSAZ)**2 + (YY41*SINAZ)**2 ) GC(4,2) = DSQRT( (XX42*COSAZ)**2 + (YY42*SINAZ)**2 ) GC(4,3) = DSQRT( (XX43*COSAZ)**2 + (YY43*SINAZ)**2 ) CCC GC(4,4) = DSQRT( (XX44*COSAZ)**2 + (YY44*SINAZ)**2 ) C CCC IF( COSAZ.NE.0.D0 ) CCC * PHCORR(2,0) = DATAN2( YY20*SINAZ, XX20*COSAZ ) * DEG PHCORR(2,1) = DATAN2( YY21*SINAZ, XX21*COSAZ ) * DEG PHCORR(2,2) = DATAN2( YY22*SINAZ, XX22*COSAZ ) * DEG CCC IF( COSAZ.NE.0.D0 ) CCC * PHCORR(3,0) = DATAN2( YY30*SINAZ, XX30*COSAZ ) * DEG PHCORR(3,1) = DATAN2( YY31*SINAZ, XX31*COSAZ ) * DEG PHCORR(3,2) = DATAN2( YY32*SINAZ, XX32*COSAZ ) * DEG PHCORR(3,3) = DATAN2( YY33*SINAZ, XX33*COSAZ ) * DEG CCC IF( COSAZ.NE.0.D0 ) CCC * PHCORR(4,0) = DATAN2( YY40*SINAZ, XX40*COSAZ ) * DEG PHCORR(4,1) = DATAN2( YY41*SINAZ, XX41*COSAZ ) * DEG PHCORR(4,2) = DATAN2( YY42*SINAZ, XX42*COSAZ ) * DEG PHCORR(4,3) = DATAN2( YY43*SINAZ, XX43*COSAZ ) * DEG CCC PHCORR(4,4) = DATAN2( YY44*SINAZ, XX44*COSAZ ) * DEG C DFACT2 = 1.D0 + K2 - H2 DFACT3 = 1.D0 + K3 - H3 DFACT4 = 1.D0 + K4 - H4 C ELSE IF( KIND.EQ.4 .OR. KIND.EQ.5 .OR. KIND.EQ.6 ) THEN C C STRAIN TIDE WITH LOVE'S NUMBER H AND L C IF( KIND.EQ.5 ) AZ = 90.D0 IF( KIND.EQ.6 ) AZ = 45.D0 COSAZ = DCOS(AZ*RAD) SINAZ = DSIN(AZ*RAD) COSAZ2 = COSAZ*COSAZ SINAZ2 = SINAZ*SINAZ COSIAZ = COSAZ*SINAZ C PSI = PHI - EP THETA = 90.D0*RAD - PSI CTH = DCOS(THETA) CTH2 = CTH*CTH STH = DSIN(THETA) STH2 = STH*STH C G = CV * 1.D9/GRAV/6371023.6D2 * RR**2 / RE RC = RR/RE RC4 = RC*RC CCC GC20H = G*0.5D0*(1.D0-3.D0*CTH2) * H2 GC21H = G*2.D0*CTH*STH * H2 GC22H = G*STH2 * H2 CCC GC30H = G*1.11803D0*CTH*(3.D0-5.D0*CTH2) * RC * H3 GC31H = G*0.72618D0*STH*(1.D0-5.D0*CTH2) * RC * H3 GC32H = G*2.59808D0*CTH*STH2 * RC * H3 GC33H = G*STH2*STH * RC * H3 CCC GC40H = G*0.12500D0*(3.D0-30.D0*CTH2+35.D0*CTH2*CTH2) CCC * * RC4 * H4 GC41H = G*0.47347D0*CTH*STH*(6.D0-14.D0*CTH2)* RC4 * H4 GC42H = G*0.77778D0*(1.D0-8.D0*CTH2+7.D0*CTH2*CTH2) * * RC4 * H4 GC43H = G*3.07920D0*CTH*STH2*STH * RC4 * H4 CCC GC44H = G*STH2*STH2 * RC4 * H4 C CCC GC20N = G*3.D0*(CTH2-STH2) * L2 GC21N = -G*8.D0*CTH*STH * L2 GC22N = G*2.D0*(CTH2-STH2) * L2 CCC GC30N = G*1.11803D0*CTH*(45.D0*CTH2-33.D0) * RC * L3 GC31N = G*0.72618D0*STH*(34.D0-45.D0*STH2) * RC * L3 GC32N = G*2.59808D0*CTH*(9.D0*CTH2-7.D0) * RC * L3 GC33N = G*STH*(6.D0-9.D0*STH2) * RC * L3 CCC GC40N = -G*0.12500D0*(60.D0-540.D0*CTH2+560.D0*CTH2*CTH2) CCC * * RC4 * L4 GC41N = G*0.47347D0*CTH*STH*(224.D0*CTH2-108.D0) * * RC4 * L4 GC42N = -G*0.77778D0*(16.D0-116.D0*CTH2+112.D0*CTH2*CTH2) * * RC4 * L4 GC43N = G*3.07920D0*CTH*STH*(6.D0-16.D0*STH2)* RC4 * L4 CCC GC44N = G*STH2*(12.D0-16.D0*STH2) * RC4 * L4 C CCC GC20E = G*3.D0*CTH2 * L2 GC21E = -G*4.D0*CTH*STH * L2 GC22E = G*(2.D0*CTH2-4.D0) * L2 CCC GC30E = G*1.11803D0*3.D0*(5.D0*CTH2-1.D0) * RC * L3 GC31E = G*0.72618D0*STH*(14.D0-15.D0*STH2) * RC * L3 GC32E = G*2.59808D0*CTH*(3.D0*CTH2-5.D0) * RC * L3 GC33E = -G*STH*(6.D0+3.D0*STH2) * RC * L3 CCC GC40E = G*0.12500D0*CTH2*(60.D0-140.D0*CTH2) * RC4 * L4 GC41E = G*0.47347D0*CTH*STH*(56.D0*CTH2-12.D0)*RC4 * L4 GC42E = -G*0.77778D0*(4.D0-44.D0*CTH2+28.D0*CTH2*CTH2) * * RC4 * L4 GC43E = -G*3.07920D0*CTH*STH*(6.D0+4.D0*STH2) * RC4 * L4 CCC GC44E = -G*STH2*(12.D0+4.D0*STH2) * RC4 * L4 C CCC GC20X = 0.D0 GC21X = G*4.D0*STH * L2 GC22X = -G*4.D0*CTH * L2 CCC GC30X = 0.D0 GC31X = -G*0.72618D0*20.D0*CTH*STH * RC * L3 GC32X = G*2.59808D0*4.D0*(STH2-CTH2) * RC * L3 GC33X = -G*12.D0*CTH*STH * RC * L3 CCC GC40X = 0.D0 GC41X = G*0.47347D0*STH*(84.D0*STH2-72.D0) * RC4 * L4 GC42X = G*0.77778D0*CTH*(84.D0*CTH2-60.D0) * RC4 * L4 GC43X = G*3.07920D0*STH*(18.D0*STH2-12.D0) * RC4 * L4 CCC GC44X = G*CTH*(24.D0*CTH2-24.D0) * RC4 * L4 C CCC PIN20 = ( GC20N + GC20H )*COSAZ2 + ( GC20E + GC20H )*SINAZ2 PIN21 = ( GC21N + GC21H )*COSAZ2 + ( GC21E + GC21H )*SINAZ2 PIN22 = ( GC22N + GC22H )*COSAZ2 + ( GC22E + GC22H )*SINAZ2 CCC PIN30 = ( GC30N + GC30H )*COSAZ2 + ( GC30E + GC30H )*SINAZ2 PIN31 = ( GC31N + GC31H )*COSAZ2 + ( GC31E + GC31H )*SINAZ2 PIN32 = ( GC32N + GC32H )*COSAZ2 + ( GC32E + GC32H )*SINAZ2 PIN33 = ( GC33N + GC33H )*COSAZ2 + ( GC33E + GC33H )*SINAZ2 CCC PIN40 = ( GC40N + GC40H )*COSAZ2 + ( GC40E + GC40H )*SINAZ2 PIN41 = ( GC41N + GC41H )*COSAZ2 + ( GC41E + GC41H )*SINAZ2 PIN42 = ( GC42N + GC42H )*COSAZ2 + ( GC42E + GC42H )*SINAZ2 PIN43 = ( GC43N + GC43H )*COSAZ2 + ( GC43E + GC43H )*SINAZ2 CCC PIN44 = ( GC44N + GC44H )*COSAZ2 + ( GC44E + GC44H )*SINAZ2 C POUT21 = - GC21X*COSIAZ POUT22 = - GC22X*COSIAZ POUT31 = - GC31X*COSIAZ POUT32 = - GC32X*COSIAZ POUT33 = - GC33X*COSIAZ POUT41 = - GC41X*COSIAZ POUT42 = - GC42X*COSIAZ POUT43 = - GC43X*COSIAZ CCC POUT44 = - GC44X*COSIAZ C CCC GC(2,0) = DSQRT( PIN20*PIN20 ) GC(2,1) = DSQRT( PIN21*PIN21 + POUT21*POUT21 ) GC(2,2) = DSQRT( PIN22*PIN22 + POUT22*POUT22 ) CCC GC(3,0) = DSQRT( PIN30*PIN30 ) GC(3,1) = DSQRT( PIN31*PIN31 + POUT31*POUT31 ) GC(3,2) = DSQRT( PIN32*PIN32 + POUT32*POUT32 ) GC(3,3) = DSQRT( PIN33*PIN33 + POUT33*POUT33 ) CCC GC(4,0) = DSQRT( PIN40*PIN40 ) GC(4,1) = DSQRT( PIN41*PIN41 + POUT41*POUT41 ) GC(4,2) = DSQRT( PIN42*PIN42 + POUT42*POUT42 ) GC(4,3) = DSQRT( PIN43*PIN43 + POUT43*POUT43 ) CCC GC(4,4) = DSQRT( PIN44*PIN44 + POUT44*POUT44 ) C CCC PHCORR(2,0) = - DATAN2( 0.0D0, PIN20 ) * DEG PHCORR(2,1) = - DATAN2( POUT21, PIN21 ) * DEG PHCORR(2,2) = - DATAN2( POUT22, PIN22 ) * DEG CCC PHCORR(3,0) = - DATAN2( 0.0D0, PIN30 ) * DEG PHCORR(3,1) = - DATAN2( POUT31, PIN31 ) * DEG PHCORR(3,2) = - DATAN2( POUT32, PIN32 ) * DEG PHCORR(3,3) = - DATAN2( POUT33, PIN33 ) * DEG CCC PHCORR(4,0) = - DATAN2( 0.0D0, PIN40 ) * DEG PHCORR(4,1) = - DATAN2( POUT41, PIN41 ) * DEG PHCORR(4,2) = - DATAN2( POUT42, PIN42 ) * DEG PHCORR(4,3) = - DATAN2( POUT43, PIN43 ) * DEG CCC PHCORR(4,4) = - DATAN2( POUT44, PIN44 ) * DEG C ELSE IF( KIND.EQ.7 .OR. KIND.EQ.8 ) THEN C C KIND = 7 : POTENTIAL/GRAV/RADIUS 10**(-9) STRAIN C KIND = 8 : POTENTIAL/GRAV OCEAN TIDE (CM) C CCC IF( KIND.EQ.7 ) FACT = 1.D9/GRAV/RR IF( KIND.EQ.7 ) FACT = 1.D9/GRAV/6371023.6D2 IF( KIND.EQ.8 ) FACT = 1.D0/GRAV CV = CV*FACT C PSI = PHI - EP COSPSI = DCOS( PSI ) SINPSI = DSIN( PSI ) CVR2 = CV * RR**2 / RE CCC GC(2,0) = CVR2 * ( 1.D0 - 3.D0*SINPSI**2 )/2.D0 GC(2,1) = CVR2 * DSIN( 2.D0*PSI ) GC(2,2) = CVR2 * COSPSI**2 CVR3 = CVR2 * RR/RE CCC GC(3,0) = DSQRT(5.D0)/2.D0 * CVR3 CCC * * SINPSI * ( 3.D0 - 5.D0*SINPSI**2 ) GC(3,1) = 3.D0*DSQRT(15.D0)/16.D0 * CVR3 * * COSPSI * ( 1.D0 - 5.D0*SINPSI**2 ) GC(3,2) = 3.D0*DSQRT(3.D0)/2.D0 * CVR3 * * SINPSI * COSPSI**2 GC(3,3) = CVR3 * COSPSI**3 CVR4 = CVR3 * RR/RE CCC GC(4,0) = 0.125D0 * CVR4 CCC * * ( 3.D0 - 30.D0*SINPSI**2 + 35.D0*SINPSI**4 ) GC(4,1) = 224.D0 / ( 3.D0+DSQRT(393.D0) ) * / DSQRT( 390.D0+2.D0*DSQRT(393.D0) ) * CVR4 * * DSIN( 2.D0*PSI )*( 3.D0 - 7.D0*SINPSI**2 ) GC(4,2) = 7.D0/9.D0 * CVR4 * * COSPSI**2*( 1.D0 - 7.D0*SINPSI**2 ) GC(4,3) = 16.D0/3.D0/DSQRT(3.D0) * CVR4 * * SINPSI*COSPSI**3 CCC GC(4,4) = CVR4 * COSPSI**4 C IF( KIND.EQ.7 ) THEN DFACT2 = 2.D0*H2 - 6.D0*L2 DFACT3 = 2.D0*H3 - 12.D0*L3 DFACT4 = 2.D0*H4 - 20.D0*L4 ELSE IF( KIND.EQ.8 ) THEN DFACT2 = 1.D0 + K2 - H2 DFACT3 = 1.D0 + K3 - H3 DFACT4 = 1.D0 + K4 - H4 END IF C ELSE C WRITE(6,200) KIND 200 FORMAT(1H /' INVALID KIND IN GEOGR2 KIND =',I5) STOP C END IF C IF( KIND.LE.3 .OR. KIND.GE.7 ) THEN IF( LOVENM.EQ.1 ) THEN CCC GC(2,0) = GC(2,0)*DFACT2 GC(2,1) = GC(2,1)*DFACT2 GC(2,2) = GC(2,2)*DFACT2 CCC GC(3,0) = GC(3,0)*DFACT3 GC(3,1) = GC(3,1)*DFACT3 GC(3,2) = GC(3,2)*DFACT3 GC(3,3) = GC(3,3)*DFACT3 CCC GC(4,0) = GC(4,0)*DFACT4 GC(4,1) = GC(4,1)*DFACT4 GC(4,2) = GC(4,2)*DFACT4 GC(4,3) = GC(4,3)*DFACT4 CCC GC(4,4) = GC(4,4)*DFACT4 WRITE(6,220) DFACT2 WRITE(6,230) DFACT3 WRITE(6,240) DFACT4 ELSE IF( LOVENM.EQ.2 ) THEN DFACTA = DFACT3/DFACT2 DFACTB = DFACT4/DFACT2 DFACTC = DFACT4/DFACT3 CCC GC(3,0) = GC(3,0)*DFACTA GC(3,1) = GC(3,1)*DFACTA GC(3,2) = GC(3,2)*DFACTA CCC GC(4,0) = GC(4,0)*DFACTB GC(4,1) = GC(4,1)*DFACTB GC(4,2) = GC(4,2)*DFACTB GC(4,3) = GC(4,3)*DFACTC WRITE(6,250) DFACTA WRITE(6,260) DFACTB WRITE(6,270) DFACTC END IF END IF 220 FORMAT(//' DELTA FACTOR FOR P20, P21, P22 :', F9.5 ) 230 FORMAT( ' DELTA FACTOR FOR P30, P31, P32, P33 :', F9.5 ) 240 FORMAT( ' DELTA FACTOR FOR P40, P41, P42, P43, P44 :', F9.5 ) 250 FORMAT(//' NORMALIZATION FACTOR FOR P30, P31, P32 :', F9.5 ) 260 FORMAT( ' NORMALIZATION FACTOR FOR P40, P41, P42 :', F9.5 ) 270 FORMAT( ' NORMALIZATION FACTOR FOR P43 :', F9.5 ) C RETURN C END SUBROUTINE SPHAS2( TJ, HD, ALON, TC0 ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION T(8) C COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /PHASE0/ PHCORR(4,4) COMMON /JSTUT/ TIMSYS COMMON /DTDMUT/ TDMUT(130) C DATA PI / 3.14159 26535 89793 D0 / C RAD = PI/180.D0 C TMID = TC0/36525.D0 - 1.D0 C TU19 = TJ + TIMSYS/24.D0 TU = TU19/36525.D0 - 1.0D0 TU2 = TU*TU C MM = ( TMID*100.0D0 + 100.0D0 ) IF( MM.LT. 1 ) MM = 1 IF( MM.GT.130 ) MM = 130 DETMUT = TDMUT( MM ) CCC DETMUT = 51.D0 C TD19 = TU19 + DETMUT/86400.D0 TD = TD19/36525.D0 - 1.0D0 TD2 = TD*TD C EL1 = 280.4606184D0 + 36000.7700536D0*TU * + 0.00038793D0*TU2 * - 0.0000000258D0*TU**3 T(2) = 218.316656D0 + 481267.881342D0*TD - 0.001330D0 *TD2 * + 0.0040D0*DCOS( ( 29.D0+133.D0*TMID)*RAD ) T(1) = EL1 - T(2) + (HD+TIMSYS)*15.D0 + ALON T(3) = 280.466449D0 + 36000.769822D0*TD + 0.0003036D0*TD2 * + 0.0018D0*DCOS( (159.D0+ 19.D0*TMID)*RAD ) T(4) = 83.353243D0 + 4069.013711D0*TD - 0.010324D0 *TD2 T(5) = 234.955444D0 + 1934.136185D0*TD - 0.002076D0 *TD2 T(6) = 282.937348D0 + 1.719533D0*TD + 0.0004597D0*TD2 T(7) = 248.1D0 + 32964.47D0 *TD T(8) = 81.5D0 + 22518.44D0 *TD C C DO 5 I = 1, 909 R = 0.D0 DO 1 J = 1, 8 F = IDN(I,J) R = R + T(J)*F 1 CONTINUE IS = IDN(I,1) NP = IP(I) R = R + PHCORR( NP, IS ) IF( MOD(NP+IS,2).EQ.1 ) R = R - 90.D0 K = R/360.D0 U = K PHC(I) = R - 360.D0*U 5 CONTINUE C RETURN C END SUBROUTINE SETCOS( NCOSGR, COMG, COMG0, TCOS, START, IPOTEN ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION COMG(*), COMG0(*), TCOS(*) C COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /STNGRP/ IG1(36),IG2(36),NGROUP,IMWAVE(36) C IWN3 = 377 IF( IPOTEN.EQ.2 ) IWN3 = 909 C WRITE(6,200) DO 10 I = 1, NCOSGR AVC(I+IWN3) = COMG(I) HH = ( START - TCOS(I) ) * 24.D0 PHC(I+IWN3) = DMOD( COMG0(I) + HH*COMG(I), 360.D0 ) IG1(I+NGROUP) = I + IWN3 IG2(I+NGROUP) = I + IWN3 WRITE(6,210) I, COMG(I), TCOS(I), COMG0(I), PHC(I+IWN3) 10 CONTINUE WRITE(6,220) C 200 FORMAT(1H ///' /* COSINE WAVE FITTING */' * / 5X,'I OMEGA(DEG/HOUR)', * 8X,'EPOCH',7X,'OMEGA0(DEG)',3X,'INITIAL PHASE') 210 FORMAT(1H ,I5,F16.8,F17.2,F16.7,F16.7) 220 FORMAT(1H ) C RETURN END *-TIME SUBROUTINE TIDE01( N1, N2, T, CT, ST, XT, NN ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /SAVEWS/ WS(919), AMPC(919), DIFC(919), CON(919), * AMPS(919), DIFS(919) DATA ICALL / 0 / DATA ININT / 500 / DATA RAD / 0.017453292519943D0 / C IF( ICALL.EQ.0 ) THEN DO 10 I = 1, 377 J = IDN(I,1) K = IP(I) WS(I) = AMC(I)*GC(K,J) 10 CONTINUE DO 12 I = 378, 387 12 WS(I) = 1.D0 ICALL = 1 END IF C CT = 0.D0 ST = 0.D0 C IF( MOD(NN,ININT).EQ.1 ) THEN C C INITIALIZE DELTA = XT*0.5D0*RAD DO 30 I = N1, N2 TEMPB = ( AVC(I)*T + PHC(I) )*RAD TEMPC = AVC(I)*DELTA AMPC(I) = WS(I)*DCOS(TEMPB) CT = CT + AMPC(I) TEMP = 2.D0*DSIN(TEMPC) CON(I) = TEMP*TEMP DIFC(I) = WS(I)*TEMP*DSIN(TEMPC-TEMPB) AMPS(I) = WS(I)*DSIN(TEMPB) ST = ST + AMPS(I) DIFS(I) = WS(I)*TEMP*DCOS(TEMPB-TEMPC) 30 CONTINUE C ELSE C C RECURRENCE FORMULA DO 50 LOOP IS VECTORIZABLE C AMC2(I) IS USED AS AN IAP WORK AREA DO 50 I = N1, N2 CCC AMC2(I) = AMPC(I)*CON(I) CCC DIFC(I) = DIFC(I) - AMC2(I) DIFC(I) = DIFC(I) - AMPC(I)*CON(I) AMPC(I) = AMPC(I) + DIFC(I) CT = CT + AMPC(I) CCC AMC2(I) = AMPS(I)*CON(I) CCC DIFS(I) = DIFS(I) - AMC2(I) DIFS(I) = DIFS(I) - AMPS(I)*CON(I) AMPS(I) = AMPS(I) + DIFS(I) ST = ST + AMPS(I) 50 CONTINUE C END IF C RETURN C END *-TIME SUBROUTINE TIDE02( N1, N2, T, CT, ST, XT, NN, P4FLAG, TC0 ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C INTEGER P4FLAG COMMON /CTETBL/ AVC(919),AMC(919),AMC2(919),PHC(919), * GC(4,4),IDN(919,8),IP(919) COMMON /SAVEWS/ WS(919), AMPC(919), DIFC(919), CON(919), * AMPS(919), DIFS(919) DATA ICALL / 0 / DATA ININT / 500 / DATA RAD / 0.017453292519943D0 / C IF( ICALL.EQ.0 ) THEN TMID = TC0/36525.D0 - 1.0D0 DO 10 I = 1, 909 J = IDN(I,1) K = IP(I) WS(I) = ( AMC(I) + AMC2(I)*TMID )*GC(K,J) IF( P4FLAG.EQ.0 .AND. K.EQ.4 ) WS(I) = 0.D0 10 CONTINUE DO 12 I = 910, 919 12 WS(I) = 1.D0 ICALL = 1 END IF C CT = 0.D0 ST = 0.D0 C IF( MOD(NN,ININT).EQ.1 ) THEN C C INITIALIZE DELTA = XT*0.5D0*RAD DO 30 I = N1, N2 TEMPB = ( AVC(I)*T + PHC(I) )*RAD TEMPC = AVC(I)*DELTA AMPC(I) = WS(I)*DCOS(TEMPB) CT = CT + AMPC(I) TEMP = 2.D0*DSIN(TEMPC) CON(I) = TEMP*TEMP DIFC(I) = WS(I)*TEMP*DSIN(TEMPC-TEMPB) AMPS(I) = WS(I)*DSIN(TEMPB) ST = ST + AMPS(I) DIFS(I) = WS(I)*TEMP*DCOS(TEMPB-TEMPC) 30 CONTINUE C ELSE C C RECURRENCE FORMULA DO 50 LOOP IS VECTORIZABLE C AMC2(I) IS USED AS AN IAP WORK AREA DO 50 I = N1, N2 CCC AMC2(I) = AMPC(I)*CON(I) CCC DIFC(I) = DIFC(I) - AMC2(I) DIFC(I) = DIFC(I) - AMPC(I)*CON(I) AMPC(I) = AMPC(I) + DIFC(I) CT = CT + AMPC(I) CCC AMC2(I) = AMPS(I)*CON(I) CCC DIFS(I) = DIFS(I) - AMC2(I) DIFS(I) = DIFS(I) - AMPS(I)*CON(I) AMPS(I) = AMPS(I) + DIFS(I) ST = ST + AMPS(I) 50 CONTINUE C END IF C RETURN C END SUBROUTINE RTITLE( NUMTIT, NUMINS, CTITLE, CSTNAM, * CINSTR, CUNIT, CASSOC, IAUG ) C C READING TITLS, STATION NAME, INSTRUMENT NAME, UNIT AND C ASSOCIATED DATA NAME C IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*40 CSTNAM, CINSTR(10), CUNIT, CASSOC(*), BLK40 CHARACTER*80 CTITLE(10), CBUFF DATA BLK40 / ' ' / C C SKIP C 10 READ(5,2000) CBUFF(1:4) IF( CBUFF(1:4) .NE. '9999' ) GO TO 10 2000 FORMAT(A4) C C TITLES C NUMTIT = 0 20 CONTINUE C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) CBUFF( 1:40) = BLK40 CBUFF(41:80) = BLK40 READ(5,2010,ERR=28) (CBUFF(I:I),I=1,80) 2010 FORMAT(80A1) 28 IF( CBUFF(1:4) .EQ. '9999' ) GO TO 50 NUMTIT = NUMTIT + 1 IF( NUMTIT.GT.10 ) THEN WRITE(6,2020) 2020 FORMAT(1H /' TOO MANY TITLES') STOP END IF CTITLE(NUMTIT) = CBUFF GO TO 20 C C STATION NAME C 50 CSTNAM = BLK40 C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) READ(5,2030,ERR=58) (CSTNAM(I:I),I=1,40) 2030 FORMAT(40A1) C C INSTRUMENT NAME C 58 NUMINS = 0 60 CBUFF(1:40) = BLK40 C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) READ(5,2030,ERR=68) (CBUFF(I:I),I=1,40) 68 IF( CBUFF(1:4) .EQ. '9999' ) GO TO 90 NUMINS = NUMINS + 1 IF( NUMINS.GT.10 ) THEN WRITE(6,2040) 2040 FORMAT(1H /' TOO MANY INSTRUMENT NAME') STOP END IF CINSTR(NUMINS) = CBUFF(1:40) GO TO 60 C C UNIT C 90 CUNIT = BLK40 C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) READ(5,2030,ERR=98) (CUNIT(I:I),I=1,40) C C ASSOCIATED DATA NAME C 98 IF( IAUG.EQ.0 ) RETURN DO 100 II = 1, IAUG CASSOC(II) = BLK40 C RECOVER IN CASE OF VARIABLE RECORD LENGTH (RECFM=VB) READ(5,2030,ERR=100) (CASSOC(II)(I:I),I=1,40) 100 CONTINUE C RETURN END SUBROUTINE CALEND(JA0,JB0,JC0,HR0,N,DELT,JA,JB,JC,HR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION MONTH(12) DATA MONTH /31,28,31,30,31,30,31,31,30,31,30,31/ C ILL = 0 JAX = MOD(JA0,1900) IF( JAX.LT.1 .OR. JAX.GT.199 ) ILL = 1 IF( JB0.LT.1 .OR. JB0.GT.12 ) ILL = 1 IF( JC0.LT.1 .OR. JC0.GT.31 ) ILL = 1 IF( HR0.LT.0.D0 .OR. HR0.GT.24.D0 ) ILL = 1 IF( ILL.EQ.0 ) GO TO 10 WRITE(6,200) JA0, JB0, JC0, HR0 200 FORMAT(1H /' INVALID DATE Y,M,D,H =', 3I9, F9.1, * ' STOPPED AT CALEND' ) STOP 10 CONTINUE JA = JAX JB = JB0 HR = HR0 + (N-1)*DELT JD = HR/24.D0 HR = HR - 24.D0*JD JC = JC0 + JD 1 M = MONTH(JB) IF( JB .EQ. 2 .AND. MOD(JA,4) .EQ. 0 ) M = M + 1 IF( JC .LE. M ) GO TO 9 JC = JC - M JB = JB + 1 IF( JB .LE. 12 ) GO TO 1 JB = 1 JA = JA + 1 IF( JA .GE. 100 ) JA = MOD( JA, 100 ) GO TO 1 9 CONTINUE RETURN END SUBROUTINE SUBSEA(DMIN,RESPC,ITHAUG,ABICM,TREND,RESPNS,IRREG, * TIDAL,FTRN,N,ISTART,RLIM,DC,IDC,H2,N2,N2MAX, * DATA,INTEST,LIMIT,ISHIFT,ERR,SQEOPT,CENTER,SD0, * MAXDC,MAXH2) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DOUBLE PRECISION IRREG INTEGER SPAN DIMENSION FTRN(*), DC(IDC,*), H2(N2MAX,*), ERR(*), * TREND(*), RESPNS(*), IRREG(*), * DATA(*), TIDAL(*), RESPC(*), CENTER(*) COMMON /IDATA/ IORD,LAGDBL,LAGP,NJUMP,SPAN,LAG,NCOEF,IDUMMY(4) COMMON /RDATA/ WEIGHT,DMIN0,MAXITR DATA EPS12 / 1.D-12 / C C SAVE AREA CHECK ISAVE = 0 ITEM1 = IDC*LIMIT ITEM2 = N2MAX*LIMIT IF( ITEM1*2.LE.MAXDC .AND. ITEM2*2.LE.MAXH2 ) ISAVE = 1 IF( MAXITR.EQ.1 ) ISAVE = 0 C N200 = N2 NJUMP0 = NJUMP ISTEP = 1 IFLAG = 0 DMAX0 = 1024.D0 MODE = 0 RO = DSQRT(2.0D0) ABICM = 1.D25 C ALNDT0 = 0.D0 IF( INTEST.NE.0 ) THEN DD = DMIN0*4.D0 ELSE DD = DMIN END IF C C DO 900 IIII = 1, MAXITR NJUMP=NJUMP0 N2=N200 C C CONSTRUCTION AND HOUSEHOLDER TRANSFORMATION OF MATRIX DC C WT=DD C C CALCULATION OF LOG(DET(DC'DC))*0.5 C C -------------------- ALNDTD=N*DLOG(WT)+ALNDT0 C -------------------- C C CONSTRUCTION AND HOUSEHOLDER TRANSFORMATION OF MATRIX DCX C ALNDN=0.D0 DO 100 I=1,LIMIT 100 DC(1,I)=1.D0 CALL SETX( ALNDN,DC,IDC,H2,N2,N2MAX,M1,ICOUNT,FTRN,N,WT, * N,ISTART,ISTEP,RLIM, * ITHAUG,ALNDTD,INTEST,LIMIT,TIDAL, * CENTER ) C C LEAST SQUARES COMPUTATION C K=M1+N2-1 K=MOD(K,LIMIT) + 1 SQE=H2(N2,K)**2 C C INTERPRETATION OF THE SOLUTION C C **** C C ABIC COMPUTATION C AN = ICOUNT DO 200 I=1,LIMIT TEM=DABS(DC(1,I)) 200 ALNDN=ALNDN + DLOG(TEM) IF( N2.NE.1 .AND. WEIGHT.GT.0.D0 ) THEN N2M1 = N2 - 1 NTEM = N2M1 - NJUMP DO 300 I = 1, NTEM M1I = M1 + I - 1 M1I = MOD(M1I,LIMIT) + 1 TEM = DABS(H2(I,M1I)) 300 ALNDN = ALNDN + DLOG(TEM) END IF SD=0.D0 IF( ICOUNT.LE.0 ) GO TO 500 SQE0=SQE/AN ALSQE=AN*DLOG(SQE0) SD=DSQRT(SQE0) ABIC=ALSQE + 2.D0*(ALNDN-ALNDTD) ABIC=ABIC+NJUMP*2.D0 IF(WEIGHT .LE. 0.D0) ABIC=ABIC+2.D0*NCOEF WRITE(6,1300) DD, ABIC, ALSQE, ALNDN, ALNDTD C C END OF BASIC ROUTINE C C C MINIMUM ABIC PROCEDURE C IF( ABIC.GE.ABICM ) GO TO 500 IF( ABICM-ABIC.LT.0.0001D0 ) GO TO 500 ABICM = ABIC SQEOPT=SQE0 DMIN = DD C C SAVING IF( ISAVE.EQ.1 ) THEN DO 10010 K2 = 1, LIMIT DO 10010 K1 = 1, IDC 10010 DC(K1,LIMIT+K2) = DC(K1,K2) DO 10020 K2 = 1, LIMIT DO 10020 K1 = 1, N2MAX 10020 H2(K1,LIMIT+K2) = H2(K1,K2) CC WITHOUT SUBSCRIPT RANGE CHECK CC DO 10010 KK = 1, ITEM1 CC10010 DC(KK,LIMIT+1) = DC(KK,1) CC DO 10020 KK = 1, ITEM2 CC10020 H2(KK,LIMIT+1) = H2(KK,1) END IF N2SAVE = N2 NJSAVE = NJUMP C SD0 = SD DD = RO*DD IF( IIII.EQ.2 ) MODE = 1 GO TO 700 500 IF( MODE.EQ.1 ) GO TO 1000 MODE = 1 RO = 1.D0/RO DD = DD*RO*RO 700 CONTINUE IF( DD .LE. DMAX0+EPS12 ) GO TO 800 IF( MODE.EQ.0 ) THEN MODE = 1 RO = 1.D0/RO DD = DD*RO*RO GO TO 700 ELSE DD = DMAX0 IFLAG = 1 GO TO 1000 END IF 800 CONTINUE IF( DD .LT. DMIN0-EPS12 ) THEN DD = DMIN0 IFLAG = -1 GO TO 1000 END IF 900 CONTINUE C C 1000 CONTINUE C C C CONSTRUCTION AND HOUSEHOLDER TRANSFORMATION OF MATRIX DC C WT=DMIN C C CALCULATION OF LOG(DET(DC'DC))*0.5 C C -------------------- ALNDTD=N*DLOG(WT)+ALNDT0 C -------------------- C C CONSTRUCTION AND HOUSEHOLDER TRANSFORMATION OF MATRIX DCX C IF( ISAVE.EQ.1 ) THEN DO 10030 K2 = 1, LIMIT DO 10030 K1 = 1, IDC 10030 DC(K1,K2) = DC(K1,LIMIT+K2) DO 10040 K2 = 1, LIMIT DO 10040 K1 = 1, N2MAX 10040 H2(K1,K2) = H2(K1,LIMIT+K2) CC WITHOUT SUBSCRIPT RANGE CHECK CC DO 10030 KK = 1, ITEM1 CC10030 DC(KK,1) = DC(KK,LIMIT+1) CC DO 10040 KK = 1, ITEM2 CC10040 H2(KK,1) = H2(KK,LIMIT+1) N2 = N2SAVE NJUMP = NJSAVE ELSE IF( MAXITR.EQ.1 .OR. IFLAG.NE.0 ) THEN N2 = N2SAVE NJUMP = NJSAVE ELSE NJUMP = NJUMP0 N2 = N200 ALNDN = 0.D0 DO 1050 I=1,LIMIT 1050 DC(1,I)=1.D0 CALL SETX( ALNDN,DC,IDC,H2,N2,N2MAX,M1,ICOUNT,FTRN,N,WT, * N,ISTART,ISTEP,RLIM, * ITHAUG,ALNDTD,INTEST,LIMIT,TIDAL, * CENTER ) END IF C C LEAST SQUARES COMPUTATION C ND=SPAN+N2 IF(ISHIFT .EQ. 0) ND=N2 DO 1060 I=1,ND 1060 ERR(I)=0.D0 CALL SOLVE(DC,IDC,H2,N2,N2MAX,RESPNS,M1,SQE,ND,LIMIT,ERR) C WRITE(6,1400) ABICM, DMIN, SD0 IF( ISHIFT .EQ. 0 ) THEN DO 1100 J = 1, NCOEF 1100 RESPC(J) = RESPNS(J) ELSE CALL DECODE(TREND,RESPNS,IRREG,TIDAL,RESPC,ERR, * SPAN,ISTART,ITHAUG,INTEST,DATA) END IF IF( IFLAG .EQ. 1 ) WRITE(6,1500) DMAX0 IF( IFLAG .EQ. -1 ) WRITE(6,1600) DMIN0 RETURN 1300 FORMAT( 1H ,' ABIC(',1PD12.4,' ) = ',0PF10.2, *3X,'ALSQE =',1PD14.5,4X,'ALNDN =',D13.5,4X,'ALNDTD =',D13.5) 1400 FORMAT(1H /' MINIMUM ABIC =',F10.2,' ATTAINED AT VMIN =', * 1PD12.4,' SD =',D12.4) 1500 FORMAT(1H ,'**** D IS HITTING THE UPPER BOUND ',1PD12.4,' -----', * ' TRY LOWER VALUES OF ORDER') 1600 FORMAT(1H ,'**** D IS HITTING THE LOWER BOUND ',1PD12.4,' -----', * ' TRY HIGHER VALUES OF ORDER') C END SUBROUTINE SETX( ALNDN,H1,N1,H2,N2,N2MAX,M1,ICOUNT,FTRN,N,WT, * NDATA,ISTART,ISTEP,RLIM, * ITHAUG,ALNDTD,INTEST,LIMIT,D, * CENTER ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CCC EXTERNAL YREADO,YREADA C DIMENSION H1(N1,*),H2(N2MAX,*),H3(500),H4(500),D(*) DIMENSION FTRN(*),TI(10),T0(10),DJUMP(300),WORK(3),ERR(3) DIMENSION CENTER(*) COMMON /IDATA/ ID,LAGDBL,LAGP,NJUMP,ISPAN,LAG, * NCOEF,ITH,ITHP1,LAGINT,IFCR COMMON /RDATA/ WEIGHT,DMIN0,MAXITR COMMON /LPCNT/ LPOUT COMMON /VALIDD/ NUMUSE, MAXJMP C LAGPP1=LAGP+1 IIII=ISTART-ISTEP NTEM=ITH + 2+(ITHAUG-ITH)*LAGPP1+NJUMP IF(NTEM .EQ. 0) NTEM=1 DO 100 I=1,100 100 DJUMP(I)=0.D0 DO 200 I=1,LIMIT DO 200 J=1,N2 200 H2(J,I)=0.D0 LAGP1=NCOEF+1 N2M1=N2-1 J0 = 0 IF( ITH.GT.0 ) THEN CCC STABIL = WEIGHT*0.001D0 STABIL = WEIGHT*0.01D0 IF( STABIL.GT.0.01D0 ) STABIL = 0.01D0 H2(1,1) = STABIL H2(2,2) = STABIL J0 = 2 NGROUP = ITH/2 - 1 - IFCR DO 400 K = 1, NGROUP DO 400 J = 1, 2 J0 = J0 + 1 H2(J0-2,J0) = -WEIGHT*CENTER(K) H2(J0, J0) = WEIGHT*CENTER(K) 400 CONTINUE H2(J0-1,J0+1) = STABIL H2(J0, J0+2) = STABIL IF( IFCR.EQ.1 ) THEN J0 = J0 + 1 H2(J0,J0+2) = STABIL J0 = J0 + 1 H2(J0,J0+2) = STABIL END IF END IF IF( ITHP1.LE.ITHAUG ) THEN DO 700 K = ITHP1, ITHAUG DO 700 J = 1, LAGPP1 J0 = J0 + 1 IF( J.LT.LAGPP1 ) H2(J0,J0+3) = -WEIGHT H2(J0,J0+2) = WEIGHT 700 CONTINUE END IF CALL HUSHLD(H2,NTEM,N2,N2MAX,D) IF(N2M1 .LE. 0) GO TO 1300 IF(WEIGHT .LE. 0.D0) GO TO 1300 DO 1200 J=1,N2M1 TEM=DABS(H2(J,J)) IF(J .LE. NCOEF)ALNDTD=ALNDTD+DLOG(TEM) 1200 CONTINUE 1300 CONTINUE CALL SETD(T0,ID,WT) DO 1400 I=1,ID 1400 TI(I) = FTRN(I)*WT CALL INIT(TI,ID) ICOUNT=0 M1=0 ID0=ID+1 IPOS=0 K=0 JMPCHK=0 JUMPNN = 0 NUMUSE = 0 C DO 3300 I=1,N IPOS = IPOS + 1 N3 = IPOS IF(N3 .GT. ID0) N3 = ID0 JTEM = ID0 - N3 DO 1500 J=1,N3 JTEM = JTEM + 1 1500 H3(J) = T0(JTEM) DO 1600 J=1,N2 1600 H4(J) = 0.D0 IF(I .GT. ID) GO TO 1800 H4(N2)=TI(I) IF(INTEST.EQ.0) GO TO 1800 DO 1700 J=I,ID JTEM=J-I+1 KTEM=N2+JTEM-1-ID 1700 H4(KTEM)=T0(JTEM) 1800 CONTINUE CALL EXQRCN(ALNDN,H1,N1,H2,N2,N2MAX,H3,N3,H4,M1,IPOS,LIMIT) IIII=IIII+ISTEP YIIII=YREADO(IIII) N3 = -1 IF(I .GT. NDATA) GO TO 3200 IF(YIIII .GT. RLIM .AND. RLIM .GT. 0.D0) GO TO 3200 IF(YIIII .GE. -RLIM .OR. RLIM .LE. 0.D0) GO TO 1900 JMPCHK = 1 JUMPNN = JUMPNN + 1 IF( JUMPNN.LE.MAXJMP ) GO TO 3200 WRITE(6,1890) MAXJMP 1890 FORMAT(1H /' ERROR TOO MANY STEPS MAXJMP =',I3) STOP 1900 NUMUSE = NUMUSE + 1 C************************** N3=-1 H3(1) = 1.D0 J0=0 IF(ITH .LE. 0) GO TO 2100 CALL YREADT(IIII,H4,ITH) J0=ITH C 2100 CONTINUE IF(ITHP1 .GT. ITHAUG) GO TO 2300 DO 2210 LL=ITHP1,ITHAUG NF0=LL-ITH L=IIII-LAGP*LAGINT IF(L .LE. 0) GO TO 3200 DO 2200 J=1,LAGPP1 WEEKLL=YREADA(NF0,L) IF(WEEKLL .GT. RLIM .AND. RLIM .GT. 0.D0) GO TO 3200 L=L+LAGINT J0=J0+1 H4(J0)=WEEKLL 2200 CONTINUE 2210 CONTINUE C************************** 2300 N3=1 IF(JMPCHK .EQ. 0) GO TO 2900 KTEM=NJUMP+1-INTEST IF(N2.LT.N2MAX) GO TO 2500 CALL SOLVE(H1,N1,H2,N2,N2MAX,WORK,M1,SQE,3,LIMIT,ERR) DJUMP(NJUMP)=WORK(1)+WORK(2)*DJUMP(NJUMP) NJUMPM=NJUMP-INTEST IF(LPOUT.GE.2) WRITE(6,3400) NJUMPM IF(LPOUT.GE.2) WRITE(6,3500) WORK(1),WORK(2),DJUMP(NJUMP) L=N2M1 LM1=L-1 DO 2400 J=1,LIMIT H2(L,J)=WORK(2)*H2(L,J)+WORK(1)*H2(LM1,J) 2400 H2(LM1,J)=0.D0 KTEM=KTEM-3 GO TO 2700 2500 CONTINUE NJUMP=NJUMP+1 DJUMP(NJUMP)=1.D0 LENGTH=M1+N2 LENGTH=MOD(LENGTH,LIMIT)+1 N2=N2+1 N2M1=N2M1+1 L=N2+1 LM1=N2 DO 2600 J=1,N2 2600 H2(J,LENGTH)=0.D0 2700 CONTINUE DO 2800 K=1,KTEM L=L-1 LM1=LM1-1 DO 2800 J=1,LIMIT H2(L,J)=H2(LM1,J) 2800 H2(LM1,J)=0.D0 2900 CONTINUE H4(N2)=YIIII IF(LAGP1 .GT. N2M1) GO TO 3100 L=0 DO 3000 J=LAGP1,N2M1 L=L+1 3000 H4(J)=DJUMP(L) 3100 ICOUNT=ICOUNT+1 JMPCHK=0 3200 CALL EXQRCN(ALNDN,H1,N1,H2,N2,N2MAX,H3,N3,H4,M1,IPOS,LIMIT) 3300 CONTINUE C M1=M1-1 M1=MOD(M1,LIMIT) + 1 RETURN 3400 FORMAT(1H /' ***** MORE THAN',I4,' JUMPS ****') 3500 FORMAT(1H ,' EARLIEST TWO JUMPS',1P2D11.3,' ARE AGGREGATED', * ' TO',D11.3) C END *-TIME DOUBLE PRECISION FUNCTION YREADO(I) C C READING OBSERVED TIDAL DATA C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( LBUFF = 4500 ) PARAMETER ( MAXASC = 3 ) C DIMENSION YO(LBUFF,2), YA(LBUFF,2,MAXASC) REAL YT(LBUFF,2,MAXGR2) COMMON /BUFOTA/ YO, YT, YA, INCORE DATA NBLOCK /0/ C IF( INCORE.EQ.1 ) GO TO 300 C KL=(I-1)/LBUFF K1=KL+1 K2=I-KL*LBUFF C 70 CONTINUE IF(K1 .EQ. NBLOCK) GO TO 100 IF(K1 .EQ. NBLOCK-1) GO TO 100 IF(K1 .GT. NBLOCK) GO TO 50 REWIND 30 NBLOCK=0 50 CONTINUE NBLOCK=NBLOCK+1 IBUFF=MOD(NBLOCK,2)+1 READ(30) ( YO(J,IBUFF), J=1,LBUFF ) GO TO 70 C 100 IBUFF=MOD(K1,2)+1 YREADO=YO(K2,IBUFF) RETURN C CCC DO NOT CHECK SUBSCRIPT RANGE 300 YREADO = YO(I,1) RETURN C END *-TIME SUBROUTINE YREADT(I,BUF,ITH) C C READING THEORETICAL TIDE DATA C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( LBUFF = 4500 ) PARAMETER ( MAXASC = 3 ) C DIMENSION BUF(*) DIMENSION YO(LBUFF,2), YA(LBUFF,2,MAXASC) REAL YT(LBUFF,2,MAXGR2) COMMON /BUFOTA/ YO, YT, YA, INCORE DATA NBLOCK /0/ C IF( INCORE.EQ.1 ) GO TO 300 C KL=(I-1)/LBUFF K1=KL+1 K2=I-KL*LBUFF C 70 CONTINUE IF(K1 .EQ. NBLOCK) GO TO 100 IF(K1 .EQ. NBLOCK-1) GO TO 100 IF(K1 .GT. NBLOCK) GO TO 50 REWIND 31 NBLOCK=0 50 CONTINUE NBLOCK=NBLOCK+1 IBUFF=MOD(NBLOCK,2)+1 DO 58 KTEM = 1, ITH, 2 READ(31) ( YT(J,IBUFF,KTEM), J=1,LBUFF ), * ( YT(J,IBUFF,KTEM+1), J=1,LBUFF ) 58 CONTINUE GO TO 70 C 100 IBUFF=MOD(K1,2)+1 DO 200 K = 1, ITH 200 BUF(K) = YT(K2,IBUFF,K) RETURN C CCC DO NOT CHECK SUBSCRIPT RANGE 300 DO 400 K = 1, ITH 400 BUF(K) = YT(I,1,K) RETURN C END *-TIME DOUBLE PRECISION FUNCTION YREADA(NF0,I) C C READING ASSOCIATED DATA C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) PARAMETER ( LBUFF = 4500 ) PARAMETER ( MAXASC = 3 ) C DIMENSION YO(LBUFF,2), YA(LBUFF,2,MAXASC) REAL YT(LBUFF,2,MAXGR2) COMMON /BUFOTA/ YO, YT, YA, INCORE DIMENSION NBLOCK(MAXASC) DATA NBLOCK /MAXASC*0/ C IF( INCORE.EQ.1 ) GO TO 300 C KL=(I-1)/LBUFF K1=KL+1 K2=I-KL*LBUFF C 70 CONTINUE IF(K1 .EQ. NBLOCK(NF0)) GO TO 100 IF(K1 .EQ. NBLOCK(NF0)-1) GO TO 100 NFILE=NF0+31 IF(K1 .GT. NBLOCK(NF0)) GO TO 50 REWIND NFILE NBLOCK(NF0)=0 50 CONTINUE NBLOCK(NF0)=NBLOCK(NF0)+1 IBUFF=MOD(NBLOCK(NF0),2)+1 READ(NFILE) (YA(J,IBUFF,NF0),J=1,LBUFF) GO TO 70 C 100 IBUFF=MOD(K1,2)+1 YREADA=YA(K2,IBUFF,NF0) RETURN C CCC DO NOT CHECK SUBSCRIPT RANGE 300 YREADA = YA(I,1,NF0) RETURN C END SUBROUTINE HUSHLD( X,N,K,MJ1,D) C C C THIS SUBROUTINE TRANSFORMS MATRIX X INTO AN UPPER TRIANGULAR FORM C BY HOUSEHOLDER TRANSFORMATION. C C INPUTS: C X: ORIGINAL N*K DATA MATRIX C MJ1: ABSOLUTE DIMENSION OF X C N: NUMBER OF ROWS OF X, NOT GREATER THAN MJ1 C K: NUMBER OF COLUMNS OF X. NOT GREATER THAN MJ1 C OUTPUT: C X: IN UPPER TRIANGULAR FORM C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(MJ1,*) , D(*) C TOL=1.0D-38 C IF( MJ1 .GT. 100) STOP MNK=K IF(N.LE.K) MNK=N-1 DO 100 II=1,MNK H = 0.0D00 DO 10 I=II,N D(I) = X(II,I) ABSLD=DABS(D(I)) IF(ABSLD.LE.TOL) D(I)=0.0D-00 10 H = H + D(I)*D(I) IF( H .GT. TOL ) GO TO 20 G = 0.0D00 GO TO 100 20 G = DSQRT( H ) F=X(II,II) IF( F .GE. 0.0D00 ) G = -G D(II) = F-G H = H - F*G C C FORM (I - D*D'/H) * X, WHERE H = D'D/2 C II1 = II+1 DO 30 I=II1,N 30 X(II,I) = 0.D0 IF( II .EQ. K ) GO TO 100 DO 60 J=II1,K S = 0.0D00 DO 40 I=II,N 40 S = S + D(I)*X(J,I) S = S/H DO 50 I=II,N 50 X(J,I) = X(J,I) - D(I)*S 60 CONTINUE 100 X(II,II) = G C RETURN C END SUBROUTINE SETD(W,ID,C) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION W(*) IDP1 = ID + 1 W(IDP1) = C DO 10 J=1,ID 10 W(J) = 0.D0 DO 50 J=1,ID I = IDP1 - J - 1 DO 50 K=1,J I = I + 1 50 W(I) = W(I) - W(I+1) RETURN END SUBROUTINE INIT(W,ID) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION W(*) IDP1 = ID + 1 W(IDP1) = 0.D0 DO 5 J=1,ID DO 5 K=1,ID 5 W(K) = W(K+1) - W(K) DO 20 J=1,ID 20 W(J) = -W(J) RETURN END SUBROUTINE EXQRCN(ALNDN,H1,N1,H2,N2,N2MAX,H3,N3,H4,MACC,IPOS, * LIMIT) C C FOR FORTRAN77 & NOIAP VERSION C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION H1(N1,*), H2(N2MAX,*), H3(*), H4(*) DATA EPS30 /1.D-30/ C IF( IPOS .GT. MACC ) THEN MACC = IPOS MTEM = IPOS - 1 M1 = MOD( MTEM, LIMIT ) + 1 ALNDN = ALNDN + DLOG( DABS( H1(1,M1) ) ) DO 10 J = 1, N1 10 H1(J,M1) = 0.D0 M1N2 = M1 + N2 IF( M1N2 .GT. LIMIT ) M1N2 = M1N2 - LIMIT DO 20 J = 1, N2 20 H2(J,M1N2) = 0.D0 END IF IF( N3 .LT. 0 ) RETURN C M0 = IPOS - N3 - 1 MM = MOD(M0,LIMIT) + 1 DO 70 J = 1, N3 MM = MM + 1 IF( MM .GT. LIMIT ) MM = 1 IF( DABS(H3(J)) .LT. EPS30 ) GO TO 70 SQD = DSQRT( H1(1,MM)**2 + H3(J)**2 ) C = H1(1,MM)/SQD S = H3(J)/SQD H1(1,MM) = SQD M = 1 DO 40 K = J+1, N3 M = M + 1 IF(M .GT. N1) GO TO 50 H1M = H1(M,MM) H1(M,MM) = C*H1M + S*H3(K) 40 H3(K) = C*H3(K) - S*H1M 50 CONTINUE DO 60 K = 1, N2 H2K = H2(K,MM) H2(K,MM) = C*H2K + S*H4(K) H4(K) = C*H4(K) - S*H2K 60 CONTINUE 70 CONTINUE C MM = M1 DO 90 J = 1, N2 MM = MM + 1 IF(MM .GT. LIMIT) MM = 1 IF( DABS(H4(J)) .LT. EPS30 ) GO TO 90 SQD = DSQRT( H2(J,MM)**2 + H4(J)**2 ) C = H2(J,MM)/SQD S = H4(J)/SQD H2(J,MM) = SQD DO 80 K = J+1, N2 H2K = H2(K,MM) H2(K,MM) = C*H2K + S*H4(K) H4(K) = C*H4(K) - S*H2K 80 CONTINUE 90 CONTINUE C RETURN END SUBROUTINE SOLVE(H1,N1,H2,N2,N2MAX,A,M1,SQE,NANS,LIMIT,ERR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION H1(N1,*), H2(N2MAX,*), A(*), ERR(*) C DO 100 LER = 1, NANS K = M1 + N2 - 1 K = MOD(K,LIMIT) + 1 KA = NANS JJ = NANS - 1 IF( LER.NE.NANS ) THEN SQE = 0.D0 KAM1 = KA - 1 DO 10 KKA = 1, KAM1 10 A(KKA) = 0.D0 A(LER) = 1.D0 ELSE SQE = H2(N2,K)**2 KKA = KA KK = K DO 20 J = 1, JJ KK = KK - 1 IF( KK.EQ.0 ) KK = LIMIT KKA = KKA - 1 A(KKA) = H2(N2,KK) 20 CONTINUE END IF C KCOPY = K DO 80 J = 1, JJ KA = KA - 1 IF( A(KA).EQ.0.D0 ) GO TO 80 K = KCOPY - J IF( K.LE.0 ) K = K + LIMIT*( (-K)/LIMIT + 1 ) C IF( J.GE.N2 ) THEN IF( LER.NE.NANS ) GO TO 100 A(KA) = A(KA)/H1(1,K) IF( LER.LT.NANS ) ERR(KA) = ERR(KA) + A(KA)**2 LTEM = K L = KA DO 30 I = 2, N1 L = L - 1 LTEM = LTEM - 1 IF( LTEM.EQ.0 ) LTEM = LIMIT IF( L.LE.0 ) GO TO 80 A(L) = A(L) - A(KA)*H1(I,LTEM) 30 CONTINUE ELSE A(KA) = A(KA)/H2(N2-J,K) IF( LER.LT.NANS ) ERR(KA) = ERR(KA) + A(KA)**2 KAM1 = KA - 1 IF( KAM1.LE.0 ) GO TO 80 LTEM = K KATEM = KA AKA = A(KA) N2MJ = N2-J IF( LTEM .GT. KAM1 ) THEN LOFFS = LTEM - KAM1 - 1 DO 40 L = 1, KAM1 40 A(L) = A(L) - AKA*H2(N2MJ,L+LOFFS) ELSE DO 50 L = 1, KAM1 KATEM = KATEM - 1 LTEM = LTEM - 1 IF( LTEM.EQ.0 ) LTEM = LIMIT A(KATEM) = A(KATEM) - AKA*H2(N2MJ,LTEM) 50 CONTINUE END IF END IF 80 CONTINUE 100 CONTINUE C RETURN END SUBROUTINE DECODE(TREND,RESPNS,IRREG,TIDAL,W,ERR, * NN,ISTART,ITHAUG,INTEST,YS) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER ( MAXGRP = 31, MAXGR2 = MAXGRP*2 ) C C THIS SUBROUTINE COMPUTES C TREND C SEAS0 C RESPNS=Y - SEAS0 C DOUBLE PRECISION IRREG CCC EXTERNAL YREADO,YREADA C DIMENSION BUF(MAXGR2) DIMENSION TREND(*),RESPNS(*),IRREG(*), * ERR(*),YS(*),W(*),TIDAL(*) COMMON /IDATA/ IORDER,LAGDBL,LAGP,NJUMP,ISPAN,LAG, * NCOEF,ITH,ITHP1,LAGINT,IFCR C JP=NJUMP+1 AJUMP=0.D0 RLIM=1.D25 LAGPP1=LAGP+1 N=NN C NTEM = N+1 N7 = NCOEF + NJUMP CALL COPY(W,RESPNS,N7,NTEM) CALL COPY(ERR,ERR,N7,NTEM) DO 10 I=1,N TREND(I)=RESPNS(I) 10 CONTINUE C II=ISTART-1 DO 80 I=1,N II=II+1 IF(I .GT. NN) GO TO 20 YS(I)=YREADO(II) IF(JP .LE. INTEST) GO TO 20 IF(YS(I) .GT. -1.D25) GO TO 20 JP=JP-1 AJUMP=AJUMP+W(NCOEF+JP) C 20 CONTINUE SUM=0.D0 J0=0 IF( ITH.NE.0 ) THEN CALL YREADT(II,BUF,ITH) DO 30 K=1,ITH 30 SUM=SUM+W(K)*BUF(K) J0=ITH END IF C TIDAL(I)=SUM SUM=0.D0 IF(ITHP1 .GT. ITHAUG) GO TO 70 DO 62 K=ITHP1,ITHAUG NF0 = K - ITH JJ=II-LAGP*LAGINT IF(JJ .LE. 0) GO TO 64 DO 60 J=1,LAGPP1 WEEKKJ=YREADA(NF0,JJ) IF(DABS(WEEKKJ) .GE. RLIM) GO TO 64 JJ=JJ+LAGINT J0=J0+1 SUM=SUM+W(J0)*WEEKKJ 60 CONTINUE 62 CONTINUE GO TO 70 64 SUM=1.D30 C 70 CONTINUE RESPNS(I)=SUM C************************** IRREG(I) = YS(I) - RESPNS(I) - TIDAL(I) - AJUMP IF(DABS(YS(I)).GT.RLIM.OR.DABS(RESPNS(I)).GT.RLIM) IRREG(I)=1.D30 IRREG(I)=IRREG(I)-TREND(I) 80 CONTINUE C RETURN END SUBROUTINE COPY(X,Y,NY,IY) C THIS SUBROUTINE COPIES Y INTO X. C INPUTS: C X: MX*NX MATRIX C Y: MY*NY MATRIX C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(*), Y(*) IYM1 = IY - 1 DO 50 I=1,NY 50 X(I)=Y(I+IYM1) RETURN END SUBROUTINE FRQRSP( RAT, RESPC, ERR, LAGINT, LAGPP1, * IAUG, DEG0, PAGLEN, RSPERR ) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION RESPC(LAGPP1,*), FY(7), SXX(181,2), FRESP(181,2) DIMENSION AMYM(2), DMYM(2), FMYM(2), RAT(*) DIMENSION ERR(LAGPP1,*), ERRA(181), ERRP(181) INTEGER H1, H1END, PAGLEN, RSPERR CHARACTER*62 CC C DATA PI / 3.14159 26535 89793 D0 / C H1 = 61 H1END = H1 IF( PAGLEN.EQ.60 ) H1END = H1END - 5 C COEF = PI*LAGINT/(H1-1) DEG = 180.D0/DEG0/DBLE(H1-1) FNYQ = 180.D0 / ( DEG0*DBLE(LAGINT) ) C DO 90 IAUGT=1,IAUG C WRITE(6,2000) IAUGT 2000 FORMAT(1H1,' FREQUENCY RESPONSE OF ASSOCIATED DATA', * ' TO',I3,'-TH AUX. INPUT' /) C DO 30 I=1,H1 AR=RESPC(LAGPP1,IAUGT) AI=0.D0 ARERR = ERR(LAGPP1,IAUGT)**2 AIERR = 0.D0 IF(LAGPP1 .LE. 1) GO TO 20 C=DBLE(I-1)*COEF DO 10 J=2,LAGPP1 ARG=C*DBLE(J-1) JR = LAGPP1+1-J AR=AR+RESPC(JR,IAUGT)*DCOS(ARG) AI=AI-RESPC(JR,IAUGT)*DSIN(ARG) ARERR = ARERR + (ERR(JR,IAUGT)*DCOS(ARG))**2 AIERR = AIERR + (ERR(JR,IAUGT)*DSIN(ARG))**2 10 CONTINUE 20 AR = AR * RAT(IAUGT+1) AI = AI * RAT(IAUGT+1) ARERR = DSQRT(ARERR) * RAT(IAUGT+1) AIERR = DSQRT(AIERR) * RAT(IAUGT+1) FRESP(I,1)=DSQRT(AR**2 + AI**2) SXX(I,1)=DLOG10(FRESP(I,1))*2.D0 FRESP(I,2)=DATAN2(AI,AR)*180.D0/PI SXX(I,2)=FRESP(I,2) ERRA(I) = DSQRT( (AR*ARERR)**2 + (AI*AIERR)**2 ) / FRESP(I,1) ERRP(I) = DSQRT( (AR*AIERR)**2 + (AI*ARERR)**2 ) * / FRESP(I,1)**2 * 180.D0/PI 30 CONTINUE DMYM(1)=10.D0 CALL DSP3( SXX(1,1), H1, R8MIN ) AMYM(1)=R8MIN*10.D0 DO 40 I = 1, H1 40 SXX(I,1)=SXX(I,1)*10.D0 AMYM(2)=-180.D0 DMYM(2)=60.D0 FMYM(1)=1.D0 FMYM(2)=60.D0/360.D0 DO 80 KKK=1,2 IF(KKK.EQ.2) WRITE(6,2010) 2010 FORMAT(1H1) FY(1)=AMYM(KKK) DO 50 I=2,7 50 FY(I)=FY(I-1)+DMYM(KKK) WRITE(6,2020) (FY(I),I=1,7) 2020 FORMAT(1H ,24X,6(F6.1,4X),F6.1) IF(KKK.EQ.1) WRITE(6,2030) IF(KKK.EQ.2) WRITE(6,2040) 2030 FORMAT(1H ,'DEG/H AMPLITUDE FACTOR ',6('+---------'),'+') 2040 FORMAT(1H ,'DEG/H PHASE SHIFT IN DEG ',6('+---------'),'+') DO 60 J = 1, 62 60 CC(J:J) = ' ' DO 70 I = 1, H1END IT=(SXX(I,KKK)-AMYM(KKK))*FMYM(KKK)+1.0 IF( IT.LT. 1 ) IT = 1 IF( IT.GT.62 ) IT = 62 IM1=I-1 DIM1=IM1*DEG CC(1:1) = '!' CC(IT:IT) = '*' IF( RSPERR.EQ.0 ) THEN IF(KKK.EQ.1) WRITE(6,2060) DIM1, FRESP(I,1), CC IF(KKK.EQ.2) WRITE(6,2070) DIM1, FRESP(I,2), CC ELSE IF(KKK.EQ.1) WRITE(6,2060) DIM1, FRESP(I,1), CC, ERRA(I) IF(KKK.EQ.2) WRITE(6,2070) DIM1, FRESP(I,2), CC, ERRP(I) END IF CC(IT:IT) = ' ' 70 CONTINUE C FOR NORMAL VERSION 2060 FORMAT(1H ,F5.1,2X,1PD15.4,6X,A62,1PD17.4) 2070 FORMAT(1H ,F5.1,2X,F13.2, 8X,A62,F15.2) C FOR LONG PERIOD VERSION C2060 FORMAT(1H ,F6.3,1X,1PD15.4,6X,A62,1PD17.4) C2070 FORMAT(1H ,F6.3,1X,F13.2, 8X,A62,F15.2) 80 CONTINUE IF(LAGPP1.LE.1) WRITE(6,2080) IF(LAGPP1.GE.2) WRITE(6,2090) FNYQ 2080 FORMAT(1H ,' (LAG:NEGATIVE)') 2090 FORMAT(1H ,' (LAG:NEGATIVE)', * ' NYQUIST FREQUENCY =',F6.1,' (DEG./HOUR)') C 90 CONTINUE C RETURN END SUBROUTINE DSP3( Y, LAGH1, YMIN ) C C SEARCH FOR MINIMUM OF Y C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(LAGH1) C YMIN = Y(1) DO 800 I = 2, LAGH1 IF( Y(I).LT.YMIN ) YMIN = Y(I) 800 CONTINUE C IF( YMIN.GE.0.D0 ) THEN IMIN = YMIN ELSE YMIN = YMIN - 1.D0 IMIN = -YMIN IMIN = -IMIN END IF YMIN = IMIN C RETURN END SUBROUTINE OTLCHK(DATA,N,IDIF,LIST,L,W,LMAX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DATA(*), LIST(*) DIMENSION W(*) L=0 NMD=N-IDIF S=0.D0 S2=S COUNT=S DO 1 I=1,NMD IF(IDIF .EQ. 0) W(I)=DATA(I) IF(IDIF .NE. 0) W(I)=DATA(I+1)-DATA(I) IF(DABS(W(I)) .GT. 1.D25) GO TO 1 S=S+W(I) S2=S2+W(I)**2 COUNT=COUNT+1.D0 1 CONTINUE IF(COUNT .LT. 0.5D0) RETURN S=S/COUNT S2=S2/COUNT - S**2 SD4=DSQRT(S2)*4.D0 DO 2 I=1,NMD DIV=W(I)-S DIV=DABS(DIV) IF(DIV .GT. 1.D25) GO TO 2 IF(DIV .LE. SD4) GO TO 2 IF(L .GE. LMAX) GO TO 2 L=L+1 LIST(L)=I 2 CONTINUE RETURN END SUBROUTINE SAUTCO(X,CXX,N,LAGH1,LPOUT) C THE OUTPUTS ARE AUTOCOVARIANCES (CXX(I); I=0,LAGH) AND C AUTO CORRELATIONS (NORMALIZED COVARIANCES). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N), CXX(LAGH1), CN(361) CHARACTER F(21) C C INITIAL CONDITION PRINT OUT LAGH=LAGH1-1 C MEAN DELETION CALL SMEADL(X,N,XMEAN) C AUTO COVARIANCE COMPUTATION CALL CROSCO(X,X,N,CXX,LAGH1) C NORMALIZATION CX0=CXX(1) CALL CORNOM(CXX,CN,LAGH1,CX0,CX0) C AUTO COVARIANCE PRINT OUT IF(LPOUT .LT. 2) RETURN WRITE(6,162) N,LAGH,XMEAN WRITE(6,262) WRITE(6,163) DO 80 I=1,LAGH1 IM1=I-1 DO 90 J = 1, 21 90 F(J) = ' ' F(11) = '!' CST0=0.0D-00 CST10=10.0D-00 CST05=0.5D-00 IF(CN(I).LT.CST0) GO TO 300 NFC=CST10*CN(I)+CST05 GO TO 310 300 NFC=CST10*CN(I)-CST05 310 ANFC=NFC JJ=ANFC+11.D0 F(JJ) = '*' WRITE(6,199) IM1,CXX(I),CN(I),(F(J),J=1,21) 80 CONTINUE RETURN 162 FORMAT(1H /' AUTOCOVARIANCE CXX(I) N=',I5,5X,'LAGH=',I5, * 5X,'MEAN=',1PD15.4) 262 FORMAT(1H ,39X,'-1',9X,'0',9X,'1') 163 FORMAT(1H ,4X,'I CXX(I)',8X,'NORMALIZED',6X,4('+----'),'+') 199 FORMAT(1H ,I5,2X,1P2D14.4,5X,21A1) END SUBROUTINE SMEADL(X,N,XMEAN) C MEAN DELETION C SUMF: FUNCTION IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N) CCC EXTERNAL SUMF C AN=N XMEAN=SUMF(X,N)/AN DO 10 I=1,N 10 X(I)=X(I)-XMEAN RETURN END DOUBLE PRECISION FUNCTION SUMF( X, N ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N) SUMF = 0.D0 DO 10 I = 1, N 10 SUMF = SUMF + X(I) RETURN END SUBROUTINE CROSCO(X,Y,N,C,LAGH1) C THIS SUBROUTINE COMPUTES C(L)=COVARIANCE(X(S+L),Y(S)) C (L=0,1,...,LAGH1-1). IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),Y(N),C(LAGH1) AN=N BN1=1.0D-00 BN=BN1/AN CT0=0.0D-00 DO 10 II=1,LAGH1 I=II-1 T=CT0 IL=N-I DO 20 J=1,IL J1=J+I 20 T=T+X(J1)*Y(J) 10 C(II)=T*BN RETURN END SUBROUTINE CORNOM(C,CN,LAGH1,CX0,CY0) C NORMALIZATION OF COVARIANCE IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION C(LAGH1),CN(LAGH1) CST1=1.0D-00 DS=CST1/DSQRT(CX0*CY0) DO 10 I=1,LAGH1 10 CN(I)=C(I)*DS RETURN END SUBROUTINE SICP2(CYY,L1,N,COEF,MO,OSD,OAIC) C THIS SUBROUTINE FITS AUTOREGRESSIVE MODELS OF SUCCESSIVELY C INCREASING ORDER UP TO L(=L1-1). C INPUT: C CYY(I),I=0,L1; AUTOCOVARIANCE SEQUENCE C L1: L1=L+1, L IS THE UPPER LIMIT OF THE MODEL ORDER C N; LENGTH OF ORIGINAL DATA C OUT PUT: C COEF; AR-COEFFICIENTS C MO: ORDER OF AR C OSD: INNOVATION VARIANCE C OAIC: VALUE OF AIC IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CYY(L1), COEF(L1) DIMENSION A(101), B(101) CHARACTER FFFF, F(41), AMES(41) COMMON /LPCNT/ LPOUT DATA F / 41*' ' /, AMES / 41*' ' / C CST0=0.0D-00 CST1=1.0D-00 CST2=2.0D-00 CST20=20.0D-00 CST05=0.05D-00 CST01=0.00001D-00 L=L1-1 SD=CYY(1) AN=N OAIC=AN*DLOG(SD) OSD=SD MO=0 C C INITIAL CONDITION PRINT OUT C RAN=CST1/DSQRT(AN) SCALH=CST20 JJ0=SCALH+CST1 JJL=SCALH*CST2+CST1 AMES(1) = '+' AMES(11) = '+' AMES(JJ0) = '+' AMES(JJ0+10) = '+' AMES(JJL) = '+' IAN=SCALH*(RAN+CST05) IAN1=IAN+JJ0 IAN2=2*IAN+JJ0 LAN1=-IAN+JJ0 LAN2=-2*IAN+JJ0 IF(LPOUT.GE.2) WRITE(6,26101) IF(LPOUT.GE.2) WRITE(6,261) IF(LPOUT.GE.2) WRITE(6,26102) IF(LPOUT.GE.2) WRITE(6,262) IF(LPOUT.GE.2) WRITE(6,264) (AMES(J),J=1,JJL) IF(LPOUT.GE.2) WRITE(6,859) MO,OSD,OAIC F(JJ0) = '!' F(IAN1) = '!' F(IAN2) = '!' F(LAN1) = '!' F(LAN2) = '!' IF(LPOUT.GE.2) WRITE(6,861) (F(J),J=1,JJL) SE=CYY(2) C ITERATION START DO 400 M=1,L SDR=SD/CYY(1) IF(SDR.GE.CST01) GO TO 399 IF(LPOUT.GE.2) WRITE(6,2600) GO TO 402 399 MP1=M+1 D=SE/SD A(M)=D D2=D*D SD=(CST1-D2)*SD AM=M AIC=AN*DLOG(SD)+CST2*AM IF(M.EQ.1) GO TO 410 C A(I) COMPUTATION LM=M-1 DO 420 I=1,LM A(I)=A(I)-D*B(I) 420 CONTINUE 410 CONTINUE DO 421 I=1,M IM=MP1-I 421 B(I)=A(IM) C M,SD,AIC PRINT OUT IF(A(M).LT.CST0) GO TO 300 NFC=SCALH*(A(M)+CST05) GO TO 310 300 NFC=SCALH*(A(M)-CST05) 310 ANFC=NFC JJ=ANFC+SCALH+CST1 FFFF = F(JJ) F(JJ) = '*' IF(LPOUT.GE.2) WRITE(6,860) M,SD,AIC,A(M),(F(J),J=1,JJL) F(JJ) = FFFF C C C ----- 5/15/80 ----- C IF(OAIC.LT.AIC) GO TO 440 OAIC=AIC OSD=SD MO=M GO TO 440 C DO 430 I=1,M C 430 COEF(I)=-A(I) 440 IF(M.EQ.L) GO TO 400 SE=CYY(M+2) DO 441 I=1,M 441 SE=SE-B(I)*CYY(I+1) 400 CONTINUE 402 CONTINUE IF(LPOUT.GE.2) WRITE(6,870) OAIC,MO C ----- 5/15/80 ----- OAIC=AIC OSD=SD MO=L DO 5100 I=1,L 5100 COEF(I)=-A(I) C F(JJ0) = ' ' F(IAN1) = ' ' F(IAN2) = ' ' F(LAN1) = ' ' F(LAN2) = ' ' AMES(JJ0) = '-' AMES(JJ0+10) = '-' AMES(JJL) = '-' RETURN 26101 FORMAT(1H /////4X,'M',11X,'INNOVATION',10X,'AIC(M)= ',8X, A 'PARTIAL AUTO-') 261 FORMAT(1H ,16X,' VARIANCE',9X,'N*DLOG(SD(M))+2*M', A ' CORRELATION ', A 'PARTIAL CORRELATION (LINES SHOW +SD AND +2SD)') 26102 FORMAT(1H+,102X,'_',7X,'_') 262 FORMAT(1H ,69X,'-1',19X,'0',19X,'1') 264 FORMAT(1H ,70X,41A1) 859 FORMAT(1H ,I5,2X,2D20.5) 860 FORMAT(1H ,I5,2X,3D20.5,3X,41A1) 861 FORMAT(1H+,70X,41A1) 870 FORMAT(1H /' MINIMUM AIC(M)=',1PD12.4,2X,'ATTAINED AT M=',I5) 2600 FORMAT(1H /' ACCURACY OF COMPUTATION LOST') END SUBROUTINE SNRASP(A,B,SXX,SGME2,L,K,H1,DEG0,PAGLEN) C THIS PROGRAM COMPUTES POWER SPECTRUM OF AN AR-MA PROCESS DEFINED B C X(N)=A(1)X(N-1)+...+A(L)X(N-L)+E(N)+B(1)E(N-1)+...+B(K)E(N-K), C WHERE E(N) IS A WHITE NOISE WITH ZERO MEAN AND VARIANCE EQUAL TO C SGME2. OUTPUTS PXX(I) ARE GIVEN AT FREQUENCIES I/(2*H) C I=0,1,...,H. C REQUIRED INPUTS ARE: C L,K,H,SGME2,(A(I),I=1,L), AND (B(I),I=1,K). C 0 IS ALLOWABLE AS L AND/OR K. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C INTEGER H1, H1END, PAGLEN DIMENSION A(101),B(101) DIMENSION G(361),GR1(361),GI1(361),GR2(361),GI2(361) DIMENSION PXX(361),SXX(361),FY(10) CHARACTER*72 CC C CST0=0.0D-00 CST1=1.0D-00 DEG=180.D0/DEG0/(H1-1.D0) IF(L.LE.0) GO TO 320 DO 330 I=1,L 330 A(I)=-A(I) 320 L1=L+1 K1=K+1 G(1)=CST1 IF(L.LE.0) GO TO 400 DO 10 I=1,L I1=I+1 10 G(I1)=-A(I) C 400 CALL FOUGER(G,L1,GR1,GI1,H1) G(1)=CST1 IF(K.LE.0) GO TO 410 DO 20 I=1,K I1=I+1 20 G(I1)=B(I) C 410 CALL FOUGER(G,K1,GR2,GI2,H1) DO 30 I=1,H1 30 PXX(I)=(GR2(I)**2+GI2(I)**2)/(GR1(I)**2+GI1(I)**2)*SGME2 IF(L.LE.0) GO TO 500 DO 340 I=1,L 340 A(I)=-A(I) 500 IF(K.LE.0) GO TO 510 WRITE(6,63) C CALL SUBVCP(B,K) 510 T0=CST0 DO 520 I=1,H1 T=PXX(I) IF(T.LT.T0) T=-T 520 SXX(I)=DLOG10(T) C SEARCH FOR THE MINIMUM OF (SXX(I),I=0,LAGH) CALL DSP3( SXX, H1, R8MIN ) WRITE(6,266) 266 FORMAT(1H ,47X,'(DENSITY IN DB)') FY(1)=R8MIN*10.D0 DO 2031 I=2,7 2031 FY(I)=FY(I-1)+10.D0 WRITE(6,2033) (FY(I),I=1,7) 2033 FORMAT(1H ,24X,6(F6.1,4X),F6.1) WRITE(6,2035) 2035 FORMAT(1H ,'DEG/H RATIONAL SPECTRUM ',6('+---------'),'+') DO 2050 J = 1, 72 2050 CC(J:J) = ' ' C H1END = H1 IF( PAGLEN.EQ.60 ) H1END = H1END - 5 DO 2040 I = 1, H1END IT=(SXX(I)-R8MIN)*10.D0+1.D0 IF( IT.LT. 1 ) IT = 1 IF( IT.GT.72 ) IT = 72 IM1=I-1 DIM1=IM1*DEG CC(1:1) = '!' CC(IT:IT) = '*' WRITE(6,2060) DIM1, PXX(I), CC CC(IT:IT) = ' ' 2040 CONTINUE C C FOR NORMAL VERSION 2060 FORMAT(1H ,F5.1,2X,1PD14.4,7X,A72) C FOR LONG PERIOD VERSION C2060 FORMAT(1H ,F6.3,1X,1PD14.4,7X,A72) RETURN 63 FORMAT(1H /6X,'MA(I)') END SUBROUTINE FOUGER(G,LGP1,FC,FS,LF1) C FOURIER TRANSFORM (GOERTZEL METHOD) C THIS SUBROUTINE COMPUTES FOURIER TRANSFORM OF G(I),I=0,1,...,LG AT C FREQUENCIES K/(2*LF),K=0,1,...,LF AND RETURNS COSIN TRANSFORM IN C FC(K) AND SIN TRANSFORM IN FS(K). IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION G(LGP1),FC(LF1),FS(LF1) LG=LGP1-1 LF=LF1-1 CST0=0.0D-00 C REVERSAL OF G(I),I=1,...,LGP1 INTO G(LG3-I) LG3=LGP1+1 IF(LGP1.LE.1) GO TO 110 LG3=LGP1+1 LG4=LGP1/2 DO 100 I=1,LG4 I2=LG3-I T=G(I) G(I)=G(I2) 100 G(I2)=T 110 PI = 3.141592653589793D0 ALF=LF T=PI/ALF DO 10 K=1,LF1 AK=K-1 TK=T*AK CK=DCOS(TK) SK=DSIN(TK) CK2=CK+CK UM2=CST0 UM1=CST0 IF(LG.EQ.0) GO TO 12 DO 11 I=1,LG UM0=CK2*UM1-UM2+G(I) UM2=UM1 11 UM1=UM0 12 FC(K)=CK*UM1-UM2+G(LGP1) FS(K)=SK*UM1 10 CONTINUE RETURN END SUBROUTINE SUBVCP(HH,LC) C THIS SUBROUTINE PRINTS OUT VECTOR (HH(I),I=1,LC). IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION HH(LC) WRITE(6,39) LC WRITE(6,20) J1=0 J2=0 I=1 JC=0 16 J1=J2+1 J2=J1+9 IF(J2.LE.LC) GO TO 13 J2=LC 13 JC=JC+1 WRITE(6,21) I,(HH(J),J=J1,J2) IF(J2.GE.LC) GO TO 10 I=I+10 GO TO 16 10 RETURN 39 FORMAT(1H /' VECTOR LC=',I5) 20 FORMAT(1H /16X,'1',11X,'2',11X,'3',11X,'4',11X,'5',11X,'6',11X, * '7',11X,'8',11X,'9',10X,'10') 21 FORMAT(1H ,I5,4X,1P10D12.4) END SUBROUTINE GRAPH3( X, Y, Z, RS, TD, N, PREPRO, * RLIM, JY, JM, JD, HR, DELT ) C C THIS SUBROUTINE PRINTS OUT THE GRAPH OF TIME SERIES X. C DATA WITH VALUES GREATER THAN OR EQUAL TO RLIM C ARE IGNORED. C INPUTS: C X : N-VECTOR ( ORIGINAL * FACTOR ) C Y : " ( TREND ) C Z : " ( IRREGULER PART ) C RS: " ( RESPONSE PART ) C TD: " ( TIDAL COMPONENT ) C RLIM: RANGE OF THE VALUE OF X C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION X(*), Y(*), Z(*), RS(*), TD(*), AXIS(11) INTEGER PREPRO CHARACTER C101, CTEM CHARACTER*11 CZ CHARACTER*25 LFORM CHARACTER*101 CC, C0, C1 COMMON /PRSTEP/ STEP00 C DATA LFORM / '(1H ,18X,1PD10.3,10D10.3)' / C XMAX = -1.D25 ZMAX = XMAX XMIN = 1.D25 ZMIN = XMIN DO 10 I = 1, N IF(RLIM .LE. 0.D0) GO TO 5 IF(DABS(X(I)) .GE. RLIM) GO TO 10 IF(DABS(Z(I)) .GE. RLIM) GO TO 5 IF(Z(I) .GT. ZMAX) ZMAX = Z(I) IF(Z(I) .LT. ZMIN) ZMIN = Z(I) 5 XTEM = X(I) IF( XTEM .GT. XMAX ) XMAX = XTEM IF( XTEM .LT. XMIN ) XMIN = XTEM 10 CONTINUE ZD = (ZMAX-ZMIN)*0.1D0 ZMIN = ZMIN - 1.5D0*ZD XDIF = XMAX - XMIN IX = IDINT( DLOG10(XDIF) ) POWER = DBLE( 10**IX ) XDIFN = XDIF/POWER U = 0.2D0 IF(XDIFN .GE. 1.5D0) U = 0.25D0 IF(XDIFN .GE. 2.0D0) U = 0.5D0 IF(XDIFN .GE. 4.0D0) U = 1.D0 IF(XDIFN .GE. 9.0D0) U = 2.D0 UNIT = U * POWER NUM = IDINT( XDIF/UNIT ) + 2 LP = 100/NUM MESH = LP*NUM LPMAX = MESH + 1 IXMIN = IDINT( XMIN/UNIT ) IF(XMIN .LE. 0.D0) IXMIN = IXMIN - 1 XMIN = DBLE(IXMIN)*UNIT XMAX = XMIN + DBLE(NUM)*UNIT XD = (XMAX-XMIN)/DBLE(MESH) XD2 = XD*0.5D0 XD10 = XD*DBLE(LP) C NUMP1 = NUM + 1 REL = DABS( XMAX + XMIN )*0.5D0/( XMAX - XMIN ) XBASE = XMIN IF(REL .GT. 50.D0) THEN WRITE(6,1600) XMIN 1600 FORMAT(100X,'( +',1PD13.5,' )') XBASE = 0.D0 END IF DO 40 I = 1, NUMP1 AXIS(I) = XBASE + XD10*DBLE(I-1) IF(DABS(AXIS(I)) .LT. XD10*1.D-4) AXIS(I) = 0.D0 40 CONTINUE C WRITE(6,500) WRITE(LFORM(21:22),'(I2)') LP WRITE(6,LFORM) (AXIS(I),I=1,NUMP1) 500 FORMAT(1H1/1H ,'* : DATA'/ * 1H ,'XXXX : CONNECTING LINE'/ * 1H ,'++++ : TREND+STEP'/1H ,'I : IRREGULAR'//) C DO 50 J = 1, 101 CC(J:J) = ' ' IF( MOD(J,LP) .EQ. 1 ) THEN C1(J:J) = '!' C0(J:J) = '!' ELSE C1(J:J) = '-' C0(J:J) = ' ' END IF 50 CONTINUE XMIN2 = XMIN - XD2 C C C FOR LONG PERIOD VERSION C JMLAST = 0 C IF( JD.NE.1 ) JMLAST = JM C DO 100 I = 1, N IOBS = 0 C FOR NOMAL VERSION IF( HR.GE.0.001D0 ) THEN C FOR LONG PERIOD VERSION C IF( JM.EQ.JMLAST ) THEN CC(1:LPMAX) = C0(1:LPMAX) ELSE CC(1:LPMAX) = C1(1:LPMAX) END IF CCC ITEM = IDINT( (X(I)-XMIN2)/XD ) + 1 XPOSIT = (X(I)-XMIN2)/XD + 1.D0 IF( XPOSIT.GT.1.D4 ) THEN ITEM = 10000 ELSE IF( XPOSIT.LT.0.D0 ) THEN ITEM = 0 ELSE ITEM = IDINT( XPOSIT ) END IF ITEM0 = ITEM ITEM1 = ITEM CCC IF(I .GT. 1) ITEM0 = IDINT( (X(I-1) - XMIN2)/XD ) + 1 IF(I .GT. 1) THEN XPOSI0 = (X(I-1) - XMIN2)/XD + 1.D0 IF( DABS(XPOSI0).GT.1.D3 ) THEN ITEM0 = ITEM ELSE ITEM0 = IDINT( XPOSI0 ) END IF END IF CCC IF(I .LT. N) ITEM1 = IDINT( (X(I+1) - XMIN2)/XD ) + 1 IF(I .LT. N) THEN XPOSI1 = (X(I+1) - XMIN2)/XD + 1.D0 IF( DABS(XPOSI1).GT.1.D3 ) THEN ITEM1 = ITEM ELSE ITEM1 = IDINT( XPOSI1 ) END IF END IF IF(ITEM .GT. LPMAX) GO TO 90 IF(ITEM .LT. 1) GO TO 90 IF(ITEM0 .LT. ITEM) ITEM0 = ITEM - (ITEM-ITEM0-1)/2 IF(ITEM0 .GT. ITEM) ITEM0 = ITEM + (ITEM0-ITEM-1)/2 IF(ITEM1 .LT. ITEM) ITEM1 = ITEM - (ITEM-ITEM1-1)/2 IF(ITEM1 .GT. ITEM) ITEM1 = ITEM + (ITEM1-ITEM-1)/2 IF(ITEM1 .GT. LPMAX) ITEM1 = ITEM IF(ITEM1 .LT. 1) ITEM1 = ITEM IF(ITEM0 .GT. LPMAX) ITEM0 = ITEM IF(ITEM0 .LT. 1) ITEM0 = ITEM C JTEM0 = ITEM0 IF(JTEM0 .GT. ITEM) JTEM0 = ITEM IF(JTEM0 .GT. ITEM1) JTEM0 = ITEM1 JTEM1 = ITEM0 IF(JTEM1 .LT. ITEM) JTEM1 = ITEM IF(JTEM1 .LT. ITEM1) JTEM1 = ITEM1 DO 80 J = JTEM0, JTEM1 80 CC(J:J) = 'X' CC(ITEM:ITEM) = '*' IOBS = 1 90 ITEM = IDINT( (Y(I)-XMIN2)/XD ) + 1 IF(ITEM .GT. LPMAX) GO TO 91 IF(ITEM .LT. 1) GO TO 91 CTEM = '+' IF(CC(ITEM:ITEM) .EQ. '*') CTEM = 'C' CC(ITEM:ITEM) = CTEM 91 CONTINUE CZ = ' ' IF( Z(I).LT.1.D25 ) THEN KTEM = IDINT( (Z(I)-ZMIN)/ZD ) IF( KTEM.GE.1 .AND. KTEM.LE.11 ) CZ(KTEM:KTEM) = 'I' END IF IF( PREPRO.GE.1 ) THEN CALL KANDO( COEF, JY, JM, JD, HR, 0 ) IF( RS(I).LT.1.D25 ) THEN PRED = ( Y(I) + RS(I) + TD(I) - STEP00 ) / COEF ELSE PRED = ( Y(I) + TD(I) - STEP00 ) / COEF END IF C101 = CC(101:101) IF( IOBS.EQ.0 ) CC(101:101) = '*' IF( PREPRO.EQ.1 ) THEN WRITE(6,601) JY,JM,JD,HR, CZ, CC, PRED ELSE IF( PREPRO.EQ.2 ) THEN WRITE(6,602) JY,JM,JD,HR, CZ, CC, PRED ELSE WRITE(6,603) JY,JM,JD,HR, CZ, CC, PRED END IF CC(101:101) = C101 ELSE WRITE(6,600) JY,JM,JD,HR, CZ, CC END IF C C FOR LONG PERIOD VERSION C JMLAST = JM C CALL UPDATE(JY,JM,JD,HR,DELT) 600 FORMAT(1H ,3I2,F4.1,1X,A11,1X,A101) 601 FORMAT(1H ,3I2,F4.1,1X,A11,1X,A101,F8.1) 602 FORMAT(1H ,3I2,F4.1,1X,A11,1X,A101,F8.2) 603 FORMAT(1H ,3I2,F4.1,1X,A11,1X,A101,F8.3) C 100 CONTINUE C WRITE(6,LFORM) (AXIS(I),I=1,NUMP1) RETURN END *-TIME SUBROUTINE KANDO( COEF, IY, IM, ID, HR, IRESET ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER ( MAXCLB = 100 ) COMMON /INFORM/ IFORMT, NETALT, ETALT(MAXCLB), TJET(MAXCLB) DATA KSTART / 1 / CCC EXTERNAL TJUL C IF( IFORMT.EQ.1 ) THEN COEF = 1.D0 ELSE IF( IRESET.EQ.1 ) KSTART = 1 T = TJUL( IY, IM, ID, HR ) DO 10 K = KSTART, NETALT-1 IF( T.GE.TJET(K) .AND. T.LT.TJET(K+1) ) GO TO 20 10 CONTINUE WRITE(6,1000) NETALT, IY, IM, ID, HR 1000 FORMAT(1H /' ERROR AT KANDO NETALT =', I4, * ' IY,IM,ID,HR =', 3I4, F6.1) STOP C 20 KSTART = K COEF = ETALT(K) + ( ETALT(K+1) - ETALT(K) ) * * ( T - TJET(K) ) / ( TJET(K+1) - TJET(K) ) END IF C RETURN END SUBROUTINE BLOCKD C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C CHARACTER*8 CMODEL, CVERSI CHARACTER*60 CCALCU C COMMON /DTDMUT/ TDMUT(130) COMMON /VERNAM/ CMODEL, CVERSI, CCALCU C C VERSION NAME CMODEL = 'BAYTAP-G' CVERSI = '97-03-20' CCALCU = * 'NATIONAL ASTRONOMICAL OBSERVATORY, MIZUSAWA ' C C READING TD - UT TABLE OPEN( 95, FILE = '/home/tamura/baytap/TDMUT.DAT', * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home1/users/tamura/baytap/TDMUT.DAT', C * STATUS = 'OLD' ) C OPEN( 95, FILE = '/home/scg/u101/common/baytap/TDMUT.DAT', C * STATUS = 'OLD' ) C READ(95,*) ( TDMUT(I), I = 1, 130 ) C CLOSE( 95 ) C RETURN END