!This program calculates statistics from information in wind_gen.wdb (histogram version): !percentage erosive winds !wind power density (wpd) !wind energy !'wind erosion force' !preponderance !The latter two are calculated after: !E.L. Skidmore. 1965. Assessing wind erosion forces: Directions and relative magnitudes. !Soil Sci. Soc. Am. Proc. 29(5): 587-590. !Monthly values may be printed. !Example: !stats_from_histo -t8 < wind_gen.wdb.his > result.out !stats_from_histo -t10 -w0.2 < wind_gen.wdb.his > result.out !mo month (1..12) !d wind direction (1..16) !n number of cutoff wind speeds (n=27) !v(27) cutoff wind speeds (cdf) v(1)=0 m/s, v(2)=0.5 m/s !v2(26) central wind speeds (pdf) !f(12,16,n) cumulative freq. < cutoff wind speed (cdf) f(1)=f(2)=0 !p(12,16,n-1) normalized freq. at central wind speed (pdf) !df(12,16) frequencies for wind direction !ero_yr % erosive winds for the year (for a station) !wpd_yr wpd for the year (for a station) !ero_mo(12) % erosive winds for each month !wpd_mo(12) wpd for each month real, parameter :: pi=3.141592654 real latitude, longitude real r(16) real rr_max(12), anglem(12), force_max(12), angle_force_max(12) real df(12,16), calm(12), v_mean_mo(12), ero_mo(12), wpd_mo(12), year_fraction(12) real v(27), v2(26), interp real f(12,16,100) !cumulative freq. < cutoff wind speed (cdf) real p(12,16,100) !freq. at central wind speed (pdf) integer monl(12), mo, d, monthly, constant_density integer lat_deg, lat_min, lon_deg, lon_min integer pwd(12) !prevailing wind direction for each month integer pwed(12) !direction (one of 16) with max. erosion force for each month integer angle character*40 pound, wmoid, name, country, state character*25 argv character*3 moname(12) !names of 12 months character*5 dname(16) !names of 16 wind directions character lat_NS, lon_EW data v /0.0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5, & & 11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5, & & 25.5,30.5,35.5,40.5,45.5/ data monl /31,28,31,30,31,30,31,31,30,31,30,31/ data moname /'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG', & & 'SEP','OCT','NOV','DEC'/ data dname /'North','NNE','NE','ENE','East','ESE','SE','SSE', & & 'South','SSW','SW','WSW','West','WNW','NW','NNW'/ n=27 vt = 10 !default threshold wind speed (m/s) width = 0.1 !default integration width (m/s) nstat = 20000 !max. number of stations to be processed density = 1.2 !air density = 1.2 kg/m3 at 20 degr. C. if (iargc() .gt. 0) then do 09 i = 1, iargc() call getarg(i,argv) !command line help prompt if ((argv(2:2).eq.'?').or.(argv(2:2).eq.'h')) then write(*,*) 'Valid command line options:' write(*,*) '-? or -h Display this help screen' write(*,*) '-t# threshold wind speed in m/s, e.g. -t8' write(*,*) ' default threshold wind speed is 10 m/s' write(*,*) '-w# integration width in m/s, e.g. -w0.2' write(*,*) ' default integration width is 0.1 m/s' write(*,*) '-d use constant air density of 1.2 kg/m3' write(*,*) '-m output by month' stop else if (argv(2:2) == 't') then read(argv(3:),*) vt !threshold wind speed (m/s) else if (argv(2:2) == 'w') then read(argv(3:),*) width !integration width (m/s) else if (argv(2:2) == 'd') then constant_density = 1 else if (argv(2:2) == 'm') then monthly = 1 end if 09 continue endif do 55 i=1,12 !year fraction ~ 1/12 for each month year_fraction(i) = monl(i)/365.0 55 continue if (monthly .ne. 1) then write(*,*) write(*,*) 'WMO ID Station Name ST lat. lon. Uavg U > Ut WPD Energy' write(*,*) ' (degr.) (degr.) (m/s) (%) (W/m2) (kJ/m2/d)' end if do 999 ns=1,nstat !do first nstat stations !initialize cumulative frequencies do 77 i = 1, 12 do 78 j = 1, 16 do 79 k = 3, n f(i,j,k)=1.0 79 continue 78 continue 77 continue do 5 i = 1, 12 do 17 j = 1, 16 f(i,j,1) = 0 !freq. corresponding to v(1) = 0.0 m/s f(i,j,2) = 0 !freq. corresponding to v(2) = 0.5 m/s 17 continue 5 continue !read everything in read (*,10,end=1000) pound, wmoid, country, state, name 10 format (a2, a6,x, a2,x,a2,x,a20) !read station latitude, longitude, and elevation (m) read (*,*) lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW, elev latitude = lat_deg + lat_min/60.0 if (lat_NS=='S') latitude = -1*latitude longitude = lon_deg + lon_min/60.0 if (lon_EW=='W') longitude = -1*longitude !read direction frequencies do 111 i=1,16 read (*,*) (df(j,i),j=1,12) 111 continue read (*,*) (calm(i),i=1,12) !read calm frequencies !read all histogram info do 21 i = 1, 12 do 22 j = 1, 16 do 23 k = 3, n read(*,49,EOR=22,ADVANCE='NO') f(i,j,k) 49 format(f6.3) 23 continue 22 continue 21 continue !skip 2 lines read (*,*) !max/min ratio read (*,*) !hour of the day when max. occurs if (constant_density .ne. 1) then !calculate air density from station elevation temp = 20.0 !assume an avg. air temp. of 20 degr. C density = 348.56 * (1.013 - 0.1183*(elev/1000) + 0.0048*(elev/1000)**2) / (temp+273.1) end if !determine prevailing wind direction for each month pwd = maxloc(df,dim=2) !calculate central wind speeds for pdf do 25 i = 1, n-1 v2(i)= (v(i)+v(i+1))/2 25 continue !calculate pdf from cdf do 26 i = 1, 12 do 27 j = 1, 16 do 28 k = 1, n-1 !need to normalize (divide by bin width), !so that area under pdf 'curve' = 1 p(i,j,k) = (f(i,j,k+1)-f(i,j,k)) / (v(k+1)-v(k)) 28 continue 27 continue 26 continue if (monthly==1) then write(*,*) write(*,24) wmoid, name, state, int(elev), lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW 24 format(x, a6, x, a20, a2, ' Elevation:', i5, ' m Latitude:',2i3,a2,' Longitude:'i4,i3,a2) write(*,*) ' Month Uavg U > Ut WPD Energy Fract Max16 PWED1 PWED2 Prepond PWED3' write(*,*) ' (m/s) (%) (W/m2) (kJ/m2/d) (%) (degr.) (degr.)' end if v_mean_yr = 0 ero_yr = 0 wpd_yr = 0 do 555 mo=1,12 !12 months v_mean_mo(mo) = 0 ero_mo(mo) = 0 wpd_mo(mo) = 0 do 444 d=1,16 !16 directions !calculate % EROSIVE WINDS do 70 i = 2, n if (v(i) > vt) then !find fraction erosive winds by linear interpolation of the cdf interp =(vt-v(i-1))/(v(i)-v(i-1)) ff = f(mo,d,i-1) + interp*(f(mo,d,i)-f(mo,d,i-1)) ero = 1 - ff goto 80 end if 70 continue 80 continue !combine all direction combinations for each month ero_mo(mo) = ero_mo(mo) + ero*df(mo,d) !calm is accounted for !combine all month and direction combinations ero_yr = ero_yr + ero*df(mo,d)*year_fraction(mo) !calm is accounted for !calculate WPD and wind erosion force !calculate for first histogram bar do 73 i = 2, n if (v(i) > vt) then wpd = wpd_bar(vt, v(i), vt, p(mo,d,i-1), width, density) !integrate by adding !calculate wind erosion force v_avg = (vt+v(i))/2 !use average speed for the bar fractie = (v(i)-vt) / (v(i)-v(i-1)) !weigh on fraction of the bar r(d) = (v_avg**3) * p(mo,d,i-1) * df(mo,d) * fractie goto 83 end if 73 continue 83 continue !calculate for all following histogram bars do 93 j = i+1, n wpd = wpd + wpd_bar(v(j-1), v(j), vt, p(mo,d,j-1), width, density) !integrate by adding r(d) = r(d) + (v2(j-1)**3) * p(mo,d,j-1) * df(mo,d) !wind erosion force if (f(mo,d,j) .ge. 1.0) goto 96 !no need to look at bars of height=0 93 continue 96 continue wpd_mo(mo) = wpd_mo(mo) + wpd*df(mo,d)/100 !calm is accounted for wpd_yr = wpd_yr + wpd*df(mo,d)/100*year_fraction(mo) !calm is accounted for !calculate mean wind speed for each month and for the year v_mean = 0 do 60 i = 1, n-1 v_mean = v_mean + v2(i) * p(mo,d,i) 60 continue !combine all direction combinations for each month v_mean_mo(mo) = v_mean_mo(mo) + v_mean*df(mo,d)/100 !combine all month and direction combinations v_mean_yr = v_mean_yr + v_mean*df(mo,d)*year_fraction(mo)/100 444 continue !do next direction !account for calm winds; the average for calm is 0.25 m/s v_mean_mo(mo) = v_mean_mo(mo) + 0.25*calm(mo)/100 v_mean_yr = v_mean_yr + 0.25*calm(mo)*year_fraction(mo)/100 !determine direction (one of 16) with max. erosion force for each month pwed(mo) = maxloc(r,dim=1) rr_max(mo) = 0 force_max(mo) = 0 do 888 angle = 0, 359, 1 par = 0 per = 0 force = 0 do 777 d = 1, 16 angle2 = (d - 1)*22.5 - angle angle2 = angle2*2*pi / 360 !degrees to radians par = par + r(d) * abs(cos(angle2)) !parallel forces, eq 7 per = per + r(d) * abs(sin(angle2)) !perpendicular forces, eq 8 if (cos(angle2) > 0) then force = force + r(d) * cos(angle2) !parallel forces, only positive ones end if 777 continue !find direction of max. erosion forces if (force > force_max(mo)) then force_max(mo) = force angle_force_max(mo) = angle end if !find direction of max. RATIO of parallel to perpendicular erosion forces rr = par / per !eq. 9 if (rr > rr_max(mo)) then rr_max(mo) = rr !find which of the opposite directions has the greatest force if (force > 0.5*par) then anglem(mo) = angle else anglem(mo) = angle + 180 if (anglem(mo) > 360) anglem(mo) = anglem(mo) - 360 end if end if ! write(*,711) angle, par, per, rr !711 format (3f7.1, f12.6) 888 continue !do next angle 555 continue !do next month if (monthly==1) then !print wind statistics by month do 666 mo = 1, 12 energy = 3.6*24*wpd_mo(mo) !convert from W/m2 to kJ/m2/day percent = 100*wpd_mo(mo) / wpd_yr / 12 !percentage of the year write(*,35) moname(mo),v_mean_mo(mo),ero_mo(mo),wpd_mo(mo),int(energy),percent,dname(pwd(mo)), & & dname(pwed(mo)),int(angle_force_max(mo)),rr_max(mo),int(anglem(mo)) 35 format(a7, 2f8.1, f10.1, i10, f9.1, a12, a9 i8, x, f11.1, i8) 666 continue end if !print wpd and wind energy by station energy = 3.6*24*wpd_yr !convert from W/m2 to kJ/m2/day if (monthly==1) then write(*,39,ADVANCE='NO') 39 format(' year') write(*,34) v_mean_yr, ero_yr, wpd_yr, int(energy) 34 format(2f8.1, f10.1, i10, ' 100.0') else write(*,37,ADVANCE='NO') wmoid, name, state, latitude, longitude 37 format(x, a6, 3x, a20, a2, 2f9.2) write(*,36) v_mean_yr, ero_yr, wpd_yr, int(energy) 36 format(x, 2f8.1, f10.1, i10) end if 999 continue !do next station 1000 continue write(*,*) write(*,*) 'Uavg Average wind speed (m/s)' write(*,*) 'U > Ut Percentage Erosive Winds (percentage of winds above the threshold wind speed)' write(*,*) 'WPD Erosive Wind Power Density (W/m2)' write(*,*) 'Energy Erosive Wind Energy (kJ/m2/d)' if (monthly==1) then write(*,*) 'Fract Fraction of Annual Erosive Wind Enery' write(*,*) 'Max16 Prevailing Wind Direction (one of 16 cardinal directions)' write(*,*) 'PWED1 Prevailing Wind Erosion Direction (one of 16 cardinal directions)' write(*,*) 'PWED2 Prevailing Wind Erosion Direction (range = 0-359 degrees)' write(*,*) ' = direction of max. erosion forces' write(*,*) 'Prepond Preponderance (max. ratio of parallel to perpendicular erosion forces)' write(*,*) 'PWED3 Prevailing Wind Erosion Direction (range = 0-359 degrees)' write(*,*) ' = direction of max. ratio of parallel to perpendicular erosion forces' end if write(*,*) write(*,*) 'WPD = 0.5 r (U-Ut) U^2' write(*,*) 'where' if (constant_density==1) then write(*,*) ' r air density (1.2 kg/m3)' else write(*,*) ' r air density (function of station elevation)' end if write(*,*) ' U wind speed (m/s)' write(*,*) ' Ut threshold wind speed' write(*,*) write(*,*) 'U > Ut: based on cumulative distribution' write(*,*) 'WPD and Energy: based on probability distribution' write(*,*) write(*,13) n 13 format (' The number of cutoff wind speeds is',i3) write(*,11) vt 11 format (' The threshold wind speed is',f6.2, ' m/s') write(*,*) end function wpd_bar(v1, v2, vt, p, width_desired, density) !Calculates the wind power density (wpd, W/m2) for one histogram bar !by numerical integration !ARGUMENT DEFINITIONS: !v1 left border of histogram bar (m/s) !v2 right border of histogram bar (m/s) !vt threshold wind speed (m/s) !p normalized freq. (pdf) of histogram bar !width_desired desired integration width !density air density (kg/m3) !LOCAL DEFINITIONS: !n number of rectangles for one bar barwidth = v2 - v1 n = nint(barwidth/width_desired) if (n < 1) n = 1 width = barwidth/n !actual width may deviate a little from desired width wpd = 0 do 10 i = 1, n v = v1 + (i-0.5)*width wpd = wpd + width*p*0.5*density*(v-vt)*v**2 !integrate by adding rectangles 10 continue wpd_bar = wpd return end