/* ============================================================================ /* @(#)solradb.aml 1.2 01/14/1999 15:33:41 /*----------------------------------------------------------------------------- /* Lab for synthetic dynamic vegephenomenology /*----------------------------------------------------------------------------- /* Program: SOLRADB.AML /* Purpose: This is the Batch-mode version of Solrad.aml. It allows to /* calculate daily or average bi-weekly direct solar flux based /* on a core program developed by Lalit Kumar (see below) and a /* suite of input data from the SAMSON network. The AML reads /* transmittivity and (optionally) extraterrestrial direct nor- /* mal radiation from a files. Seasonal lapse rates for trans- /* mittivity are derived from a second file. Thus, the present /* AML allows to calculate actual direct or potential direct /* direct solar radiation. /*----------------------------------------------------------------------------- /* Usage: &r solradb /* /* where: trn = file name for transmittance data /* lps = file name for lapse rate data /* dem = file name of DEM /* out = base name of outgrid file names (w|out index) /* lat = latitude of study site (DEM) /* time = daily time interval for solar pass (in minutes) /* base = base elevation for transmittance data /* days = # of days the lapse rate data are averaged for /* sel = method for transmittance /* 1: maximum transmittance per period /* 2: horizontal rad.-based transmittance data /* 3: direct radiation based transmittance data /* 4: corr. direct rad. based transmittance data /*----------------------------------------------------------------------------- /* Notes: The core of this program was developed and used for work publi- /* shed in: Kumar, L., Skidmore, A.K., Knowles, E., (1997). Mo- /* delling Topographic Variation in Solar Radiation in a GIS En- /* vironment. International Journal of Geographical Information /* Science, 11(5): 475-497. This version intended to generate co- /* verages of potential direct solar radiation. The present ver- /* sion allows to use input data from a file, which summarizes /* SAMSON solar radiation data as generated by SUM_TRAN.EXE, a /* ForTran-based program to analyse SAMSON stations data. This /* enables to adjust the following parameters in the calculation /* of solarradiation: (1) actual direct extraterrestrial radia- /* tion, (2) avg. actual transmittance, (3) max. actual trans- /* mittance, and (4) elevational lapse rates for transmittance. /* The last parameter is derived from a second file. Thus, the /* present version allows to calculate actual direct solar ra- /* diation, instead of potential solarflux. Also, it allows to /* incorporate actual direct extraterrestrial radiation from a /* file. /*----------------------------------------------------------------------------- /* Input: transmittance, solar constant, elev. lapse rates (daily, or /* bi-weekly) /* Output: direct actual solar radiation (daily, or bi-weekly) /*----------------------------------------------------------------------------- /* History: L. Kumar - Original coding of shortwave.aml /* L. Kumar - 5/08/1997 - Last modification of shortwave.aml /* N.E. Zimmermann - 1/13/1999 - Original coding of solrad.aml /* N.E. Zimmermann - 3/04/1999 - Last modification /* N.E. Zimmermann - 3/13/1999 - New coding for new output /*============================================================================= &args .fl1 .fl2 .dem .out .lat .int .bas .day .met &if [show program] ne GRID &then grid &sys clear &s file1 = %.fl1% &s file2 = %.fl2% &s cover = %.dem% &s outgrid = %.out% &s latdeg = %.lat% &s timeint = %.int% &s basealt = %.bas% &s dayint = %.day% &s sel = %.met% &s funit1 := [open %file1% ostat1 -read] /* check for presence &if %ostat1% ne 0 &then &return &warning Error opening %file1% &s funit2 := [open %file2% ostat2 -read] /* check for presence &if %ostat2% ne 0 &then &return &warning Error opening %file2% &s line1 := [read %funit1% rstat1] /* check for data in file1 &if %rstat1% ne 0 &then &return &warning Could not read %file1%. File may be empty. &s line2 := [read %funit2% rstat2] /* check for data is file2 &if %rstat2% ne 0 &then &return &warning Could not read %file2%. File may be empty. &if ^ [exists %cover% -grid] &then &return &error Grid %cover% does not exist. &if ^ [variable latdeg] &then &return &error Please enter a valid latitude &if %latdeg% gt 90 or %latdeg% lt -90 &then &return &error Please enter latitude between -90 and 90 degrees &if ^ [variable timeint] or %timeint% <= 0 &then &return &error Please enter a time greater than zero &if ^ [variable basealt] or %basealt% < 0 &then &return &error Please enter an elevation >= zero &if ^ [variable dayint] or %dayint% <= 0 &then &return &error Please enter a Day-interval greater than zero &if %sel% NE 1 AND %sel% NE 2 AND %sel% NE 3 &then &return &warning Please enter a valid selection (1, 2, or 3) &s line1 := [read %funit1% rstat1] &if %rstat1% ne 0 &then &return &warning No data beyond header in %file1% &s line2 := [read %funit2% rstat2] &if %rstat2% ne 0 &then &return &warning No data beyond header in %file2% &do &while %rstat1% eq 0 &if %rstat2% ne 0 &then &return &warning No corresponding lapse rate data for this day &r rdtrans [unquote %line1%] &r rdlpsrt [unquote %line2%] &if (%sel% EQ 1) &then &s attn = %.trnmx% /* Is performed uncorrected, is a semi-saccurate assessment of clear-sky transmittance &if (%sel% EQ 2) &then &s attn = %.trnhr% /* This tau is uncorrected for solar altitude &if (%sel% EQ 3) &then &s attn = %.trndr% /* This tau is uncorrected for solar altitude &if (%sel% EQ 4) &then &s attn = %.trndrc% /* Is the most accurate assessment of the actual transmittance &if [exist %outgrid%_%.id1% -grid] &then &return &error Grid %outgrid%_%.id1% already exists. &if %.id1% ne %.id2% &then &return &warning The two parameter files do not have the same ID order &s jday = [calc [calc %.id1% * %dayint%] - [truncate [calc %dayint% / 2 ]] + 1] &s daystart = %jday% &s dayend = %jday% /* Here the core program of Lalit Kumar starts ------------------------------------- &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 initialgrid1 -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 NEZ: is adjusted to read Io from the file! ----------------------- &s Io = [calc 1.367 * ( 1 + 0.0344 * [cos [calc 360 * %degtorad% * %daynumber% / 365]] )] /* &s Io = %.mxdnet% / 1000.0 /*NEZ: This can alternatively be switched on instead of /*NEZ: the line above to enable to read Io from file. /* See Eq 4 &s decl = [calc 23.45 * %degtorad% * [sin [calc %degtorad% * 360 * ~ ( 284 + %daynumber% ) / 365]]] /* See Eq 6 &s sunrise = [acos [calc -1 * [tan %latitude%] * [tan %decl%]]] &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 &s M = [calc ( 1229 + ( 614 * [sin %solaralt%] ) ** 2 ) ** 0.5 - 614 * [sin %solaralt%]] &s e = 2.7182818 /* NEZ: This function was added to simplify the calculation of tau &s base = [calc ( %e% ** ( -0.65 * %M% ) + %e% ** ( -0.095 * %M% ) ) ] /* See Eq 15 /* NEZ: This Eqn. was removed to enable the elevation-dependent calcs of transmittivity /* &s Is = %Io% * 0.56 * ( %e% ** ( -0.65 * %M% ) + %e% ** ( -0.095 * %M% ) ) /* NEZ: Insted, only the Tau is calculated and combined with %Io% and %attn% later on! /* NEZ: Kumar's Tau has a maximum of 0.8, which is adjusted afterwards by the input data &if %sel% EQ 4 &then &s tau = %.tba% * %base% &else &s tau = %attn% / .8 * .56 * %base% &s solaraltdeg = %solaralt% * 57.29578 &if ( %solarazdeg% <= 180 ) &then /* NEZ: added to cope w. N-Lat.s &s azi = [calc ( 180 - %solarazdeg% )] /* NEZ: " &else /* NEZ: " &s azi = [calc ( 180 + ( 360 - %solarazdeg% ) )] /* NEZ: " sungrid = HILLSHADE (%cover%, %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 /* NEZ: This function is replaced to additionally adjust for obstructed skies /* wattsgrid = %Is% * cosi * sungrid * shaded * 60 * %timeint% wattsgrid = %Io% * (%tau% + ((%cover% - %basealt%) * %.lpsrt%)) * 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%_%.id1% = 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 /* ---- Here the AML of Lalit Kumar ends. Below is some cleaning up --------- &s line1 := [read %funit1% rstat1] &s line2 := [read %funit2% rstat2] &end &if [close %funit1%] ne 0 &then &return &warning Error closing %file1% &if [close %funit2%] ne 0 &then &return &warning Error closing %file2% /*---- Finish up the calculations ---------------------------------------------- &ty &ty &return End of calculations &ty &ty