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

References


Downloadable working versions of this model can be found on the RReDC Solar Models: Bird Simple Spectral Model page. Versions are available in C code, compiled executable FORTRAN for the PC, and a spreadsheet implementation.
APPENDIX

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 )