C/*/WINDLOAD 00 C C C THIS PROGRAM SUPERSEDES THE PROGRAM LISTED IN 'THE BUFFETING OF TALL C STRUCTURES BY STRONG WINDS' BY E. SIMIU AND D. W. LOZIER (NBS BUILDING C SCIENCE SERIES 74, OCT. 1975). C THE PROGRAM CALCULATES THE ALONG-WIND RESPONSE OF A TALL BUILDING C SUBJECTED TO A WIND LOAD WITH MEAN DIRECTION NORMAL TO A BUILDING FACE. C THE ACROSS-WIND AND THE TORSIONAL RESPONSE ARE NOT CALCULATED BY THIS C PROGRAM C THE PROGRAM IS WRITTEN IN FORTRAN AND CONSISTS OF A MAIN PROGRAM CALLED C MAIN, AND EIGHT SUBPROGRAMS, CALLED INPUT, INIT, TRIPLE, F, UTILDA, XMU, C STILDA, AND FISTAR. C C C THE CARD DECK READ BY INPUT HAS THE FOLLOWING FORM C CARD A, FORMAT (2I5) C RLIM THE NUMBER OF MODES OF VIBRATION USED IN THE CALCULATIONS, C BETWEEN 1 AND 8 C IPRINT SELECTS ONE OF TWO OPTIONS FOR OUTPUT PRINTING BY THE MAIN C PROGRAM AS THE APPROXIMATE INTEGRATION PROGRESSES FROM C FTILDA = 0 TO FTILDA = SOME CUTOFF VALUE. OPTION 1, C OBTAINED BY SETTING IPRINT=1, CAUSES RUNNING INFORMATION TO C BE PRINTED AT EACH SAMPLE VALUE OF FTILDA IN THE QUADRATURE. C OPTION 2, OBTAINED BY LEAVING IPRINT BLANK, SUPPRESSES THIS C PRINTING C C C CARD B, FORMAT (3F10.0) C H HEIGHT OF BUILDING IN METERS C BCON WIDTH OF BUILDING (I.E., HORIZONTAL DIMENSION PERPENDICULAR C TO WIND DIRECTION) IN METERS C DCON DEPTH OF BUILDING (I.E., HORIZONTAL DIMENSION PARALLEL C TO WIND DIRECTION) IN METERS C C C CARD C, FORMAT (8F10.0) C EN(I) NATURAL FREQUENCIES OF THE BUILDING IN MODES 1,2,...IN HERTZ C C C CARD D, FORMAT (8F10.0) C ZETA(I) DAMPING RATIOS IN MODES 1,2,... C C C SETS OF CARDS E, F C C A SET OF CARDS E, F IS REQUIRED FOR EACH MODE OF VIBRATION USED IN THE C THE CALCULATIONS. THE FORMAT OF CARD E IS (8F10.0). THE FORMAT OF C CARD F IS (7F10.0). C CARD E INCLUDES THE ORDINATES XMUTAB OF THE I-TH MODAL SHAPE AT C ELEVATIONS 0, 1*H/14, 2*H/14, 3*H/14, 4*H/14, 5*H/14, 6*H/14, 7*H/14. C CARD F INCLUDES THE ORDINATES XMUTAB OF THE I-TH MODAL SHAPE AT C ELEVATIONS 8*H/14, 9*H/14, 10*H/14, 11*H/14, 12*H/14, 13*H/14, 14*H/14 C C C IF TWO OR MORE VIBRATION MODES ARE TAKEN INTO ACCOUNT IN THE C CALCULATIONS, THE SET OF CARDS E,F FOR THE FIRST MODE MUST BE FOLLOWED C BY THE SET OF CARDS E, F FOR THE SECOND MODE, AND SO FORTH. C C SET OF CARDS G, H C CARD G HAS FORMAT (8F10.0) AND INCLUDES THE WEIGHT OF THE BUILDING PER C UNIT HEIGHT, XMASS, IN NEWTONS/METER AT C ELEVATIONS 0, 1*H/14, 2*H/14, 3*H/14, 4*H/14, 5*H/14, 6*H/14, 7*H/14. C CARD H HAS FORMAT (7F10.0) AND INCLUDES THE WEIGHT OF THE BUILDING PER C UNIT HEIGHT, XMASS, IN NEWTONS/METER AT C ELEVATIONS 8*H/14, 9*H/14, 10*H/14, 11*H/14, 12*H/14, 13*H/14, 14*H/14 C C C C CARD I FORMAT (I5, 5X, 7F10.0) C ICODE SELECTS PARAMETERS CORRESPONDING TO STANDARD C MICROMETEOROLOGICAL CONDITIONS FOR VARIOUS TYPES OF EXPOSURE, C AS FOLLOWS C C C OPTION ICODE = 1 CAUSES SELECTION OF PARAMETERS CORRESPONDING C TO EXPOSURE TYPICAL OF OPEN WATER C C OPTION ICODE = 2 CAUSES SELECTION OF PARAMETERS CORRESPONDING C TO EXPOSURE TYPICAL OF OPEN TERRAIN C C OPTION ICODE = 3 CAUSES SELECTION OF PARAMETERS CORRESPONDING C TO EXPOSURE TYPICAL OF SUBURBS AT CONSIDERABLE DISTANCES C FROM TOWNS, SITES IN RURAL ZONES SPARSELY BUILT-UP, WITH C TREES, HEDGES, ETC C C OPTION ICODE = 4 CAUSES SELECTION OF PARAMETERS CORRESPONDING C TO EXPOSURE TYPICAL OF TOWNS, DENSELY BUILT-UP SUBURBS, C WOODED TERRAIN C C OPTION ICODE = 5 CAUSES SELECTION OF PARAMETERS CORRESPONDING C TO EXPOSURE TYPICAL OF CENTERS OF LARGE CITIES C C Z0 ROUGHNESS LENGTH IN METERS. IF LEFT BLANK, STANDARD VALUE C FOR EXPOSURE CORRESPONDING TO CHOSEN VALUE OF ICODE WILL BE C AUTOMATICALLY SELECTED C ZPSP ZERO PLANE DISPLACEMENT IN METERS. IF LEFT BLANK, STANDARD C ZPSP = 0 WILL BE AUTOMATICALLY SELECTED C CZ EXPONENTIAL DECAY COEFFICIENT FOR VERTICAL SEPARATION IN C EXPRESSION FOR CROSS-SPECTRA OF LONGITUDINAL VELOCITY C FLUCTUATIONS. IF LEFT BLANK, STANDARD VALUE CZ = 10 WILL BE C AUTOMATICALLY SELECTED C CY EXPONENTIAL DECAY COEFFICIENT FOR HORIZ. SEPARATION IN C EXPRESSION FOR CROSS-SPECTRA OF LONGITUDINAL VELOCITY C FLUCTUATIONS. IF LEFT BLANK, STANDARD VALUE CY = 16 WILL BE C AUTOMATICALLY SELECTED C BETACN RATIO (RMS OF LONGITUDINAL TURBULENT FLUCTUATIONS)/(USTAR**2) C IF LEFT BLANK, STANDARD VALUE FOR EXPOSURE CORRESPONDING TO C CHOSEN VALUE OF ICODE WILL BE AUTOMATICALLY SELECTED C F1 PEAK SIMILARITY COORDINATE (SEE REF. 3, EQ. 9). IF LEFT BLANK C STANDARD VALUE F1 = 0.03 WILL BE AUTOMATICALLY SELECTED C FS VALUE OF SIMILARITY COORDINATE BEYOND WHICH SIMILARITY C REPRESENTATION OF LONGITUDINAL TURBULENCE SPECTRA HOLDS. C (SEE REF. 3, EQ. 10). IF LEFT BLANK, STANDARD VALUE FS = 0.2 C WILL BE AUTOMATICALLY SELECTED C C REFERENCES PERTAINING TO MICROMETEOROLOGICAL PARAMETERS USED C IN THIS PROGRAM C REF.1 MEAN WIND PROFILES AND CHANGE OF TERRAIN ROUGHNESS, BY C J. BIETRY, C. SACRE, AND E. SIMIU, JOURNAL OF THE STRUCTURAL C DIVISION, ASCE, OCT. 1978, PP.1585-1593 C C REF. 2 ON THE RELIABILITY OF GUST LOADING FACTORS, BY C B. J. VICKERY, PROCEEDINGS OF TECHNICAL MEETING CONCERNING C WIND LOADS ON BUILDINGS AND STRUCTURES HELD AT THE C NATIONAL BUREAU OF STANDARDS, 1969, BUILDING SCIENCE SERIES C 30, NATIONAL BUREAU OF STANDARDS, NOV. 1970 C C REF. 3 WIND SPECTRA AND DYNAMIC ALONG-WIND RESPONSE, BY C E. SIMIU, JOURNAL OF THE STRUCTURAL DIVISION, ASCE, SEPT. C 1974, PP. 1897-1910 C CARD J FORMAT (I5,5X,3F10.0) C C C JCODE IF THE SPECIFIED WIND SPEED AT 10 METERS ABOVE GROUND IN OPEN C TERRAIN, U10, IS GIVEN AS A MEAN HOURLY SPEED IN METERS/SEC., C THEN SET JCODE = 1. IF U10 IS GIVEN AS A FASTEST-MILE SPEED C IN MILES/HOUR, THEN SET JCODE = 2 C U10 SPECIFIED WIND SPEED AT 10 M. ABOVE GROUND IN OPEN TERRAIN C T DURATION OF STORMS IN SECONDS. IF LEFT BLANK, STANDARD VALUE C T = 3600 SEC. WILL BE AUTOMATICALLY SELECTED C P RETARDATION FACTOR (SEE REF. 1 QUOTED ABOVE). IF LEFT BLANK, C STANDARD VALUE WILL BE AUTOMATICALLY SELECTED. C C C C CARD K FORMAT (3F10.0) C C CW MEAN PRESSURE COEFFICIENT ON WINDWARD SIDE. IF LEFT BLANK, C VALUE CW = 0.8 WILL BE AUTOMATICALLY SELECTED C CL ABSOLUTE VALUE OF MEAN SUCTION COEFFICIENT ON LEEWARD SIDE. C IF LEFT BLANK, VALUE CL = 0.5 WILL BE AUTOMATICALLY SELECTED C RHO SPECIFIC WEIGHT OF AIR. IF LEFT BLANK, VALUE RHO = 12.258 C NEWTONS/METER**3 WILL BE AUTOMATICALLY SELECTED. C C C C C IF SEVERAL SETS OF DATA ARE USED, CARD A OF THE SECOND SET FOLLOWS C IMMEDIATELY CARD K OF THE FIRST SET, CARD A OF THE THIRD SET C FOLLOWS IMMEDIATELY CARD K OF THE SECOND SET, AND SO FORTH. C A BLANK CARD MUST FOLLOW CARD K OF THE LAST SET OF DATA. C C C C C C C C C MAN01650 C MAN01660 C MAN01670 INTEGER R,RVAL,SVAL,RLIM MAN01680 DIMENSION PSUM(4),WP(4),FACTOR(4),WEIGHT(4),TERM(4),TERMV(4), MAN01690 1SIGMA(4),PSUU(4),WEIGHU(4),TERU(4),XJH(4) MAN01700 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), MAN01710 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, MAN01720 2 IPRINT,RLIM MAN01730 COMMON/BLOCK1/KOUNT MAN01740 COMMON/BLOCK3/FTILDA(200),ATILDA(200),M,NF(4),NN(4) MAN01750 COMMON/BLOCK5/RVAL,SVAL MAN01760 COMMON/BLOCK6/XMUINT(8),GINT(8) MAN01770 C C ALAN HECKERT 12/1/2005: ADD A "SAVE" TO ALL PROGRAM UNITS. C THIS SAVES THE VALUE OF ALL VARIABLES BETWEEN SUBROUTINE C CALLS (THE FORTRAN-66 DEFAULT, BUT NOT THE FORTRAN-77 OR C FORTRAN-90 DEFAULT). C SAVE C 200 CONTINUE MAN01780 KOUNT=0 MAN01790 CALL INPUT MAN01800 TFAC1=CW*CW+2.0*CW*CL+CL*CL MAN01810 TFAC2=CW*CW+XN*2.0*CW*CL+CL*CL MAN01820 TCRIT=0.9*EN(1)*H/USTAR MAN01830 DO 202 J=1,4 MAN01840 XJH(J)=0.0 MAN01850 202 CONTINUE MAN01860 GH=0.0 MAN01870 DO 250 R=1,RLIM MAN01880 RVAL=R MAN01890 SVAL=RVAL MAN01900 CALL INIT MAN01910 C PERFORM FOUR-DIMENSIONAL INTEGRATION. MAN01920 DO 205 J=1,4 MAN01930 PSUU(J)=0.0 MAN01940 SIGMA(J)=0.0 MAN01950 205 PSUM(J)=0.0 MAN01960 N=NN(1) MAN01970 C M IS NUMBER OF STEPS IN FTILDA INTEGRATION MAN01980 DO 230 I=1,M MAN01990 IF (I.GT.NF(1)) N=NN(2) MAN02000 IF (I.GT.NF(1)+NF(2)) N=NN(3) MAN02010 IF (I.GT.NF(1)+NF(2)+NF(3)) N=NN(4) MAN02020 W=FTILDA(I) MAN02030 CALL TRIPLE(N,W,Z,E) MAN02040 IF (IPRINT.EQ.2) GO TO 212 MAN02050 IF (MOD(I-1,5).EQ.0) WRITE (6,210) RVAL,SVAL MAN02060 WRITE (6,211) Z,E MAN02070 210 FORMAT(19H1INTEGRAL FOR R,S =,I2,1H,,I2) MAN02080 211 FORMAT(//26H0AVG OF TRIPLE INTEGRALS =,E12.5,4X,9HERR EST =,E12.5)MAN02090 GO TO 214 MAN02100 212 IF (I.LT.M) GO TO 214 MAN02110 WRITE (6,210) RVAL,SVAL MAN02120 WRITE (6,211) Z,E MAN02130 214 CONTINUE MAN02140 WP(1)=1.0 MAN02150 WP(2)=W*W MAN02160 WP(3)=WP(2)*WP(2) MAN02170 WP(4)=WP(3)*WP(2) MAN02180 PHISTR=FISTAR(W) MAN02190 TFAC=TFAC1 MAN02200 IF (W.GE.TCRIT) TFAC=TFAC2 MAN02210 DO 215 J=1,4 MAN02220 FACTOR(J)=WP(J)*PHISTR MAN02230 WEIGHT(J)=ATILDA(I)*FACTOR(J) MAN02240 TERM(J)=Z*WEIGHT(J) MAN02250 PSUM(J)=PSUM(J)+TFAC*TERM(J) MAN02260 C PSUM IS I SUB RRL MAN02270 TERMV(J)=(WEIGHT(J)*E)*(WEIGHT(J)*E) MAN02280 SIGMA(J)=SQRT(SIGMA(J)*SIGMA(J)+TERMV(J)) MAN02290 WEIGHU(J)=ATILDA(I)*WP(J) MAN02300 TERU(J)=Z*WEIGHU(J) MAN02310 PSUU(J)=PSUU(J)+TERU(J) MAN02320 C PSUU IS I SUB RRL WITH FISTAR AND TFAC2 EQUAL UNITY. 00000670 MAN02330 215 CONTINUE MAN02340 IF (IPRINT.EQ.2 .AND. I.LT.M) GO TO 230 MAN02350 WRITE (6,220) MAN02360 220 FORMAT(//19X,12HMECH. ADMITT,9X,6HWEIGHT,11X,4HTERM,4X,11HPARTIAL MAN02370 1SUM,6X,9HTERM VAR.,10X,5HSIGMA,9X,4HTERU,1X,12HPARTIAL SUMU ) MAN02380 WP(1)=FTILDA(I) MAN02390 DO 225 J=1,4 MAN02400 WRITE (6,221) WP(J),FACTOR(J),WEIGHT(J),TERM(J),PSUM(J),TERMV(J), MAN02410 1SIGMA(J),TERU(J),PSUU(J) MAN02420 221 FORMAT(1X,7E15.8,2E13.8) MAN02430 225 CONTINUE MAN02440 230 CONTINUE MAN02450 DO 240 J=1,4 MAN02460 XJH(J)=XJH(J)*XJH(J) MAN02470 XJH(J)=XJH(J)+XMU(RVAL,1.0)*XMU(SVAL,1.0)*PSUM(J)/ MAN02480 1(XMUINT(RVAL)*XMUINT(SVAL)*(EN(RVAL)*EN(SVAL)*(H/USTAR)**2)**2) MAN02490 XJH(J)=SQRT(XJH(J)) MAN02500 240 CONTINUE MAN02510 GH=GH+XMU(R,1.0)*GINT(R)/(XMUINT(R)*(EN(R)*H/USTAR)**2) MAN02520 DEF=0.5*(CW+CL)*RHO*BCON*H*H*GH/(39.478418*XMASS(1)) MAN02530 RMSDEF=RHO*BCON*H*H*XJH(1)/(39.478418*XMASS(1)) MAN02540 RMSACC=RHO*BCON*USTAR*USTAR*XJH(3)/XMASS(1) MAN02550 XX=SQRT(2.0*ALOG(USTAR*XJH(2)*T/(H*XJH(1)))) MAN02560 PEAKDF=XX+0.577/XX MAN02570 XXX=SQRT(2.0*ALOG(USTAR*XJH(4)*T/(H*XJH(3)))) MAN02580 PEAKAC=XXX+0.577/XXX MAN02590 IF (R.EQ.1) WRITE (6,278) MAN02600 IF (R.GT.1) WRITE (6,279) R MAN02610 278 FORMAT(42H1RESPONSE BASED ON FIRST MODE OF VIBRATION) MAN02620 279 FORMAT(24H1RESPONSE BASED ON FIRST,I3, MAN02630 119H MODES OF VIBRATION) MAN02640 WRITE (6,280) DEF,RMSDEF,RMSACC,PEAKDF,PEAKAC MAN02650 280 FORMAT(29H0MEAN ALONGWIND DEFLECTION...,5X,E16.8,' METERS'/ MAN02660 134H0RMS OF FLUCTUATING DEFLECTIONS...,E16.8,' METERS'/ MAN02670 224H0RMS OF ACCELERATIONS...,10X,E16.8,' METERS PER SEC**2'/ MAN02680 330H0PEAK FACTOR FOR DEFLECTION...,4X,E16.8/ MAN02690 432H0PEAK FACTOR FOR ACCELERATION...,2X,E16.8) MAN02700 250 CONTINUE MAN02710 GO TO 200 MAN02720 END MAN02730 C/*/INPUT INP00010 SUBROUTINE INPUT INP00020 INTEGER RLIM INP00030 DIMENSION DFLTAB(5,3),CTABL(10,2) INP00040 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), INP00050 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, INP00060 2 IPRINT,RLIM INP00070 C SAVE C DATA DFLTAB/0.005,0.07,0.30,1.00,2.50,6.00,6.00,5.25,4.85, INP00080 1 4.00,0.83,1.00,1.15,1.33,1.46/ INP00090 DATA CTABL/2.,5.,10.,30.,60.,100.,200.,500.,1000.,3600.,0.65, INP00100 1 0.68,0.70,0.78,0.81,0.85,0.88,0.93,0.97,1.00/ INP00110 C INP00120 C INPUT PROGRAM CONTROL PARAMETERS INP00130 C INP00140 READ (5,1) RLIM,IPRINT INP00150 CCCC1 FORMAT(2I5) INP00910 IF (RLIM .EQ. 0) GO TO 30 INP00160 IF (IPRINT .EQ. 0) IPRINT = 2 INP00170 C INP00180 C INPUT STRUCTURAL PARAMETERS INP00190 C INP00200 READ (5,2) H,BCON,DCON INP00210 CCCC2 FORMAT(8F10.0) INP00920 READ (5,2) (EN(I),I=1,RLIM) INP00220 READ (5,2) (ZETA(I),I=1,RLIM) INP00230 DO 5 I=1,RLIM INP00240 5 READ (5,2) (XMUTAB(J,I),J=1,15) INP00250 READ (5,2) (XMASS(I),I=1,15) INP00260 C INP00270 C INPUT MICROMETEROLOGICAL PARAMETERS INP00280 C INP00290 READ (5,3) ICODE,Z0,ZPDSP,CZ,CY,BETACN,F1,FS INP00300 CCCC3 FORMAT(I5,5X,7F10.0) INP00930 IF (Z0 .EQ. 0.) Z0 = DFLTAB(ICODE,1) INP00310 IF (CZ .EQ. 0.) CZ = 10. INP00320 IF (CY .EQ.0.) CY = 16. INP00330 IF (BETACN .EQ. 0.) BETACN = DFLTAB(ICODE,2) INP00340 IF (F1 .EQ. 0.) F1 = 0.03 INP00350 IF (FS .EQ. 0.) FS = 0.2 INP00360 C INP00370 C INPUT CLIMATOLOGICAL PARAMETERS INP00380 C INP00390 READ (5,3) JCODE,U10,T,P INP00400 IF (T .EQ. 0.) T = 3600. INP00410 IF (P .EQ. 0.) P = DFLTAB(ICODE,3) INP00420 C INP00430 C INPUT AERODYNAMIC PARAMETERS INP00440 C INP00450 READ (5,2) CW,CL,RHO IF (CW .EQ.0.) CW = 0.8 INP00470 IF (CL .EQ. 0.) CL = 0.5 INP00480 IF (RHO .EQ. 0.) RHO = 12.2583 INP00490 C INP00500 C COMPUTE FRICTION VELOCITY ... USTAR INP00510 C INP00520 IF (JCODE .EQ. 1) GO TO 10 INP00530 TI = 3600./U10 INP00540 DO 15 I=1,10 INP00550 IF (TI .GT. CTABL(I,1)) GO TO 15 INP00560 CT = CTABL(I-1,2) + (( CTABL(I,2)-CTABL(I-1,2) ) INP00570 * / ( CTABL(I,1)-CTABL(I-1,1) )) * ( TI - CTABL(I-1,1) ) INP00580 U0 = 0.447 * CT * U10 INP00590 GO TO 20 INP00600 15 CONTINUE INP00610 10 U0 = U10 INP00620 20 USTAR = 0.0806 * P * U0 INP00630 C INP00640 C COMPUTE ALONGWIND CROSS-CORRELATION COEF ... XN INP00650 C INP00660 X = 6.15*EN(1)*(DCON/USTAR)/ALOG(0.67*H/Z0) INP00670 XN = 1./X - (1.-EXP(-2.*X))/(2.*X**2) INP00680 C INP00690 C OUTPUT ALL INPUT PARAMETERS INP00700 C INP00710 WRITE (6,4) RLIM,IPRINT INP00720 WRITE (6,6) H,BCON,DCON INP00730 WRITE (6,7) (EN(I),I=1,RLIM) INP00740 WRITE (6,8) (ZETA(I),I=1,RLIM) INP00750 WRITE (6,12) INP00760 DO 45 I=1,RLIM INP00770 45 WRITE (6,14) I,(XMUTAB(J,I),J=1,15) INP00780 WRITE (6,16) INP00790 WRITE (6,19) (XMASS(I),I=1,15) INP00800 WRITE (6,9) Z0,ZPDSP,CZ,CY,BETACN,F1,FS INP00810 IF (JCODE .EQ. 2) GO TO 35 INP00820 WRITE (6,17) U10,T GO TO 40 INP00840 35 WRITE (6,11) U10,T 40 WRITE (6,18) CW,CL,RHO RETURN INP00870 30 WRITE (6,13) INP00880 STOP INP00890 13 FORMAT(///47H0NORMAL EXIT, END-OF-FILE ENCOUNTERED ON INPUT.) INP00900 1 FORMAT(2I5) INP00910 2 FORMAT(8F10.0) INP00920 3 FORMAT(I5,5X,7F10.0) INP00930 4 FORMAT(1H1,35X,26HVALUES OF INPUT PARAMETERS// INP00940 119H0CONTROL PARAMETERS/ INP00950 114H0 RLIM =,I3/ INP00960 214H IPRINT =,I3) INP00970 6 FORMAT(/22H0STRUCTURAL PARAMETERS/ INP00980 114H0 H =,F12.3,' METERS'/ INP00990 214H BCON =,F12.3,' METERS'/ INP01000 314H DCON =,F12.3,' METERS') INP01010 7 FORMAT(33H EN(I),I=1,RLIM ... (IN HERTZ),4F10.3,33X,4F10.3) INP01020 8 FORMAT(33H ZETA(I),I=1,RLIM ... ,4F10.3,33X,4F10.3) INP01030 12 FORMAT(33H XMUTAB(J,I),J=1,15,I=1,RLIM ... ) INP01040 14 FORMAT(7H MODE,I3,7F10.5/10X,7F10.5/10X,F10.5) INP01050 16 FORMAT(40H XMASS(I),I=1,15 ... (IN NEWTONS/METER)) INP01060 19 FORMAT(8F10.1) INP01070 9 FORMAT(/31H0MICROMETEOROLOGICAL PARAMETERS/ INP01080 114H0 Z0 =,F12.3,' METERS'/ INP01090 114H ZPDSP =,F12.3,' METERS'/ INP01100 214H CZ =,F12.3/ INP01110 214H CY =,F12.3/ INP01120 314H BETACN =,F12.3/ INP01130 414H F1 =,F12.3/ INP01140 514H FS =,F12.3) INP01150 17 FORMAT(/26H0CLIMATOLOGICAL PARAMETERS/ INP01160 114H0 U10 =,F12.3,' MEAN HOURLY SPEED (METERS/SECOND)'/ INP01170 214H T =,F12.3,' SECONDS') 11 FORMAT(/26H0CLIMATOLOGICAL PARAMETERS/ INP01200 114H0 U10 =,F12.3,' FASTEST-MILE WIND SPEED (MILES/HOUR)'/ INP01210 214H T =,F12.3,' SECONDS') 18 FORMAT(/23H0AERODYNAMIC PARAMETERS/ INP01240 114H0 CW =,F12.3/ INP01250 114H CL =,F12.3/ INP01260 214H RHO =,F12.3,' NEWTONS/M**3') END INP01290 C/*/INIT INT00010 SUBROUTINE INIT INT00020 C THIS SUBROUTINE SETS UP THE INTEGRATION PARAMETERS FOR THE CURRENT INT00030 C VALUE OF RVAL. NOTE THIS INTEGRATION SCHEME IS VALID ONLY FOR THE INT00040 C CASE RVAL=SVAL. INT00050 INTEGER RLIM,RVAL,SVAL INT00060 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), INT00070 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, INT00080 2 IPRINT,RLIM INT00090 COMMON/BLOCK3/FTILDA(200),ATILDA(200),M,NF(4),NN(4) INT00100 COMMON/BLOCK5/RVAL,SVAL INT00110 C SAVE C NF(1)=10 INT00120 NF(2)=40 INT00130 NF(3)=20 INT00140 NF(4)=10 INT00150 NN(1)=5 INT00160 NN(2)=6 INT00170 NN(3)=6 INT00180 NN(4)=6 INT00190 FPEAK=EN(RVAL)*H/USTAR INT00200 M=1 INT00210 STEPOL=0. INT00220 DO 110 I=1,4 INT00230 IF (I.EQ.1) FTILDA(M)=0. INT00240 IF (I.EQ.2) FTILDA(M) = (FPEAK-2.)/30. INT00250 IF (I.EQ.3) FTILDA(M)=FPEAK-2. INT00260 IF (I.EQ.4) FTILDA(M)=FPEAK+2. INT00270 IF (I.EQ.1) STEPNU=((FPEAK-2.)/30.)/FLOAT(NF(1)) INT00280 IF (I.EQ.2) STEPNU=((29.*FPEAK-58.)/30.)/FLOAT(NF(2)) INT00290 IF (I.EQ.3) STEPNU=4./FLOAT(NF(3)) INT00300 IF (I.EQ.4) STEPNU=(2.*FPEAK )/FLOAT(NF(4)) INT00310 ATILDA(M)=(STEPOL+STEPNU)/2. INT00320 MM=NF(I)-1 INT00330 DO 101 J=1,MM INT00340 MJ=M+J INT00350 FTILDA(MJ)=FTILDA(M)+FLOAT(J)*STEPNU INT00360 101 ATILDA(MJ)=STEPNU INT00370 M=M+MM+1 INT00380 110 STEPOL=STEPNU INT00390 M=M-1 INT00400 ATILDA(M)=ATILDA(M)+(2.*FPEAK-STEPOL)*3./20. INT00410 WRITE (6,115) RVAL,SVAL,M,NF,NN,FPEAK INT00420 115 FORMAT(60H1VALUES OF INTEGRATION PARAMETERS FOR VIBRATION MODES R,INT00430 1S = ,I2,1H,,I2/ INT00440 19H0 M =,I4/ INT00450 217H0(NF(I),I=1,4)...,4I4/ INT00460 317H0(NN(I),I=1,4)...,4I4/ INT00470 49H0 FPEAK =,F12.3) INT00480 WRITE (6,116) (FTILDA(I),I=1,M) INT00490 116 FORMAT(21H0(FTILDA(I),I=1,M)...,10F8.3/(21X,10F8.3)) INT00500 WRITE (6,117) (ATILDA(I),I=1,M) INT00510 117 FORMAT(21H0(ATILDA(I),I=1,M)...,10F8.3/(21X,10F8.3)) INT00520 RETURN INT00530 END INT00540 C/*/TRIPLE TRP00010 SUBROUTINE TRIPLE(N,W,Z,E) TRP00020 C TRIPLE INTEGRAL EVALUATED BY A MONTE CARLO METHOD. TRP00030 C INTEGRAL EVALUATED IS I SUB R,S OF FTILDA. TRP00040 C INPUT...W, THE CURRENT VALUE OF FTILDA TRP00050 C N, THE NUMBER OF SUBDIVISIONS TO TAKE ON EACH EDGE OF TRP00060 C THE UNIT CUBE FOR THE MONTE CARLO QUADRATURE TRP00070 C R AND S ARE IN LABELLED COMMON AND ARE NEEDED ONLY IN TRP00080 C THE FUNCTION SUBROUTINE F THAT EVALUATES THE INTEGRAND TRP00090 C OF THE TRIPLE INTEGRAL TRP00100 C OUTPUT...Z, THE APPROXIMATION TO THE TRIPLE INTEGRAL TRP00110 C E, AN ERROR ESTIMATE FOR Z TRP00120 DOUBLE PRECISION D,DFRA,DK,DKSQ,DPRIME TRP00130 COMMON/BLOCK1/KOUNT TRP00140 DIMENSION D(3),DPRIME(3),Q(3),X(3),XPRIME(3),P(3),R(3),RPRIME(3), TRP00150 1S(3),SPRIME(3) TRP00160 C SAVE C C TRP00170 C SIX SEQUENCES OF PSEUDO-RANDOM NUMBERS ARE GENERATED TRP00180 C FOR THE MONTE CARLO QUADRATURE. THEY ARE OF THE FORM TRP00190 C TRP00200 C X SUB K = FRACT. PART OF (K**2)*SQRT(P) TRP00210 C TRP00220 C WHERE P TAKES THE PRIME VALUES 2,3,5,7,11 AND 13. TRP00230 C VALUES OF K AND K**2 ARE STORED IN DOUBLE PRECISION TRP00240 C TO ALLOW THE LARGEST POSSIBLE VALUE OF K**2. THE TRP00250 C FOLLOWING ARITHMETIC STATEMENT FUNCTION CONTAINS TRP00260 C A NON-STANDARD INTRINSIC FUNCTION, DMOD, TO PERMIT TRP00270 C THE DOUBLE PRECISION FRACTIONAL PART TO BE TAKEN. TRP00280 C C NOTE: ALAN HECKERT 10/2004: CONVERT STATEMENT FUNCTION TO C EXTERNAL FUNCTION TO AVOID COMPILATION PROBLEMS ON C FORTRAN 77/90 COMPILERS. CCCCC DFRA(D)=DMOD(D,1.0D0) TRP00290 EXTERNAL DFRA C IF (KOUNT.GT.0) GO TO 10 TRP00300 DK=0.0D0 TRP00310 D(1)=DFRA(DSQRT(2.0D0)) TRP00320 D(2)=DFRA(DSQRT(3.0D0)) TRP00330 D(3)=DFRA(DSQRT(5.0D0)) TRP00340 DPRIME(1)=DFRA(DSQRT(7.0D0)) TRP00350 DPRIME(2)=DFRA(DSQRT(11.0D0)) TRP00360 DPRIME(3)=DFRA(DSQRT(13.0D0)) TRP00370 10 CONTINUE TRP00380 XN=N TRP00390 RXN=1.0/XN TRP00400 RXNH=0.5*RXN TRP00410 N3=N*N*N TRP00420 XN3=N3 TRP00430 RXN3=1.0/XN3 TRP00440 RXN3H=0.5*RXN3 TRP00450 Y=0.0 TRP00460 YPRIME=0.0 TRP00470 E=0.0 TRP00480 DO 30 K3=1,N TRP00490 DO 30 K2=1,N TRP00500 DO 30 K1=1,N TRP00510 DK=DK+1.0D0 TRP00520 DKSQ=DK*DK TRP00530 Q(1)=K1-1 TRP00540 Q(2)=K2-1 TRP00550 Q(3)=K3-1 TRP00560 DO 15 I=1,3 TRP00570 X(I)=DFRA(DKSQ*D(I)) TRP00580 XPRIME(I)=DFRA(DKSQ*DPRIME(I)) TRP00590 15 CONTINUE TRP00600 DO 20 I=1,3 TRP00610 Q(I)=Q(I)*RXN TRP00620 P(I)=Q(I)+RXNH TRP00630 R(I)=Q(I)+RXN*X(I) TRP00640 RPRIME(I)=Q(I)+RXN*XPRIME(I) TRP00650 PI2=2.0*P(I) TRP00660 S(I)=PI2-R(I) TRP00670 SPRIME(I)=PI2-RPRIME(I) TRP00680 20 CONTINUE TRP00690 FR=F(R,W) TRP00700 FS=F(S,W) TRP00710 G=FR+FS TRP00720 FRPRIM=F(RPRIME,W) TRP00730 FSPRIM=F(SPRIME,W) TRP00740 GPRIME=FRPRIM+FSPRIM TRP00750 Y=Y+G TRP00760 YPRIME=YPRIME+GPRIME TRP00770 E=E+(G-GPRIME)**2 TRP00780 30 CONTINUE TRP00790 Y=Y*RXN3H TRP00800 YPRIME=YPRIME*RXN3H TRP00810 Z=0.5*(Y+YPRIME) TRP00820 E=0.25*RXN3*SQRT(E) TRP00830 RETURN TRP00840 END TRP00850 C/*/F F0000010 FUNCTION F(V,W) F0000020 C CALCULATION OF INTEGRAND FOR TRIPLE INTEGRAL F0000030 C Y SUB R,S OF FTILDA F0000040 C THE PROGRAM VARIABLE W STANDS FOR FTILDA. THE VALUES F0000050 C OF R AND S ARE STORED IN COMMON BLOCK /BLOCK5/. F0000060 C THE VARIABLE KOUNT IN /BLOCK1/ IS SET TO ZERO WHEN F0000070 C NEW INPUT PARAMETERS ARE READ BY SUBROUTINE INPUT. IT IS INCREASED F0000080 C BY ONE EACH TIME THIS SUBROUTINE IS EXECUTED, SO AT THE END OF THE F0000090 C WHOLE FOUR-DIMENSIONAL INTEGRATION, IT CONTAINS THE NUMBER OF TIMES F0000100 C THE INTEGRAND FOR THE TRIPLE INTEGRAL WAS CALCULATED. F0000110 INTEGER RLIM,RVAL,SVAL F0000120 DIMENSION V(3) F0000130 COMMON/BLOCK1/KOUNT F0000140 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), F0000150 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, F0000160 2 IPRINT,RLIM F0000170 COMMON/BLOCK5/RVAL,SVAL F0000180 C SAVE C IF (KOUNT.GT.0) GO TO 2000 F0000190 C INITIALIZATION. F0000200 CON1=2.0*CZ F0000210 CON2=CY*BCON/(CZ*H) F0000220 2000 CONTINUE F0000230 Z1=V(1) F0000240 Z2=V(2) F0000250 TT=V(3) F0000260 UT1=UTILDA(Z1) F0000270 UT2=UTILDA(Z2) F0000280 ALPHA=(1.0-TT)*XMU(RVAL,Z1)*XMU(SVAL,Z2)*UT1*UT2 F0000290 ALPHA=ALPHA*SQRT(STILDA(W,Z1,UT1)*STILDA(W,Z2,UT2)) F0000300 BETA=CON1*SQRT((Z1-Z2)**2+(CON2*TT)**2)/(UT1+UT2) F0000310 F=2.0*ALPHA*EXP(-BETA*W) F0000320 KOUNT=KOUNT+1 F0000330 RETURN F0000340 END F0000350 C/*/STILDA STL00010 FUNCTION STILDA(W,Z,UT) STL00020 C THIS FUNCTION IS RELATED TO THE ASSUMED REPRESENTATION STL00030 C OF THE SPECTRUM OF LONGITUDINAL WIND FLUCTUATIONS. STILDA STL00040 C IS THE EXPRESSION STL00050 C N*S(Z,N)/(FTILDA*USTAR**2) STL00060 C STL00070 INTEGER RLIM STL00080 COMMON/BLOCK1/KOUNT STL00090 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), STL00100 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, STL00110 2 IPRINT,RLIM STL00120 C SAVE C IF (KOUNT.GT.0) GO TO 20 STL00130 C STL00140 C INITIALIZATION FOR FORMULA DEPENDENT ON THE PARAMETER F1. STL00150 12 BETA1=0.39*FS**(-2.0/3.0) STL00160 BETA2=(2.0/3.0)*BETA1 STL00170 BETA3=BETACN-BETA1 STL00180 BETA4=FS-2.0*F1 STL00190 BETA5=0.5*FS-2.0*F1 STL00200 C2=(BETA4*BETA3-BETA2*BETA5)/((1.5+ALOG(FS/F1))*BETA4-BETA5) STL00210 IF (BETA4) 16,15,16 STL00220 15 B2=2.0*(BETA2*(1.5+ALOG(2.0))-BETA3)/(FS*FS) STL00230 GO TO 17 STL00240 16 B2=(BETA2-C2)/(FS*BETA4) STL00250 17 B1=B2-C2/(F1*F1) STL00260 DTILDA=ZPDSP/H STL00270 ZA=(10.0+ZPDSP)/H STL00280 C STL00290 C FORMULA DEPENDENT ON THE PARAMETER F1. STL00300 20 CONTINUE STL00310 ZZ=AMAX1(Z,ZA) STL00320 ZUT=(ZZ-DTILDA)/UT STL00330 A=W*ZUT STL00340 IF (A.GE.FS) GO TO 202 STL00350 IF (A.GE.F1) GO TO 201 STL00360 STILDA=B1*ZUT*(A-2.0*F1) STL00370 GO TO 205 STL00380 201 STILDA=B2*ZUT*(A-2.0*F1)+C2/W STL00390 GO TO 205 STL00400 202 STILDA=0.26/(W*A**0.666667) STL00410 205 RETURN STL00420 END STL00430 C/*/XMU XMU00010 FUNCTION XMU(I,Z) XMU00020 C COMPUTATION OF MODAL SHAPES. XMU00030 INTEGER RLIM XMU00040 COMMON/BLOCK1/KOUNT XMU00050 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), XMU00060 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,TT,CW,CL,RHO,XN, XMU00070 2 IPRINT,RLIM XMU00080 COMMON/BLOCK6/XMUINT(8),GINT(8) XMU00090 C SAVE C C THE NONDIMENSIONAL COORDINATES 1,2,3,...15, REPRESENT POINTS WITHXMU00100 C ELEVATIONS 0, 1*H/14, 2*H/14,......14*H/14, RESPECTIVELY. XMU00110 IF (KOUNT .GT. 0) GO TO 5 XMU00120 C INITIALIZATION. XMU00130 DO 3 J=1,8 XMU00140 GINT(J)=0.0 XMU00150 3 XMUINT(J)=0.0 XMU00160 C COMPUTE INTEGRALS OF XMU(I,Z)**2 AND OF XMU(I,Z)*UTILDA(Z)**2. XMU00170 C THESE ARE DENOTED BY XMUINT AND GINT, RESPECTIVELY, AND THE XMU00180 C VALUES ARE RETURNED IN COMMON BLOCK /BLOCK6/. XMU00190 DO 2 J=1,RLIM XMU00200 S=0.0 XMU00210 T=0.0 XMU00220 DO 1 K=2,14 XMU00230 S=S+XMUTAB(K,J)*UTILDA(FLOAT(K-1)/14.0)**2 XMU00240 1 T=T+(XMUTAB(K,J)**2 * (XMASS(K)/XMASS(1))) XMU00250 S=S+0.5*XMUTAB(15,J)*UTILDA(1.0)**2 XMU00260 T=T+(0.5*XMUTAB(15,J)**2 * (XMASS(15)/XMASS(1))) XMU00270 GINT(J)=S/14.0 XMU00280 2 XMUINT(J)=T/14.0 XMU00290 5 CONTINUE XMU00300 Y=14.0*Z XMU00310 J=Y XMU00320 IF (J.EQ.14) GO TO 30 XMU00330 XMU=XMUTAB(J+1,I)+(Y-FLOAT(J))*(XMUTAB(J+2,I)-XMUTAB(J+1,I)) XMU00340 RETURN XMU00350 30 XMU=XMUTAB(J+1,I) XMU00360 RETURN XMU00370 END XMU00380 C/*/FISTAR FTR00010 FUNCTION FISTAR(W) FTR00020 C CALCULATION OF THE FUNCTION PHI STAR OF FTILDA. FTR00030 C W STANDS FOR FTILDA. FTR00040 INTEGER RLIM,RVAL,SVAL FTR00050 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), FTR00060 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, FTR00070 2 IPRINT,RLIM FTR00080 COMMON/BLOCK5/RVAL,SVAL FTR00090 C SAVE C A=(W/EN(RVAL))*(USTAR/H) FTR00100 E=(W/EN(SVAL))*(USTAR/H) FTR00110 G=2.*ZETA(RVAL)*A FTR00120 D=2.*ZETA(SVAL)*E FTR00130 A=1.-A*A FTR00140 E=1.-E*E FTR00150 FISTAR=(A*E+G*D)/((A*A+G*G)*(E*E+D*D)) FTR00160 RETURN FTR00170 END FTR00180 C/*/UTILDA UTL00010 FUNCTION UTILDA(Z) UTL00020 C NONDIMENSIONALIZED MEAN WIND SPEED AT ELEVATION Z. UTL00030 INTEGER RLIM UTL00040 COMMON/BLOCK1/KOUNT UTL00050 COMMON/BLOCK2/H,BCON,DCON,XMASS(15),EN(8),ZETA(8),XMUTAB(15,8), UTL00060 1 Z0,ZPDSP,CZ,CY,BETACN,F1,FS,USTAR,T,CW,CL,RHO,XN, UTL00070 2 IPRINT,RLIM UTL00080 C SAVE C IF (KOUNT.GT.0) GO TO 1 UTL00090 C INITIALIZATION. UTL00100 ZA=(10.0+ZPDSP)/H UTL00110 UCONST=2.5*ALOG(10.0/Z0) UTL00120 DTILDA=ZPDSP/H UTL00130 HOZ0=H/Z0 UTL00140 1 CONTINUE UTL00150 UTILDA=UCONST UTL00160 IF (Z.GT.ZA) UTILDA=2.5*ALOG((Z-DTILDA)*HOZ0) UTL00170 RETURN UTL00180 END UTL00190 DOUBLE PRECISION FUNCTION DFRA(D) C C CONVERT FROM STATEMENT FUNCTION TO EXTERNAL FUNCTION. C WOULD NOT COMPILE ON FORTRAN 77/90 COMPILERS. C DOUBLE PRECISION D C DFRA=DMOD(D,1.0D0) RETURN END