subroutine sunmap(sdate,r2,r5smap,halfdy,dsunmp) c c +++PURPOSE+++ c This routine determines the amount of heat unit calories c available for melting snow. The logic as well as the variable c comments are taken from "Computational Algorithm for c Solar Radiation on Mountain Slopes" paper written by c Lloyd W. Swift Jr. and R.J. Luxmoore, 06-73. The value is then c distributed over the daily radiation curve and passed to the c weather data structure as HRADMJ. c c Author(s): Cully Hession and Bruce Lucord, USDA-ARS-NCSRL c Revised by John Witte, UofMn WCES @ USDA-ARS-NCSRL c Date: 03/24/93 c Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051 c August 1994 c c +++ARGUMENT DECLARATIONS+++ c real r2,r3,r5smap,halfdy c integer sdate,hour,irdflg real r2,r3,halfdy,r5smap,dsunmp integer sdate c c +++ARGUMENT DEFINITIONS+++ c r2 - Measured daily solar radiation from cligen file (Ly/day). c r3 - Potential solar beam on horizontal surface (MJ/m^2/day). c halfdy - Number of hours from sunrise until noon. c sdate - Julian date that is being run. c c + + + PARAMETERS + + + c include 'pmxhil.inc' include 'pmxpln.inc' c c +++COMMON BLOCKS+++ include 'cangie.inc' c read: eqlat,delong,radlat,radpot c include 'cclim.inc' c read: rpoth,hradmj c c +++LOCAL VARIABLES+++ real d,e,r1,solcon,t,t7,t6,x,t1,t0,t3,t2,r4,pi,t4, 1 f,cratio,langmj c c +++LOCAL DEFINITIONS+++ c d - Declination of the sun in radian units,(+)North (-)South. c e - Radius vertor of the sun in radian units. c r1 - Solar constant factor. c solcon - Solar constant (Ly/min). c t - Temporary variable. c t7 - Hour Angle of sunset on horizontal surface (degrees). c t6 - Hour Angle of sunrise on the equivalent slope surface. c x - Temp variable used to determine HA on the equiv slope. c t1 - Hour Angle of sunset on the horizontal surface (deg). c t0 - Hour Angle of sunrise on the horizontal surface (deg). c t3 - Hour Angle of sunset on the sloping surface (deg). c t2 - Hour Angle of sunrise on the sloping surface (deg). c r4 - Potential solar beam on the sloping surface (MJ/m^2/day) c pi - Estimated value for 22/7. c t4 - Number of hours before solar noon of sunrise on slope. c langmj - Conversion factor from Cal/cm^2 = Ly to MJ/M2. c f - Slope factor (r4/r3). c r6 - Est solar rad on sloping surface map area (MJ/m^2/day). c cratio - Hourly radiation adjustment ratio calculated in RADCURV. c c +++OUTPUT FORMATS+++ c2000 format(' Inside SunMap Routine ',/,' d = ',f10.3,/,' e = ', c 1 f10.3,/,' r1 = ', c 2 f10.3,/,' t7 = ',f10.3,/,' t6 = ',f10.3,/,' t3 = ',f10.3, c 3 /,' t2 = ',f10.3,/,' r4 = ',f10.3,/,' t4 = ',f10.3,/) c c +++END SPECIFICATIONS+++ c pi = 3.141593 solcon = 1.94 langmj = 0.04184 c c -- begin calculations -- d = 0.00698 - 0.4067 * cos((sdate + 10.0) * 0.0172) e = 1.00000 - 0.0167 * cos((sdate - 3.0) * 0.0172) r1 = (60.0 * solcon) / (e * e) x = -((sin(eqlat)) / (cos(eqlat))) * ((sin(d)) / (cos(d))) if (x .gt. 1.0) x = 1.0 if (x .lt.-1.0) x = -1.0 c t = acos(x) t7 = t - delong t6 = -t - delong c -- NOTE that the "x" calculations are written out the long way rather c -- than by simply using the "tan" function. This is do to a bug in c -- Lahey compiler with tangent that Dennis F. found on AT&T 6300's. c c x = -((sin(radlat))/(cos(radlat)))*((sin(d))/(cos(d))) c c Use tangent function - as 6300's obsolete and now using newer c Lahey compiler. dcf 2/22/94 c x = - (tan(radlat) * tan(d)) if (x .gt. 1.0) x = 1.0 if (x .lt.-1.0) x = -1.0 t = acos(x) t1 = t t0 = -t t3 = t7 if (t7 .ge. t1) t3 = t1 t2 = t6 if (t6 .le. t0) t2 = t0 r4 = r1 * ((sin(d) * sin(eqlat) * (t3 - t2) * 12. / pi) 1+ (cos(d) * cos(eqlat) * (sin(t3 + delong) - sin(t2 + delong)) 2* 12. / pi)) radpot=r4 r4 = r4 * langmj c reza 3/2/94 radpot = r4 c -- Now the potential radiation on slope is in units of MJ/m^2. t4 = t2 * 12.0 / pi halfdy = abs(t4) r3 = r1 * ((sin(d) * sin(radlat) * (t1 - t0) * 12.0 / pi) 1 + (cos(d) * cos(radlat) * (sin(t1) - sin(t0)) * 12.0 / pi)) r3 = r3 * langmj c -- Now the potential radiation on horiz is also in MJ/m^2. f = r4 / r3 c Reza, the next line will make r5smap to be in LY c because r2=radly, we need to convert r5smap to MJ r5smap = f * r2 c Reza changed the next line, 3/9/94 c -- r5smap = r5smap * langmj r5smap = r5smap * langmj c -- The above was commented out because f and r2 are already in c -- the desired units. c -- Make sure that r5smap is in MJ/m*m units. c r6 = r5smap / cos(radinc) c -- The above line is commented out due to Bob's suggestion c -- we do so. It was used to project sloping surface to the c -- horizontal surface (slightly smaller area) resulting in c -- slightly higher radiation. We feel that when dealing with c -- non-mountain slope segments, we don't need such an adjustment. rpoth = r3 dsunmp = d c creza 11/16/94 c if(irdflg .eq. 1) then c -- Now we call the RADCUR routine to calc hourly radiation. c call radcur(sdate,hour,radlat,d,cratio) c cratio = cratio / r3 c -- Now we know the ratio of daily radiation to hourly radiation. c reza, in the next line hradmj is in MJ c hradmj = r5smap * cratio c endif creza change 11/16/94 return end