Qual2K: a modeling Framework for Simulating River and Stream Water Quality (Version 11)
APPENDIX B: SOLAR POSITION, SUNRISE, AND SUNSET CALCULATIONS
Yüklə 1,65 Mb.
|
- Bu səhifədəki naviqasiya:
- APPENDIX C: SEDIMENT-WATER HEAT EXCHANGE
APPENDIX B: SOLAR POSITION, SUNRISE, AND SUNSET CALCULATIONSThe sunrise/sunset and solar position functions are a VBA translation of NOAA's sunrise/sunset calculator and NOAA's solar position calculator at the following web pages:
The calculations in the NOAA Sunrise/Sunset and Solar Position Calculators are based on equations from Astronomical Algorithms, by Jean Meeus. The sunrise and sunset results have been verified by NOAA to be accurate to within a minute for locations between +/- 72° latitude, and within 10 minutes outside of those latitudes. Five main functions are included for use from Excel worksheets or VBA programs:
A subroutine is also provided that calculates solar azimuth (az), solar elevation (el):
The sign convention for the main functions and subroutine is:
The Excel/VBA functions and subroutines for solar position, sunrise, and sunset times are as follows: Option Explicit
'// Convert radian angle to degrees radToDeg = (180# * angleRad / Application.WorksheetFunction.Pi()) End Function Function degToRad(angleDeg) '// Convert degree angle to radians degToRad = (Application.WorksheetFunction.Pi() * angleDeg / 180#) End Function Function calcJD(year, month, day) '***********************************************************************/ '* Name: calcJD '* Type: Function '* Purpose: Julian day from calendar day '* Arguments: '* year : 4 digit year '* month: January = 1 '* day : 1 - 31 '* Return value: '* The Julian day corresponding to the date '* Note: '* Number is returned for start of day. Fractional days should be '* added later. '***********************************************************************/ Dim A As Double, B As Double, jd As Double If (month <= 2) Then year = year - 1 month = month + 12 End If
A = Application.WorksheetFunction.Floor(year / 100, 1) B = 2 - A + Application.WorksheetFunction.Floor(A / 4, 1) jd = Application.WorksheetFunction.Floor(365.25 * (year + 4716), 1) + _ Application.WorksheetFunction.Floor(30.6001 * (month + 1), 1) + day + B - 1524.5 calcJD = jd 'gp put the year and month back where they belong If month = 13 Then month = 1 year = year + 1 End If
If month = 14 Then month = 2 year = year + 1 End If
End Function Function calcTimeJulianCent(jd) '***********************************************************************/ '* Name: calcTimeJulianCent '* Type: Function '* Purpose: convert Julian Day to centuries since J2000.0. '* Arguments: '* jd : the Julian Day to convert '* Return value: '* the T value corresponding to the Julian Day '***********************************************************************/
calcTimeJulianCent = t End Function Function calcJDFromJulianCent(t) '***********************************************************************/ '* Name: calcJDFromJulianCent '* Type: Function '* Purpose: convert centuries since J2000.0 to Julian Day. '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* the Julian Day corresponding to the t value '***********************************************************************/
calcJDFromJulianCent = jd End Function Function calcGeomMeanLongSun(t) '***********************************************************************/ '* Name: calGeomMeanLongSun '* Type: Function '* Purpose: calculate the Geometric Mean Longitude of the Sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* the Geometric Mean Longitude of the Sun in degrees '***********************************************************************/
Do If (l0 <= 360) And (l0 >= 0) Then Exit Do If l0 > 360 Then l0 = l0 - 360 If l0 < 0 Then l0 = l0 + 360 Loop calcGeomMeanLongSun = l0 End Function Function calcGeomMeanAnomalySun(t) '***********************************************************************/ '* Name: calGeomAnomalySun '* Type: Function '* Purpose: calculate the Geometric Mean Anomaly of the Sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* the Geometric Mean Anomaly of the Sun in degrees '***********************************************************************/ Dim m As Double m = 357.52911 + t * (35999.05029 - 0.0001537 * t) calcGeomMeanAnomalySun = m End Function
'* Name: calcEccentricityEarthOrbit '* Type: Function '* Purpose: calculate the eccentricity of earth's orbit '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* the unitless eccentricity '***********************************************************************/
calcEccentricityEarthOrbit = e End Function
'* Name: calcSunEqOfCenter '* Type: Function '* Purpose: calculate the equation of center for the sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* in degrees '***********************************************************************/
Dim c As Double m = calcGeomMeanAnomalySun(t) mrad = degToRad(m) sinm = Sin(mrad) sin2m = Sin(mrad + mrad) sin3m = Sin(mrad + mrad + mrad) c = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) _ + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289 calcSunEqOfCenter = c End Function Function calcSunTrueLong(t) '***********************************************************************/ '* Name: calcSunTrueLong '* Type: Function '* Purpose: calculate the true longitude of the sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* sun's true longitude in degrees '***********************************************************************/
c = calcSunEqOfCenter(t) O = l0 + c calcSunTrueLong = O End Function
'* Name: calcSunTrueAnomaly (not used by sunrise, solarnoon, sunset) '* Type: Function '* Purpose: calculate the true anamoly of the sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* sun's true anamoly in degrees '***********************************************************************/
c = calcSunEqOfCenter(t) v = m + c calcSunTrueAnomaly = v End Function
'* Name: calcSunRadVector (not used by sunrise, solarnoon, sunset) '* Type: Function '* Purpose: calculate the distance to the sun in AU '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* sun radius vector in AUs '***********************************************************************/
e = calcEccentricityEarthOrbit(t) r = (1.000001018 * (1 - e * e)) / (1 + e * Cos(degToRad(v))) calcSunRadVector = r End Function
'* Name: calcSunApparentLong (not used by sunrise, solarnoon, sunset) '* Type: Function '* Purpose: calculate the apparent longitude of the sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* sun's apparent longitude in degrees '***********************************************************************/
lambda = O - 0.00569 - 0.00478 * Sin(degToRad(omega)) calcSunApparentLong = lambda
'* Name: calcMeanObliquityOfEcliptic '* Type: Function '* Purpose: calculate the mean obliquity of the ecliptic '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* mean obliquity in degrees '***********************************************************************/
e0 = 23# + (26# + (seconds / 60#)) / 60# calcMeanObliquityOfEcliptic = e0 End Function Function calcObliquityCorrection(t) '***********************************************************************/ '* Name: calcObliquityCorrection '* Type: Function '* Purpose: calculate the corrected obliquity of the ecliptic '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* corrected obliquity in degrees '***********************************************************************/
e = e0 + 0.00256 * Cos(degToRad(omega)) calcObliquityCorrection = e End Function Function calcSunRtAscension(t) '***********************************************************************/ '* Name: calcSunRtAscension (not used by sunrise, solarnoon, sunset) '* Type: Function '* Purpose: calculate the right ascension of the sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* sun's right ascension in degrees '***********************************************************************/
Dim alpha As Double e = calcObliquityCorrection(t) lambda = calcSunApparentLong(t) tananum = (Cos(degToRad(e)) * Sin(degToRad(lambda))) tanadenom = (Cos(degToRad(lambda))) 'original NOAA code using javascript Math.Atan2(y,x) convention: ' var alpha = radToDeg(Math.atan2(tananum, tanadenom)); ' alpha = radToDeg(Application.WorksheetFunction.Atan2(tananum, tanadenom))
alpha = radToDeg(Application.WorksheetFunction.Atan2(tanadenom, tananum)) calcSunRtAscension = alpha
'* Type: Function '* Purpose: calculate the declination of the sun '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* sun's declination in degrees '***********************************************************************/ Dim e As Double, lambda As Double, sint As Double, theta As Double e = calcObliquityCorrection(t) lambda = calcSunApparentLong(t) sint = Sin(degToRad(e)) * Sin(degToRad(lambda)) theta = radToDeg(Application.WorksheetFunction.Asin(sint)) calcSunDeclination = theta
'* Name: calcEquationOfTime '* Type: Function '* Purpose: calculate the difference between true solar time and mean '* solar time '* Arguments: '* t : number of Julian centuries since J2000.0 '* Return value: '* equation of time in minutes of time '***********************************************************************/ Dim epsilon As Double, l0 As Double, e As Double, m As Double Dim y As Double, sin2l0 As Double, sinm As Double Dim cos2l0 As Double, sin4l0 As Double, sin2m As Double, Etime As Double
l0 = calcGeomMeanLongSun(t) e = calcEccentricityEarthOrbit(t) m = calcGeomMeanAnomalySun(t) y = Tan(degToRad(epsilon) / 2#) y = y ^ 2 sin2l0 = Sin(2# * degToRad(l0)) sinm = Sin(degToRad(m)) cos2l0 = Cos(2# * degToRad(l0)) sin4l0 = Sin(4# * degToRad(l0)) sin2m = Sin(2# * degToRad(m))
- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m calcEquationOfTime = radToDeg(Etime) * 4# End Function Function calcHourAngleSunrise(lat, SolarDec)
'* Name: calcHourAngleSunrise '* Type: Function '* Purpose: calculate the hour angle of the sun at sunrise for the '* latitude '* Arguments: '* lat : latitude of observer in degrees '* solarDec : declination angle of sun in degrees '* Return value: '* hour angle of sunrise in radians '***********************************************************************/
sdRad = degToRad(SolarDec) HAarg = (Cos(degToRad(90.833)) / (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad)) HA = (Application.WorksheetFunction.Acos(Cos(degToRad(90.833)) _ / (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad))) calcHourAngleSunrise = HA End Function Function calcHourAngleSunset(lat, SolarDec) '***********************************************************************/ '* Name: calcHourAngleSunset '* Type: Function '* Purpose: calculate the hour angle of the sun at sunset for the '* latitude '* Arguments: '* lat : latitude of observer in degrees '* solarDec : declination angle of sun in degrees '* Return value: '* hour angle of sunset in radians '***********************************************************************/
sdRad = degToRad(SolarDec) HAarg = (Cos(degToRad(90.833)) / (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad)) HA = (Application.WorksheetFunction.Acos(Cos(degToRad(90.833)) _ / (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad))) calcHourAngleSunset = -HA End Function Function calcSunriseUTC(jd, Latitude, longitude) '***********************************************************************/ '* Name: calcSunriseUTC '* Type: Function '* Purpose: calculate the Universal Coordinated Time (UTC) of sunrise '* for the given day at the given location on earth '* Arguments: '* JD : julian day '* latitude : latitude of observer in degrees '* longitude : longitude of observer in degrees '* Return value: '* time in minutes from zero Z '***********************************************************************/ Dim t As Double, eqtime As Double, SolarDec As Double, hourangle As Double Dim delta As Double, timeDiff As Double, timeUTC As Double Dim newt As Double
SolarDec = calcSunDeclination(t) hourangle = calcHourAngleSunrise(Latitude, SolarDec)
timeDiff = 4 * delta ' in minutes of time timeUTC = 720 + timeDiff - eqtime ' in minutes
eqtime = calcEquationOfTime(newt) SolarDec = calcSunDeclination(newt) hourangle = calcHourAngleSunrise(Latitude, SolarDec) delta = longitude - radToDeg(hourangle) timeDiff = 4 * delta timeUTC = 720 + timeDiff - eqtime ' in minutes calcSunriseUTC = timeUTC End Function Function calcSolNoonUTC(t, longitude) '***********************************************************************/ '* Name: calcSolNoonUTC '* Type: Function '* Purpose: calculate the Universal Coordinated Time (UTC) of solar '* noon for the given day at the given location on earth '* Arguments: '* t : number of Julian centuries since J2000.0 '* longitude : longitude of observer in degrees '* Return value: '* time in minutes from zero Z '***********************************************************************/ Dim newt As Double, eqtime As Double, solarNoonDec As Double, solNoonUTC As Double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 + longitude / 360#)
solarNoonDec = calcSunDeclination(newt) solNoonUTC = 720 + (longitude * 4) - eqtime calcSolNoonUTC = solNoonUTC End Function Function calcSunsetUTC(jd, Latitude, longitude) '***********************************************************************/ '* Name: calcSunsetUTC '* Type: Function '* Purpose: calculate the Universal Coordinated Time (UTC) of sunset '* for the given day at the given location on earth '* Arguments: '* JD : julian day '* latitude : latitude of observer in degrees '* longitude : longitude of observer in degrees '* Return value: '* time in minutes from zero Z '***********************************************************************/ Dim t As Double, eqtime As Double, SolarDec As Double, hourangle As Double Dim delta As Double, timeDiff As Double, timeUTC As Double Dim newt As Double t = calcTimeJulianCent(jd) ' // First calculates sunrise and approx length of day eqtime = calcEquationOfTime(t) SolarDec = calcSunDeclination(t) hourangle = calcHourAngleSunset(Latitude, SolarDec)
timeDiff = 4 * delta timeUTC = 720 + timeDiff - eqtime
eqtime = calcEquationOfTime(newt) SolarDec = calcSunDeclination(newt) hourangle = calcHourAngleSunset(Latitude, SolarDec) delta = longitude - radToDeg(hourangle) timeDiff = 4 * delta timeUTC = 720 + timeDiff - eqtime ' // in minutes calcSunsetUTC = timeUTC End Function Function sunrise(lat, lon, year, month, day, timezone, dlstime) '***********************************************************************/ '* Name: sunrise '* Type: Main Function called by spreadsheet '* Purpose: calculate time of sunrise for the entered date '* and location. '* For latitudes greater than 72 degrees N and S, calculations are '* accurate to within 10 minutes. For latitudes less than +/- 72° '* accuracy is approximately one minute. '* Arguments: ' latitude = latitude (decimal degrees) ' longitude = longitude (decimal degrees) ' NOTE: longitude is negative for western hemisphere for input cells ' in the spreadsheet for calls to the functions named ' sunrise, solarnoon, and sunset. Those functions convert the ' longitude to positive for the western hemisphere for calls to ' other functions using the original sign convention ' from the NOAA javascript code. ' year = year ' month = month ' day = day ' timezone = time zone hours relative to GMT/UTC (hours) ' dlstime = daylight savings time (0 = no, 1 = yes) (hours) '* Return value: '* sunrise time in local time (days) '***********************************************************************/ Dim longitude As Double, Latitude As Double, jd As Double Dim riseTimeGMT As Double, riseTimeLST As Double ' change sign convention for longitude from negative to positive in western hemisphere longitude = lon * -1 Latitude = lat If (Latitude > 89.8) Then Latitude = 89.8 If (Latitude < -89.8) Then Latitude = -89.8 jd = calcJD(year, month, day) ' // Calculate sunrise for this date riseTimeGMT = calcSunriseUTC(jd, Latitude, longitude) ' // adjust for time zone and daylight savings time in minutes riseTimeLST = riseTimeGMT + (60 * timezone) + (dlstime * 60) ' // convert to days sunrise = riseTimeLST / 1440 End Function Function solarnoon(lat, lon, year, month, day, timezone, dlstime) '***********************************************************************/ '* Name: solarnoon '* Type: Main Function called by spreadsheet '* Purpose: calculate the Universal Coordinated Time (UTC) of solar '* noon for the given day at the given location on earth '* Arguments: ' year ' month
' day '* longitude : longitude of observer in degrees ' NOTE: longitude is negative for western hemisphere for input cells ' in the spreadsheet for calls to the functions named ' sunrise, solarnoon, and sunset. Those functions convert the ' longitude to positive for the western hemisphere for calls to ' other functions using the original sign convention ' from the NOAA javascript code. '* Return value: '* time of solar noon in local time days '***********************************************************************/ Dim longitude As Double, Latitude As Double, jd As Double Dim t As Double, newt As Double, eqtime As Double Dim solarNoonDec As Double, solNoonUTC As Double ' change sign convention for longitude from negative to positive in western hemisphere longitude = lon * -1 Latitude = lat If (Latitude > 89.8) Then Latitude = 89.8 If (Latitude < -89.8) Then Latitude = -89.8 jd = calcJD(year, month, day) t = calcTimeJulianCent(jd) newt = calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 + longitude / 360#) eqtime = calcEquationOfTime(newt) solarNoonDec = calcSunDeclination(newt) solNoonUTC = 720 + (longitude * 4) - eqtime ' // adjust for time zone and daylight savings time in minutes solarnoon = solNoonUTC + (60 * timezone) + (dlstime * 60)
solarnoon = solarnoon / 1440 End Function Function sunset(lat, lon, year, month, day, timezone, dlstime) '***********************************************************************/ '* Name: sunset '* Type: Main Function called by spreadsheet '* Purpose: calculate time of sunrise and sunset for the entered date '* and location. '* For latitudes greater than 72 degrees N and S, calculations are '* accurate to within 10 minutes. For latitudes less than +/- 72° '* accuracy is approximately one minute. '* Arguments: ' latitude = latitude (decimal degrees) ' longitude = longitude (decimal degrees) ' NOTE: longitude is negative for western hemisphere for input cells ' in the spreadsheet for calls to the functions named ' sunrise, solarnoon, and sunset. Those functions convert the ' longitude to positive for the western hemisphere for calls to ' other functions using the original sign convention ' from the NOAA javascript code. ' year = year ' month = month ' day = day ' timezone = time zone hours relative to GMT/UTC (hours) ' dlstime = daylight savings time (0 = no, 1 = yes) (hours) '* Return value: '* sunset time in local time (days) '***********************************************************************/ Dim longitude As Double, Latitude As Double, jd As Double Dim setTimeGMT As Double, setTimeLST As Double ' change sign convention for longitude from negative to positive in western hemisphere longitude = lon * -1 Latitude = lat If (Latitude > 89.8) Then Latitude = 89.8 If (Latitude < -89.8) Then Latitude = -89.8 jd = calcJD(year, month, day) ' // Calculate sunset for this date setTimeGMT = calcSunsetUTC(jd, Latitude, longitude) ' // adjust for time zone and daylight savings time in minutes setTimeLST = setTimeGMT + (60 * timezone) + (dlstime * 60) ' // convert to days sunset = setTimeLST / 1440 End Function Function solarazimuth(lat, lon, year, month, day, _ hours, minutes, seconds, timezone, dlstime) '***********************************************************************/ '* Name: solarazimuth '* Type: Main Function '* Purpose: calculate solar azimuth (deg from north) for the entered '* date, time and location. Returns -999999 if darker than twilight '* '* Arguments: '* latitude, longitude, year, month, day, hour, minute, second, '* timezone, daylightsavingstime '* Return value: '* solar azimuth in degrees from north '* '* Note: solarelevation and solarazimuth functions are identical '* and could be converted to a VBA subroutine that would return '* both values. '* '***********************************************************************/ Dim longitude As Double, Latitude As Double Dim Zone As Double, daySavings As Double Dim hh As Double, mm As Double, ss As Double, timenow As Double Dim jd As Double, t As Double, r As Double Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double Dim SolarDec As Double, earthRadVec As Double, solarTimeFix As Double Dim trueSolarTime As Double, hourangle As Double, harad As Double Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double Dim azimuth As Double, exoatmElevation As Double Dim step1 As Double, step2 As Double, step3 As Double Dim refractionCorrection As Double, Te As Double, solarzen As Double
Latitude = lat If (Latitude > 89.8) Then Latitude = 89.8 If (Latitude < -89.8) Then Latitude = -89.8 Zone = timezone * -1 daySavings = dlstime * 60 hh = hours - (daySavings / 60) mm = minutes ss = seconds
timenow = hh + mm / 60 + ss / 3600 + Zone jd = calcJD(year, month, day) t = calcTimeJulianCent(jd + timenow / 24#) r = calcSunRadVector(t) alpha = calcSunRtAscension(t) theta = calcSunDeclination(t) Etime = calcEquationOfTime(t) eqtime = Etime SolarDec = theta '// in degrees earthRadVec = r
trueSolarTime = hh * 60# + mm + ss / 60# + solarTimeFix '// in minutes
trueSolarTime = trueSolarTime - 1440 Loop hourangle = trueSolarTime / 4# - 180# '// Thanks to Louis Schwarzmayr for the next line: If (hourangle < -180) Then hourangle = hourangle + 360# harad = degToRad(hourangle) csz = Sin(degToRad(Latitude)) * _ Sin(degToRad(SolarDec)) + _ Cos(degToRad(Latitude)) * _ Cos(degToRad(SolarDec)) * Cos(harad) If (csz > 1#) Then csz = 1#
ElseIf (csz < -1#) Then csz = -1# End If zenith = radToDeg(Application.WorksheetFunction.Acos(csz)) azDenom = (Cos(degToRad(Latitude)) * Sin(degToRad(zenith))) If (Abs(azDenom) > 0.001) Then azRad = ((Sin(degToRad(Latitude)) * _ Cos(degToRad(zenith))) - _ Sin(degToRad(SolarDec))) / azDenom If (Abs(azRad) > 1#) Then If (azRad < 0) Then azRad = -1# Else azRad = 1# End If End If
If (hourangle > 0#) Then azimuth = -azimuth End If Else
If (Latitude > 0#) Then azimuth = 180# Else azimuth = 0# End If End If
If (azimuth < 0#) Then azimuth = azimuth + 360# End If exoatmElevation = 90# - zenith If (exoatmElevation > 85#) Then refractionCorrection = 0# Else Te = Tan(degToRad(exoatmElevation)) If (exoatmElevation > 5#) Then refractionCorrection = 58.1 / Te - 0.07 / (Te * Te * Te) + _ 0.000086 / (Te * Te * Te * Te * Te) ElseIf (exoatmElevation > -0.575) Then step1 = (-12.79 + exoatmElevation * 0.711) step2 = (103.4 + exoatmElevation * (step1)) step3 = (-518.2 + exoatmElevation * (step2)) refractionCorrection = 1735# + exoatmElevation * (step3) Else refractionCorrection = -20.774 / Te End If refractionCorrection = refractionCorrection / 3600# End If solarzen = zenith - refractionCorrection solarazimuth = azimuth End Function Function solarelevation(lat, lon, year, month, day, _ hours, minutes, seconds, timezone, dlstime) '***********************************************************************/ '* Name: solarazimuth '* Type: Main Function '* Purpose: calculate solar azimuth (deg from north) for the entered '* date, time and location. Returns -999999 if darker than twilight '* '* Arguments: '* latitude, longitude, year, month, day, hour, minute, second, '* timezone, daylightsavingstime '* Return value: '* solar azimuth in degrees from north '* '* Note: solarelevation and solarazimuth functions are identical '* and could converted to a VBA subroutine that would return '* both values. '* '***********************************************************************/ Dim longitude As Double, Latitude As Double Dim Zone As Double, daySavings As Double Dim hh As Double, mm As Double, ss As Double, timenow As Double Dim jd As Double, t As Double, r As Double Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double Dim SolarDec As Double, earthRadVec As Double, solarTimeFix As Double Dim trueSolarTime As Double, hourangle As Double, harad As Double Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double Dim azimuth As Double, exoatmElevation As Double Dim step1 As Double, step2 As Double, step3 As Double Dim refractionCorrection As Double, Te As Double, solarzen As Double
Latitude = lat If (Latitude > 89.8) Then Latitude = 89.8 If (Latitude < -89.8) Then Latitude = -89.8 Zone = timezone * -1 daySavings = dlstime * 60 hh = hours - (daySavings / 60) mm = minutes ss = seconds
timenow = hh + mm / 60 + ss / 3600 + Zone jd = calcJD(year, month, day) t = calcTimeJulianCent(jd + timenow / 24#) r = calcSunRadVector(t) alpha = calcSunRtAscension(t) theta = calcSunDeclination(t) Etime = calcEquationOfTime(t) eqtime = Etime SolarDec = theta '// in degrees earthRadVec = r
trueSolarTime = hh * 60# + mm + ss / 60# + solarTimeFix '// in minutes
trueSolarTime = trueSolarTime - 1440 Loop hourangle = trueSolarTime / 4# - 180# '// Thanks to Louis Schwarzmayr for the next line: If (hourangle < -180) Then hourangle = hourangle + 360# harad = degToRad(hourangle) csz = Sin(degToRad(Latitude)) * _ Sin(degToRad(SolarDec)) + _ Cos(degToRad(Latitude)) * _ Cos(degToRad(SolarDec)) * Cos(harad) If (csz > 1#) Then csz = 1#
ElseIf (csz < -1#) Then csz = -1# End If zenith = radToDeg(Application.WorksheetFunction.Acos(csz)) azDenom = (Cos(degToRad(Latitude)) * Sin(degToRad(zenith))) If (Abs(azDenom) > 0.001) Then azRad = ((Sin(degToRad(Latitude)) * _ Cos(degToRad(zenith))) - _ Sin(degToRad(SolarDec))) / azDenom If (Abs(azRad) > 1#) Then If (azRad < 0) Then azRad = -1# Else azRad = 1# End If End If
If (hourangle > 0#) Then azimuth = -azimuth End If Else
If (Latitude > 0#) Then azimuth = 180# Else azimuth = 0# End If End If
If (azimuth < 0#) Then azimuth = azimuth + 360# End If exoatmElevation = 90# - zenith If (exoatmElevation > 85#) Then refractionCorrection = 0# Else Te = Tan(degToRad(exoatmElevation)) If (exoatmElevation > 5#) Then refractionCorrection = 58.1 / Te - 0.07 / (Te * Te * Te) + _ 0.000086 / (Te * Te * Te * Te * Te) ElseIf (exoatmElevation > -0.575) Then step1 = (-12.79 + exoatmElevation * 0.711) step2 = (103.4 + exoatmElevation * (step1)) step3 = (-518.2 + exoatmElevation * (step2)) refractionCorrection = 1735# + exoatmElevation * (step3) Else refractionCorrection = -20.774 / Te End If refractionCorrection = refractionCorrection / 3600# End If solarzen = zenith - refractionCorrection solarelevation = 90# - solarzen End Function Sub solarposition(lat, lon, year, month, day, _ hours, minutes, seconds, timezone, dlstime, _ solarazimuth, solarelevation, earthRadVec)
'* Name: solarposition '* Type: Subroutine '* Purpose: calculate solar azimuth (deg from north) '* and elevation (deg from horizon) for the entered '* date, time and location. '* '* Arguments: '* latitude, longitude, year, month, day, hour, minute, second, '* timezone, daylightsavingstime '* Return value: '* solar azimuth in degrees from north '* solar elevation in degrees from horizon '* earth radius vector (distance to the sun in AU) '* '* Note: solarelevation and solarazimuth functions are identical '* and could converted to a VBA subroutine that would return '* both values. '* '***********************************************************************/ Dim longitude As Double, Latitude As Double Dim Zone As Double, daySavings As Double Dim hh As Double, mm As Double, ss As Double, timenow As Double Dim jd As Double, t As Double, r As Double Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double 'Dim SolarDec As Double, earthRadVec As Double, solarTimeFix As Double Dim SolarDec As Double, solarTimeFix As Double Dim trueSolarTime As Double, hourangle As Double, harad As Double Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double Dim azimuth As Double, exoatmElevation As Double Dim step1 As Double, step2 As Double, step3 As Double Dim refractionCorrection As Double, Te As Double, solarzen As Double longitude = lon * -1 Latitude = lat If (Latitude > 89.8) Then Latitude = 89.8 If (Latitude < -89.8) Then Latitude = -89.8 Zone = timezone * -1 daySavings = dlstime * 60 hh = hours - (daySavings / 60) mm = minutes ss = seconds
timenow = hh + mm / 60 + ss / 3600 + Zone jd = calcJD(year, month, day) t = calcTimeJulianCent(jd + timenow / 24#) r = calcSunRadVector(t) alpha = calcSunRtAscension(t) theta = calcSunDeclination(t) Etime = calcEquationOfTime(t) eqtime = Etime SolarDec = theta '// in degrees earthRadVec = r
trueSolarTime = hh * 60# + mm + ss / 60# + solarTimeFix '// in minutes
trueSolarTime = trueSolarTime - 1440 Loop hourangle = trueSolarTime / 4# - 180# '// Thanks to Louis Schwarzmayr for the next line: If (hourangle < -180) Then hourangle = hourangle + 360# harad = degToRad(hourangle) csz = Sin(degToRad(Latitude)) * _ Sin(degToRad(SolarDec)) + _ Cos(degToRad(Latitude)) * _ Cos(degToRad(SolarDec)) * Cos(harad) If (csz > 1#) Then csz = 1#
ElseIf (csz < -1#) Then csz = -1# End If zenith = radToDeg(Application.WorksheetFunction.Acos(csz)) azDenom = (Cos(degToRad(Latitude)) * Sin(degToRad(zenith))) If (Abs(azDenom) > 0.001) Then azRad = ((Sin(degToRad(Latitude)) * _ Cos(degToRad(zenith))) - _ Sin(degToRad(SolarDec))) / azDenom If (Abs(azRad) > 1#) Then If (azRad < 0) Then azRad = -1# Else azRad = 1# End If End If
If (hourangle > 0#) Then azimuth = -azimuth End If Else
If (Latitude > 0#) Then azimuth = 180# Else azimuth = 0# End If End If
If (azimuth < 0#) Then azimuth = azimuth + 360# End If exoatmElevation = 90# - zenith If (exoatmElevation > 85#) Then refractionCorrection = 0# Else Te = Tan(degToRad(exoatmElevation)) If (exoatmElevation > 5#) Then refractionCorrection = 58.1 / Te - 0.07 / (Te * Te * Te) + _ 0.000086 / (Te * Te * Te * Te * Te) ElseIf (exoatmElevation > -0.575) Then step1 = (-12.79 + exoatmElevation * 0.711) step2 = (103.4 + exoatmElevation * (step1)) step3 = (-518.2 + exoatmElevation * (step2)) refractionCorrection = 1735# + exoatmElevation * (step3) Else refractionCorrection = -20.774 / Te End If refractionCorrection = refractionCorrection / 3600# End If solarzen = zenith - refractionCorrection solarazimuth = azimuth solarelevation = 90# - solarzen End Sub APPENDIX C: SEDIMENT-WATER HEAT EXCHANGEAlthough the omission of sediment-water heat exchange is usually justified for deeper systems, it can have a significant impact on the heat balance for shallower streams. Consequently, sediment-water heat exchange is included in QUAL2K. A major impediment to its inclusion is that incorporating sediment heat transfer often carries a heavy computational burden. This is because the sediments are usually represented as a vertically segmented distributed system. Thus, inclusion of the mechanism results in the addition of numerous sediment segments for each overlying water element. In the present appendix, I derive a computationally-efficient lumped approach that yields comparable results to the distributed methods. The conduction equation is typically used to simulate the vertical temperature distribution in a distributed sediment (Figure a)  ()
This model can be subjected to the following boundary conditions: 
where T = sediment temperature [oC], t = time [s], ? = sediment thermal diffusivity [m2 s?1], and z = depth into the sediments [m], where z = 0 at the sediment-water interface and z increases downward,  = mean temperature of overlying water [oC], Ta = amplitude of temperature of overlying water [oC], ? = frequency [s?1] = 2?/Tp, Tp = period [s], and ? = phase lag [s]. The first boundary condition specifies a sinusoidal Dirichlet boundary condition at the sediment-water interface. The second specifies a constant temperature at infinite depth. Note that the mean of the surface sinusoid and the lower fixed temperature are identical.
Figure . Alternate representations of sediments: (a) distributed and (b) lumped. Applying these boundary conditions, Eq. (230) can be solved for (Carslaw and Jaeger 1959)  ()
where ?’ [m?1] is defined as  ()
The heat flux at the sediment water interface can then be determined by substituting the derivative of Eq. (231) into Fourier’s law and evaluating the result at the sediment-water interface (z = 0) to yield  ()
where J(0, t) = flux [W/m2]. An alternative approach can be developed using a first-order lumped model (Figure b), 
where Hsed = the thickness of the sediment layer [m], ?s = sediment density [kg/m3], and Cps = sediment specific heat [joule (kg oC)]?1]. Collecting terms gives, 
where 
After initial transient have died out, this solution to this equation is  ()
which can be used to determine the flux as  ()
It can be shown that Eqs. (231) and (234) yield identical results if the depth of the single layer is set at  ()
Water quality models typically consider annual, weekly and diel variations. Using ? = 0.0035 cm2/s (Hutchinson 1957), the single-layer depth that would capture these frequencies can be calculated as 2.2 m, 30 cm and 12 cm, respectively. Because QUAL2K resolves diel variations, a value on the order of 12 cm should be selected for the sediment thickness. We have chosen of value of 10 cm as being an adequate first estimate because of the uncertainties of the river sediment thermal properties (Table ). 1 Notice that time is measured in seconds in this and other formulas used to characterize hydraulics. This is how the computations are implemented within Q2K. However, once the hydraulic characteristics are determined they are converted to day units to be compatible with other computations. 2 The nomenclature FNU stands for Filtration, Nitrification inhibition and Ultimate 3 ly/d = langley per day. A langley is equal to a calorie per square centimeter. Note that a ly/d is related to the ?E/m2/d by the following approximation: 1 ?E/m2/s ? 0.45 Langley/day (LIC-OR, Lincoln, NE). 4 The conversion, m3 = 1000 L is included because all mass balances express volume in m3, whereas total inorganic carbon is expressed as mole/L. 5 Dr. Pieter Tans, NOAA/ESRL (www.esrl.noaa.gov/gmd/cgg/trends 6 Note that although it will almost always be negligible, Eq. (193) PO43– ? for completeness. Yüklə 1,65 Mb. Dostları ilə paylaş:
|