/* ======================================================================== /* @(#)shortwavb.aml 1.4 09/14/2000 11:22:17 /*------------------------------------------------------------------------- /* Lab for synthetic dynamic vegephenomenology /*------------------------------------------------------------------------- /* Program: SHORTWAVB.AML /* Purpose: This program calculates the shortwave radiation received /* at the surface of the earth over a period of time. For the /* given day(s), it calculates the sunset and sunrise times /* and intergrates solar radiation from sunrise to sunset each /* day. /* Reference: This AML program was used originally for work published as: /* Kumar, L., Skidmore, A.K., Knowles, E., (1997) Modelling /* Topographic Variation in Solar Radiation in a GIS Envi- /* ronment. International Journal of Geographical Informa- /* tion Science, 11(5): 475-497. /*------------------------------------------------------------------------- /* Usage: &r shortwavb /*------------------------------------------------------------------------- /* Notes: This AML was originally coded as "shortwave.aml" by Lalit /* Kumar. For northern hemispheres, the code had an error. /* When calculating the terrain-based overhsadowing, the /* sun passed always north of the target cell. However, the /* calculation of the radiation intensity was correct. To /* adjust the code for northern hemispheres, the shortwavc.aml /* was developed. /* The current version allows batch operations. It is based /* on the corrected version of "shortwavc.aml". /* In ArcInfo 8.0.2 raises to fractional powers cause errors, /* which is a bug that will be fixed in AI 8.2. CALVIN TOLKAMP /* suggested the workaround for this problem. /*------------------------------------------------------------------------- /* Addresses: Lalit Kumar /* School of Geography /* University of New South Wales /* Sydney 2052 /* Australia /* /* p2114659@geog.unsw.edu.au [nez: doesn't work anymore] /* /* Niklaus E. Zimmermann /* Swiss Federal Research Institute WSL /* Dept. of Landscape Dynamics /* CH-8903 Birmensdorf /* Switzerland /* /* niklaus.zimmermann@wsl.ch /* /*------------------------------------------------------------------------- /* Input: DEM /* Output: radiation covers /*------------------------------------------------------------------------- /* History: Lalit Kumar - ... - Original coding /* Lalit Kumar - 5/08/1997 - Last Modifications /* Niklaus E. Zimmermann - 11/21/1998 - Corrections (N-hemisph.) /* Niklaus E. Zimmermann - 9/08/2000 - Adjustments to AI 8.0.2 /* Niklaus E. Zimmermann - 9/26/2002 - Corr. latitudes > +/-66° /* Niklaus E. Zimmermann - 10/07/2003 - Corr kill initialgrid(1) /*========================================================================= /* /* &args cover outgrid latdeg daystart dayend timeint &if [show program] ne GRID &then grid &if ^ [exists %cover% -grid] &then &return &error Grid %cover% does not exist. &if [exist %outgrid% -grid] &then &return &error Grid %outgrid% already exists. &if ^ [variable latdeg] &then &return &error need a valid latitude &if ^ [variable daystart] or %daystart% < 1 or %daystart% > 365 &then &return &error please enter a start day between 1 and 365 &if ^ [variable dayend] or %dayend% < 1 or %dayend% > 365 &then &return &error please enter an end day between 1 and 365 &if ^ [variable timeint] or %timeint% <= 0 &then &return &error please enter a time greater than zero &s time = %timeint% * 0.0043633 /* convert minutes into radians &s daynumber = %daystart% &s pi = 22 / 7 &s degtorad = 2 * %pi% / 360 &s latitude = %latdeg% * %degtorad% &if [exist slopegrid -grid] &then kill slopegrid all slopegrid = slope(%cover%) &if [exist slopegrid2 -grid] &then kill slopegrid2 all slopegrid2 = slopegrid * %degtorad% &if [exist aspectgrid -grid] &then kill aspectgrid all aspectgrid = aspect(%cover%) &if [exist aspectgrid1 -grid] &then kill aspectgrid1 all /* change so that surface azimuth is zero at south if (aspectgrid == -1) aspectgrid1 = 0 else if (aspectgrid le 180) aspectgrid1 = 180 - aspectgrid else aspectgrid1 = 540 - aspectgrid &if [exist aspectgrid2 -grid] &then kill aspectgrid2 all aspectgrid2 = aspectgrid1 * %degtorad% &if [exist outgrid0 -grid] &then kill outgrid0 all outgrid0 = hillshade (%cover%, 0, 0, shadow) &if [exist initialgrid -grid] &then kill initialgrid all initialgrid = outgrid0 * 0 kill outgrid0 all kill slopegrid all kill aspectgrid all kill aspectgrid1 all &if [exist sungrid -grid] &then kill sungrid all &if [exist wattsgrid -grid] &then kill wattsgrid all &if [exist wattsgrid2 -grid] &then kill wattsgrid2 all &if [exist shaded -grid] &then kill shaded all &if [exist cosi -grid] &then kill cosi all &do &until %daynumber% > %dayend% /* See Eq 7 &s Io = [calc 1.367 * ( 1 + 0.034 * [cos [calc 360 * %degtorad% * ~ %daynumber% / 365]] )] /* See Eq 4 &s decl = [calc 23.45 * %degtorad% * [sin [calc %degtorad% * 360 * ~ ( 284 + %daynumber% ) / 365]]] /* See Eq 6 /* This line needed to be adjusted to latitudes exceeding 66.5° N/S &if [calc -1 * [tan %latitude%] * [tan %decl%]] lt -1.0 &then /*nez: check that x is not >1.0 for acos(x) &s sunrise = [acos -1.0] /*nez: " &else &if [calc -1 * [tan %latitude%] * [tan %decl%]] gt 1.0 &then /*nez: " &s sunrise = [acos 1.0] /*nez: " &else /*nez: " &s sunrise = [acos [calc -1 * [tan %latitude%] * [tan %decl%]]] /*orig " &s sunset = -1 * %sunrise% /* To ensure calculations are at half the time interval &s hourangle = %sunrise% - ( %time% / 2 ) &s pass = 1 &do &until %hourangle% < %sunset% &ty day : %daynumber% &ty pass: %pass% /* See Eq 2 &s solaralt = [asin [calc [sin %latitude%] * [sin %decl%] + [cos %latitude%] * ~ [cos %decl%] * [cos %hourangle%]]] &s test = [tan %decl%] / [tan %latitude%] /* See Eq 3 &if ( [cos %hourangle%] > %test% ) &then &s solaraz = [asin [calc [cos %decl%] * [sin %hourangle%] / [cos %solaralt%]]] &else &if ( [cos %hourangle%] < %test% ) &then &s solaraz = %pi% - [asin [calc [cos %decl%] * [sin %hourangle%] / [cos %solaralt%]]] &else &if ( %test% = [cos %hourangle%] and %hourangle% >= 0 ) &then &s solaraz = %pi% / 2 &else &if ( %test% = [cos %hourangle%] and %hourangle% < 0 ) &then &s solaraz = -1 * %pi% / 2 &end /* don't ask :)) &if ( %solaraz% >= 0 ) &then &s solarazdeg = %solaraz% * 57.29578 &else &s solarazdeg = 360 - ( [abs %solaraz%] * 57.29578 ) /* See Eq 11 /* Below is the original code, which doesn't run under AI 8.0.2 (no raise to fractional powers) /*&s M = [calc ( 1229 + ( 614 * [sin %solaralt%] ) ** 2 ) ** 0.5 - 614 * [sin %solaralt%]] /* Calvin Tolkamp proposed the workaround below: &s M = [calc [sqrt [calc 1229 + ( 614 * [sin %solaralt%] ) ** 2 ]] - 614 * [sin %solaralt%]] &s e = 2.7182818 /*This is to back-calculate the tau-value /*nez /*&s tau = [calc ( 0.56 * ( %e% ** ( -0.65 * %M% ) + %e% ** ( -0.095 * %M% ) ) )] /*nez /*&ty Solaralt: %solaralt%, Air mass: %M%, Tau: %tau% /*nez /* See Eq 15 /* Below is the original code, which doesn't run under AI 8.0.2 (no raise to fractional powers) /*&s Is = %Io% * 0.56 * ( %e% ** ( -0.65 * %M% ) + %e% ** ( -0.095 * %M% ) ) /* Calvin Tolkamp proposed the workaround below: &s Is = %Io% * 0.56 * ( [exp [calc -0.65 * %M% ]] + [exp [calc -0.095 * %M% ]] ) &s solaraltdeg = %solaralt% * 57.29578 /*sungrid = HILLSHADE (%cover%, %solarazdeg%, %solaraltdeg%, shadow) /*nez: was disabled because it is wrong in N-lat. &if ( %solarazdeg% <= 180 ) &then /*nez: added to correct for northern lat.s &s azi = [calc ( 180 - %solarazdeg% )] /*nez: " &else /*nez: " &s azi = [calc ( 180 + ( 360 - %solarazdeg% ) )] /*nez: " /*&type %azi%, %solarazdeg% /*nez: " sungrid = HILLSHADE (%cover%, %azi%, %solaraltdeg%, shadow) /*nez: " /*&ty outgrid = hillshade(DEM, %solarazdeg%, %solaraltdeg%, SHADOW) /*nez: test for corrections -> are OK /*&ty outgrid = hillshade(DEM, %azi%, %solaraltdeg%, SHADOW) /*nez: " /* See Eq 16 cosi = sin(%decl%) * (sin(%latitude%) * cos(slopegrid2) - cos(%latitude%) * ~ sin(slopegrid2) * cos(aspectgrid2)) + cos(%decl%) * cos(%hourangle%) * ~ (cos(%latitude%) * cos(slopegrid2) + sin(%latitude%) * sin(slopegrid2) * ~ cos(aspectgrid2) ) + cos(%decl%) * sin(slopegrid2) * sin(aspectgrid2) * ~ sin(%hourangle%) /* Trust me on the following two lines (You always need to fudge a little). if (cosi lt 0) shaded = 0 else shaded = 1 wattsgrid = %Is% * cosi * sungrid * shaded * 60 * %timeint% wattsgrid2 = wattsgrid + initialgrid kill wattsgrid all kill initialgrid all kill sungrid all kill shaded all kill cosi all rename wattsgrid2 initialgrid &s hourangle = %hourangle% - %time% &s pass = %pass% + 1 &end &s daynumber = %daynumber% + 1 &end %outgrid% = int(initialgrid) /* NOTE : THE UNITS FOR RADIATION IN THE OUTGRID ARE kJ/m^2/timeperiod /* If these are to be converted to MJ/m^2/day, then they should be divided by 1000 * no. of days. kill initialgrid all kill aspectgrid2 all kill slopegrid2 all