!*************************************************************** ! This program adjusts a file of observation wind values to ! 1-min averaged and 10-meter height. ! ! Written by Yee Lau (lau@gri.msstate.edu) 9/16/15 ! and Dr. Pat Fitzpatrick (fitz@gri.msstate.edu) ! ! Affiliation: ! ! Dr. Pat Fitzpatrick, Associate Research Professor of Meteorology ! GeoSystems Research Institute ! Mississippi State University ! ! References: ! ! Vickery, P. J. and Skerlj, P. F., "Hurricane Gust Factors Revisited", ! J. Structural Engineering - ASCE, Vol. 131, No. 5, May 2005, pp. 825-832. ! ! Krayer, W. R. and Marshall, R. D., "Gust Factors Applied to Hurricane Winds", ! Bull. Am. Meteorol. Soc., 74(5), 1992, p. 613 - 617. ! ! Also see: http://ijsetr.com/uploads/453126Ma%20Tin%20Nilar%20Tun.doc ! ! To convert to 1-min averaged wind, the following conversion factors are used: ! Marine ! 1 hr - 1.17 ! 10 minute - 1.11 ! 8 minute - 1.08 ! 3 minute - 1.05 ! 2 minute - 1.03 ! Land ! 1 hr - 1.22 ! 10 minute - 1.16 ! 3 minute - 1.10 ! 2 minute - 1.08 ! ! To convert to 10-meter height, the following formula is used: ! wten = wref*(zten/zref)^0.143 ! ! According to NDBC, moored buoy wind is 8-minute averages, and ! cman buoy wind is 2-minute averages. ! ! RAWS's wind is 10-minute averages measured at a height of 20 feet (6.096m). ! ASOS wind is 2-minute averages measured at a height of mostly 33 feet (10.0584m). ! ! Usage: ! adjustObsWind1minAvg10mHeight ndbcInfoFile missingValInNDBCInfoFile asosInfoFile missingValInASOSInfoFile miscInfoFile missingValInMiscInfoFile inputDataFile missingValInDataFile outputDataFile ! ! Input Arguments: ! 1. NDBC's station information file (both NDBC and non-NDBC stations are included) ! This file must have 1 row of header followed by multiple rows ! of station information arranged in 6 columns as follows: ! i. station id (max 5 characters located at position 1 thru 5) ! ii. station owner (max 12 characters located at position 7 thru 18) ! iii station elevation ! iv. air temperature height ! v. anemometer height ! vi. barometer height ! ! All measurements are in meters. Missing value is -999. ! ! 2. Missing value in the NDBC station information file ! ! 3. ASOS tower information file ! This file must have 1 row of header followed by multiple rows ! of station information arranged in 4 columns as follows: ! i. station id (max 5 characters located at position 1 thru 5) ! ii. city (max 25 characters located at position 7 thru 31) ! iii state ! iv. tower height in feet, missing value is -999 ! ! 4. Missing value in the ASOS tower information file ! ! 5. Misc. tower information file ! This file must have 1 row of header followed by multiple rows ! of station information arranged in 3 columns as follows: ! i. station id (max 5 characters located at position 1 thru 5) ! ii. station owner (max 12 characters located at position 7 thru 18) ! iii. tower height in feet, missing value is -999 ! ! 6. Missing value in the misc. tower information file ! ! 7. input data file ! This file must have one or more rows of station wind information ! arranged in at least 6 columns as follows: ! i. station id (max 5 characters located at postiion 1 thru 5) ! ii statoin owner (max 12 characters located at position 7 thru 18) ! iii. wind direction in degrees ! iv. wind speed in meter per second ! v. latitude ! vi. longitude ! ! 8. Missing value in the input data file ! ! 9. output data file ! This file will have at least 7 columns as follows: ! i. station id ! ii. station owner ! iii. wind direction in degrees ! iv. 1-minute wind speed at 10 meters height in meter per second ! v. latitude ! vi. longitude ! vii. wind speed adjustment flag ! 0 if data is not adjusted because there's no height and time-averaged information, ! 1 if data is adjusted to 10-meter height only, and there's no time-averaged information, ! 2 if data is adjusted to 1-minute average only, and there's no height information, ! 3 if data is adjusted to 10-meter height and 1-minute average, ! 4 if data is already in 10-meter height, and there's no time-averaged information, ! 5 if data is already in 1-minute average, and there's no height information, ! 6 if data is already in 10-meter height, only time is adjusted to 1-minute average, ! 7 if data is already in 1-mnute average, only height is ajusted to 10-meter, ! 8 if data is already in 10-meter height and 1-minute average, and no adjustment is needed. ! ! !*************************************************************** Program adjustObsWind1minAvg10mHeight implicit none character*100 ndbcInfoFile, asosInfoFile, miscInfoFile character*100 inputDataFile character*100 outputDataFile character*200 line character*10 curDate character*20 missNDBCInfoStr, missASOSInfoStr, missMiscInfoStr, missDataStr character*5 stationName, curTime, curId, curState character*25 curCity character*12 ownerName, curOwner character*12, allocatable :: owner(:), miscOwner(:) character*5, allocatable :: id(:), asosId(:), miscId(:) character*4, asosIdNew character(len=53), parameter :: outfmt = "(A5,1x,A12,1x,F14.6,1x,F14.6,1x,F12.6,1x,F12.6,1x,I1)" real, allocatable :: aneHt(:), siteHt(:), airTempHt(:), baroHt(:) real, allocatable :: asosHt(:), miscHt(:) real cp, rho, cH, von, g, pi, zten, alpha character*4 yyStr character*2 mmStr, ddStr, hhStr, minStr integer yy,mm,dd,hh,min,numStations,numTowers,numMisc,cnt,ios,foundIndex,foundType integer iargc,i,pos,flag real z, wdir, wspd, adjwspd, missNDBCInfo, missASOSInfo, missMiscInfo, lat, lon, missData, height real xleft, xmid, xright, ustar, zo real MARINE_1hr_to_1min, MARINE_10min_to_1min, MARINE_8min_to_1min, MARINE_3min_to_1min, MARINE_2min_to_1min real LAND_1hr_to_1min, LAND_10min_to_1min, LAND_3min_to_1min, LAND_2min_to_1min, RAWS_Height, ASOS_Height real FEET_to_METER parameter(cp=1004,rho=1.275,cH=0.003,von=0.4) parameter(g=9.8,pi=3.14159,zten=10.0,alpha=0.015) parameter(MARINE_1hr_to_1min=1.17,MARINE_10min_to_1min=1.11,MARINE_8min_to_1min=1.08) parameter(MARINE_3min_to_1min=1.05,MARINE_2min_to_1min=1.03) parameter(LAND_1hr_to_1min=1.22,LAND_10min_to_1min=1.16,LAND_3min_to_1min=1.10,LAND_2min_to_1min=1.08) parameter(RAWS_Height=6.096,ASOS_Height=10.0) parameter(FEET_to_METER=0.3048) external f INTERFACE FUNCTION is_numeric (string) logical is_numeric CHARACTER(len=*), INTENT(IN) :: string END FUNCTION is_numeric END INTERFACE INTERFACE FUNCTION to_upper (str) result (string) Character(*), Intent(In) :: str Character(LEN(str)) :: string END FUNCTION to_upper END INTERFACE !*************************************************************** ! check the number of input arguments and assign them to variables !*************************************************************** if (iargc() .ne. 9) then print *, 'usage:adjustObsWind1minAvg10mHeight ndbcInfoFile missingValInNDBCInfoFile ' print *, 'asosInofFile missingValInASOSInfoFile miscInfoFile missingValInMiscInfoFile ' print *, 'inputDataFile missingValInDataFile outputDataFile' stop endif call getarg(1,ndbcInfoFile) call getarg(2,missNDBCInfoStr) call getarg(3,asosInfoFile) call getarg(4,missASOSInfoStr) call getarg(5,miscInfoFile) call getarg(6,missMiscInfoStr) call getarg(7,inputDataFile) call getarg(8,missDataStr) call getarg(9,outputDataFile) read(missNDBCInfoStr, *) missNDBCInfo print *, 'missingValInNDBCInfoFile=', missNDBCInfo read(missASOSInfoStr, *) missASOSInfo print *, 'missingValInASOSInfoFile=', missASOSInfo read(missMiscInfoStr, *) missMiscInfo print *, 'missingValInMiscInfoFile=', missMiscInfo read(missDataStr, *) missData print *, 'missingValInDataFile=', missData !******************************************************** ! read all NDBC station information !******************************************************** open(unit=10,file=ndbcInfoFile,status='unknown') numStations=0 ! ignore header line read(10, *)line do read(10,*,iostat=ios)line if (ios .ne. 0) exit numStations=numStations+1 enddo rewind(10) print *, 'numStationsInNDBCInfoFile=', numStations ! ignore header line read(10, *)line allocate(id(numStations), owner(numStations)) allocate(siteHt(numStations), aneHt(numStations), airTempHt(numStations), baroHt(numStations)) cnt=0 do i=1,numStations read(10,*,iostat=ios)curId,curOwner,siteHt(i),airTempHt(i),aneHt(i),baroHt(i) id(i) = adjustl(to_upper(curId)) owner(i) = adjustl(to_upper(curOwner)) if (ios .ne. 0) exit cnt=cnt+1 enddo close(10) if (cnt .ne. numStations) then print *, "cnt=", cnt print *, "Program Error! Number of NDBC stations not matching!" stop endif !******************************************************** ! read all ASOS tower information !******************************************************** open(unit=15,file=asosInfoFile,status='unknown') numTowers=0 ! ignore header line read(15, *)line do read(15,*,iostat=ios)line if (ios .ne. 0) exit numTowers=numTowers+1 enddo rewind(15) print *, 'numTowersInASOSInfoFile=', numTowers ! ignore header line read(15, *)line allocate(asosId(numTowers), asosHt(numTowers)) cnt=0 do i=1,numTowers read(15,*,iostat=ios)curId,curCity,curState,height asosId(i) = adjustl(to_upper(curId)) asosHt(i) = height * FEET_to_METER if (ios .ne. 0) exit cnt=cnt+1 enddo close(15) if (cnt .ne. numTowers) then print *, "cnt=", cnt print *, "Program Error! Number of ASOS towers not matching!" stop endif !******************************************************** ! read all misc. tower information !******************************************************** open(unit=18,file=miscInfoFile,status='unknown') numMisc=0 ! ignore header line read(18, *)line do read(18,*,iostat=ios)line if (ios .ne. 0) exit numMisc=numMisc+1 enddo rewind(18) print *, 'numTowersInMiscInfoFile=', numMisc ! ignore header line read(18, *)line allocate(miscId(numMisc), miscOwner(numMisc), miscHt(numMisc)) cnt=0 do i=1,numMisc read(18,*,iostat=ios)curId,curOwner,height miscId(i) = adjustl(to_upper(curId)) miscOwner(i) = adjustl(to_upper(curOwner)) miscHt(i) = height * FEET_to_METER if (ios .ne. 0) exit cnt=cnt+1 enddo close(18) if (cnt .ne. numMisc) then print *, "cnt=", cnt print *, "Program Error! Number of misc. towers not matching!" stop endif !******************************************************** ! Open input data file to read and output file to write !******************************************************** open(unit=20,file=inputDataFile,status='unknown') open(unit=30,file=outputDataFile,status='unknown') !******************************************************** ! Loop to check if id exists in NDBC station information file ! or ASOS tower information file. Also, check owner for RAWS, ASOS etc. ! If height or time-average information exist, then adjust wind to ! 10-meter height and 1-minute average accordingly, and set adjustment flag. !******************************************************** cnt=0 do read(20,*,iostat=ios) curId, curOwner, wdir, wspd, lat, lon if (ios .ne. 0) exit if (wspd == missData) then cycle endif stationName = adjustl(to_upper(curId)) ownerName = adjustl(to_upper(curOwner)) cnt = cnt + 1 print * write(*,"(I4.4,A)") cnt, '--------------------------------------------------' write(*,outfmt) stationName, ownerName, wdir, wspd, lat, lon, 0 foundIndex = 0 foundType = 0 ! check station info from NDBC do i=1,numStations if (trim(id(i)) .eq. trim(stationName)) then foundIndex = i foundType = 1 exit endif end do ! check station info from ASOS if (foundIndex .eq. 0) then do i=1, numTowers if (trim(asosId(i)) .eq. trim(stationName)) then foundIndex = i foundType = 2 exit else ! handle special cases in ASOS file ! case 1: for example, GPT in ASOS instead of KGPT if (len_trim(asosId(i)) .eq. 3) then asosIdNew = 'K'//trim(asosId(i)) if (asosIdNew .eq. trim(stationName)) then foundIndex = i foundType = 2 exit endif ! case 2: ASOS file has 3 KBIX towers (KBIX1, KBIX2 or KBIX3). Fortunately, they are all at the same height. else if (trim(stationName) .eq. 'KBIX' .and. asosId(i)(1:4) .eq. 'KBIX') then foundIndex = i foundType = 2 exit endif endif end do endif ! check station info from misc. if (foundIndex .eq. 0) then do i=1, numMisc if (trim(miscId(i)) .eq. trim(stationName)) then foundIndex = i foundType = 3 exit endif end do endif !******************************************************** ! if id is not found, check for RAWS or ASOS, and adjust accordingly. ! Otherwise, ajust the height and time-average based on available information. ! Also, check for MARITIME data and adjust to 1-minute. !******************************************************** flag = 0 print *, 'foundIndex=', foundIndex ! id is not found, check for RAWS or ASOS, and adjust if (foundIndex == 0) then if (trim(ownerName) .eq. 'RAWS') then adjwspd = wspd*(zten/RAWS_Height)**0.143 adjwspd = adjwspd * LAND_10min_to_1min print *, 'Is RAWS, adjusted to 1min and 10m = ',adjwspd flag = 3 else if (trim(ownerName) .eq. 'ASOS') then adjwspd = wspd * LAND_2min_to_1min print *, 'Is ASOS, already in 10m, adjusted to 1min = ',adjwspd flag = 6 else adjwspd = wspd print *, 'No station height or average information, no adjustment = ',wspd flag = 0 endif write(30,outfmt) adjustl(stationName), adjustl(ownerName), wdir, adjwspd, lat, lon, flag cycle endif ! id is found in NDBC station information file if (foundType == 1) then print *, 'aneHt=', aneHt(foundIndex) if (aneHt(foundIndex) .eq. missNDBCInfo) then adjwspd = wspd print *, 'No station height information = ',wspd flag = 0 else if (aneHt(foundIndex) .ne. 10) then height = aneHt(foundIndex) if (siteHt(foundIndex) .ne. missNDBCInfo) then height = height + siteHt(foundIndex) endif adjwspd = wspd*(zten/height)**0.143 flag = 1 print *, 'Found aneHt, adjusted 10m = ', adjwspd else print *, 'Found aneHt at 10m = ', wspd adjwspd = wspd flag = 4 endif if (trim(ownerName) .eq. 'MARITIME') then if (is_numeric(stationName)) then adjwspd = adjwspd * MARINE_8min_to_1min print *, 'MARITIME, is numeric, adjusted 8min to 1min = ', adjwspd else adjwspd = adjwspd * MARINE_2min_to_1min print *, 'MARITIME, not numeric, adjusted 2min to 1min = ', adjwspd endif flag = flag + 2 endif else if (foundType == 2) then ! id is found in ASOS tower information file print *, 'asosHt=', asosHt(foundIndex) if (asosHt(foundIndex) .eq. missASOSInfo) then adjwspd = wspd print *, 'No tower height information = ',wspd flag = 0 else if (asosHt(foundIndex) .ne. 10) then adjwspd = wspd*(zten/asosHt(foundIndex))**0.143 flag = 1 print *, 'Found asosHt, adjusted 10m = ', adjwspd else print *, 'Found asosHt at 10m = ', wspd adjwspd = wspd flag = 4 endif adjwspd = adjwspd * LAND_2min_to_1min print *, 'ASOS, adjusted from 2min to 1min = ', adjwspd flag = flag + 2 else if (foundType == 3) then ! id is found in misc. station information file print *, 'miscHt=', miscHt(foundIndex) if (miscHt(foundIndex) .eq. missMiscInfo) then adjwspd = wspd print *, 'No misc. tower height information = ',wspd flag = 0 else if (miscHt(foundIndex) .ne. 10) then adjwspd = wspd*(zten/miscHt(foundIndex))**0.143 flag = 1 print *, 'Found miscHt, adjusted 10m = ', adjwspd else print *, 'Found miscHt at 10m = ', wspd adjwspd = wspd flag = 4 endif endif write(30,outfmt) adjustl(stationName), adjustl(ownerName), wdir, adjwspd, lat, lon, flag if (ios < 0) exit enddo close(20) close(30) deallocate(id, asosId, owner, aneHt, airTempHt, siteHt, baroHt, asosHt, miscId, miscOwner, miscHt) print * print *, 'numObsInDataFile=', cnt end subroutine bisect(f,xleft,xright,xmid,a,b) !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! Bisection program ! ! Written by Pat Fitzpatrick, Jackson State University ! ! Originally written February 14, 1996 ! Updated May 14, 1998 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! This program calculates the root of a function by solving f(x)=0 ! with an error tolerance of 0.001 (this can be modified by changing ! eps). For example, if the user wants to find the cube root of 10, ! f(x) = x**3 - 10 ! ! User must define the function f(x) and the interval for within which ! the root exists (xleft & xright) in a function subprogram. ! !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! Variables used: ! ! f(x) --- function to be solved for f(x)=0 ! eps --- error tolerance ! xleft --- left side of interval ! xright --- right side of interval ! xmid --- midpoint of interval ! !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !********************************************************************* ! Set error tolerance !********************************************************************* eps=0.001 !********************************************************************* ! Check to see if interval properly chosen. If ! f(xright,a,b)*f(xleft,a,b) > 0 , ! then interval improperly chosen because either xleft and xright are ! both left of the root or both right of the root ! ! If interval incorrect, program stops, and a message prints to ! screen notifying the user of the problem. !********************************************************************* if( f(xright,a,b)*f(xleft,a,b).gt.0.0) then print*,f(xright,a,b),f(xleft,a,b) print*,a,b print*,'Improper interval chosen; program stops' stop endif !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! Main bisection routine exists here. There are two separate ! sections, since either f(xleft,a,b) < 0 and f(xright,a,b) > 0 , ! or f(xleft,a,b)) > 0 and f(xright,a,b)) > 0 . ! If after 40 iterations no ! answer is found, then either eps has been set too small or ! more loops are required, and the program stops with a warning ! message. ! ! The steps are as follows: ! ! 1) The midpoint value of the interval "xmid" is computed. ! ! 2) f(xmid,a,b) is then computed. ! ! 3) If f(xmid,a,b) < absolute value of error tolerance, ! calculation is completed. Goto line 20 and print the ! answer. Otherwise, continue to 4) ! ! 4) Insert xmid into the interval limit which has the same function ! sign as f(xmid,a,b) . Goto part 1. !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! This section is used if f(xleft,a,b)<0 and f(xright,a,b)>0 !********************************************************************* if(f(xleft,a,b).lt.0.0) then do i=1,40 xmid=(xright+xleft)/2.0 if(f(xmid,a,b).lt.eps.and.f(xmid,a,b).gt.(-1.0*eps))goto 20 if(f(xmid,a,b).lt.0.0) then xleft=xmid else xright=xmid endif enddo else !********************************************************************* ! Otherwise this section is used since f(xleft,p,thetae)>0 and ! f(xright,p,thetae)<0 !********************************************************************* do i=1,40 xmid=(xright+xleft)/2.0 if(f(xmid,a,b).lt.eps.and.f(xmid,a,b).gt.(-1.0*eps)) goto 20 if(f(xmid,a,b).gt.0.0) then xleft=xmid else xright=xmid endif enddo endif print*,'Iterated solution not found in subroutine bisect' stop 20 return end function f(ustar,z,wspd) von=0.4 g=9.8 alpha=0.015 dum=z*g/(alpha*ustar**2) f=(von*wspd)**2/(log(dum))**2-ustar**2 return end FUNCTION is_numeric(string) Implicit none CHARACTER(len=*), INTENT(IN) :: string LOGICAL :: is_numeric REAL :: x INTEGER :: e READ(string,*,IOSTAT=e) x is_numeric = e == 0 END FUNCTION is_numeric Function to_upper (str) Result (string) ! ============================== ! Changes a string to upper case ! ============================== Implicit None Character(*), Intent(In) :: str Character(LEN(str)) :: string Integer :: ic, i Character(26), Parameter :: cap = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' Character(26), Parameter :: low = 'abcdefghijklmnopqrstuvwxyz' ! Capitalize each letter if it is lowecase string = str do i = 1, LEN_TRIM(str) ic = INDEX(low, str(i:i)) if (ic > 0) string(i:i) = cap(ic:ic) end do End Function to_upper