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 NOTE 12/2007: MODIFIED BY ALAN HECKERT TO USE LESS RIGID C FORMATTING FOR THE INPUT. IT WILL PRINT A C PROMPT AND THEN ACCEPT THE INPUT IN FREE-FORMAT C MODE. 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 C C NOTE 12/2007 (ALAN HECKERT): C MODIFY THIS ROUTINE TO SIMPLIFY THE INPUT. C C 1) GENERATE A PROMPT FOR EACH INPUT LINE C 2) IF ALL FIELDS EXPECTED, THEN JUST READ USING FREE-FORMAT C OPTION. IF SOME FIELDS CAN BE LEFT BLANK, THEN PARSE THE C LINE MANUALLY TO DETERMINE HOW MANY FIELDS ACTUALLY C ENTERED. ALL INPUTS ARE EITHER INTEGER OR REAL, SO THIS C WILL SIMPLIFY THE CONVERSION. C C NOTE 9/2009 (ALAN HECKERT): ECHO INPUT TO "windload.echo" TO SAVE C PARAMETERS ENTERED MANUALLY C PARAMETER (MAXLIM=50) INTEGER LOWLIM(MAXLIM) INTEGER UPPLIM(MAXLIM) CHARACTER*132 IATEMP LOGICAL IBLANK LOGICAL IERROR INTEGER ICOUNT C 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 LOGICAL IECHO 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 DATA ICOUNT /0/ C DATA IECHO/.FALSE./ C C INPUT PROGRAM CONTROL PARAMETERS INP00130 C INP00140 ICOUNT=ICOUNT+1 IBLANK=.FALSE. IERROR=.FALSE. C IF(.NOT.IECHO)THEN OPEN(10,FILE='windload.echo') IECHO=.TRUE. ENDIF C IF (ICOUNT.EQ.1)THEN WRITE(6,*) 'Enter Inputs for Windload Program:' WRITE(6,*) ' ' WRITE(6,*) 'The original version of this program required a' WRITE(6,*) 'rigid input format. This updated version now allows' WRITE(6,*) 'the data to be entered in free-format. A few ', 1 'comments.' WRITE(6,*) ' ' WRITE(6,*) ' 1. Data values should be separated by at least one' WRITE(6,*) ' space, comma, or tab.' WRITE(6,*) ' ' WRITE(6,*) ' 2. Some prompts will have a "(same line)" at the' WRITE(6,*) ' end. This means that all values should be' WRITE(6,*) ' entered on the same line (although the spacing' WRITE(6,*) ' between data values is arbitrary). If the' WRITE(6,*) ' "(same line)" is not specified, multiple values' WRITE(6,*) ' can be split over multiple lines.' WRITE(6,*) ' ' WRITE(6,*) ' 3. The original version of this program allowed' WRITE(6,*) ' default values to be entered by leaving the' WRITE(6,*) ' appropriate field blank. This still works for' WRITE(6,*) ' blank fields at the END of the line. However,' WRITE(6,*) ' all fields preceding a field for which you' WRITE(6,*) ' enter a non-default value must also have a' WRITE(6,*) ' value entered. Enter a 0 to use the default.' WRITE(6,*) ' ' ENDIF C WRITE(6,*) 'Section 1: Input Program Control Parameters' WRITE(6,*) 'Enter Number of Vibration Modes (1 - 8)' WRITE(6,*) ' Print Option for each Sample (1 = Yes, ', 1 '0 = Suppress) (0)' CCCCC READ (5,1) RLIM,IPRINT INP00150 IATEMP=' ' READ (5,'(A132)') IATEMP WRITE(10,'(A132)') IATEMP NCHAR=132 CALL EXTLIM(IATEMP,NCHAR,LOWLIM,UPPLIM,MAXLIM, 1 NUMLIM,ILAST,IBLANK) IF(IBLANK) THEN WRITE(6,*) 'Normal Exit to Windload Program.' WRITE(6,*) 'First input line was blank.' STOP ENDIF IFRST=LOWLIM(1) ILAST=UPPLIM(1) NCHAR=ILAST-IFRST+1 CALL EXTINT(IATEMP(IFRST:ILAST),NCHAR,IVALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the number of vibration modes.' WRITE(6,*) 'Windload program exited.' STOP ELSE RLIM=IVALUE ENDIF IF(NUMLIM.GE.2)THEN IFRST=LOWLIM(1) ILAST=UPPLIM(1) NCHAR=ILAST-IFRST+1 CALL EXTINT(IATEMP(IFRST:ILAST),NCHAR,IVALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading print control parameter.' WRITE(6,*) 'Windload program exited.' STOP ENDIF IPRINT=IVALUE ELSE IPRINT=0 ENDIF C IF (RLIM .EQ. 0) GO TO 30 INP00160 CCCCC IF (IPRINT .EQ. 0) IPRINT = 2 INP00170 IF (IPRINT .NE. 1) IPRINT = 2 INP00170 C INP00180 C INPUT STRUCTURAL PARAMETERS INP00190 C INP00200 WRITE(6,*) 'Section 2: Input Structural Parameters' WRITE(6,*) 'Enter Height, Width, and Depth of Building (meters):' READ (5,*,ERR=9010,END=9010) H,BCON,DCON WRITE (10,'(G15.7)') H,BCON,DCON C WRITE(6,*) 'Enter Natural Frequencies of Building for each ', 1 ' vibration mode (Hertz):' READ (5,*,ERR=9010,END=9010) (EN(I),I=1,RLIM) WRITE (10,*) (EN(I),I=1,RLIM) C WRITE(6,*) 'Enter Damping Ratios for each vibration mode:' READ (5,*,ERR=9010,END=9010) (ZETA(I),I=1,RLIM) WRITE (10,'(8G15.7)') (ZETA(I),I=1,RLIM) C DO 5 I=1,RLIM INP00240 WRITE(6,*) 'For Vibration Mode ',I,', enter the ordinates of ', 1 'XMUTAB at' WRITE(6,*) 'elevations 0, 1*H/14, 2*H/14, ..., 14*H/14:' READ (5,*,ERR=9010,END=9010) (XMUTAB(J,I),J=1,15) WRITE (10,'(8G15.7)') (XMUTAB(J,I),J=1,15) 5 CONTINUE C WRITE(6,*) 'Enter the weight of the building per unit height, ', 1 'XMASS, in Newtown/meter at' WRITE(6,*) 'elevations 0, 1*H/14, 2*H/14, ..., 14*H/14:' READ (5,*,ERR=9010,END=9010) (XMASS(I),I=1,15) WRITE (10,'(8G15.7)') (XMASS(I),I=1,15) C INP00270 C INPUT MICROMETEROLOGICAL PARAMETERS INP00280 C INP00290 WRITE(6,*) 'Section 3: Input Micrometerological Parameters ', 1 '(same line)' WRITE(6,*) 'Enter Option Code (1 = Open Water, 2 = Open Terrain,', 1 '3 = Rural, 4 = Towns, 5 = Large Cities)' WRITE(6,*) ' Roughness length (meters, 0 for default)' WRITE(6,*) ' (1 => 0.005, 2 => 0.07, 3 => 0.3, ', 1 '4 => 1.0, 5=> 2.5)' WRITE(6,*) ' Zero plane dispalcement (meters, 0 for default)' WRITE(6,*) ' Exponential decay coefficient for vertical ', 1 'separation (0 for default (10.0))' WRITE(6,*) ' Exponential decay coefficient for horizontal ', 1 'separation (0 for default (16.0))' WRITE(6,*) ' RMS of Longitudinal turbulent fluctuations/', 1 'USTAR**2 (0 for default)' WRITE(6,*) ' (1 => 6.0, 2 => 6.0, 3 => 5.25, ', 1 '4 => 4.85, 5=> 4.0)' WRITE(6,*) ' Peak similarity coordinate (0 for default ', 1 '(0.03))' WRITE(6,*) ' Similarity coordinate beyond which similarity ', 1 'representation of' WRITE(6,*) ' longitudinal turbulence spectra holds ', 1 '(0 for default (0.2))' CCCCC READ (5,3) ICODE,Z0,ZPDSP,CZ,CY,BETACN,F1,FS INP00300 IATEMP=' ' READ (5,'(A132)') IATEMP WRITE (10,'(A132)') IATEMP NCHAR=132 CALL EXTLIM(IATEMP,NCHAR,LOWLIM,UPPLIM,MAXLIM, 1 NUMLIM,ILAST,IBLANK) IF(IBLANK) THEN WRITE(6,*) 'No micrometeorological paramaeters entered.' WRITE(6,*) 'Windload program exited.' STOP ENDIF IFRST=LOWLIM(1) ILAST=UPPLIM(1) NCHAR=ILAST-IFRST+1 CALL EXTINT(IATEMP(IFRST:ILAST),NCHAR,IVALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE ICODE=IVALUE ENDIF C IF(NUMLIM.GE.2)THEN IFRST=LOWLIM(2) ILAST=UPPLIM(2) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE Z0=VALUE ENDIF ELSE Z0=0. ENDIF C IF(NUMLIM.GE.3)THEN IFRST=LOWLIM(3) ILAST=UPPLIM(3) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE ZPDSP=VALUE ENDIF ELSE ZPDSP=0. ENDIF C IF(NUMLIM.GE.4)THEN IFRST=LOWLIM(4) ILAST=UPPLIM(4) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE CZ=VALUE ENDIF ELSE CZ=0. ENDIF C IF(NUMLIM.GE.5)THEN IFRST=LOWLIM(5) ILAST=UPPLIM(5) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE CY=VALUE ENDIF ELSE CY=0. ENDIF C IF(NUMLIM.GE.6)THEN IFRST=LOWLIM(6) ILAST=UPPLIM(6) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE BETACN=VALUE ENDIF ELSE BETACN=0. ENDIF C IF(NUMLIM.GE.7)THEN IFRST=LOWLIM(7) ILAST=UPPLIM(7) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE F1=VALUE ENDIF ELSE F1=0. ENDIF C IF(NUMLIM.GE.8)THEN IFRST=LOWLIM(8) ILAST=UPPLIM(8) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the micrometeorological ', 1 'parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE FS=VALUE ENDIF ELSE FS=0. ENDIF C 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 WRITE(6,*) 'Section 4: Input Climatological Parameters ', 1 '(same line)' WRITE(6,*) 'Enter Option Code (1 = speed meters/second, 2 = ', 1 'speed miles/hour)' WRITE(6,*) ' Specified wind speed at 10 meters' WRITE(6,*) ' Duration of storm (seconds, 0 for default ', 1 '(3600))' WRITE(6,*) ' Retardation factor (0 for default)' WRITE(6,*) ' (1 => 0.83, 2 => 1.0, 3 => 1.15, ', 1 '4 => 1.33, 5=> 1.46)' CCCCC READ (5,3) JCODE,U10,T,P INP00400 IATEMP=' ' READ (5,'(A132)') IATEMP WRITE (10,'(A132)') IATEMP NCHAR=132 CALL EXTLIM(IATEMP,NCHAR,LOWLIM,UPPLIM,MAXLIM, 1 NUMLIM,ILAST,IBLANK) IF(IBLANK) THEN WRITE(6,*) 'No climatological paramaeters entered.' WRITE(6,*) 'Windload program exited.' STOP ENDIF IFRST=LOWLIM(1) ILAST=UPPLIM(1) NCHAR=ILAST-IFRST+1 CALL EXTINT(IATEMP(IFRST:ILAST),NCHAR,IVALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the climatological parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE JCODE=IVALUE ENDIF C IFRST=LOWLIM(2) ILAST=UPPLIM(2) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the climatological parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE U10=VALUE ENDIF C IF(NUMLIM.GE.3)THEN IFRST=LOWLIM(3) ILAST=UPPLIM(3) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the climatological parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE T=VALUE ENDIF ELSE T=0. ENDIF C IF(NUMLIM.GE.4)THEN IFRST=LOWLIM(4) ILAST=UPPLIM(4) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the climatological parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE P=VALUE ENDIF ELSE P=0. ENDIF C 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 WRITE(6,*) 'Section 5: Input Aerodynamic Parameters (same line)' WRITE(6,*) 'Enter Mean pressure coefficient on windward side ', 1 '(0 for default (0.8))' WRITE(6,*) ' Absolute value of mean suction coefficient ', 1 'on leeward side (0 for default (0.5))' WRITE(6,*) ' Specific weight of air (0 for default) ', 1 '(12.2583)' CCCCC READ (5,2) CW,CL,RHO IATEMP=' ' READ (5,'(A132)') IATEMP WRITE (10,'(A132)') IATEMP NCHAR=132 CALL EXTLIM(IATEMP,NCHAR,LOWLIM,UPPLIM,MAXLIM, 1 NUMLIM,ILAST,IBLANK) C IF(NUMLIM.GE.1)THEN IFRST=LOWLIM(1) ILAST=UPPLIM(1) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the aerodynamic parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE CW=VALUE ENDIF ELSE CW=0. ENDIF C IF(NUMLIM.GE.2)THEN IFRST=LOWLIM(2) ILAST=UPPLIM(2) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the aerodynamic parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE CL=VALUE ENDIF ELSE CL=0. ENDIF C IF(NUMLIM.GE.3)THEN IFRST=LOWLIM(3) ILAST=UPPLIM(3) NCHAR=ILAST-IFRST+1 CALL EXTREA(IATEMP(IFRST:ILAST),NCHAR,VALUE,IERROR) IF(IERROR)THEN WRITE(6,*) 'Error reading the aerodynamic parameters.' WRITE(6,*) 'Windload program exited.' STOP ELSE RHO=VALUE ENDIF ELSE RHO=0. ENDIF C 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 C 9010 CONTINUE WRITE(6,*)'An unexpected error or end of file encountered' WRITE(6,*)'on the input. Windload program exited.' STOP C 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 C C SEPTEMBER 2009. ECHO INPUT TO FILE "windload.echo" C 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 CCCCC IF (BETA4) 16,15,16 STL00220 IF (BETA4.EQ.0.) THEN B2=2.0*(BETA2*(1.5+ALOG(2.0))-BETA3)/(FS*FS) ELSE B2=(BETA2-C2)/(FS*BETA4) ENDIF 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 subroutine extlim(iatemp,nchar,lowlim,upplim,maxlim, 1 numlim,ilast,iblank) c c purpose--This subroutine searches a record for tabs, spaces, or c commas. It returns the first (lowlim) and last character for c each field in the record. c character*(*) iatemp integer lowlim(*) integer upplim(*) logical iblank c character*1 itab c iblank=.false. c c Initialize c do 10 i=1,maxlim lowlim(i)=0 upplim(i)=0 10 continue c c Find first non-blank character c ifrst=nchar do 100 i=1,nchar iindx=ichar(iatemp(i:i)) if (iatemp(i:i).eq.' ' .or. iatemp(i:i).eq.',' .or. 1 iindx.le.32)goto100 ifrst=i goto199 100 continue iblank=.true. goto9000 199 continue c c Find last non-blank character (ilast indicates a blank line) c ilast=nchar do 200 i=nchar,1,-1 iindx=ichar(iatemp(i:i)) if (iatemp(i:i).eq.' ' .or. iatemp(i:i).eq.',' .or. 1 iindx.le.32)goto200 ilast=i goto299 200 continue iblank=.true. goto9000 299 continue c c Search record for separators: space, commas, or non-printing c character (i.e., ASCII index <= 32). Be sure that multiple c spaces count as single separator. c nlim=1 lowlim(nlim)=ifrst iflag=0 do 300 i=ifrst,ilast iindx=ichar(iatemp(i:i)) if (iatemp(i:i).eq.' ' .or. iatemp(i:i).eq.',' .or. 1 iindx.le.32)then if (iflag.eq.0)then iflag=1 upplim(nlim)=i-1 else goto300 end if else if (iflag.eq.1) then nlim=nlim+1 iflag=0 if(nlim.le.maxlim)then lowlim(nlim)=i else write(6,*)'Error: maximum number of fields (',maxlim, 1 ') exceeded' write(6,*)'Windload program exited.' stop endif endif end if 300 continue upplim(nlim)=ilast c 9000 continue c c Check for empty fields. c numlim=0 do 9010 i=1,maxlim if (lowlim(i).gt.upplim(i))then lowlim(i)=0 upplim(i)=0 elseif (lowlim(i).le.0 .or. upplim(i).le.0)then lowlim(i)=0 upplim(i)=0 else numlim=numlim+1 end if 9010 continue 9019 continue if (numlim.eq.0)then iblank=.true. endif c return end subroutine extint(iatemp,nchar,ivalue,ierror) c c purpose--This subroutine extracts an integer value from c a character string. Look for first and last c non-blank characters (if entire field is blank, c return a 0). If a non-integer value encountered, c print an error message. c character*(*) iatemp logical ierror character*10 iformt c ierror=.false. c c Find first non-blank character c istrt=0 do 100 i=1,nchar if (iatemp(i:i).eq.' ')goto100 istrt=i goto199 100 continue ivalue=0 goto9090 199 continue c c Find last non-blank character c ilast=nchar do 200 i=nchar,istrt,-1 if (iatemp(i:i).eq.' ')goto200 ilast=i goto299 200 continue ivalue=0 goto9090 299 continue c c Truncated if a decimal point found c do 2200 i=istrt,ilast if (iatemp(i:i).eq.'.')then ilast=i-1 goto2299 endif 2200 continue 2299 continue c c Convert text to integer c do 300 i=istrt,ilast if (iatemp(i:i).eq.' ')then iatemp(i:i)='0' else itemp=ichar(iatemp(i:i)) if (itemp.lt.48 .or. itemp.gt.57)then write(6,*)'Error: non-integer value encountered when' write(6,*)' extracting integer value.' write(6,*)' Character is: ',char(itemp) write(6,*)' Windload program exited.' stop end if end if 300 continue ntemp=ilast-istrt+1 iformt='(i )' if (ntemp.le.9)then write(iformt(3:3),'(i1)')ntemp else if (ntemp.le.99)then write(iformt(3:4),'(i2)')ntemp else if (ntemp.le.999)then write(iformt(3:5),'(i3)')ntemp else if (ntemp.le.9999)then write(iformt(3:6),'(i4)')ntemp else if (ntemp.le.99999)then write(iformt(3:7),'(i5)')ntemp else if (ntemp.le.999999)then write(iformt(3:8),'(i6)')ntemp end if read(iatemp(istrt:ilast),iformt,err=9010)ivalue goto9090 c 9010 continue write(6,*)'Error: Unable to convert integer field' write(6,*)'Windload program exited.' stop c 9090 continue return end subroutine extrea(iatemp,nchar,avalue,ierror) c c purpose--This subroutine extracts a real value from c a character string. Look for first and last c non-blank characters (if entire field is blank, c return a 0.). The input string should contain only c "+", "-", ".", or a number. If any other character c encountered, the rest of the string will be c ignored. c c note --Add support for exponential notation: c c 2.60E+05 c c character*(*) iatemp logical ierror character*10 iformt c avalue=-9999.0 ntot=0 nleft=0 nright=0 c c Find first non-blank character c if (nchar.le.0) goto9090 if (iatemp(1:nchar).eq.' ') goto9090 if (nchar.eq.1) then ival=ichar(iatemp(1:1)) if (ival.le.47 .or. ival.gt.127) goto9090 endif c istrt=0 do 100 i=1,nchar if (iatemp(i:i).eq.' ')goto100 istrt=i goto199 100 continue goto9090 199 continue c c Find last non-blank character c ilast=nchar do 200 i=nchar,istrt,-1 if (iatemp(i:i).eq.' ')goto200 ilast=i goto299 200 continue goto9090 299 continue c c Check for "E" format. Handle this by extracting c the real number before "E" and then extracting c the power after the "E". c ipose=0 ipower=0 do 400 i=istrt,ilast if (iatemp(i:i).eq.'e' .or. iatemp(i:i).eq.'E')then icurr=i+1 ilast=i-1 ipose=i if (iatemp(icurr:icurr).eq.'+') then read(iatemp(i+2:i+3),'(I2)')ipower goto409 else if (iatemp(icurr:icurr).eq.'-') then read(iatemp(i+2:i+3),'(I2)')ipower ipower=-ipower goto409 else read(iatemp(i+1:i+2),'(I2)')ipower goto409 endif endif 400 continue 409 continue c c Convert text to real number c c Look for "+", "-", ".", or number. First character that c is not one of these will cause the rest of the field to c be ignored (this basically should strip off units). c Also, "+" and "-" only valid if the first element of the c field. Only one "." will be recognized. c iflag=0 do 300 i=istrt,ilast if (iatemp(i:i).eq.' ')then ilast=i-1 goto399 elseif (iatemp(i:i).eq.'.')then if (iflag.eq.0)then iflag=1 ntot=ntot+1 else ilast=i-1 goto399 endif elseif (iatemp(i:i).eq.'+')then if (ntot.eq.0)then afact=1.0 ntot=ntot+1 else ilast=i-1 goto399 endif elseif (iatemp(i:i).eq.'-')then if (ntot.eq.0)then afact=-1.0 ntot=ntot+1 else ilast=i-1 goto399 endif else itemp=ichar(iatemp(i:i)) if (itemp.lt.48 .or. itemp.gt.57)then ilast=i-1 goto399 else if (iflag.eq.0)then ntot=ntot+1 nleft=nleft+1 else ntot=ntot+1 nright=nright+1 endif end if end if 300 continue 399 continue iformt='(f . )' nmaxd=9 if (nright.gt.nmaxd)then nsubt=nright-nmaxd ilast=ilast-nsubt ntot=ntot-nmaxd nright=nmaxd endif write(iformt(3:4),'(i2)')ntot write(iformt(6:6),'(i1)')nright c read(iatemp(istrt:ilast),iformt,err=9010)avalue c c Now check for exponential power c if (ipose.gt.0 .and. ipower.ne.0) then avalue=avalue*10**ipower endif c goto9090 c 9010 continue write(6,*)'Error: Unable to convert real field' goto9090 c 9090 continue return end