C C*********************************************************************** C* * C* PROGRAM 'T A P W I N D': * C* --------------------- * C* * C* WRITTEN BY: M.M. NASSAR, JULY 1975. * C* * C* FOR PREDICTING THE TEMPERATURE AND PRESSURE VARIATIONS AT A * C* TIDE-GUAGE, FROM THE AVAILABLE INFORMATION AT ANY THREE POINTS * C* CLOSELY SURROUNDING THE TIDE-GUAGE: BY APPLYING THE LEAST-SQUARES * C* APPROXIMATION (INTERPOLATION) TECHNIQUES. * C*-*-*-*-*-*-**-*-*-*-*-*-*-*-*--**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*--**-*-*- C* MODIFIED IN THE FALL OF 1976 * C* TO INCORPORATE THE INVESTIGATION OF THE EFFECT OF GEOSTROPHIC WIND * C* ON THE ANALYSIS OF 'MEAN SEA LEVEL RECORDS'. * C*-*-*-*-*-*-**-*-*-*-*-*-*-*-*--**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*--**-*-*- C* * C*********************************************************************** C IMPLICIT REAL*8(A-H,O-Z) DIMENSION SLG(700),DTG(700),DPG(700),RD(4,700) DIMENSION GNAME(5) DIMENSION F1(700),F2(700),F3(700),FC(700) DIMENSION DUMVEC(700),IRD(12) DIMENSION RDBAR(3) DIMENSION VWE(700),VWN(700),MILES(8) DATA IAST/'*'/ C C SET-UP SOME CONVERSION CONSTANTS. C PI=3.141592653589793D 0 RHODEG=180.D 0/PI RHOSEC=RHODEG*3600.0D 0 C C READ-IN THE TOTAL NUMBER OF YEARS IN WHICH DATA IS AVAILABLE, C THE TOTAL NUMBER OF INPUTTED RIVER DISCHARGE DATA SETS, AND THE DATA REFERENCE C NUMBER THAT CONTAINS THE RELEVANT INFORMATION AT THE TIDE GUAGE. C WHICH IS STORED ON DISC. C READ (5,1) NYEARS,NRIVTL,IDG,IWIN,TGAP,DTIM,EH,EK 1 FORMAT (4I4,4X,4F10.3) C C COMPUTE THE SIZE OF THE VECTORS OF DATA (MSL, TEMP., ..ETC.). C NSIZE=NYEARS*12 C C READ-IN THE HEADER CARD OF THE TIDE GUAGE STATION: NUMBER; NAME; C LAT.; LONG.; STARTING YEAR OF OBS.; AND STARTING MONTH IN THIS YEAR. C READ(5,2) NUMG,(GNAME(I),I=1,5),PHIG,ALONGG,IYST,IMST 2 FORMAT(I6,4X,5A4,2F5.2,4X,I4,2X,I2) C C READ-IN THE LATITUDE (IN DEG.) OF THE TIDE-GAUGE LOCATION IN QUESTION. C IN FULL--- FOR THE PURPOSE OF GEOSTROPHIC WIND COMPUTATIONS. C ALSO** READ-IN THE AZIMUTH (IN DEGREES) OF THE SHORELINE AT THE TIDE- C GAUGE W.R.T. A A RIGHT-HANDED SYSTEM OF TWO-AXES: THE TANGENTIAL AND C THE NORMAL TO THE SHORELINE*** FKR THE PURPOSE OF GEOSTROPHIC WIND. C READ(5,570) ALATG,SHORAZ 570 FORMAT(F7.4,2X,F3.0) ALATGR=ALATG/RHODEG SHAZRD=SHORAZ/RHODEG C C READ-IN THE COORDINATES (LAT. & LONG.) OF THE THREE POINTS CLOSELY C SURROUNDING THE TIDE GUAGE, AND CONTAINING TEMPERATURE AND PRESSURE C DATA, TO BE USED FOR INTERPOLATION AT THE TIDE-GUAGE ITSELF. C C NOTE THAT IF ONE DATA SET (TEMP., PRESSURE) FROM ONLY ONE STATION C WILL BE CONSIDERED, A NUMBER GREATER THAN '10' SHOULD BE INPUT ON C THIS CARD ,ANYWHERE BETWEEN COLS. 1-5, AND THE REST OF THE CARD BLANK. C READ(5,3) PHI1,ALONG1,PHI2,ALONG2,PHI3,ALONG3 3 FORMAT(6F6.2) IF(PHI1.GT.10.0) GO TO 999 C C COMPUTE THE MEAN LAT. AND LONG. OF THE THREE POINTS FORMING THE TRIANGLE C PHI0=(PHI1+PHI2+PHI3)/3.D 0 ALONG0=(ALONG1+ALONG2+ALONG3)/3.D 0 C C COMPUTE THE CONVERSION FACTORS FROM CURVILINEAR (LAT. & LONG.) C TO LOCAL RECTANGULAR COORDINATES (X, Y). C R=6378.D 0 FACT1=R/RHODEG PHIRAD=PHI0/RHODEG FACT2=R*DCOS(PHIRAD)/RHODEG C C COMPUTE THE X (+VE NORTH) AND Y (+VE WEST) RECTANGULAR LOCAL COORDINATES C RELATIVE TO THE CENTRE OF THE TRIANGLE AS AN ORIGIN. C X1=(PHI1-PHI0)*FACT1 X2=(PHI2-PHI0)*FACT1 X3=(PHI3-PHI0)*FACT1 XG=(PHIG-PHI0)*FACT1 Y1=(ALONG1-ALONG0)*FACT2 Y2=(ALONG2-ALONG0)*FACT2 Y3=(ALONG3-ALONG0)*FACT2 YG=(ALONGG-ALONG0)*FACT2 C C COMPUTE THE DETERMINANT OF THE LEAST-SQUARES APPROXIMATION NORMAL C EQUATIONS MATRIX. C XISQ=X1**2+X2**2+X3**2 YISQ=Y1**2+Y2**2+Y3**2 XIYI=X1*Y1+X2*Y2+X3*Y3 DET=XISQ*YISQ-XIYI**2 C C COMPUTE THE THREE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL FOR C (TEMP., PRESSURE) AT THE TIDE-GUAGE LOCATION. C NOTE THAT THESE COEFFICIENTS ARE NOT DEPENDENT ON TIME ,AND ARE C CONSTANTS FOR EACH TRIANGLE UNDER CONSIDERATION. C D1=1.D 0/3.D 0+(YG*Y1*XISQ+XG*X1*YISQ-(XG*Y1+YG*X1)*XIYI)/DET D2=1.D 0/3.D 0+(YG*Y2*XISQ+XG*X2*YISQ-(XG*Y2+YG*X2)*XIYI)/DET D3=1.D 0/3.D 0+(YG*Y3*XISQ+XG*X3*YISQ-(XG*Y3+YG*X3)*XIYI)/DET C 999 CONTINUE C C******************************************************************R************ C SET-UP SOME CONSTANTS FOR THE PURPOSE OF GEOSTROPHIC-WIND COMPUTATIONS C C SET-UP A VALUE FOR THE DENSITY OF THE LOWER ATMOSPHERE ( G/CM**3 ) C DENS=0.001D 0 C C COMPUTE AN APPROXIMATE VALUE FOR THE ANGULAR VELOCITY OF THE EARTH.. C .. IN : RAD./SEC. ) C OMEGA=(2.0D 0*PI)/86400.0D 0 C C COMPUTE SOME FACTOR ,CONSTANT AT A PARTICULAR LOCATION, TO BE USED C IN COMPUTING THE GEOSTROPHIC WIND SPEED COMPONENTS. FROM THE PRESSURE C VFACT=2.0D 0*DENS*OMEGA*DSIN(ALATGR) C C SET-UP THE CONVERSION FACTOR TO CONVERT THE PRESSURE GRADIENTS FROM : C INCHES(HG)/KM. TO G*CM**-2*S**-2 C GRADNT=0.33863788D 0 C C COMPUTE THE CONVERSION FACTOR TO CONVERT THE WIND SPEED FROM:- C CM./SEC. TO M./SEC. C SPEEDF=1.0D 00/100.0D 00 C******************************************************************R************ C C INITIALIZE ZERO VALUES FOR THE ELEMENTS OF DATA VECTORS AT THE TIDE-GUAGE. C DO 10 I=1,NSIZE SLG(I)=0.0D 0 DTG(I)=0.0D 0 DPG(I)=0.0D 0 F1(I)=0.0D 0 F2(I)=0.0D 0 F3(I)=0.0D 0 DUMVEC(I)=0.0D 0 VWN(I)=0.0D 0 VWE(I)=0.0D 0 10 CONTINUE DO 20 I=1,4 DO 20 J=1,NSIZE RD(I,J)=0.0D 0 20 CONTINUE C C READ-IN AND STORE THE MONTHLY MEANS OF M.S.L. DATA AT THE TIDE- C GUAGE, FOR THE WHOLE PERIOD OF INVESTIGATION. C K=1 DO 30 I=1,NYEARS JACK=K+11 READ(5,4)ICODE,PHIG,ALONGG,IYEAR,(SLG(J),J=K,JACK) 4 FORMAT(A1,2F3.2,I4,12F5.2) K=K+12 30 CONTINUE C C READ-IN AND STORE THE MONTHLY MEANS OF THE RIVERS DISCHARGE DATA C SETS (UP TO THREE), NEAR THE TIDE-GUAGE. C C IT SHOULD BE NOTED THAT HERE THERE IS A POSSIBILITY TO ACCEPT MORE THAN C THREE SETS OF RIVER DISCHARGES- IN WHICH CASE TWO OR MORE SETS OF THEM C HAVE TO BE ADDED (LUMPED) TOGETHER. THE CHOSEN CRITERION HERE IS :- C EACH SET OF RIVER DISCHARGE DATA DECK HAS TO BE ENDED-UP WITH A BLANK C CARD, UNLESS THE CURRENT SET IN THE INPUT SEQUENCE IS REQUIRED TO BE C ADDED TO THE PREVIOUS ONE, WHERE THIS ENDING BLANK CARD OF THE CURRENT C SET MUST INCLUDE A NON-ZERO INTEGER IN COLUMN NUBER ONE. C IF(NRIVTL.EQ.0) GO TO 33 ILUMP=0 DO 41 II=1,NRIVTL KK=1 IROW=II-ILUMP DO 40 I=1,NYEARS JIM=KK+11 READ(5,5)ICODE,PHI,ALONG,IYEAR,(RD(IROW,J),J=KK,JIM) 5 FORMAT(A1,2F3.2,I4,12G5.0) KK=KK+12 40 CONTINUE READ(5,300) IADD 300 FORMAT(I1) IF(IADD.EQ.0) GO TO 41 ISTORE=IROW-1 DO 401 J=1,NSIZE 401 RD(ISTORE,J)=RD(ISTORE,J)+RD(IROW,J) ILUMP=ILUMP+1 41 CONTINUE NRIVER=NRIVTL-ILUMP C C CLEANING THE DUMMY VALUES STORED INTO THE RIVER DISCHARGES ARRAY , C WHICH EXCEED THE ACTUAL NUMBER OF DESIRED RIVER DISCHARGES, AND THEN C RE-INITIALIZE THEM BACK TO ZEROS. C IF(NRIVER.LT.3) GO TO 625 GO TO 650 625 ISAMY=NRIVER+1 DO 675 I=ISAMY,4 DO 675 J=1,NSIZE 675 RD(I,J)=0.0D 0 C 650 CONTINUE 33 CONTINUE C C PRINT-OUT THE REQUIRED INFORMATION AT THE TIDE-GUAGE STATION. C WRITE(6,15) IDG 15 FORMAT('1',//,5X,'GENERAL INFORMATION ABOUT THE TIDE-GUAGE :-',/,5 1X,'----------------------------------------',//,5X,'DATA SET REFER 2ENCE NUMBER ON DISC IS =',I4/) WRITE(6,16) NUMG,(GNAME(I),I=1,5) 16 FORMAT(5X,'TIDE-GUAGE NUMBER IS =',I8//,5X,'TIDE-GUAGE NAME IS = 1',5A4//) WRITE(6,17) IYST,IMST,NSIZE 17 FORMAT(5X,'DATA STARTING YEAR IS =',I5//,5X,'DATA STARTING MONTH O 1F THE ABOVE YEAR IS =',I3//,5X,'SIZE OF ARRAYS OF DIFFERENT TYPES 2OF DATA IS =',I4//) WRITE(6,501) PHI1,ALONG1,PHI2,ALONG2,PHI3,ALONG3 501 FORMAT(//,5X,'LAT. & LONG. OF STATION NO. 1 ARE =',2F6.2,/,5X,'LAT 1. & LONG. OF STATION NO. 2 ARE =',2F6.2,/,5X,'LAT. & LONG. OF STAT 2ION NO. 3 ARE =',2F6.2/) WRITE (6,601) TGAP,DTIM,EH,EK 601 FORMAT (//5X,'GAP IN DATA SIGNALLED BY:',F7.0//5X,'TEMP. DATA PHAS @E SHIFT =',F8.4,/5X,'TEMP. DATA LINEAR MODEL PARAM. =',2F10.4) IF(PHI1.GT.10.0) WRITE(6,502) 502 FORMAT(///,10X,'ONLY ONE SET OF TEMP. & PRESSURE IS INPUTTED'/,13X 1,'**NO INTERPOLATION IS REQUIRED**'//) WRITE(6,18) 18 FORMAT('1',//,5X,'VECTOR OF MONTHLY MEAN SEA LEVEL IS:') K=1 DO 90 I=1,NYEARS IYR=IYST+I-1 KL=K+11 WRITE(6,19) IYR,(SLG(J),J=K,KL) 19 FORMAT (3X,I4,3X,12F8.2) K=K+12 90 CONTINUE IF (NRIVTL.EQ.0) GO TO 593 WRITE(6,21) 21 FORMAT(/////,5X,'FIRST SET OF MONTHLY MEAN RIVER DISCHARGE IS:') K=0 DO 100 I=1,NYEARS IYR=IYST+I-1 DO 101 J=1,12 K=K+1 IRD(J)=RD(1,K) 101 CONTINUE WRITE(6,22) IYR,(IRD(JOE),JOE=1,12) 22 FORMAT (3X,I4,3X,12I8) 100 CONTINUE WRITE(6,23) 23 FORMAT(//,5X,'SECOND SET OF MONTHLY MEAN RIVER DISCHARGE IS:') K=0 DO 110 I=1,NYEARS IYR=IYST+I-1 DO 111 J=1,12 K=K+1 IRD(J)=RD(2,K) 111 CONTINUE WRITE(6,22) IYR,(IRD(JOE),JOE=1,12) 110 CONTINUE WRITE(6,24) 24 FORMAT(//,5X,'THIRD SET OF MONTHLY MEAN RIVER DISCHARGE IS:') K=0 DO 120 I=1,NYEARS IYR=IYST+I-1 DO 121 J=1,12 K=K+1 IRD(J)=RD(3,K) 121 CONTINUE WRITE(6,22) IYR,(IRD(JOE),JOE=1,12) 120 CONTINUE C 593 IF(PHI1.GT.10.0) GO TO 888 C C CALL THE SUBROUTINE 'F U N I N T', FOR THE C PREDICTION OF THE MONTHLY MEANS OF TEMPERATURE AT THE C TIDE-GUAGE, FROM THE INFORMATION AVAILABLE AT THE SURROUNDING THREE C POINTS, BY THE LEAST-SQUARES APPROXIMATION TECHNIQUES. C WRITE(6,70) 70 FORMAT('1'//,' THE FOLLOWING PLOT AND ITS ASSOCIATED LEGEND 1 PERTAIN TO THE ***TEMPERATURE***') C CALL FUNINT(D1,D2,D3,NYEARS,NSIZE,F1,F2,F3,DTG) C C CALL THE SUBROUTINE 'F U N I N T', FOR THE C PREDICTION OF THE MONTHLY MEANS OF THE PRESSURE AT THE TIDE- C GUAGE STATION FROM THE AVAILABLE INFORMATION AT THE SURROUNDING THREE C STATIONS, BY THE LEAST-SQUARES APPROXIMATION TECHNIQUES. C WRITE(6,80) 80 FORMAT('1'//,' THE FOLLOWING PLOT AND ITS ASSOCIATED LEGEND 1PERTAINS TO THE ***PRESSURE***') C CALL FUNINT(D1,D2,D3,NYEARS,NSIZE,F1,F2,F3,DPG) C GO TO 777 888 CONTINUE C C READ-IN THE ONLY ONE AVAILABLE SET OF OBSERVED TEMPERATURE, AT ONE C STATION, AND HENCE STORE IT AS IF IT WAS TAKEN AT THE TIDE-GUAGE. C K=1 DO 210 I=1,NYEARS JOE=K+11 READ(5,211)ICODE,PHI,ALONG,IYEAR,(DTG(J),J=K,JOE) 211 FORMAT(A1,2F3.2,I4,12F5.1) K=K+12 210 CONTINUE C C READ-IN THE ONLY ONE AVAILABLE SET OF OBSERVED PRESSURE AT ONE C STATION, AND HENCE STORE IT AS IF IT WAS TAKEN AT THE TIDE-GUAGE. C K=1 DO 212 I=1,NYEARS JOE=K+11 READ(5,211)ICODE,PHI,ALONG,IYEAR,(DPG(J),J=K,JOE) K=K+12 212 CONTINUE C 777 CONTINUE IF (IWIN.EQ.0) GO TO 599 C C **** READ IN OBSERVED MONTHLY TOTAL WIND MILEAGES FOR 8 OR 16 DIRECTIONS C ND=8 ANG=PI*0.25D0 CF=5280.D0*0.3048D0/864.D2 DO 499 I=1,NYEARS DO 498 J=1,12 KD=1 ITOT=0 595 READ (5,299) IYEAR,MON,KARD,MILES,MTOT 299 FORMAT (I4,I3,A1,2X,8I5,I10) WRITE (6,602) IYEAR,MON,KARD,MILES,MTOT,ND,KD 602 FORMAT (10X,I4,I3,A1,2X,8I5,I10,2I5) IF (KARD.NE.IAST) GOTO 598 ND=16 ANG=ANG/2.D0 598 K=(IYEAR-IYST)*12+MON-IMST+1 DO 497 ID=1,8 AZ=(ID-1)*ANG+(KD-1)*PI VWN(K)=VWN(K)+MILES(ID)*DCOS(AZ) VWE(K)=VWE(K)+MILES(ID)*DSIN(AZ) ITOT=ITOT+MILES(ID) 497 CONTINUE IF (ND.EQ.8) GO TO 597 IF (KD.EQ.2) GO TO 596 KD=2 GO TO 595 597 IF (ITOT.EQ.MTOT) GO TO 596 WRITE (6,298) IYEAR,MON,MILES,MTOT 298 FORMAT (//5X,'ERROR IN WIND MILEAGES:'/5X,I4,I3,3X,8I5,I10) C **** FIND AV. MILES PER HOUR AND CONVERT TO M/S 596 CONV=CF/NDAYS(IYEAR,MON) VWN(K)=VWN(K)*CONV VWE(K)=VWE(K)*CONV 498 CONTINUE 499 CONTINUE GO TO 594 C C READ-IN AND STORE THE PREDICTED MONTHLY MEANS OF THE NORTH AND EAST C GRADIENTS OF THE ATMOSPHERIC PRESSURE IN THE REGION SURROUNDING THE C TIDE-GUAGE OF INTEREST- AS COMPUTED FROM A SEPARATE PROGRAMM. C 599 K=0 DO 503 I=1,NYEARS READ(5,505) ICODE,IYEAR,(IRD(JIM),JIM=1,12) 505 FORMAT(A2,1X,I4,1X,12I6) DO 504 J=1,12 K=K+1 504 VWN(K)=IRD(J)*1.0D-07 503 CONTINUE K=0 DO 506 I=1,NYEARS READ(5,505) ICODE,IYEAR,(IRD(JIM),JIM=1,12) DO 507 J=1,12 K=K+1 507 VWE(K)=IRD(J)*1.0D-07 506 CONTINUE C C COMPUTE AND STORE THE NORTH AND EAST COMPONENTS OF GEOSTROPHIC WIND SPEED C FROM THE ESTIMATED PRESSURE GRADIENTS IN THOSE DIRECTIONS ,RESPECTIVELY. C ** IN M./SEC. UNITS ** C DO 508 J=1,NSIZE DPDX=VWE(J) DPDY=VWN(J) VN=+(DPDX/VFACT)*GRADNT*SPEEDF VE=-(DPDY/VFACT)*GRADNT*SPEEDF VWN(J)=VN VWE(J)=VE 508 CONTINUE C C PRINT-OUT THE ESTIMATED COMPONENTS OF GEOSTROPHIC WIND SPEED ** IN M./SEC. ** C 594 WRITE(6,509) 509 FORMAT('1'//,5X,'VECTOR OF MONTHLY MEANS OF PREDICTED NORTH COMPON 1ENT OF THE GEOSTROPHIC-WIND SPEED IN M./SEC. :') K=1 DO 510 I=1,NYEARS IYR=IYST+I-1 KL=K+11 WRITE(6,19) IYR,(VWN(J),J=K,KL) 510 K=K+12 WRITE(6,511) 511 FORMAT(////,5X,'VECTOR OF MONTHLY MEANS OF PREDICTED EAST COMPONEN 1T OF THE GEOSTROPHIC-WIND SPEED IN M./SEC. :') K=1 DO 512 I=1,NYEARS IYR=IYST+I-1 KL=K+11 WRITE(6,19) IYR,(VWE(J),J=K,KL) 512 K=K+12 C C CALL THE SUBROUTINE 'T S P L O T' ,TO COMPUTE AND PRINT-OUT THE MEAN C MAXIMUM AND MINIMUM VALUES OF THE PREDICTED TEMPERATURE AT THE TIDE- C GUAGE, AS WELL AS PLOTTING THE CORRESPONDING TIME-SERIES. C WRITE(6,460) 460 FORMAT('1',//,5X,'***THE FOLLOWING PLOT IS FOR THE PREDICTED TEMPE 1RATURE SERIES AT THE TIDE-GUAGE***'/) C CALL TSPLOT(DTG,NSIZE,TGBAR) C C COMPUTE AND PLOT THE ''NEW'' TEMPERATURE SERIES:- C WRITE(6,541) 541 FORMAT('1'///,5X,'THE FOLLOWING PLOT AND RESULTS PERTAIN TO THE NE 1W GENERATED TEMPERATURE SERIES'/) NN=0 DO 411 I=1,NSIZE IF (DTG(I).EQ.TGAP) GO TO 191 NN=NN+1 IF (I.NE.NSIZE) GO TO 411 191 IF (NN.NE.0) GO TO 195 DTG(I)=0.0D0 GO TO 411 195 IS=I-NN IF (NN.LT.12) GO TO 192 IF (MOD(NN,2).NE.0) NN=NN-1 CALL FFS (DTG(IS),NN,FC) CALL FINT (FC,NN,DTIM,DTG(IS)) DTG(IS)=0.D0 GO TO 193 192 DO 412 J=IS,I DTG(J)=0.D0 412 CONTINUE 193 NN=0 411 CONTINUE DO 413 I=1,NSIZE IF (DTG(I).EQ.TGAP) GO TO 413 DTG(I)=DTG(I)*EH+EK 413 CONTINUE C CALL TSPLOT(DTG,NSIZE,TGBAR) C C C C COMPUTE THE VARIATIONS IN TEMPERATURE VALUES (FROM THEIR MEAN VALUE) C FOR THE PREDICTED VALUES AT THE TIDE-GUAGE. C C 0 D0.0=RABGT DO 216 I=1,NSIZE IF(DTG(I).EQ.0.0) GO TO 216 DTG(I)=DTG(I)-TGBAR 216 CONTINUE C C PRINT-OUT THE 'PREDICTED' VARIATIONS IN TEMPERATURE AT THE TIDE-GUAGE. C WRITE(6,25) 25 FORMAT('1'//,5X,'PREDICTED MONTHLY MEAN TEMPERATURE VARIATIONS AT 1THE TIDE-GUAGE ARE:') K=1 DO 130 I=1,NYEARS KL=K+11 WRITE(6,26) (DTG(J),J=K,KL) 26 FORMAT(10X,12F8.2) K=K+12 130 CONTINUE C C CALL THE SUBROUTINE 'T S P L O T' ,TO COMPUTE AND PRINT-OUT THE MEAN C MAXIMUM AND MINIMUM VALUES OF THE PREDICTED PRESSURE AT THE TIDE- C GUAGE, AS WELL AS PLOTTING THE CORRESPONDING TIME-SERIES. C WRITE(6,465) 465 FORMAT('1',//,5X,'***THE FOLLOWING PLOT IS FOR THE PREDICTED PRESS 1URE SERIES AT THE TIDE-GUAGE***'/) C CALL TSPLOT(DPG,NSIZE,PGBAR) C C C COMPUTE THE VARIATIONS IN THE PRESSURE VALUES (FROM THEIR MEAN VALUE) C FOR THE PREDICTED VALUES AT THE TIDE-GUAGE. C DO 221 I=1,NSIZE IF(DPG(I).EQ.0.0) GO TO 221 DPG(I)=DPG(I)-PGBAR 221 CONTINUE C C C PRINT-OUT THE 'PREDICTED' VARIATIONS IN PRESSURE AT THE TIDE-GUAGE. C WRITE(6,27) 27 FORMAT('1'//,5X,'PREDICTED MONTHLY MEAN PRESSURE VARIATIONS AT THE 1 TIDE-GUAGE ARE:') K=1 DO 140 I=1,NYEARS KL=K+11 WRITE(6,26) (DPG(J),J=K,KL) K=K+12 140 CONTINUE IF(NRIVTL.EQ.0) GO TO 32 C C CALL THE SUBROUTINE 'T S P L O T' ,TO COMPUTE AND PRINT-OUT THE MEAN C MAXIMUM AND MINIMUM VALUES OF THE SETS OF ACCEPTED RIVER DISCHARGE C SERIES TO BE UNDER CONSIDERATION (UP TO 3 SETS), AS WELL AS PLOTTING C THE CORRESPONDING TIME SERIES FOR EACH SET. C DO 201 I=1,NRIVER WRITE(6,470) I 470 FORMAT('1',//,5X,'***THE FOLLOWING PLOT IS FOR SET NO.',I2,' OF RI 1VER DISCHARGE SERIES***'/) DO 202 J=1,NSIZE 202 DUMVEC(J)=RD(I,J) C CALL TSPLOT(DUMVEC,NSIZE,DISBAR) RDBAR(I)=DISBAR C 201 CONTINUE C C COMPUTE THE VARIATIONS IN THE MONTHLY RIVER DISCHARGE (FROM THE CORRESPONDING C MEAN VALUE) AT EVERY AVAILABLE RIVER LOCATION. C DO 203 I=1,NRIVER DO 203 J=1,NSIZE IF(RD(I,J).EQ.0.0) GO TO 203 RD(I,J)=RD(I,J)-RDBAR(I) 203 CONTINUE C C PRINT-OUT THE COMPUTED VARIATIONS IN RIVER DISCHARGE VALUES AT EACH C GIVEN RIVER LOCATION. C WRITE(6,207) 207 FORMAT('1') DO 204 I=1,NRIVER WRITE(6,205) I 205 FORMAT(/,/,/,5X,'COMPUTED VARIATIONS IN RIVER DISCHARGE VALUES FOR 1 **SET NO.',I2,'** ARE:'/) K=0 DO 204 JJ=1,NYEARS DO 206 J=1,12 K=K+1 IRD(J)=RD(I,K) 206 CONTINUE WRITE(6,22) (IRD(JOE),JOE=1,12) 204 CONTINUE 32 CONTINUE C C CALL THE SUBROUTINE 'T S P L O T' TO COMPUTE AND PRINT-OUT THE MEAN, C MAX. & MIN. VALUES OF THE PREDICTED NORTH COMPONENT OF THE GEOSTROPHIC- C WIND SPEED AT THE TIDE-GUAGE LOCATION- AND TO PLOT THE CORRESPONDING C TIME SERIES. C WRITE(6,561) 561 FORMAT('1'///,5X,'THE FOLLOWING PLOT IS FOR THE PREDICTED NORTH CO 1MPONENT OF THE GEOSTROPHIC-WIND SPEED'/) C CALL TSPLOT(VWN,NSIZE,VWNBAR) C C CALL THE SUBROUTINE 'T S P L O T' TO COMPUTE AND PRINT-OUT THE MEAN, C MAX. & MIN. VALUES OF THE PREDICTED EAST COMPONENT OF THE GEOSTROPHIC- C WIND SPEED AT THE TIDE-GUAGE LOCATION- AND TO PLOT THE CORRESPONDING C TIME SERIES. C WRITE(6,562) 562 FORMAT('1'///,5X,'THE FOLLOWING PLOT IS FOR THE PREDICTED EAST COM 1PONENT OF THE GEOSTROPHIC-WIND SPEED'/) C CALL TSPLOT(VWE,NSIZE,VWEBAR) C C C COMPUTE THE TANGENTIAL AND THE NORMAL 'TO THE SHORELINE' COMPONENT- C SQUARED OF THE GEOSTROPHIC WIND - FROM THE PREVIOUSLY COMPUTED NORTH C AND EAST COMPONENTS,***, AND STORE THEM IN PLACE OF THE N & E COMPONENTS. C DO 543 I=1,NSIZE TRICK1=VWN(I) TRICK2=VWE(I) VWN(I)=TRICK1*DCOS(SHAZRD)+TRICK2*DSIN(SHAZRD) VWE(I)=-TRICK1*DSIN(SHAZRD)+TRICK2*DCOS(SHAZRD) C SQUARE THE TANGENTIAL COMPONENT: WITHOUT THE SIGN BEEN TAKEN INTO ACCOUNT C AND SQURE THE NORMAL COMPONENT WITH THE SIGN OF IT APPEARS INTO THE SQUARE. VWN(I)=VWN(I)*VWN(I) VWE(I)=VWE(I)*DABS(VWE(I)) 543 CONTINUE C C PLOT THE TANGENTIAL AND THE NORMAL COMPONENTS TIME-SERIES. C WRITE(6,544) 544 FORMAT('1'//,5X,'THE FOLLOWING PLOT AND RESULTS PERTAIN TO THE TAN 1GENTIAL COMPONENT-SQUARED OF THE GEOSTROPHIC-WIND'/) C CALL TSPLOT(VWN,NSIZE,VWNBAR) C C C COMPUTE THE VARIATIONS OF WIND-SPEED TANG.-COMPONENT (FROM ITS MEAN VALUE). C VWNBAR=0.0D 0 DO 513 I=1,NSIZE IF(VWN(I).EQ.0.0) GO TO 513 VWN(I)=VWN(I)-VWNBAR 513 CONTINUE C C PRINT-OUT THE PREDICTED VARIATIONS OF WIND-SPEED TANG.-COMPONENT. C WRITE(6,514) 514 FORMAT('1',//,5X,'PREDICTED MONTHLY MEAN VARIATIONS OF TANG.-COMPO 1NENT OF GEOSTROPHIC-WIND SPEED :') K=1 DO 515 I=1,NYEARS IYR=IYST+I-1 KL=K+11 WRITE(6,19) IYR,(VWN(J),J=K,KL) 515 K=K+12 C WRITE(6,545) 545 FORMAT('1',///,5X,'THE FOLLOWING PLOT AND RESULTS PERTAIN TO THE N 1ORMAL-COMPONENT-SQUARED OF THE GEOSTROPHIC WIND'/) C CALL TSPLOT(VWE,NSIZE,VWEBAR) C C C COMPUTE THE VARIATIONS OF WIND-SPEED NORM-COMPONENT (FROM ITS MEAN VALUE). C VWEBAR=0.0D 0 DO 516 I=1,NSIZE IF(VWE(I).EQ.0.0) GO TO 516 VWE(I)=VWE(I)-VWEBAR 516 CONTINUE C C PRINT-OUT THE VARIATIONS OF NORM-COMPONENT OF WIND-SPEED C WRITE(6,517) 517 FORMAT(///,5X,'PREDICTED MONTHLY MEAN VARIATIONS OF NORM-COMPONENT 1 OF GEOSTROPHIC-WIND SPEED :') K=1 DO 518 I=1,NYEARS IYR=IYST+I-1 KL=K+11 WRITE(6,19) IYR,(VWE(J),J=K,KL) 518 K=K+12 C C******************************************************************************* C C TRANSFER THE TANGENTIAL AND NORMAL COMPONENTS OF THE GEOSTROPHIC WIND TO BE C STORED WITHIN THE RIVER-DISCHARGE MATRIX***FOR THE PURPOSE OF PUTTING C THEM O NDISC FOR THE M.S.L. ANALYSIS PROGRAM. C NOTE: THIS IS POSSIBLE HERE SINCE ALL THE PORTS DEALT WITH IN EASTERN C CANADA HAVE AT THE MOST ONLY ONE RIVER DISCHARGE CONSIDERED. C KIM=2 IF(NRIVTL.EQ.0) KIM=1 KISS=KIM+1 DO 547 I=1,NSIZE RD(KIM,I)=VWN(I) RD(KISS,I)=VWE(I) 547 CONTINUE IF (IDG.EQ.0) GO TO 194 C C******************************************************************************* C C WRITE THE NECESSARY INFORMATION AT THE TIDE GUAGE ON THE DISC-- C TO BE STORED FOR USE BY ANOTHER PROGRAM. C REWIND IDG WRITE(IDG,301) NUMG 301 FORMAT(I6) DO 302 I=1,5 WRITE(IDG,303) GNAME(I) 303 FORMAT(A4) 302 CONTINUE WRITE(IDG,304) IYST 304 FORMAT(I4) WRITE(IDG,305) IMST 305 FORMAT(I2) WRITE(IDG,306) NSIZE 306 FORMAT(I3) C DO 307 I=1,NSIZE WRITE(IDG,308) SLG(I) 308 FORMAT(F10.5) 307 CONTINUE DO 309 I=1,NSIZE 309 WRITE(IDG,308) DTG(I) DO 310 I=1,NSIZE 310 WRITE(IDG,308) DPG(I) C DO 311 I=1,3 DO 311 J=1,NSIZE WRITE(IDG,312) RD(I,J) 312 FORMAT(F10.5) 311 CONTINUE C 194 STOP END C C*************************************************************************** SUBROUTINE FUNINT(D1,D2,D3,NY,NS,F1,F2,F3,FG) IMPLICIT REAL*8(A-H,O-Z) DIMENSION F1(NS),F2(NS),F3(NS),FG(NS) C C READ-IN AND STORE THE VALUES OF THE OBSERVED FUNCTION AT THE THREE C STATIONS FORMING A TRIANGLE AROUND THE TIDE-GUAGE , TO BE USED FOR C INTERPOLATION. C K=1 DO 100 I=1,NY JOE=K+11 READ(5,6)ICODE,PHI,ALONG,IYEAR,(F1(J),J=K,JOE) 6 FORMAT(A1,2F3.2,I4,12F5.1) K=K+12 100 CONTINUE K=1 DO 200 I=1,NY JOE=K+11 READ(5,6)ICODE,PHI,ALONG,IYEAR,(F2(J),J=K,JOE) K=K+12 200 CONTINUE K=1 DO 300 I=1,NY JOE=K+11 READ(5,6)ICODE,PHI,ALONG,IYEAR,(F3(J),J=K,JOE) K=K+12 300 CONTINUE C C INITIALIZE ZERO VALUES FOR THE FUNCTIONAL VALUES AT THE TIDE-GUAGE. C DO 75 I=1,NS FG(I)=0.0D 0 75 CONTINUE C C COMPUTE THE VALUES OF THE INTERPOLATED FUNCTION AT THE TIDE-GUAGE. C DO 400 I=1,NS FG(I)=D1*F1(I)+D2*F2(I)+D3*F3(I) 400 CONTINUE C C CALLING THE SUBROUTINE 'F P L O T', TO PLOT THE THREE TIME SERIES C OF OBSERVED FUNCTIONAL VALUES AT THE THREE STATIONS ,AND ALSO THE C PREDICTED VALUES AT THE TIDE-GUAGE- TO VISUALIZE ITS VARIATIONS. C CALL FPLOT(F1,F2,F3,FG,NS) C RETURN END C C*************************************************************************** SUBROUTINE FPLOT(F1,F2,F3,FG,NS) IMPLICIT REAL*8(A-H,O-Z) DIMENSION F1(NS),F2(NS),F3(NS),FG(NS) DIMENSION IPLOT(30),IFRAME(30) DOUBLE PRECISION DABS,DMAX1,DMIN1,DFLOAT DATA IBLANK/' '/,IDASH/'|'/,IPLUS/'M'/ DATA IF1/'.'/,IF2/'+'/,IF3/'X'/,IFG/'*'/ C C SPECIFY THE DESIRED LATERAL SPACE FOR THE PLOT ,IN TERMS OF THE C EQUIVALENT NUMBER OF CHARACTERS ON THE LINE PRINTER 'NCOL'. C NCOL=30 C C SET-UP VALUES OF 'M' FOR THE ELEMENTS OF THE ARRAY 'IFRAME', C TO BE USED IN FRAMING THE PLOT. C DO 9 I=1,NCOL IFRAME(I)=IPLUS 9 CONTINUE C C INITIALIZE SOME REASONABLE VALUES - TO HELP EXTRACTING THE MAXIMUM C AND MINIMUM OBSERVED (OR PREDICTED) VALUES OF THE FUNCTION UNDER C CONSIDERATION- THIS IS REQUIRED MAINLY FOR PLOTTING ROUTINES. C FMAX=0.0D 0 FMIN=1.0D 06 C C EXTRACT THE MAX. & MIN. VALUES OF THE FUNCTION THAT OCCURRED AT ANY C OF THE THREE STATIONS OR THE TIDE-GUAGE. C DO 500 I=1,NS FMAX=DMAX1(FMAX,F1(I)) FMAX=DMAX1(FMAX,F2(I)) FMAX=DMAX1(FMAX,F3(I)) FMAX=DMAX1(FMAX,FG(I)) IF(F1(I).EQ.0.0) GO TO 10 FMIN=DMIN1(FMIN,F1(I)) 10 IF(F2(I).EQ.0.0) GO TO 20 FMIN=DMIN1(FMIN,F2(I)) 20 IF(F3(I).EQ.0.0) GO TO 30 FMIN=DMIN1(FMIN,F3(I)) 30 IF(FG(I).EQ.0.0) GO TO 500 FMIN=DMIN1(FMIN,FG(I)) 500 CONTINUE C C COMPUTE THE MEAN VALUE OF THE INTERPOLATED (PREDICTED) FUNCTION C AT THE TIDE-GUAGE. C C NOTE:- THE ZERO VALUES WHICH REPRESENT MISSING DATA ARE EXCLUDED C FROM THE COMPUTATIONS. C NZEROG=0 SUMFG=0.0D 0 DO 600 I=1,NS IF(FG(I).EQ.0.0) NZEROG=NZEROG+1 SUMFG=SUMFG+FG(I) 600 CONTINUE NN=NS-NZEROG FGBAR=SUMFG/DFLOAT(NN) C C INITIALIZE THE PLOTTING ARRAY. C DO 60 I=1,NCOL IPLOT(I)=IBLANK 60 CONTINUE IPLOT(1)=IPLUS IPLOT(NCOL)=IPLUS C X=DFLOAT(NCOL-1)/(FMAX-FMIN) Y=FMIN*X-1.0 KBAR=FGBAR*X-Y IPLOT(KBAR)=IDASH C C PRINT-OUT A LEGEND FOR THE DIFFERENT FUNCTIONS IN THE PLOT C WRITE(6,5) 5 FORMAT(/,5X,'LEGEND :',2X,'.-REPRESENTS THE FUNCTIONAL VALUES AS O 1BSERVED AT STATION NO. 1 ',/,15X,'+-REPRESENTS THE FUNCTIONAL VALU 2ES AS OBSERVED AT STATION NO. 2',/,15X,'X-REPRESENTS THE FUNCTIONA 3L VALUES AS OBSERVED AT STATION NO. 3',/,15X,'*-REPRESENTS THE FUN 4CTIONAL VALUES AS PREDICTED AT THE TIDE-GUAGE',/,15X,'M-REPRESENTS 5 THE FRAME AROUND THE PLOT') C C PRINT-OUT THE MEAN, MAX. & MIN. VALUES OF THE TIME SERIES. C WRITE(6,15) NS,FMAX,FMIN,FGBAR 15 FORMAT(/,/,10X,'TOTAL NUMBER OF OBSERVED FUNCTIONAL VALUES IS =',I 15,/,10X,'MAXIMUM OBSERVED FUNCTIONAL VALUE IS =',F8.3,/,10X,'MINIM 1UM OBSERVED FUNCTIONAL VALUE IS =',F8.3,/,10X,'MEAN PREDICTED FUNC 2TIONAL VALUE AT THE TIDE-GUAGE IS =',F8.3//) C C COMPUTE AND PRINT-OUT THE SCALE OF THE PLOT, BY CONSIDERING THE LINE C PRINTER # 1403 (AT U.N.B.), WHICH PRINTS 10 CHARACTERS PER INCH. C SINCH=10.0 SCALE=DABS(FMAX-FMIN)/(DFLOAT(NCOL)/SINCH) WRITE(6,70) SCALE 70 FORMAT(7X,'NO.',5X,'SCALE OF THE PLOT IS: 1 INCH REPRESENTS',E10.3 1,' UNITS') C C PLOT THE GIVEN TIME SERIES, WITH ITS FUNCTIONAL VALUES. C WRITE(6,40) (IFRAME(J),J=1,NCOL) C DO 72 I=1,NS KF1=F1(I)*X-Y KF2=F2(I)*X-Y KF3=F3(I)*X-Y KFG=FG(I)*X-Y IF(KF1.LT.1) KF1=1 IF(KF2.LT.1) KF2=1 IF(KF3.LT.1) KF3=1 IF(KFG.LT.1) KFG=1 IF(KF1.GT.NCOL) KF1=NCOL IF(KF2.GT.NCOL) KF2=NCOL IF(KF3.GT.NCOL) KF3=NCOL IF(KFG.GT.NCOL) KFG=NCOL C IPLOT(KF1)=IF1 IPLOT(KF2)=IF2 IPLOT(KF3)=IF3 IPLOT(KFG)=IFG C WRITE(6,80) I,(IPLOT(J),J=1,NCOL) C C RE-INITIALIZE THE PLOTTING ARRAY. C IPLOT(KF1)=IBLANK IPLOT(KF2)=IBLANK IPLOT(KF3)=IBLANK IPLOT(KFG)=IBLANK IPLOT(1)=IPLUS IPLOT(NCOL)=IPLUS IPLOT(KBAR)=IDASH 72 CONTINUE WRITE(6,40) (IFRAME(J),J=1,NCOL) 40 FORMAT(15X,30A1) 80 FORMAT(6X,I4,5X,30A1) RETURN SUBROUTINE FFS (F,NN,A) C C******************************************************************************* END SUBROUTINE TSPLOT(F,N,FBAR) C C******************************************************************************* C C SUBROUTINE 'T S P L O T': C FOR COMPUTING AND PRINTING THE MAXIMUM, C MINIMUM AND MEAN VALUES OF A GIVEN TIME SERIES. ALSO IT PLOTS THE C CORRESPONDING TIME SERIES ALONG WITH THE FUNCTIONAL VALUES AND C THE SCALE OF THE PLOT. C 8791 BEF NOSREDNA -- NOISICERP ELGNIS ROF DEIFIDOM C C******************************************************************************* C IMPLICIT REAL*8(A-H,O-Z) C DOUBLE PRECISION DABS,DMAX1,DMIN1,DFLOAT C DIMENSION F(N) DIMENSION IPLOT(30),IFRAME(30) DATA IBLANK/' '/,ISTAR/'*'/,IPLUS/'+'/,IDASH/'|'/ C C SPECIFY THE DESIRED LATERAL SPACE FOR THE PLOT ,IN TERMS OF THE C EQUIVALENT NUMBER OF CHARACTERS ON THE LINE PRINTER 'NCOL'. C NCOL=30 C C SET-UP VALUES OF '+' FOR THE ELEMENTS OF THE ARRAY 'IFRAME', C TO BE USED IN FRAMING THE PLOT. C DO 9 I=1,NCOL IFRAME(I)=IPLUS 9 CONTINUE C C COMPUTE THE MEAN, MAXIMUM AND MINIMUM FUNCTIONAL VALUES FOR THE C GIVEN TIME-SERIES (REQUIRED TO BE PLOTTED). C C NOTE:- THE ZERO VALUES WHICH REPRESENT MISSING DATA ARE EXCLUDED C FROM THE COMPUTATIONS. C NZEROF=0 SUMF=0.0 DO 1 I=1,N IF(F(I).EQ.0.0) NZEROF=NZEROF+1 C C 1+FOREZN=FOREZN )9.99.QE.)I(F(FI SUMF=SUMF+F(I) 1 CONTINUE NSUMF=N-NZEROF C )FMUSN(TAOLF /FMUS=RABF FBAR=SUMF/DFLOAT(NSUMF) C C FMAX=0.0D 0 C C 03E0.1-=XAMF FMIN=1.0D 10 C C 03E0.1=NIMF DO 2 I=1,N IF(F(I).EQ.0.0) GO TO 2 C C 2 OT OG )9.99.QE.)I(F(FI FMAX=DMAX1(FMAX,F(I)) C C ))I(F,XAMF(1XAMA=XAMF FMIN=DMIN1(FMIN,F(I)) C C ))I(F,NIMF(1NIMA=NIMF 2 CONTINUE C C INITIALIZE THE PLOTTING ARRAY. C DO 10 I=1,NCOL IPLOT(I)=IBLANK 10 CONTINUE IPLOT(1)=IPLUS IPLOT(NCOL)=IPLUS C X=DFLOAT(NCOL-1)/(FMAX-FMIN) C C )NIMF-XAMF(/)1-LOCN(TAOLF =X Y=FMIN*X-1.0 KBAR=FBAR*X-Y IPLOT(KBAR)=IDASH C C PRINT-OUT THE MEAN, MAX. & MIN. VALUES OF THE TIME SERIES. C WRITE(6,60) FMAX,FMIN,FBAR 60 FORMAT(5X,'MAX. FUNCTIONAL VALUE IN THE SERIES = ',E11.4,/,5X,'MIN 1. FUNCTIONAL VALUE IN THE SERIES = ',E11.4,/,5X,'MEAN FUNCTIONAL V 2ALUE IN THE SERIES = ',E11.4/) C C COMPUTE AND PRINT-OUT THE SCALE OF THE PLOT, BY CONSIDERING THE LINE C PRINTER # 1403 (AT U.N.B.), WHICH PRINTS 10 CHARACTERS PER INCH. C SINCH=10.0 SCALE=DABS(FMAX-FMIN)/(DFLOAT(NCOL)/SINCH) C C )HCNIS/)LOCN(TAOLF (/)NIMF-XAMF(SBA =ELACS WRITE(6,70) SCALE 70 FORMAT(6X,'NO.',5X,'VALUE',5X,'SCALE OF THE PLOT IS: 1 INCH REPRES 1ENTS',E10.3,' UNITS') C C PLOT THE GIVEN TIME SERIES, WITH ITS FUNCTIONAL VALUES. C WRITE(6,40) (IFRAME(J),J=1,NCOL) C INITIALIZE ANY ARBITRARY INTEGER 'KF' FOR THE POSITION OF THE '*' KF=3 DO 30 I=1,N IF(F(I).EQ.0.0) GO TO 27 C C 72 OT OG )9.99.QE.)I(F(FI KF=F(I)*X-Y IF(KF.LT.1) KF=1 IF(KF.GT.NCOL) KF=NCOL IPLOT(KF)=ISTAR 27 CONTINUE WRITE(6,20) I,F(I),(IPLOT(J),J=1,NCOL) C C RE-INITIALIZE THE PLOTTING ARRAY. C IPLOT(KF)=IBLANK IPLOT(1)=IPLUS IPLOT(NCOL)=IPLUS IPLOT(KBAR)=IDASH 30 CONTINUE WRITE(6,40) (IFRAME(J),J=1,NCOL) 20 FORMAT(5X,I4,2X,E11.4,3X,30A1) 40 FORMAT(25X,30A1) RETURN END C C SUBROUTINE FFS (FINITE FOURIER SERIES) C ********** *** C C DOUBLE PRECISION ROUTINE TO COMPUTE FOURIER COEFFICIENTS USING RECURSIVE C FORMULATIONS OF THE TRIG FUNCTIONS. (SEE "NUMERICAL METHODS FOR C SCIENTISTS AND ENGINEERS", RICHARD W. HAMMING, MCGRAW-HILL, 1962, P.71) C C INPUT C ----- C C F(NN) REAL*8 ARRAY OF TIME-SERIES FUNCTION VALUES TO BE FOURIER C TRANSFORMED. C C NN INT*4 DIMENSION OF ARRAY F C C OUTPUT C ------ C C A(NN) REAL*8 ARRAY OF FOURIER COEFFICIENTS IN ORDER: C A(0),A(1),B(1),A(2),B(2), ... A(N-1),B(N-1),A(N) C C DIMENSION OF A MUST BE AT LEAST NN. C C******************************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(1),A(1),SINK(300),COSK(300) DATA PI/3.141592653589793D0/,LN/-1/ N=NN/2 NN=N*2 A(1)=0.D0 NNA=NN-2 DO 401 I=1,NN A(1)=A(1)+F(I) 401 CONTINUE A(1)=A(1)/N NEW=0 IF (LN.EQ.NN) GO TO 101 LN=NN PIN=PI/N NEW=1 INIT=1 101 DO 402 K=1,N I=K*2 IF (NEW.EQ.1) CALL MAF (PIN,SINK(K),COSK(K),INIT) TCOSK=2.D0*COSK(K) UB=0 UA=F(NN) DO 403 M=1,NNA U=TCOSK*UA-UB+F(NN-M) UB=UA UA=U 403 CONTINUE A(I)=(COSK(K)*UA-UB+F(1))/N IF (K.EQ.N) GO TO 102 A(I+1)=SINK(K)*UA/N 402 CONTINUE 102 A(NN)=A(NN)/2.D0 RETURN END SUBROUTINE FINT (C,NN,DT,F) C C******************************************************************************* C C SUBROUTINE FINT (FOURIER INTERPOLATION) C ********** **** C C DOUBLE PRECISION ROUTINE TO INTERPOLATE A REQUIRED VALUE OF A FUNCTION C GIVEN ITS FOURIER TRANSFORM COEFFICIENTS C C******************************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION C(1),F(1) DATA PI/3.141592653589793D0/ N=NN/2 OM=PI/N DO 401 I=1,NN IF (I.EQ.1.AND.DT.LT.0.D0) GO TO 401 IF (I.EQ.NN.AND.DT.GT.0.D0) GO TO 401 X=(I-1+DT)*OM F(I)=C(1)/2.D0 INIT=1 DO 402 J=1,N K=J*2 CALL MAF (X,SINJX,COSJX,INIT) F(I)=F(I)+C(K)*COSJX IF (J.LT.N) F(I)=F(I)+C(K+1)*SINJX 402 CONTINUE 401 CONTINUE RETURN END SUBROUTINE MAF (T,SINMT,COSMT,INIT) C C******************************************************************************* C C SUBROUTINE MAF (MULTIPLE ANGLE FORMULA) C ********** *** C C DOUBLE PRECISION ROUTINE TO COMPUTE SIN AND COS OF MULTIPLE ANGLES, C FOR USE IN FOURIER TRANSFORM. C BASED ON RECURSIVE FORMULAE: C C SIN(MT)=V(M)*SIN(T) C COS(MT)=V(M)*COS(T)-V(M-1) C C WHERE V(0)=0 C V(1)=1 C V(M)=2*COS(T)*V(M-1)-V(M-2) (M=2,3, ... ) C C INPUT C ----- C C T REAL*8 T IS THE BASIC ANGLE IN RADIANS. C C INIT INT*4 IF INIT=1 A NEW SEQUENCE OF MULTIPLE ANGLES IS C INITIALIZED AND MAF RESETS INIT=0. C IF INIT=0 THE MULTIPLE ANGLE SEQUENCE IS CONTINUED C USING THE SAME BASIC ANGLE T. MAF C AUTOMATICALLY INCREMENTS M ON EACH CALL. C C OUTPUT C ------ C C SINMT ) REAL*8 RETURN THE VALUES OF SIN(MT) AND COS(MT). ON THE C COSMT ) FIRST CALL, M=1. ON SUBSEQUENT CALLS M=M+1 AS C LONG AS INIT=0. C C INIT INT*4 MAF AUTOMATICALLY RESETS INIT=0 WHEN A NEW C SEQUENCE IS INITIALIZED. C C C******************************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DATA PI2/6.283185307179586D0/ IF (INIT.NE.1) GO TO 501 INIT=0 TR=DMOD(T,PI2) SINT=DSIN(TR) COST=DCOS(TR) VM1=1.D0 VM2=0.D0 SINMT=SINT COSMT=COST RETURN 501 VM=2.D0*COST*VM1-VM2 SINMT=VM*SINT COSMT=VM*COST-VM1 VM2=VM1 VM1=VM RETURN END FUNCTION NDAYS (IYEAR,MON) DIMENSION ND(12) DATA ND/31,28,31,30,31,30,31,31,30,31,30,31/ NDAYS=ND(MON) IF (MON.EQ.2.AND.MOD(IYEAR,4).EQ.0) NDAYS=29 RETURN END