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
