| Simple Solar Spectral Model for Direct and Diffuse Irradiance on Horizontal and Tilted Planes at the Earth's Surface for Cloudless Atmospheres |
Table of Contents
PROGRAM LISTING
The FORTRAN listing of the original computer program SPCTRAL2, used on a CONTROL Data Corporation Cyber 750 computer, is presented here.
PROGRAM SPCTRAL2 ( INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT, +TAPE1,TAPE2,TAPE3) C C...TAPE1 = INPUT PARAMETERS C...TAPE2 = EXTRATERRESTRIAL SPECTRAL IRRADIANCE, ABSORPTION COEFFICIENTS. C...TAPE3 = OUTPUT FILE C...TAPE5 = TERMINAL SCREEN C C...THIS CODE WAS WRITTEN BY RICHARD BIRD AND CALCULATES SPECTRAL C IRRADIANCE OR PHOTON FLUX ON A TILTED OR HORIZONTAL SURFACE. C MODIFICATIONS TO THE CODE TO CONVERT SPECTRAL IRRADIANCE TO PHOTON C FLUX PER WAVELENGTH OR ELECTRON VOLT WERE MADE BY CAROL RIORDAN. C COMMENTS AND PROGRAM INSTRUCTIONS WERE ADDED BY RIORDAN 8/84. C C ******** INPUT VARIABLES ******** C AI = ANGLE OF INCIDENCE OF DIRECT BEAM ON FLAT SURFACE (DEG) C ALPHA= POWER ON ANGSTRON TURBIDITY EXPRESSION (1.14 FOR RURAL) C NDAY = JULIAN DAY C NW = NUMBER OF WAVELENGTHS FOR THIS RUN C O3 = OZONE AMOUNT (ATM CM) C RHO = GROUND ALBEDO C TILT = TILT ANGLE OF SURFACE FROM THE HORIZONTAL (DEG) C SPR = SURFACE PRESSURE (MILLIBARS) C TAU5 = THE AEROSOL OPTICAL DEPTH AT 0.5 MICRONS (BASE E) C W = PRECIPITABLE WATER VAPOR (CM) C WV(I) & R(I) GROUND REFLECTIVITY AT SPECIFIC WAVELENGTHS C Z = SOLAR ZENITH ANGLE (DEG) C C ***** TO RUN THIS PROGRAM ***** C C... IN THE FIRST LINE OF INPUT ON TAPE1, CHOOSE FROM THE C FOLLOWING OPTIONS: C COLLECTOR MODE C (MODE) 1=GLOBAL NORMAL AND DIRECT NORMAL C (TRACKING FLAT PLATE/CONCENTRATOR) C 2=GLOBAL TILT (FIXED FLAT PLATE) C 3=GLOBAL HORIZONTAL C C OUTPUT UNITS (NUNITS) 1=IRRADIANCE PER WAVELENGTH C 2=PHOTON FLUX PER WAVELENGTH C 3=PHOTON FLUX PER ELECTRON VOLT C C NUMBER OF SPECTRA (NUMSPT) C C CHECK THE INPUT PARAMETERS ON THE SECOND LINE OF TAPE1, ESPECIALLY C TAU5, W, TILT, SPR, NDAY. Z AND AI ARE READ IN THE "DO 15 LOOP". C C ******* THE PROGRAM ***** DIMENSION WV(6),R(6) C C...READ INPUT READ(1,1,) MODE,NUNITS,NUMSPT READ(1,2) TAU5,ALPHA,O3,W,TILT,SPR,NDAY,NW READ(1,3) WV(1),R(1),WV(2),R(2),WV(3),R(3),WV(4),R(4),WV(5),R(5), + WV(6),R(6) C C...THESE CONSTANTS ARE NEEDED TO CALCULATE PHOTON FLUX. C XH=6.626176*(10.**(-34.)) XC=2.99792458*(10.**8.) EVOLT=1.6021892*(10.**(-19.)) CONST=(1./(XH*XC))*(10.**(-10.)) XC=XC*(10.**6.) C C...RADIANS PER DEGREE RPD=.0174533 C C...ERV IS THE EARTH RADIUS VECTOR-CORRECTION FOR SUN DISTANCE TT=6.283185*(NDAY-1)/365. TT2=2.*TT ERV=1.00011+.034221*COS(TT)+.00128*SIN(TT) + + .000719*COS(TT2)+.000077*SIN(TT2) C C...WRITE OUT INITIAL CONDITIONS IF(MODE .EQ. 1) WRITE(5,20) IF(MODE .EQ. 2) WRITE(5,21) IF(MODE .EQ. 3) WRITE(5,22) WRITE(5,4) NUNITS,NUMSPT WRITE(5,5) NDAY,ERV WRITE(5,7) TAU5,ALPHA,W,O3 WRITE(5,8) WV(1),R(1),WV(2),R(2),WV(3),R(3),WV(4),R(4), + WV(5),R(5),WV(6),R(6) C C...OMEG AND OMEGP ARE USED IN THE SINGLE SCATTERING ALBEDO C CALCULATIONS. JUSTUS HAS A FORM OF OMEGL THAT VARIES WITH C RELATIVE HUMIDITY AS WELL AS WAVELENGTH. C OMEG=0.945 OMEGP=0.095 C C...FS AND FSP ARE THE FORWARD/TOTAL SCATTERING RATIOS AS C A FUNCTION OF ZENITH ANGLE. (AIR MASS FOR FSP IS FIXED IN THIS CODE, C FS IS IN THE "DO 15 LOOP". GG IS THE AEROSOL ASSYMETRY C FACTOR (0.65 USED FOR RURAL). C GG=0.65 ALG=ALOG(1.-GG) AFS=ALG*(1.459+ALG*(0.1595+ALG*0.4129)) BFS=ALG*(0.0783+ALG*(-0.3824-ALG*0.5874)) FSP=1.-0.5*EXP((AFS+BFS/1.8)/1.8) C C...RR IS USED IN THE OZONE MASS EXPRESSION IN THE "DO 15 LOOP". RR=22./6370. C C...THESE EXPRESSIONS DEPEND ON ZENITH AND INCIDENCE ANGLE AND C SHOULD BE INSIDE THE LOOP THAT CHANGES Z AND AI. C C...FOR GLOBAL NORMAL SPECTRA (MODE=1), INCIDENT ANGLE(AI)=0 AND C TILT= Z. FOR GLOBAL HORIZONTAL SPECTRA C (MODE=3), INCIDENCE ANGLE (AI)=Z AND TILT=0. C C...DO LOOP FOR SEVERAL SPECTRA (NUMSPT) D0 15 KK=1, NUMSPT READ(1,10) Z, AI IF (MODE .EQ. 1) AI=0. IF (MODE .EQ. 1) TILT=Z IF (MODE .EQ. 3) TILT=0. IF (MODE .EQ. 3) AI=Z COSTLT=COS(TILT*RPD) CI=COS(AI*RPD) ZCOS=COS(Z*RPD) ZSIN=SIN(Z*RPD) FS=1.-0.5*EXP((AFS+BFS*ZCOS)*ZCOS) C C...RELATIVE OPTICAL AIR MASS. AM=1./(ZCOS+.15*(93.885-Z)**(-1.253)) C...PRESSURE CORRECTED AIR MASS. AMP=AM*SPR/1013. C C...OZONE MASS. AMO=(1.+RR)/(ZCOS**2.+2.*RR)**.5 C C...WRITE OUT INITAL CONDITIONS TO TTERMINAL SCREEN C WRITE(5,17) Z,TILT,AI WRITE(5,6) AM,SPR,AMP,AMO C C...INITIALIZE THE INTERPOLATION COUNTER FOR RHO. NR=2 C...REWIND THE TAPE WITH ET SPECTRUM AND ABSORPTION COEFFICIENTS FOR C EACH SPECTRUM. REWIND 2 C C...DO LOOP FOR NUMBER OF WAVELENGTHS (NW) DO 14 I=1,NW READ(2,9) WVL,H0,AW,AO,AU C C...CORRECT EXTRATERRESTRIAL IRRADIANCE FOR EARTH RADIUS VECTOR. H0=H0*ERV C C C...OMEGL IS THE WAVELENGTH-DEPENDENT SINGLE SCATTERING ALBEDO C OMEGL=OMEG*EXP(-OMEGP*(ALOG(WVL/0.4))**2.) C C...THE FOLLOWING STATEMENTS PRODUCE THE WAVELENGTH-DEPENDENT ALBEDO C IF(WVL .GT. WV(NR)) NR=NR+1 SLP=(R(NR)-R(NR-1))/(WV(NR)-WV(NR-1)) RHO=SLP*(WVL-WV(NR-1))+R(NR-1) C C...CALCULATE TRANSMITTANCE TR=EXP(-AMP/(WVL**4.*(115.6406-1.335/WVL**2.))) TO=EXP(-AO*O3*AMO) TW=EXP(-.2385*AW*W*AM/(1.+20.07*AW*W*AM)**.45) TU=EXP(-1.41*AU*AMP/((1.+118.93*AU*AMP)**.45)) DELA=TAU5*(WVL/.5)**(-ALPHA) TAS=EXP(-OMEGL*DELA*AM) TAA=EXP(-(1.-OMEGL)*DELA*AM) TA=EXP(-DELA*AM) C DIR=H0*TR*TO*TW*TU*TA C...CALCULATE DIRECT COMPONENT OF IRRADIANCE ONTO SURFACE DIRSUR=DIR*CI C C...DRAY AND DAER HAVE BEEN MODIFIED BY BIRD. NOTE POWER ON TR TERM. C DRAY=H0*ZCOS*TO*TW*TU*TAA*(1.-TR**.95)*.5 DAER=H0*ZCOS*TO*TW*TU*TAA*TR**1.5*(1.-TAS)*FS TRP=EXP(-1.8/(WVL**4.*(115.6406-1.335/WVL**2))) TWP=EXP(-.2385*AW*W*1.8/((1.+20.07*AW*W*1.8)**.45)) TUP=EXP(-1.41*AU*1.8/((1.+118.93*AU*1.8)**.45)) TASP=EXP(-OMEGL*DELA*1.8) TAAP=EXP(-(1.-OMEGL)*DELA*1.8) RHOA=TUP*TWP*TAAP*(.5*(1.-TRP)+(1.-FSP)*TRP*(1.-TASP)) DRGD=(DIR*ZCOS+(DRAY+DAER))*RHO*RHOA/(1.-RHO*RHOA) DIF=DRAY+DAER+DRGD C C...CRC IS A UV CORRECTION FACTOR C CRC=1.0 IF (WVL .LE. .45) CRC=(WVL+.55)**1.8 C C...DIFFUSE ON A HORIZONTAL SURFACE DIF=DIF*CRC C C...TOTAL ON A HORIZONTAL SURFACE DTOT=DIR*ZCOS+DIF C...MAKE DIFS=DIF TO USE FORMAT 12 FOR GLOBAL HORIZONTAL CASE. DIFS=DIF IF (MODE .EQ. 3) GOTO 79 C C...GROUND REFLECTED COMPONENT REFS=DTOT*RHO*(1.0-COSTLT)/2.0 C C...THE THREE FOLLOWING STATEMENTS ARE THE HAY TILT ALGORITHM C C...ANISOTROPY INDEX AII=DIR/H0 C C...CIRCUMSOLAR AND ISOTROPIC COMPONENT (WEIGHTED BY DIF*AII AND C DIF*(1-AII). DIFSC=DIF*AII*CI/ZCOS DIFSI=DIF*(1.0-AII)*(1.0+COSTLT)/2.0 C C...DIFFUSE ON TILTED SURFACE. DIFS=DIFSC+DIFSI+REFS C C...TOTAL ON TILTED SURFACE. DTOT=DIR*CI+DIFS C C...PHOTON FLUX CALCULATIONS. C C...WRITE SPECTRAL IRRADIANCE OR PHOTON FLUX OUTPUT TO TAPE3. 79 IF (NUNITS .EQ.1) GOTO 12 PFWVGL=DTOT*WVL*CONST ENERGY=(XH*XC)/WVL E=ENERGY/EVOLT PFEVGL=(PFWVGL*WVL)/E PFWVDN=DIR*WVL*CONST PFEVDN=(PFWVDN*WVL)/E C PFWVDS=DIRSUR*WV*CONST C PFEVDS=(PFWVDS*WVL)/E C.... DECIDE IF USING DIF OR DIFS C PFWVDF=DIF*WVL*CONST C PFEVDF=(PFWVDF*WVL)/E IF (NUNITS .EQ. 2) GOTO 11 WRITE (3,18) E,PFEVGL,PFEVDN GOTO 13 11 WRITE (3,18) WVL,PFWVGL,PFWVDN GOTO 13 12 WRITE (3,16) WVL,DTOT,DIR,DIFS 13 CONTINUE C 14 CONTINUE 15 CONTINUE C C****** FORMAT STEMENTS ***** C 1 FORMAT(3I5) 2 FORMAT(6F6.3,2I5) 3 FORMAT(12F6.2) 4 FORMAT(2X,"NUNITS=",I6,3X,"NUMSPT=",I6/) 5 FORMAT(2X,"NDAY =",I7,3X,ERV =",F7.4) 6 FORMAT(2X,"AM =",F7.2,3X,"SPR =",F7.2,3X,"AMP =",F7.2,3X, +"AMO =",F7.2) 7 FORMAT(2X,"TAU5 =",F7.2,3X,"ALPHA=",F7.2,3X,"W =",F7.2, +3X,"O3 =",F7.3) 8 FORMAT(2X,"W1=",F6.2,2X,"R1=",F6.2/2CX,"W2=",F6.2,2X, +"R2=",F6.2/2X,"W3=",F6.2,2X,"R3=",F6.2/2X,"W4=",F6.2,2X, +"R4=",F6.2/2X,"W5=",F6.2,2X,"R5=",F6.2/2X,"W6=",F6.2,2X, +"R6=",F6.2/) 9 FORMAT(5F10.4) 10 FORMAT(2F10.4) 16 FORMAT(F7.4,3F10.4) 17 FORMAT(2X,"Z =",F7.2,3X,"TILT =",F7.2,3X,"AI =",F7.2) 18 FORMAT(F10.5,2E10.4) 20 FORMAT(1X,"GLOBAL NORMAL AND DIRECT NORMAL") 21 FORMAT(1X,"GLOBAL TILT") 22 FORMAT(1X,"GLOBAL HORIZONTAL") STOP END
Table of Contents
Return to RReDC Homepage ( http://www.nrel.gov/rredc )
