c file: newdb.for c - program to calculate weibull parameters c and extract direction distribution, c and weather parameters. c - data is read in from WERIS tables c initialize character line1*23, larr(23), tab*13, id*9, runfil*35, c line2*117, l5arr(117), unit*6, line3*56, spgr(27)*9 equivalence (line1,larr(1)), (unit,l5arr(63)), e (tab,larr(11)), (line3,l5arr(1)), (id,larr(2)) integer i, ii(17), ipwed(12), i j, i k, kk, i l, i m, mxhr(12) logical fexist real f(16,19), ht, d(20), r pwed(12), prepon(12), pospar(12) double precision a(27,27), r b, r c(25), cc(17,12), ck(17,12), r dir(18,12), den(12), r jj(17), r norm(27,27), r rat(12), r snx, snx2, sny, snxy, sum(18), wt(18,28) c definitions c a - frequency of occurrence of wind speed by wind direction (%) c b - intercept in regression c c - speed group c cc - weibull scale parameter c ck - weibull shape parameter c dir - frequency of wind direction c den - air density for the month c f - frequency of occurrence of speeds over 7 m/s c i - loop counter on k - direction c j - loop counter on l - speed group c jj - flag for speed group where only one group has frequency c k - number of directions c l - number of speed groups c m - number of months c ht - annemometer height c mxhr - hour at which the max wind speed occurrs c norm - normalized frequencies within a direction c rat - ratio of max to min wind speed c sn* - regression variables c sum - sum of frequencies within a direction c wt - weighting factor c input formats 10 format (a117) 11 format (a9,18f6.1) 12 format (27f5.1) 13 format (a23) 14 format (a56,f6.1,a6) 205 format (a61,i2,a21,f4.1) 305 format (a35,12f7.2) 315 format (a35,12i7) c output formats 61 format (12f6.2) 62 format (a9) 63 format (12f5.1) 68 format (12i5) 69 format (12f5.2) 71 format (16f4.1) 920 FORMAT(12i4) c open files c open (unit = 11, access = 'append', file = 'new.db') c open (unit = 7, access = 'append', file = 'rweq.db') open (unit = 11, access = 'append', file = 'jnk.db') open (unit = 7, access = 'append', file = 'rjnk.db') write(*,*) 'please enter input file name' 5 read (*,'(a)') runfil inquire(file = runfil, exist = fexist) if(.not.fexist) then open (unit = 9, access = 'append', file = 'error.db') write(9,*) runfil, '- file not found' write(11,*) ' newdb error:',runfil,' - WERIS file not found' write(7,*) ' rweqdb error:',runfil,' - WERIS file not found' goto 999 endif open (unit = 10,file = runfil) 1 read (10,13,end=999) line1 if (tab .eq. '.05 ') goto 200 if (tab .eq. '.10 ') goto 300 if (tab .ne. '.12A ') goto 1 c do table # 12 do 100 m = 1,1 if (m .eq. 1) goto 16 read (10,13) line1 read (10,13) line1 16 do 50 i=1,3 read (10,10) line2 50 continue read (10,14) line3,ht,unit do 51 i = 4,13 read (10,10) line2 51 continue c convert feet to meters if needed if (unit .ne. 'METERS') ht = ht * 0.3048 l = 25 k = 17 do 20 j = 1,26 read(10,11) spgr(j), (a(i,j),i = 1,18) 20 continue c get frequencies for speed >= 7 c for prepon4.for do 42 j = 7,25 do 43 i = 1,16 kk=j-6 f(i,kk) = a(i,j) 43 continue 42 continue do 6 j = 1,25 a(17,j) = a(18,j) 6 continue do 7 i = 1,k dir(i,m) = a(i,26) 7 continue c if only one speed group c has freq, then flag do 21 i = 1,k ii(i) = 0 do 26 j = 1,l if (a(i,j) .ne. 0.) then ii(i) = ii(i) + 1 if (j .le. 20) then jj(i) = j else if (j .eq. 21) then jj(i) = j + 2 else if (j .eq. 22) then jj(i) = j + 6 else if (j .eq. 23) then jj(i) = j + 10 else if (j .eq. 24) then jj(i) = j + 14 else if (j .eq. 25) then jj(i) = j + 18 end if end if 26 continue 21 continue c end points of speed groups do 8 j = 1,20 8 c(j) = j + 0.5 c(21) = 25.5 do 9 j = 22,25 9 c(j) = c(j-1) + 5.0 c adjust for anamometer height do 18 j = 1,25 c(j) = c(j) * (10./ht)**0.142857 c write(11,*) 'x(',j,')=',c(j) c log of speed group end points c(j) = dlog(c(j)) 18 continue c data for yi - ? do 24 i = 1,k c sum sum(i) = 0. do 22 j = 1,l sum(i) = sum(i) + a(i,j) c d(j) = dlog(j*.1) c write(11,*) 'dlog(',(j*.1),')=',d(j) 22 continue c if no winds for this direction if (ii(i) .eq. 0) sum(i) = .01 c normalize for each direction do 23 j = 1,l norm(i,j) = a(i,j) / sum(i) 23 continue 24 continue c data for yi and ni - ? do 32 i = 1,k c cumulate frequencies 25 do 30 j = 1,l wt(i,j) = norm(i,j) if(j .gt. 1) norm(i,j) = wt(i,j) + norm(i,j-1) c write(11,9998)i,j,wt(i,j),i,j,norm(i,j),j,c(j) c9998 format('wt(',i2,',',i2,')=',f5.2,' norm(',i2,',',i2,')=',f5.2, &' logc(',i2,')=',f5.2) 30 continue c double logarithm sum(i) = norm(i,25) do 31 j = 1,l if (norm(i,j) .lt. 0.000001) norm(i,j) = 0.000001 if (norm(i,j) .gt. 0.999999) norm(i,j) = 0.999999 norm(i,j) = dlog ( - (dlog (1.0 - norm(i,j) ) ) ) write(11,9999) i,j,norm(i,j) 9999 format(' dblognorm(',i2,',',i2,')=',f5.2) 31 continue 32 continue c regression do 40 i = 1,k c if no winds for this direction if (ii(i) .eq. 0) go to 47 snx2 = 0.0 snx = 0.0 sny = 0.0 snxy = 0.0 do 39 j = 1,l sny = sny + norm(i,j) * wt(i,j) snx = snx + c(j) * wt(i,j) snx2 = snx2 + c(j) * c(j) * wt(i,j) 39 snxy = snxy + norm(i,j) * c(j) * wt(i,j) c numerator for a snxy = snxy - snx * sny / sum(i) c denomonator for a snx2 = snx2 - snx * snx / sum(i) ck(i,m) = snxy / snx2 b = (sny - ck(i,m) * snx) / sum(i) cc(i,m) = dexp (-b / ck(i,m)) c if only one sp group for this c direction - assign another cc if (ii(i) .eq. 1) goto 46 go to 45 46 cc(i,m) = 1.12 * jj(i) c and assign ck accordingly ck(i,m) = 2. * cc(i,m) go to 45 c if all frequencies = 0. for this c direction then c = 0. and k = 1. 47 cc(i,m) = 0. ck(i,m) = 1. 45 continue 40 continue c get rweq database wind info. call prepon4 (m, f, ht, pwed, prepon, pospar) 100 continue c output to WERM db c write (11,62) id do 80 i = 1,k-1 80 write (11,63) (dir(i,m),m = 1,12) write (11,63) (dir(17,m),m = 1,12) do 70 i = 1,k-1 70 write (11,61) (cc(i,m),m = 1,12) do 75 i = 1,k-1 75 write (11,61) (ck(i,m),m = 1,12) write (11,63) (rat(i),i = 1,12) write (11,68) (mxhr(i),i = 1,12) c output to RWEQ db c write (7,62) id write (7,61) (cc(17,m),m = 1,12) write (7,61) (ck(17,m),m = 1,12) write (7,69) (den(i),i = 1,12) do 85 kk=1,12 ipwed(kk) = pwed(kk) c write(*,*) ipwed(i) 85 continue WRITE(7,920)(iPWED(KK),KK=1,12) WRITE(7,63)(PREPON(KK),KK=1,12) WRITE(7,69)(POSPAR(KK),KK=1,12) go to 1 c read table # 05 200 do 210 i = 1,8 210 read (10,10) line2 do 220 i = 1,12 220 read (10,205) line3,mxhr(i),line3,rat(i) goto 1 c read table # 10 300 do 310 i = 1,7 310 read (10,10) line2 read (10,305) line3,(den(i),i = 1,12) goto1 999 stop end