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, it real WindDir(2), WindSpeed(2), PeakWind(2), BP, RH, Rain15, - TempC(2), TempF(2), DwPt integer TowHeight, iTime_Now character*1 Stability(2) character*3 WindDirCh(2) character*4 Tow(2) character*5 cTime character*8 cDate character*16 EPZ(2) common WindSpeed, PeakWind, BP, RH, Rain15, TempC, TempF, WindDir, - DwPt, TowHeight, Stability, WindDirCh, cTime, cDate, EPZ Tow(1) = '1208' Tow(2) = '1209' c raw data file do it = 1, 2 if (it.eq. 1) then open(10,file='TOWK.txt') else open(10,file='TOWL.txt') end if 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) if (it .eq. 1) then BP = rData(9) RH = rData(8) DwPt = rData(6) Solar = rData(11) Rain15 = rdata(13) else BP = rData(6) Solar = rData(8) Rain15 = rdata(10) end if do i = 1, 2 TempC(i) = rData(i) TempF(i) = rData(i+2) if (it .eq. 1) then WindDir(i)= rData(i+14) WindSpeed(i) = rData(i+20) PeakWind(i) = rData(i+24) else WindDir(i)= rData(i+11) WindSpeed(i) = rData(i+17) PeakWind(i) = rData(i+21) end if if (WindDir(i) .ne. -99) then WindDirCH(i) = GetCH(WindDir(i)) EPZ(i) = get_planning_sectors(WindDir(i)) else WindDirCH(i) = '?' EPZ(i) = 'X' endif end do if (Solar .ne. -99 .and. TempC(2).ne. -99 .and. TempC(1) .ne. -99) - then deltat = 1.1 * (TempC(2) - TempC(1)) if (WindSpeed(1).ne. -99) Then Stability(1) =stabcalc(day_night, WindSpeed(1), solar, deltat) else Stability(1) = '?' endif if (WindSpeed(2).ne. -99) Then Stability(2) =stabcalc(day_night, WindSpeed(2), solar, deltat) else Stability(2) = '?' endif else do i = 1, 2 Stability(i) = '?' end do endif call GetPrecip(Rain15, iHr, iMin, it) if (it .eq. 1) then call writeTow(Tow(it), 10) call writeTow(Tow(it), 60) else call writeTow(Tow(it), 10) call writeTow(Tow(it), 30) end if end do 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(Rain15, Hr, Mn, i) implicit none real PrcTot, Rain15 integer Hr, Mn, i character*11 FN(2) common/PTotal/ PrcTot FN(1) = 'PrcTotK.txt' FN(2) = 'PrcTotL.txt' open(15,file=FN(i)) read(15,*) PrcTot close(15) open(15,file=FN(i)) if (Hr .eq. 0 .and. Mn .lt. 18) then PrcTot = 0. write(15,'(f5.2)') PrcTot endif if (Rain15 .gt. 0.0) then PrcTot = PrcTot + Rain15 write(15,'(f5.2)') PrcTot endif close(15) return end c------------------------------------------ subroutine writeTow(Tow, Hgt) implicit none integer Hgt, i real PrcTot common/PTotal/ PrcTot character*5 cRain15 character*100 fname character*180 outLIne real WindDir(2), WindSpeed(2), PeakWind(2), BP, RH, Rain15, - TempC(2), TempF(2), DwPt integer TowHeight, iTime_Now character*1 Stability(2) character*3 WindDirCh(2) character*4 Tow character*5 cTime character*8 cDate character*16 EPZ(2) common WindSpeed, PeakWind, BP, RH, Rain15, TempC, TempF, WindDir, - DwPt, TowHeight, Stability, WindDirCh, cTime, cDate, EPZ if (hgt .eq. 10) then i = 1 else if (hgt .eq. 30) then i = 2 else if (hgt .eq. 60) 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)),'F ', int(TempC(i)),'C DwPt ', - int(DwPt),'C','WFr', int(WindDir(i)), WindDirCh(i), - EPZ(i),'WSp ', WindSpeed(i), ' MPH PkWd ', PeakWind(i), - 'BP ', BP, ' in Stability ',Stability(i), - 'RH', int(RH), '% 24h Prc', PrcTot, ' in' 100 format(a8,1x,a4,2a4,i3, - a5,i3,a2,i3,a7,i3,a1, - a3,i4,1x,a3,1x,a12, - 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. 295.0 .or. wd .le. 5.0) then get_planning_sectors = myTrim(get_planning_sectors, 'A') end if if (wd .ge. 320.0 .or. wd .le. 55.0) then get_planning_sectors = myTrim(get_planning_sectors, 'B') end if if (wd .ge. 15.0 .And. wd .le. 100.0) then get_planning_sectors = myTrim(get_planning_sectors, 'C') end if if (wd .ge. 40.0 .And. wd .le. 145.0) then get_planning_sectors = myTrim(get_planning_sectors, 'D') end if if (wd .ge. 95.0 .And. wd .le. 180.0) then get_planning_sectors = myTrim(get_planning_sectors, 'E') end if if (wd .ge. 150.0 .And. wd .le. 255.0) then get_planning_sectors = myTrim(get_planning_sectors, 'F') end if if (wd .ge. 210.0 .And. wd .le. 285.0) then get_planning_sectors = myTrim(get_planning_sectors, 'G') end if if (wd .ge. 195.0 .And. wd .le. 260.0) then get_planning_sectors = myTrim(get_planning_sectors, 'H') end if if (wd .ge. 175.0 .And. wd .le. 245.0) then get_planning_sectors = myTrim(get_planning_sectors, 'J') end if if (wd .ge. 0.0 .And. wd .le. 360.0) then get_planning_sectors = myTrim(get_planning_sectors, 'K') end if if (wd .ge. 190.0 .And. wd .le. 250.0) then get_planning_sectors = myTrim(get_planning_sectors, 'L') end if if (wd .ge. 210.0 .And. wd .le. 255.0) then get_planning_sectors = myTrim(get_planning_sectors, 'M') end if if (wd .ge. 215.0 .And. wd .le. 260.0) then get_planning_sectors = myTrim(get_planning_sectors, 'N') end if if (wd .ge. 220.0 .And. wd .le. 280.0) then get_planning_sectors = myTrim(get_planning_sectors, 'P') end if if (wd .ge. 220.0 .And. wd .le. 290.0) then get_planning_sectors = myTrim(get_planning_sectors, 'Q') end if if (wd .ge. 245.0 .And. wd .le. 310.0) then get_planning_sectors = myTrim(get_planning_sectors, 'R') end if if (wd .ge. 265.0 .And. wd .le. 345.0) then get_planning_sectors = myTrim(get_planning_sectors, 'S') end if if (wd .ge. 230.0 .And. wd .le. 320.0) then get_planning_sectors = myTrim(get_planning_sectors, 'X') end if if (wd .ge. 215.0 .And. wd .le. 290.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