program WriteMetBoard implicit none character*1 stabcalc character*3 GetCh character*10 cDum character*16 get_planning_Sectors character*250 inLine real rData(100), deltat, Solar integer i, j, Jday, IY, IM, ID, day_night, Daylite, ihr, imin, - offset real WindDir(2,2), WindSpeed(2,2), PeakWind(2,2), BP, RH(2), - Rain15,TempC(3,2), TempF(2,2), Soil integer TowHeight, iTime_Now character*1 Stability(2) character*3 WindDirCh(2,2) character*4 Tow character*5 cTime character*8 cDate character*16 EPZ(2,2) common WindSpeed, PeakWind, BP, RH, Rain15, TempC, TempF, - WindDir, Soil, TowHeight, Stability, WindDirCh, Tow, - cTime, cDate, EPZ c read raw data file for tower C open(10,file='edas1.txt') read(10,'(a)')cDum read(10,'(a8,2x,a)')cDate, Inline close(10) i = 1 do while(inLine(i:i) .ne. ',') i = i+ 1 end do read(inline(1:i-1),'(a)') cTime Call readL(inline(i+1:), rData) read(cDate,'(i2,2(1x,i2))')IM, ID, IY i = 1 do while(cTime(i:i) .ne. ':') i = i+ 1 end do read(cTime(1:i-1),'(I2)') ihr read(cTime(i+1:),'(I2)') imin c add 17 mins to time imin = imin + 17 if (imin .gt. 59) then imin = imin - 60 ihr = ihr + 1 if (ihr .gt. 23) then ihr = ihr - 24 id = id + 1 call chkDMY(ID, IM, IY) end if endif c write(cDate,'(i2.2,2(a1, i2.2) )')im, '/', id, '/', iy write(cTime, '(2i2.2)') ihr, imin call JULIAN(IY,IM,ID,Jday) iTime_Now = ihr*100 + imin day_night = daylite(JDay, iTime_Now) BP = rData(12) Soil = rData(41) Solar = rData(14) c use all three tower heights to calculate stablity at two tower A/B heights do i = 1,3 TempC(i,1) = rData(i) WindSpeed(i,1) = rData(i+23) end do if (Solar .ne. -99 .and. TempC(3,1).ne. -99 .and. - TempC(1,1) .ne. -99) then deltat = 1.1 * (TempC(3,1) - TempC(1,1)) if (WindSpeed(1,1).ne. -99) Then Stability(1) = - stabcalc(day_night, WindSpeed(1,1), solar, deltat) else Stability(1) = '?' endif if (WindSpeed(2,1).ne. -99) Then Stability(2) = - stabcalc(day_night, WindSpeed(2,1), solar, deltat) else Stability(2) = '?' endif else do i = 1, 2 Stability(i) = '?' end do endif Rain15 = rdata(17) close(10) call GetPrecip c read tower A and B data open(10,file='TowAB_15.txt') read(10,'(a)')cDum read(10,'(a8,2x,a)')cDate, Inline close(10) i = 1 do while(inLine(i:i) .ne. ',') i = i+ 1 end do Call readL(inline(i+1:), rData) close(10) offset = 0 do i = 1, 2 TempC(i,1) = rData(1+offset) TempF(i,1) = rData(3+offset) WindDir(i,1) = rData(7+offset) WindSpeed(i,1) = rData(13+offset) PeakWind(i,1) = rData(17+offset) RH(i) = rData(6+offset) TempC(i,2) = rData(2+offset) TempF(i,2) = rData(4+offset) WindDir(i,2) = rData(8+offset) WindSpeed(i,2) = rData(14+offset) PeakWind(i,2) = rData(18+offset) do j = 1, 2 if (WindDir(i,j) .ne. -99) then WindDirCH(i,j) = GetCH(WindDir(i,j)) EPZ(i,j) = get_planning_sectors(WindDir(i,j)) else WindDirCH(i,j) = '?' EPZ(i,j) = 'X' endif end do offset = 20 end do Tow = 'TOWA' call writeTow(10) call writeTow(30) Tow = 'TOWB' call writeTow(15) call writeTow(30) stop end c------------------------------------------ subroutine readL(inline, rData) implicit none integer i, l, iStart, iEnd, ios real rData(100) character*250 inline character*10 cDum i = 1 do while (i .le. 250) if (inline(i:i) .ne. ' ') then l = i endif i = i+1 end do i = 1 iStart = 1 iEnd = 2 do while (iEnd .lt. l .and. i .le. 100) if (inline(iend:iend) .eq. ',') then iEnd = iEnd - 1 cDum = inline(iStart:iStart+2) if (cDum .eq. 'N/A' .or. - cDum .eq. 'INV') then rData(i) = -99. else cDum = inline(iStart:iEnd) read(cDum,'(f10.0)', iostat=ios) rdata(i) if (ios .ne. 0) then rdata(i) = -99. Endif End if i = i + 1 iStart = iEnd + 2 iEnd = iStart else iEnd = iEnd + 1 Endif end do return end c------------------------------------------ character*3 function GetCH(dir) implicit none Real dir, delta, lolim, hilim integer i, quad delta = 22.5 hilim = 11.25 lolim = -11.25 quad = -1 c do i = 1, 17 c if ((dir >= lolim) .and. (dir < hilim)) then c quad = i c endif c lolim = lolim + delta c hilim = hilim + delta c end do if (dir < 0) then GetCH = '??' else if ((dir >= 0) .and. (dir < 11.25)) then GetCH = 'N ' else if ((dir >= 11.25) .and. (dir < 33.75)) then GetCH = 'NNE' else if ((dir >= 33.75) .and. (dir < 56.25)) then GetCH = 'NE ' else if ((dir >= 56.75) .and. (dir < 78.75)) then GetCH = 'ENE' else if ((dir >= 78.75) .and. (dir < 101.25)) then GetCH = 'E ' else if ((dir >= 101.25) .and. (dir < 123.75)) then GetCH = 'ESE' else if ((dir >= 123.75) .and. (dir < 146.25)) then GetCH = 'SE ' else if ((dir >= 146.25) .and. (dir < 168.75)) then GetCH = 'SSE' else if ((dir >= 168.75) .and. (dir < 191.25)) then GetCH = 'S ' else if ((dir >= 191.25) .and. (dir < 213.75)) then GetCH = 'SSW' else if ((dir >= 213.75) .and. (dir < 236.25)) then GetCH = 'SW ' else if ((dir >= 236.25) .and. (dir < 258.75)) then GetCH = 'WSW' else if ((dir >= 258.75) .and. (dir < 281.25)) then GetCH = 'W ' else if ((dir >= 281.25) .and. (dir < 303.75)) then GetCH = 'WNW' else if ((dir >= 303.75) .and. (dir < 326.25)) then GetCH = 'NW ' else if ((dir >= 326.25) .and. (dir < 348.75)) then GetCH = 'NNW' else if ((dir >= 348.75) .and. (dir < 361.00)) then GetCH = 'N ' end if return end c------------------------------------------ subroutine getPrecip implicit none real PrcTot common/PTotal/ PrcTot open(15,file='PrcTotA.txt') read(15,*) PrcTot close(15) return end c------------------------------------------ subroutine writeTow(Hgt) implicit none integer Hgt, i, j real PrcTot common/PTotal/ PrcTot character*5 cRain15 character*100 fname character*180 outLIne real WindDir(2,2), WindSpeed(2,2), PeakWind(2,2), BP, RH(2), - Rain15,TempC(3,2), TempF(2,2), Soil integer TowHeight, iTime_Now character*1 Stability(2) character*3 WindDirCh(2,2) character*4 Tow character*5 cTime character*8 cDate character*16 EPZ(2,2) common WindSpeed, PeakWind, BP, RH, Rain15, TempC, TempF, WindDir, - Soil, TowHeight, Stability, WindDirCh, Tow, cTime, cDate, - EPZ if (hgt .eq. 10 .or. hgt .eq. 15) then j = 1 else if (hgt .eq. 30) then j = 2 endif if (tow .eq. 'TOWA') then i = 1 else if (tow .eq. 'TOWB') then i = 2 endif if (Rain15 .lt. 0.) then cRain15 = ' N/A ' else write(cRain15,'(f5.2)') Rain15 endif write(OutLine,100)cDate, ctime(1:4), ' hr ', tow, hgt, -'Temp ',int(TempF(i,j)),'F ', int(TempC(i,j)),'C Soil ', - int(Soil),'C','WFr', int(WindDir(i,j)), WindDirCh(i,j), - 'EPZ-', EPZ(i,j),'WSp ', WindSpeed(i,j), ' MPH PkWd ', - PeakWind(i,j), - 'BP ', BP, ' in Stability ',Stability(i), - 'RH', int(RH(i)), '% 24h Prc', PrcTot, ' in' 100 format(a8,1x,a4,2a4,i3, - a5,i3,a2,i3,a8,i2,a1, - a3,i4,1x,a3,1x,a4,a8, - a4,f4.1, a12, f4.1, - a3, f5.2, a15,a1, - a2, i4, a10, f5.2, a3) do i = 1, 178 if (outLine(i:i+1) .eq. '-9') then if (outLine(i:i+5) .eq. '-99.00') then outLine(i:i+5) = 'N/A ' else outLine(i:i+2) = 'N/A' endif endif end do write(fname,'(3a,i3.3,a)') 'ELED', Tow, '_', Hgt, '.txt' open (90,file=fname) write(90,'(a)') outLine close(90) return end c------------------------------------------ c get_planning_sectors() returns a string of emergency planning sectors for c ORNL, based on the direction the wind is FROM. The direction-to-sectors c conversion was derived by Kevin R Birdwell. c c The sectors string is returned in alphabetical order and is padded with c spaces to ten characters (maximum number of sectors for any direction). c character*16 function get_planning_sectors(wD) real wD character*16 temp, myTrim c unknown wind direction: 10 spaces if (wd .lt. 0.0) then get_planning_sectors = ' ' return endif C normalize direction to be 0-360 degrees do while (wd .gt. 360.0) wd = wd - 360.0 end do get_planning_sectors = '' C The following if statements build up the string one character at a time. C A and S are OR'ed because each straddles 0 degrees. if (wd .ge. 345.0 .or. wd .le. 105.0) then get_planning_sectors = myTrim(get_planning_sectors, 'A') end if if (wd .ge. 20.0 .And. wd .le. 90.0) then get_planning_sectors = myTrim(get_planning_sectors, 'B') end if if (wd .ge. 45.0 .And. wd .le. 115.0) then get_planning_sectors = myTrim(get_planning_sectors, 'C') end if if (wd .ge. 60.0 .And. wd .le. 125.0) then get_planning_sectors = myTrim(get_planning_sectors, 'D') end if if (wd .ge. 85.0 .And. wd .le. 150.0) then get_planning_sectors = myTrim(get_planning_sectors, 'E') end if if (wd .ge. 105.0 .And. wd .le. 175.0) then get_planning_sectors = myTrim(get_planning_sectors, 'F') end if if (wd .ge. 110.0 .And. wd .le. 195.0) then get_planning_sectors = myTrim(get_planning_sectors, 'G') end if if (wd .ge. 135.0 .And. wd .le. 205.0) then get_planning_sectors = myTrim(get_planning_sectors, 'H') end if if (wd .ge. 135.0 .And. wd .le. 195.0) then get_planning_sectors = myTrim(get_planning_sectors, 'J') end if if (wd .ge. 45.0 .And. wd .le. 155.0) then get_planning_sectors = myTrim(get_planning_sectors, 'K') end if if (wd .ge. 155.0 .And. wd .le. 215.0) then get_planning_sectors = myTrim(get_planning_sectors, 'L') end if if (wd .ge. 165.0 .And. wd .le. 220.0) then get_planning_sectors = myTrim(get_planning_sectors, 'M') end if if (wd .ge. 185.0 .And. wd .le. 240.0) then get_planning_sectors = myTrim(get_planning_sectors, 'N') end if if (wd .ge. 200.0 .And. wd .le. 240.0) then get_planning_sectors = myTrim(get_planning_sectors, 'P') end if if (wd .ge. 205.0 .And. wd .le. 285.0) then get_planning_sectors = myTrim(get_planning_sectors, 'Q') end if if (wd .ge. 235.0 .And. wd .le. 315.0) then get_planning_sectors = myTrim(get_planning_sectors, 'R') end if if (wd .ge. 285.0 .or. wd .le. 45.0) then get_planning_sectors = myTrim(get_planning_sectors, 'S') end if if (wd .ge. 0.0 .And. wd .le. 360.0) then get_planning_sectors = myTrim(get_planning_sectors, 'X') end if if (wd .ge. 155.0 .And. wd .le. 250.0) Then get_planning_sectors = myTrim(get_planning_sectors, 'Y') end if return end c------------------------------------------ character*16 function myTrim(sStr, cChr) character*16 sStr character*1 cChr integer i i = 1 do while (i .le. 16) if (sStr(i:i) .eq. ' ') then sStr(i:i) = cChr i = 17 end if i = i+1 end do myTrim = sStr return end c------------------------------------------ integer function daylite(JDay, Time_now) implicit none integer Jday real Latitude, Longitude, Radians real*8 AngleODay, DoubleDayAngle, SinDay, CosDay, SinDoubleDay, - CosDoubleDay, Sigma, DSin, DCos, MeridianTime, HCos, - SolarHourAngle, TimeOSun real Sunrise, Sunset integer HourOSun, MinOSun, day_start, day_end, TimeZone, DSTStrt, - DSTEnd, Time_Now Latitude = 35.56 Longitude = 84.23 TimeZone = 5 DSTStrt = 110 DSTEnd = 293 Radians = 57.29578 if (JDay .gt. DSTStrt .And. JDay .lt. DSTEnd) then TimeZone = TimeZone - 1 end if C Determine the angular (radian) fraction of year for date C .0172028 = 360 (degrees in circle) / 365.242 (days in year) * 57.29578 AngleODay = (JDay - 1) * .0172028 DoubleDayAngle = AngleODay * 2 SinDay = sin(AngleODay) CosDay = cos(AngleODay) SinDoubleDay = sin(DoubleDayAngle) CosDoubleDay = cos(DoubleDayAngle) C Allow for ellipticity of earth's orbit Sigma = 279.9384 + (AngleODay * Radians) + 1.914827 * SinDay - - .079525 * CosDay + .019938 * SinDoubleDay - .00162 * - CosDoubleDay C sine of solar declination DSin = .39785 * sin(Sigma / Radians) DCos = sqrt(1 - DSin * DSin) C Time (hrs) of meridian passage MeridianTime = 12 + .12357 * SinDay - .004289 * CosDay + .153809 * - SinDoubleDay + .060783 * CosDoubleDay HCos = (-sin(Latitude / Radians) * DSin) / - (cos(Latitude / Radians) * DCos) C Solar hour angle of sunrise-sunset SolarHourAngle = (atan2(sqrt(1 - HCos * HCos), HCos) / 15) * - Radians TimeOSun = MeridianTime - SolarHourAngle + - (Longitude / 15 - TimeZone) HourOSun = int(TimeOSun) MinOSun = (TimeOSun - HourOSun) * 60 Sunrise = HourOSun * 100 + MinOSun TimeOSun = MeridianTime + SolarHourAngle + - (Longitude / 15 - TimeZone) HourOSun = int(TimeOSun ) MinOSun = (TimeOSun - HourOSun) * 60 Sunset = HourOSun * 100 + MinOSun day_start = Sunrise + 100 day_end = Sunset - 100 if (time_now .gt. day_start .And.time_now .lt. day_end) then daylite = 1 else daylite = 0 endif return end c------------------------------------------ C Degrees C or degrees F will work for deltat in stabcalc() because it C is compared with zero: If deltaF is 0 then so is deltaC, and vice versa. Character*1 function stabcalc(day_night, wspeed, solar, deltat) implicit none integer day_night real wspeed, solar, deltat if (day_night .eq. 1) then if (solar .lt. 0.25) then stabcalc = 'D' else if (solar .le. 0.96) then if (wspeed .ge. 11.2) then stabcalc = 'D' else if (wspeed .ge. 4.5) then stabcalc = 'C' else stabcalc = 'B' endif endif else if (solar .le. 1.32) then if (wspeed .gt. 13.4) then stabcalc = 'D' else if (wspeed .ge. 11.2) then stabcalc = 'C' else if (wspeed .ge. 4.5) then stabcalc = 'B' else stabcalc = 'A' endif else if (solar .gt. 1.32) then if (wspeed .ge. 11.2) then stabcalc = 'C' else if (wspeed .ge. 6.8) then stabcalc = 'B' else stabcalc = 'A' endif endif else if (deltat .lt. 0) then if (wspeed .ge. 4.5) then stabcalc = 'D' else stabcalc = 'E' endif elseif (deltat .lt. 4) then if (wspeed .gt. 5.5) then stabcalc = 'D' else if (wspeed .ge. 4.5) then stabcalc = 'E' else stabcalc = 'F' endif elseif (deltat .ge. 4) then if (wspeed .ge. 3.4) then stabcalc = 'F' else stabcalc = 'G' endif endif endif end c------------------------------------------ SUBROUTINE JULIAN(IY,IM,ID,JDAY) C THIS SUBROUTINE DETERMINES THE JULIAN DAY NUMBER GO TO(10,20,30,40,50,60,70,80,90,100,110,120)IM 10 JDAY = ID RETURN 20 JDAY = ID + 31 GO TO 130 30 JDAY = ID + 59 GO TO 130 40 JDAY = ID + 90 GO TO 130 50 JDAY = ID + 120 GO TO 130 60 JDAY = ID + 151 GO TO 130 70 JDAY = ID + 181 GO TO 130 80 JDAY = ID + 212 GO TO 130 90 JDAY = ID + 243 GO TO 130 100 JDAY = ID + 273 GO TO 130 110 JDAY = ID + 304 GO TO 130 120 JDAY = ID + 334 C CHECK FOR LEAP YEAR 130 CHECK = (FLOAT(IY)/4.0) - (IY/4) IF (CHECK.EQ.0.AND.IM.GE.3)JDAY = JDAY + 1 RETURN END c------------------------------------------ SUBROUTINE chkDMY(id, im, iy) GO TO(10,20,30,40,50,60,70,80,90,100,110,120)IM 10 if (id.gt. 31) then im = im + 1 RETURN end if C CHECK FOR LEAP YEAR 20 CHECK = (FLOAT(IY)/4.0) - (IY/4) IF (CHECK.EQ.0.)then if (id.gt. 29) then im = im + 1 endif else if (id.gt. 28) then im = im + 1 end if end if return 30 if (id.gt. 31) then im = im + 1 RETURN end if 40 if (id.gt. 30) then im = im + 1 RETURN end if 50 if (id.gt. 31) then im = im + 1 RETURN end if 60 if (id.gt. 30) then im = im + 1 RETURN end if 70 if (id.gt. 31) then im = im + 1 RETURN end if 80 if (id.gt. 31) then im = im + 1 RETURN end if 90 if (id.gt. 30) then im = im + 1 RETURN end if 100 if (id.gt. 31) then im = im + 1 RETURN end if 110 if (id.gt. 30) then im = im + 1 RETURN end if 120 if (id.gt. 31) then im = 1 iy = iy + 1 RETURN end if END