      PROGRAM LITE
      REAL X0(4,30,7),BETA(4,30,2),PSI(4,30,2),RHO(4,30),XN(30,400,3),
     1     DAYF(3,400),DELF(400),F(4,30,30),WSOL(3,400),DSOL(4,30), 
     2     TAU(10),DW(5),LL(4,30,9),A(4,30),COEF(30,30),LS(3),Y(3), 
     3     YM(3),XY(3,4,2),LUM(4,30),AK2(10),X(3),COFB(30,1),FS(30),
     4     ALPHA(30),FR(30),AIK(30,400),LUMF(30,400),LUMI(30,400),
     5     XYW(5,2,2),XSUN(5,4,2),ANG(5),OHW(5),OHREF(5), 
     6     RING(5),OHTAU(5),RS(4),WW(5),WH(5),WCX(5),WCY(5),
     8     W(12),BW(12),
     9     PLTTMP(4,3)
C 
      INTEGER IDS(30),IOB(4,30,30),JOB(30,30,5),KSUN(4,30),NW(6),TZN, 
     1     KOB(4,30,30),KD(30,11,11),ICO(30),JCO(30,5),NM(30),NN(30), 
     2     ISUN(5),NUM(30),IJOB(26),KJOB(4,30,30),KKOB(4,30,30),
     3     IOH(5),KOH(5),ISTART(5),ITYPW(5),ISD(5), 
     4     IMG(5,20,20),ICNR(5,4) 
      INTEGER*4 ISEED
C 
      COMMON /BLK1/ X0,LL,XDP,YDP 
      COMMON /BLK2/ NW,PI,IOB,ICLD,ISEED,NBUNDL,A,KJOB
      COMMON /BLK3/ AK2,TAU,NC
      COMMON /BLK4/ BETAS,LS,AUX(3),GRNDL 
      COMMON /BLK5/ ICO,JCO,NSW 
      COMMON /BLK6/ XN,AIK,LUMF,LUMI,WSOL,DAYF,DELF,IMG 
     2              ,ICNR 
C     ISNGRD - DIMENSION OF IMG(), DIRECT SUN IMAGE ARRAY 
      ISNGRD=20 
      PI=ATAN(1.)*4.
      RAD=PI/180. 
      NBUNDL=2000 
      ISEED=123456 
C INID=1.. SIMPLIFIED INPUT FOR RECTANGULAR ROOMS 
C          IWRITE=1.. ILLUMINATION DATA FOR WORKING SURFACE ARE PRINTED 
C                =2.. DAYLIGHT FACTORS ARE PRINTED
C                =3.. ILLUMINATION AND DAYLIGHT FACTORS ARE PRINTED 
C INID=2.. FULL INPUT FOR ARBITRARY ROOMS,
C          IWRITE=1.. BRIEF OUTPUT (DATA FOR WORKING SURFACE ONLY)
C                =2.. DETAILED OUTPUT (DATA FOR ALL SURFACES) 
C IPLOT=0.. NO PLOT IS BEING GENERATED
C      =1.. ILLUMINATION LEVELS ARE PLOTTED (EQUALLY-SPACED LINES)
C      =2.. ILLUMINATION LEVELS ARE PLOTTED (EQUALLY-SPACED LEVELS) 
C      =3.. DAYLIGHT FACTORS ARE PLOTTED (EQUALLY-SPACED LEVELS)
C NUMIT= # OF ITERATIONS
C IOH=0.. NO OVERHANGS, 
C    =1.. OVERHANG ON TOP,
C    =2.. OVERHANGS ON BOTH SIDES,
C    =3.. OVERHANG ON BOTTOM ONLY,
C    =4.. OVERHANGS ON TOP AND BOTH SIDES,
C    =5.. OVERHANGS ON ALL FOUR SIDES.
C IWMOD=1.. INPUT IS SUN POSITION, SOLAR RADIATION, ETC., 
C        =2.. INPUT IS LATITUDE,LONGITUDE,TIME,ETC. 
C ITYPW=1.. CLEAR WINDOW, 
C      =2.. WINDOW WITH SHADING DEVICE, 
C      =3.. DIFFUSING WINDOW. 
C ISD  =1.. SHEER CURTAIN 
C 
C INPUT ROUTINES
      OPEN(UNIT=5,FILE='SLITE.IN',STATUS='OLD')
      OPEN(UNIT=7,FILE='SLITE.OUT',STATUS='NEW')
C    *     RECL=200)
      OPEN(UNIT=8,FILE='SLITE.PLT',STATUS='NEW')
C 
      CALL LOGO(6)
      READ(5,*) INID,IWMOD,NUMIT,IWRITE,IPLOT 
      IF(INID.NE.1) GOTO 31 
C 
C *** SIMPLIFIED INPUT - LIMITED TO RECTANGULAR ROOM WITH 
C *** RECTANGULAR WINDOWS, AND SIMPLE OUTSIDE OBSTRUCTIONS
C 
C ROOM SIZE AND LOCATION
C 
      READ(5,*) RW,RD,RH
      READ(5,*) RE,RO 
C 
C WINDOW DATA 
C NWINDO = NUMBER OF WINDOWS, 
C LOCW=1.. WINDOWS ARE IN THE FRONT WALL, 
C     =2.. WINDOWS ARE IN THE CEILING (SKYLIGHTS) 
      READ(5,*) NWINDO,LOCW 
C INITIALIZATION OF AUXILIARY VALUES
      IF (NWINDO.GE.1.AND.NWINDO.LE.5) GOTO 4010
      WRITE(7,8002) 
      STOP
4010  NSS=0 
      NC=0
      NSW=0 
      ND=0
      DO 4000 I=1,NWINDO
      IDS(I)=0
      OHW(I)=0. 
      RING(I)=0.
      OHREF(I)=0. 
      KOH(I)=0
 4000 OHTAU(I)=0. 
C DATA FOR EACH WINDOW
      DO 4100 I=1,NWINDO
      READ(5,*) WW(I),WH(I),WCX(I),WCY(I),NM(I),NN(I),RHO(1,I)
      READ(5,*) ITYPW(I),AK2(I),TAU(I),DW(I),IOH(I) 
      NC=NC+1/ITYPW(I)
      NSW=NSW+ITYPW(I)/2-ITYPW(I)/3 
      ND=ND+ITYPW(I)/3
      IF (LOCW.EQ.1) IDS(I)=2 
      ALPHA(I)=1. 
      ISD(I)=1
      IF(ITYPW(I).NE.2) GOTO 4050 
C SHEER CURTAIN 
      READ(5,*) ALPHA(I)
 4050 IF(ITYPW(I).EQ.3) ALPHA(I)=0. 
      IF(IOH(I).EQ.0) GOTO 4100 
C OVERHANG DATA 
      READ(5,*) OHW(I),RING(I),OHREF(I),KOH(I),OHTAU(I) 
      IF (KOH(I).EQ.0) OHTAU(I)=0.
 4100 CONTINUE
C WINDOW ORDERING (PUT IN SEQUENCE  CLEAR--SHADED--DIFFUSING) 
      IF(ITYPW(NWINDO).GE.ITYPW(1).OR.NWINDO.GT.2) GOTO 4110
      AAA=AK2(1)
      AK2(1)=AK2(2) 
      AK2(2)=AAA
      AAA=TAU(1)
      TAU(1)=TAU(2) 
      TAU(2)=AAA
      AAA=DW(1) 
      DW(1)=DW(2) 
      DW(2)=AAA 
      AAA=ALPHA(1)
      ALPHA(1)=ALPHA(2) 
      ALPHA(2)=AAA
      AAA=OHW(1)
      OHW(1)=OHW(2) 
      OHW(2)=AAA
      AAA=RING(1) 
      RING(1)=RING(2) 
      RING(2)=AAA 
      AAA=OHREF(1)
      OHREF(1)=OHREF(2) 
      OHREF(2)=AAA
      AAA=OHTAU(1)
      OHTAU(1)=OHTAU(2) 
      OHTAU(2)=AAA
      AAA=WW(1) 
      WW(1)=WW(2) 
      WW(2)=AAA 
      AAA=WH(1) 
      WH(1)=WH(2) 
      WH(2)=AAA 
      AAA=WCX(1)
      WCX(1)=WCX(2) 
      WCX(2)=AAA
      AAA=WCY(1)
      WCY(1)=WCY(2) 
      WCY(2)=AAA
      AAA=RHO(1,1)
      RHO(1,1)=RHO(1,2) 
      RHO(1,2)=AAA
      III=ITYPW(1)
      ITYPW(1)=ITYPW(2) 
      ITYPW(2)=III
      III=IOH(1)
      IOH(1)=IOH(2) 
      IOH(2)=III
      III=KOH(1)
      KOH(1)=KOH(2) 
      KOH(2)=III
      III=NM(1) 
      NM(1)=NM(2) 
      NM(2)=III 
      III=NN(1) 
      NN(1)=NN(2) 
      NN(2)=III 
C 
C REFLECTIVITIES AND DISCRETIZATION FOR WALLS 
C 
C IF LOCW=1, SEQUENCE IS WINDOWS-LEFT WALL-REAR WALL-RIGHT WALL-CEILING-
C FLOOR-FRONT WALL; IF LOCW=2, FRONT WALL AND CEILING ARE INTERCHANGED
C TO PLACE SURFACE WITH CUT-OUT (WINDOWS) LAST
 4110 DO 4120 I=1,6 
      IC=I+NWINDO-1 
      IF(I.EQ.1) IC=NWINDO+2*(4-LOCW) 
      IF(I.EQ.5) IC=NWINDO+2*(1+LOCW) 
 4120 READ(5,*) NM(IC),NN(IC),RHO(1,IC) 
C 
C WORKING SURFACE 
C 
      READ(5,*) NM(NWINDO+7),NN(NWINDO+7),WSE 
C 
C DESCRIPTION OF OUTSIDES 
C 
       READ(5,*) IOPP 
C 
C OPPOSING BUILDING 
C 
      IF(IOPP.EQ.1) READ(5,*) HO,WO,XO,DO,RHO(2,3)
C 
C SUBJECT BUILDING
C 
      READ(5,*) HS,WS,XS,RHO(2,1) 
C SET AUXILIARY VALUES
  
      NW(2)=IOPP+2
      NWC=5 
      NWO=1 
      NWS=1 
      IF (LOCW.EQ.1) NSS=1
      N0=NC+NSW 
      N01=NC+1
      N1=N0+ND
      N11=N1+1
      N2=N1+NWC 
      NT=N2+NWO 
      NT1=NT+1
      NTS=NT1 
      NSP1=NSS+1
C CALCULATE VALUES FOR CUT-OUT AND OBSTRUCTION IDENTIFIERS
      IMAX=NWINDO+6 
      IMAX1=IMAX+1
      ICO(IMAX)=NWINDO
      DO 4170 I=1,IMAX1 
      DO 4160 J=1,IMAX
 4160 IOB(1,I,J)=0
 4170 IOB(1,I,I)=-1 
      DO 4175 I=1,NWINDO
      JCO(IMAX,I)=I 
      IOB(1,I,IMAX)=-1
 4175 IOB(1,IMAX,I)=-1
      IOB(1,IMAX1,NWINDO+5)=-1
C 
      DO 4180 I=1,IMAX1 
      DO 4180 J=1,3 
      IOB(2,I,J)=0
      IOB(2,J,J)=-1 
 4180 KOB(2,I,J)=-1 
      DO 4190 I=1,5 
      DO 4190 J=2,3 
 4190 KOB(2,I+NWINDO,J)=0 
      KOB(2,NWINDO+5,2)=-1
      KOB(2,NWINDO+7,3)=0 
      RO=RO*RAD 
      SRO=SIN(RO) 
      CRO=COS(RO) 
C 
C DESCRIPTION OF WINDOWS
C 
      DO 4200 I=1,NWINDO
      X0(1,I,1)=(WCX(I)-WW(I)/2.)*SRO 
      X0(1,I,2)=-(WCX(I)-WW(I)/2.)*CRO
      X0(1,I,3)=RE+WCY(I)-WH(I)/2.
      X0(1,I,4)=WW(I) 
      X0(1,I,5)=0.
      X0(1,I,6)=WW(I) 
      X0(1,I,7)=WH(I) 
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=0.
      PSI(1,I,1)=RO-PI/2. 
      PSI(1,I,2)=0. 
      IF (LOCW.NE.2) GOTO 4200
      X0(1,I,1)=X0(1,I,1) - (WCY(I)-WH(I)/2.)*CRO 
      X0(1,I,2)=-WCX(I)*CRO-(WCY(I)-WH(I)/2.)*SRO+WW(I)*.5*SRO
      X0(1,I,3)=RE+RH 
      BETA(1,I,2)=PI/2. 
      PSI(1,I,2)=RO-PI
 4200 CONTINUE
C 
C FRONT WALL
C 
      I=NWINDO+2*(4-LOCW) 
      X0(1,I,1)=-RW/2.*SRO
      X0(1,I,2)=RW/2.*CRO 
      X0(1,I,3)=RE
      X0(1,I,4)=RW
      X0(1,I,5)=0.
      X0(1,I,6)=RW
      X0(1,I,7)=RH
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=0.
      PSI(1,I,1)=RO-PI/2. 
      PSI(1,I,2)=0. 
C 
C LEFT WALL 
C 
      I=NWINDO+1
      X0(1,I,1)=-RD*CRO-RW/2.*SRO 
      X0(1,I,2)=-RD*SRO+RW/2.*CRO 
      X0(1,I,3)=RE
      X0(1,I,4)=RD
      X0(1,I,5)=0.
      X0(1,I,6)=RD
      X0(1,I,7)=RH
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=0.
      PSI(1,I,1)=RO 
      PSI(1,I,2)=0. 
C 
C REAR WALL 
C 
      I=NWINDO+2
      X0(1,I,1)=-RD*CRO+RW/2.*SRO 
      X0(1,I,2)=-RD*SRO-RW/2.*CRO 
      X0(1,I,3)=RE
      X0(1,I,4)=RW
      X0(1,I,5)=0.
      X0(1,I,6)=RW
      X0(1,I,7)=RH
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=0.
      PSI(1,I,1)=RO+PI/2. 
      PSI(1,I,2)=0. 
C 
C RIGHT SIDE
C 
      I=NWINDO+3
      X0(1,I,1)=RW/2.*SRO 
      X0(1,I,2)=-RW/2.*CRO
      X0(1,I,3)=RE
      X0(1,I,4)=RD
      X0(1,I,5)=0.
      X0(1,I,6)=RD
      X0(1,I,7)=RH
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=0.
      PSI(1,I,1)=RO+PI
      PSI(1,I,2)=0. 
C 
C CEILING 
C 
      I=NWINDO+2*(1+LOCW) 
      X0(1,I,1)=RW/2.*SRO-RD*CRO
      X0(1,I,2)=-RW/2.*CRO-RD*SRO 
      X0(1,I,3)=RE+RH 
      X0(1,I,4)=RW
      X0(1,I,5)=0.
      X0(1,I,6)=RW
      X0(1,I,7)=RD
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=PI/2. 
      PSI(1,I,1)=RO+PI/2. 
      PSI(1,I,2)=RO 
C 
C FLOOR 
C 
      I=NWINDO+5
      X0(1,I,1)=-RW/2.*SRO-RD*CRO 
      X0(1,I,2)=RW/2.*CRO-RD*SRO
      X0(1,I,3)=RE
      X0(1,I,4)=RW
      X0(1,I,5)=0.
      X0(1,I,6)=RW
      X0(1,I,7)=RD
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=PI/2. 
      PSI(1,I,1)=RO-PI/2. 
      PSI(1,I,2)=RO 
C 
C WORKING SURFACE 
C 
      I=NWINDO+7
      X0(1,I,1)=-RW/2.*SRO-RD*CRO 
      X0(1,I,2)=RW/2.*CRO-RD*SRO
      X0(1,I,3)=RE+WSE
      X0(1,I,4)=RW
      X0(1,I,5)=0.
      X0(1,I,6)=RW
      X0(1,I,7)=RD
      BETA(1,I,1)=PI/2. 
      BETA(1,I,2)=PI/2. 
      PSI(1,I,1)=RO-PI/2. 
      PSI(1,I,2)=RO 
      RHO(1,I)=1. 
C 
C OUTSIDE WINDOW WALL 
C 
      X0(2,1,1)=(.5*WS+XS)*SRO
      X0(2,1,2)=-(.5*WS+XS)*CRO 
      X0(2,1,3)=0.
      X0(2,1,4)=WS
      X0(2,1,5)=0.
      X0(2,1,6)=WS
      X0(2,1,7)=HS
      BETA(2,1,1)=PI/2. 
      BETA(2,1,2)=0.
      PSI(2,1,1)=RO+PI/2. 
      PSI(2,1,2)=0. 
C 
C OUTSIDE GROUND
C 
      IF (IOPP.EQ.-1) GOTO 4210 
      HGR=WS
      IF(IOPP.EQ.1) HGR=DO
      X0(2,2,1)=-(0.75*WS-XS)*SRO 
      X0(2,2,2)= (0.75*WS-XS)*CRO 
      X0(2,2,3)=0.
      X0(2,2,4)=1.5*WS
      X0(2,2,5)=0.
      X0(2,2,6)=1.5*WS
      X0(2,2,7)=HGR 
      BETA(2,2,1)=PI/2. 
      BETA(2,2,2)=PI/2. 
      PSI(2,2,1)=RO-PI/2. 
      PSI(2,2,2)=RO 
      IF(IOPP.LE.0) GOTO 4210 
C 
C OPPOSING BUILDING 
C 
      X0(2,3,1)=DO*CRO-(.5*WO-XO)*SRO 
      X0(2,3,2)=DO*SRO+(.5*WO-XO)*CRO 
      X0(2,3,3)=0.
      X0(2,3,4)=WO
      X0(2,3,5)=0.
      X0(2,3,6)=WO
      X0(2,3,7)=HO
      BETA(2,3,1)=PI/2. 
      BETA(2,3,2)=0.
      PSI(2,3,1)=RO-PI/2. 
      PSI(2,3,2)=0. 
C 
C INPUT DATA PRINT-OUT
C 
 4210 RO=RO/RAD 
      WRITE(7,6300) RW,RD,RH,RO 
      IF(LOCW.EQ.1) WRITE(7,6302) NWINDO
      IF(LOCW.EQ.2) WRITE(7,6303) NWINDO
      DO 4220 I=1,NWINDO
      WRITE(7,6310) I,WW(I),WH(I),WCX(I),WCY(I),NM(I),NN(I) 
      IF(ITYPW(I).EQ.1) WRITE(7,6311) ITYPW(I)
      IF(ITYPW(I).NE.2) GOTO 4215 
      IF(ISD(I).EQ.1) WRITE(7,6312) ITYPW(I),ALPHA(I) 
 4215 IF(ITYPW(I).EQ.3) WRITE(7,6313) ITYPW(I)
      WRITE(7,6320) AK2(I),DW(I),TAU(I),RHO(1,I)
      IF(IOH(I).EQ.0) GOTO 4220 
      IF(IOH(I).EQ.1) WRITE(7,6031) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.2) WRITE(7,6032) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.3) WRITE(7,6033) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.4) WRITE(7,6034) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.5) WRITE(7,6035) OHW(I),RING(I),OHREF(I) 
      IF(KOH(I).EQ.0) WRITE(7,6036) 
      IF(KOH(I).EQ.1) WRITE(7,6037) OHTAU(I)
      IF(KOH(I).EQ.2) WRITE(7,6038) OHTAU(I)
 4220 CONTINUE
      NWP6=NWINDO+2*(4-LOCW)
      NWP4=NWINDO+2*(1+LOCW)
      WRITE(7,6331) NM(NWP6),NN(NWP6),RHO(1,NWP6) 
      WRITE(7,6332) NM(NWINDO+1),NN(NWINDO+1),RHO(1,NWINDO+1) 
      WRITE(7,6333) NM(NWINDO+2),NN(NWINDO+2),RHO(1,NWINDO+2) 
      WRITE(7,6334) NM(NWINDO+3),NN(NWINDO+3),RHO(1,NWINDO+3) 
      WRITE(7,6335) NM(NWP4),NN(NWP4),RHO(1,NWP4) 
      WRITE(7,6336) NM(NWINDO+5),NN(NWINDO+5),RHO(1,NWINDO+5) 
      WRITE(7,6337) NM(NWINDO+7),NN(NWINDO+7),WSE 
      WRITE(7,6340) WS,HS,XS,RHO(2,1) 
      IF(IOPP.EQ.1) WRITE(7,6350) WO,HO,XO,DO,RHO(2,3)
      GOTO 39 
C ***  COMPREHENSIVE INPUT FOR GENERAL ROOM;
C ***  NON-RECTANGULAR ROOMS, NON-RECTANGULAR SURFACES, OBSTRUCTIONS, ETC.
C 
C DATA FOR ROOM 
C 
   31 WRITE(7,6005) 
C NUMBER OF DIFFERENT TYPES OF SURFACES 
      READ(5,*) NC,NSW,ND,NWC,NWO,NWS,NSS 
      WRITE(7,6000)NC,NSW,ND,NWC,NWO,NWS,NSS
      N0=NC+NSW 
      N01=NC+1
      N1=N0+ND
      N11=N1+1
      N2=N1+NWC 
      NT=N2+NWO 
      NT1=NT+1
      NTS=NT+NWS
      NSP1=NSS+1
      IF (NWS.LE.3) GOTO 5170 
      WRITE(7,8007) 
      STOP
5170  IF (NSS.LE.2) GOTO 5180 
      WRITE(7,8008) 
      STOP
5180  IF (N1*2+NWC+NWO+NWS.LE.30) GOTO 5200 
      WRITE(7,8010) 
      STOP
C 
5200  DO 10 I=1,NT
C INTERNAL SURFACE DATA 
      READ(5,*)   (X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ=1, 
     1            2),RHO(1,I),NM(I),NN(I) 
      WRITE(7,6010)I,(X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ=1,
     1            2),RHO(1,I) 
      WRITE(7,6020)NM(I),NN(I)
      IF (NM(I)*NN(I).LE.400) GOTO 5205 
      WRITE(7,8011) I 
      STOP
5205  IF (RHO(1,I).GE.0..AND.RHO(1,I).LE.1.) GOTO 5210
      WRITE(7,8012) I 
      STOP
5210  DO 11 IZ=1,2
      BETA(1,I,IZ)=RAD*BETA(1,I,IZ) 
   11 PSI(1,I,IZ)=RAD*PSI(1,I,IZ) 
C 
C NUMBER AND NAMES OF CUT-OUTS
      IF(I.LE.N2) GO TO 13
      READ(5,*) ICO(I),(JCO(I,M),M=1,ICO(I))
      IF (ICO(I).LT.1) GOTO 13
      DO 5020 M=1,ICO(I)
      IF (JCO(I,M).GE.1.AND.JCO(I,M).LE.NTS) GOTO 5020
      WRITE(7,8014) I 
      STOP
5020  CONTINUE
C 
C NUMBER AND NAMES OF OBSTRUCTORS 
C 
 13   READ(5,*) (IOB(1,I,J),(JOB(I,J,K),K=1,IOB(1,I,J)),J=1,NT)
      IDS(I)=1
      IF(I.GT.N1) GO TO 10
C WINDOW DATA 
      READ(5,*) IS,AK2(I),TAU(I),DW(I),IOH(I) 
      IF (IS.LE.NSS) GOTO 5220
      WRITE(7,8020) I 
      STOP
5220  IF (TAU(I).GE.0..AND.TAU(I).LE.1.) GOTO 5230
      WRITE(7,8022) I 
      STOP
 5230 IF (AK2(I).GE.0..AND.AK2(I).LE.1.) GOTO 5240
      WRITE(7,8024) I 
      STOP
5240  IF(IOH(I).EQ.0) GOTO 16 
C OVERHANG DATA 
      READ(5,*) OHW(I),RING(I),OHREF(I),KOH(I),OHTAU(I) 
      IF (KOH(I).EQ.0) OHTAU(I)=0.
      IF (OHREF(I).GE.0..AND.OHREF(I).LE.1.) GOTO 5250
      WRITE(7,8034) I 
      STOP
5250  IF (OHTAU(I).GE.0..AND.OHTAU(I).LE.1.) GOTO 16
      WRITE(7,8036) I 
      STOP
   16 IDS(I)=IS+1-1/(1+IS)
      ALPHA(I)=1. 
      IF(I.GT.N0)ALPHA(I)=0.
      ISD(I)=1
      WRITE(7,6030)IS,AK2(I),TAU(I),DW(I) 
      IF(I.LE.NC.OR.I.GT.N0) GO TO 18 
C SHEER CURTAIN 
      READ(5,*) ALPHA(I)
      WRITE(7,6064) ALPHA(I)
   18 IF(IOH(I).EQ.0) GOTO 10 
      IF(IOH(I).EQ.1) WRITE(7,6031) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.2) WRITE(7,6032) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.3) WRITE(7,6033) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.4) WRITE(7,6034) OHW(I),RING(I),OHREF(I) 
      IF(IOH(I).EQ.5) WRITE(7,6035) OHW(I),RING(I),OHREF(I) 
      IF(KOH(I).EQ.0) WRITE(7,6036) 
      IF(KOH(I).EQ.1) WRITE(7,6037) OHTAU(I)
      IF(KOH(I).EQ.2) WRITE(7,6038) OHTAU(I)
   10 CONTINUE
C 
C DATA FOR WORKING SURFACES 
C 
      DO 35 I=NT1,NTS 
      READ(5,*) (X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ
     1             =1,2),NM(I),NN(I)
      WRITE(7,6070)I,(X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ=1 
     1             ,2),NM(I),NN(I)
      IF (NM(I)*NN(I).LE.400) GOTO 5270 
      WRITE(7,8011) I 
      STOP
5270  RHO(1,I)=1. 
      DO 36 IZ=1,2
      BETA(1,I,IZ)=RAD*BETA(1,I,IZ) 
   36 PSI(1,I,IZ)=RAD*PSI(1,I,IZ) 
      READ(5,*) (IOB(1,I,J),(JOB(I,J,K),K=1,IOB(1,I,J)),J=1,NT)
   35 CONTINUE
C 
C DATA FOR OUTSIDE ENCLOSURES 
C 
      NW(1)=0 
      IF(NSS.EQ.0)GO TO 39
C SCAN OVER ALL OUTSIDE ENCLOSURES
      DO 20 K=2,NSP1
C NUMBER OF SOLID WALLS IN OUTSIDE ENCLOSURE
      READ(5,*) NW(K) 
      K2=K-1
      WRITE(7,6040)K2,NW(K) 
      NTK=NW(K) 
      DO 20 I=1,NTK 
C SURFACE DATA
      READ (5,*)  (X0(K,I,IZ),IZ=1,7),(BETA(K,I,IZ),PSI(K,I,IZ),IZ=1, 
     1             2),RHO(K,I)
      WRITE(7,6010)I,(X0(K,I,IZ),IZ=1,7),(BETA(K,I,IZ),PSI(K,I,IZ),IZ=1,
     1             2),RHO(K,I)
      IF (RHO(K,I).GE.0..AND.RHO(K,I).LE.1.) GOTO 5290
      WRITE(7,8005) I,K 
      STOP
5290  DO 17 IZ=1,2
      BETA(K,I,IZ)=RAD*BETA(K,I,IZ) 
   17 PSI(K,I,IZ)=RAD*PSI(K,I,IZ) 
C NUMBER OF OBSTRUCTORS 
 20   READ(5,*) (IOB(K,I,J),(KJOB(K,I,J),L=1,IOB(K,I,J)),J=1,NTK)
C NUMBER OF OBSTRUCTORS BETWEEN INSIDE AND OUTSIDE SURFACES 
C 
      DO 30 K=2,NSP1
      NTK=NW(K) 
      DO 30 I=1,NTS 
 30   READ(5,*) (KOB(K,I,J),(KKOB(K,I,L),L=1,KOB(K,I,J)),J=1,NTK)
   39 WSOR=PSI(1,NT1,1) 
      NW(1)=0 
C 
C ILLUMINATION DATA 
C 
      MFIRST=1
      MLAST=1 
      MSTEP=1 
      NTSTEP=1
      IDAY=1
      GOTO(5400,34,5405)IWMOD 
C SPECIFIC SOLAR DATA ARE SUPPLIED
5400  READ(5,*) SOLC1,ETSOL,SKYC1,ETSKY,ICLD,RHOGR
      WRITE(7,6050)SOLC1,ETSOL,SKYC1,ETSKY,ICLD,RHOGR 
      SKYC=ETSKY*SKYC1
      IF(ICLD.EQ.1.OR.ICLD.EQ.3) GO TO 32 
C CLEAR DAY 
      READ(5,*) HSUND,PSID
      GOTO 5410 
C GEOGRAPHICAL DATA ARE SUPPLIED; SCANS THROUGH MONTHS, HOURS 
   34 READ(5,*) ALAT,ALONG,TZN,ALT,ICLD,RHOGR 
      READ(5,*) MFIRST,MLAST,MSTEP,IDAY,TFIRST,TLAST,DELTAT 
      WRITE(7,6100) ALAT,ALONG,TZN,ALT,ICLD,RHOGR 
      WRITE(7,6110) MFIRST,MLAST,MSTEP,IDAY,TFIRST,TLAST,DELTAT 
      IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 42
      READ(5,*) (W(M),BW(M),M=MFIRST,MLAST,MSTEP) 
      WRITE(7,6120) 
      DO 6125 M=MFIRST,MLAST,MSTEP
      IF (W(M).GE.0..AND.W(M).LE.10.) GOTO 6124 
      WRITE(7,8030) M 
      STOP
6124  IF (BW(M).GE.0..AND.BW(M).LE..5) GOTO 6125
      WRITE(7,8032) M 
      STOP
 6125 WRITE(7,6121) M,W(M),BW(M)
   42 NTSTEP=(TLAST-TFIRST)/DELTAT+1.1
      GOTO 32 
C 
C SUN POSITION( ALTITUDE AND AZIMUTH) DATA ARE SUPPLIED 
C 
5405  READ(5,*) HSUND,PSID,ALT,ICLD,RHOGR 
      WRITE(7,6105) HSUND,PSID,ALT,ICLD,RHOGR 
      IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 5410
      READ(5,*) W(1),BW(1)
5410  HSUN=HSUND*RAD
      PSIS=PSID*RAD 
      IISUN=1 
32    IF(INID.EQ.1.AND.IOPP.NE.-1)RHO(2,2)=RHOGR
      IF (RHOGR.GE.0..AND.RHOGR.LE.1.) GOTO 5300
      WRITE (6,8026)
      STOP
C 
C EVALUATION OF DIRECTION COSINES 
C 
 5300 DO 40 K=1,NSP1
      NTK=NTS 
      IF(K.GE.2) NTK=NW(K)
      DO 40 I=1,NTK 
      SBX=SIN(BETA(K,I,1))
      CBX=COS(BETA(K,I,1))
      SBY=SIN(BETA(K,I,2))
      CBY=COS(BETA(K,I,2))
      SPX=SIN(PSI(K,I,1)) 
      CPX=COS(PSI(K,I,1)) 
      SPY=SIN(PSI(K,I,2)) 
      CPY=COS(PSI(K,I,2)) 
      LL(K,I,1)=SBX*CPX 
      LL(K,I,2)=SBX*SPX 
      LL(K,I,3)=CBX 
      LL(K,I,4)=SBY*CPY 
      LL(K,I,5)=SBY*SPY 
      LL(K,I,6)=CBY 
      LL(K,I,7)=SBX*SPX*CBY-CBX*SBY*SPY 
      LL(K,I,8)=CBX*SBY*CPY-SBX*CPX*CBY 
      LL(K,I,9)=SBX*SBY*(SPY*CPX-CPY*SPX) 
   40 CONTINUE
C 
C CALCULATION OF NODAL COORDINATES
C 
      DO 1100 I=1,NTS 
      K=0 
      NM1=NM(I) 
      NN1=NN(I) 
      AA=X0(1,I,7)/(NM1*NN1)
      BB=X0(1,I,4)-X0(1,I,6)+X0(1,I,5)
      DO 1000 N=1,NN1 
      ETA=(NN1+.5-N)/NN1
      CC=X0(1,I,4)-BB*ETA 
      DO 1000 M=1,NM1 
      XSI=(M-.5)/NM1
      XIP=X0(1,I,5)*ETA+CC*XSI
      YIP=X0(1,I,7)*ETA 
      DO 1010 L=1,3 
 1010 XN(I,K+1,L)=X0(1,I,L)+LL(1,I,L)*XIP+LL(1,I,L+3)*YIP 
C CHECK WHETHER A NODAL POINT FALLS ONTO A CUT-OUT
      IF(I.LE.N2.OR.I.GT.NT) GOTO 1020
      IICO=ICO(I) 
      DO 1030 IC=1,IICO 
      J=JCO(I,IC) 
      XJP=0.
      YJP=0.
      DO 1040 L=1,3 
      XYP=XN(I,K+1,L)-X0(1,J,L) 
      XJP=XJP+XYP*LL(1,J,L) 
 1040 YJP=YJP+XYP*LL(1,J,L+3) 
      ETAJ=YJP/X0(1,J,7)
      IF(ABS(ETAJ-.5).GT..5) GOTO 1030
      XSIJ=(XJP-X0(1,J,5)*ETAJ)/(X0(1,J,4)-(X0(1,J,4)-X0(1,J,6)+
     1          X0(1,J,5))*ETAJ)
      IF(ABS(XSIJ-.5).GT..5) GOTO 1030
      GOTO 1000 
 1030 CONTINUE
 1020 K=K+1 
      AIK(I,K)=AA*CC
 1000 CONTINUE
C NUMBER OF NODES ON SURFACE
 1100 NUM(I)=K
C 
C CALCULATION OF COORDINATES FOR OUTSIDE SURFACE OF WINDOW OPENING
C (STORE THESE AUXILIARY SURFACES BEYOND WORKING SURFACE, I.E.@ NTS+I)
C 
      DO 45 I=1,N1
      DO 43 J=1,3 
      X0(1,I+NTS,J)=X0(1,I,J)-(DW(I)+1.E-3)*LL(1,I,J+6) 
      DO 43 K=1,3 
      JJ=J+3*(K-1)
   43 LL(1,I+NTS,JJ)=LL(1,I,JJ) 
      DO 45 J=4,7 
   45 X0(1,I+NTS,J)=X0(1,I,J) 
C 
C CALCULATION OF OUTSIDE SURFACE AREAS
C 
      IF(NSS.EQ.0) GOTO 49
      DO 44 K=2,NSP1
      NTK=NW(K) 
      DO 44 I=1,NTK 
      LUM(K,I)=0
   44 A(K,I)=(X0(K,I,4)+X0(K,I,6)-X0(K,I,5))*X0(K,I,7)/2. 
C 
C GEOMETRIC DESCRIPTION OF WINDOW OVERHANGS 
C (STORED BEYOND OUTSIDE ENCLOSURES, I.E. 1+NSS+1)
C 
   49 NOH=NSP1+1
      IOHK=0
      DO 490 I=1,N1 
      IF(IOH(I).EQ.0) GOTO 490
      IN=I+NTS
      IOHK=IOHK+1 
      ISTART(I)=IOHK
      X2Y=X0(1,IN,5)/X0(1,IN,7) 
      RTX2Y=SQRT(1.+X2Y*X2Y)
      X31Y=(X0(1,IN,4)-X0(1,IN,6))/X0(1,IN,7) 
      RTX31Y=SQRT(1.+X31Y*X31Y) 
      XLT=X0(1,IN,5)-RING(I)*(RTX2Y-X2Y)
      YLT=X0(1,IN,7)+RING(I)
      XTR=X0(1,IN,6)+RING(I)*(RTX31Y-X31Y)
      YTR=YLT 
      IF(IOH(I).EQ.2.OR.IOH(I).EQ.3) GOTO 425 
C TOP OVERHANG
      DO 410 L=1,3
  410 X0(NOH,IOHK,L)=X0(1,IN,L)+LL(1,I,L)*XLT+LL(1,I,L+3)*YLT 
      X0(NOH,IOHK,4)=OHW(I) 
      X0(NOH,IOHK,5)=0. 
      X0(NOH,IOHK,6)=OHW(I) 
      X0(NOH,IOHK,7)=XTR-XLT
C 
      DO 420 L=1,3
      LL(NOH,IOHK,L)=-LL(1,I,L+6) 
      LL(NOH,IOHK,L+3)=LL(1,I,L)
  420 LL(NOH,IOHK,L+6)=-LL(1,I,L+3) 
      RHO(NOH,IOHK)=OHREF(I)
      IF(IOH(I).LE.1) GOTO 490
C 
      IOHK=IOHK+1 
  425 XBL=-RING(I)*(RTX2Y+X2Y)
      YBL=-RING(I)
      XRB=X0(1,IN,4)+RING(I)*(RTX31Y+X31Y)
      YRB=YBL 
      IF(IOH(I).EQ.3) GOTO 455
C LEFT AND RIGHT OVERHANGS
      DO 430  L=1,3 
      X0(NOH,IOHK,L)=X0(1,IN,L)+LL(1,I,L)*XBL+LL(1,I,L+3)*YBL 
  430 X0(NOH,IOHK+1,L)=X0(1,IN,L)+LL(1,I,L)*XTR+LL(1,I,L+3)*YTR 
      X0(NOH,IOHK,4)=OHW(I) 
      X0(NOH,IOHK,5)=0. 
      X0(NOH,IOHK,6)=OHW(I) 
      X0(NOH,IOHK,7)=(X0(1,IN,7)+2.*RING(I))*RTX2Y
      X0(NOH,IOHK+1,4)=OHW(I) 
      X0(NOH,IOHK+1,5)=0. 
      X0(NOH,IOHK+1,6)=OHW(I) 
      X0(NOH,IOHK+1,7)=(X0(1,IN,7)+2.*RING(I))*RTX31Y 
C 
      DO 440 L=1,3
      LL(NOH,IOHK,L)=-LL(1,I,L+6) 
      LL(NOH,IOHK,L+3)=(LL(1,I,L)*X2Y+LL(1,I,L+3))/RTX2Y
      LL(NOH,IOHK+1,L)=-LL(1,I,L+6) 
  440 LL(NOH,IOHK+1,L+3)=(LL(1,I,L)*X31Y-LL(1,I,L+3))/RTX31Y
      DO 450 L=1,3
      L1=L+1-3*(L/3)
      L2=L+2-3*(L/2)
      LL(NOH,IOHK,L+6)=LL(NOH,IOHK,L1)*LL(NOH,IOHK,L2+3)
     1                -LL(NOH,IOHK,L2)*LL(NOH,IOHK,L1+3)
  450 LL(NOH,IOHK+1,L+6)=LL(NOH,IOHK+1,L1)*LL(NOH,IOHK+1,L2+3)
     1                  -LL(NOH,IOHK+1,L2)*LL(NOH,IOHK+1,L1+3)
      RHO(NOH,IOHK)=OHREF(I)
      IOHK=IOHK+1 
      RHO(NOH,IOHK)=OHREF(I)
      IF(IOH(I).EQ.2.OR.IOH(I).EQ.4) GOTO 490 
C 
C BOTTOM OVERHANG 
      IOHK=IOHK+1 
  455 DO 460 L=1,3
  460 X0(NOH,IOHK,L)=X0(1,IN,L)+LL(1,I,L)*XRB+LL(1,I,L+3)*YRB 
      X0(NOH,IOHK,4)=OHW(I) 
      X0(NOH,IOHK,5)=0. 
      X0(NOH,IOHK,6)=OHW(I) 
      X0(NOH,IOHK,7)=XRB-XBL
C 
      DO 470 L=1,3
      LL(NOH,IOHK,L)=-LL(1,I,L+6) 
      LL(NOH,IOHK,L+3)=-LL(1,I,L) 
  470 LL(NOH,IOHK,L+6)=LL(1,I,L+3)
      RHO(NOH,IOHK)=OHREF(I)
  490 CONTINUE
      ISTART(N1+1)=IOHK+1 
      ZENL0=0 
C 
C START CYCLE FOR MONTHS AND HOURS (IF IWMOD=2) 
C 
      DO 3400 MONTH=MFIRST,MLAST,MSTEP
      DO 3400 ITIME=1,NTSTEP
      IF(IWMOD.NE.2) GOTO 494 
      TIME=TFIRST + DELTAT*(ITIME-1)
      WRITE(7,6130) MONTH,TIME
C 
C ZENITH LUMINANCE
C 
494   DO 495 K=2,NSP1 
      NK=NW(K)
      DO 495 I=1,NK 
495   DSOL(K,I)=0.
      DO 496 I=1,NWS
      KMAX=NUM(I+NT)
      DO 496 K=1,KMAX 
496   WSOL(I,K)=0.
      DO 497 I=1,NTS
      KMAX=NUM(I) 
      DO 497 K=1,KMAX 
497   LUMI(I,K)=0.
      GO TO(50,60,50,60),ICLD 
C ICLD=1.. OVERCAST SKY, SKY LUMINANCE DISTRIBUTION FROM C.I.E. 
C ICLD=2.. CLEAR SKY WITH DIRECT SUN INCLUDED.
C ICLD=3.. UNIFORM-LUMIN@ANCE SKY, ZENITH VALUE FROM OVERCAST C.I.E.
C ICLD=4.. CLEAR SKY ONLY 
   50 AK=5.*(3.-2.*RHOGR)/7.*(1/ICLD) 
      AUX(1)=1./(1.+AK) 
      AUX(2)=AUX(1)*AK
      SOLC=0. 
      GOTO(57,58,59)IWMOD 
C LUMINANCE FROM INSOLATION DATA
57    ZENL=(1.+AK)*SKYC/((1.+.667*AK))
      GRNDL=SKYC*RHOGR/ZENL 
      GO TO 100 
C LUMINANCE FROM GEOGRAPHICAL DATA
C DETERMINE HSUN=SUN ALTITUDE ANGLE 
58    CALL SOLAR(ALAT,ALONG,TZN,MONTH,TIME,IISUN,HSUN,PSIS) 
59    ZENL=92.94*(.123+8.6*SIN(HSUN))*PI 
      SKYC=(1.+.667*AK)/(1.+AK)*ZENL
      GRNDL=SKYC*RHOGR/ZENL 
      IF(MONTH-MFIRST+ITIME.EQ.1) GOTO 100
C 
C SET UP LUMINANCES FOR OTHER TIMES IF NOT CLEAR SKY
C   (FOR NON-CLEAR SKY CALCUL@ATIONS NEED TO BE DONE ONLY ONCE, AS
C   GEOMETRY AND SKY DISTRIBUTION REMAIN UNCHANGED; ONLY THE VALUE
C   OF ZENL=ZENITH LUMINANCE CHANGES) 
C 
C RATIO OF NEW AND OLD ZENITH LUMINANCES
      ZR=ZENL/ZENL0 
      ZENL0=ZENL
C OUTSIDE LUMINANCES
      IF (NSS.LE.0) GOTO 2405 
      DO 2400 K=1,NSS 
      NKT=NW(K+1) 
      DO 2400 I=1,NKT 
 2400 LUM(K+1,I)=ZR*LUM(K+1,I)
2405  IF(IOHK.EQ.0) GOTO 2415 
C OVERHANG LUMINANCES 
      DO 2410 II=1,IOHK 
 2410 LUM(NOH,II)=ZR*LUM(NOH,II)
C INSIDE AND WORKING SURFACES 
 2415 DO 2420 I=1,NTS 
      KMAX=NUM(I) 
      DO 2420 K=1,KMAX
 2420 LUMF(I,K)=ZR*LUMF(I,K)
      GOTO 1405 
C 
C ESTABLISH SUN POSITION FOR THE GEOGRAPHICAL DATA GIVEN
   60 IF(IWMOD.EQ.2)CALL SOLAR(ALAT,ALONG,TZN,MONTH,TIME, 
     1                IISUN,HSUN,PSIS)
      BETAS=PI/2.-HSUN
      HSUND=HSUN/RAD
      PSID=PSIS/RAD 
      IF(IWMOD.EQ.2)WRITE(7,6140) HSUND,PSID
C INTEGRATE NORMALIZED LUMINANCE OVER SKY 
      P=SKYSUM(BETAS) 
      AUX(3)=.27385*(.91+10*EXP(-3*BETAS)+.45*COS(BETAS)**2)
      IF (IWMOD.GT.1) GOTO 67 
C IF SOLAR DATA KNOWN CALCULATE ZENL FROM SKY INTEGRAL
      SOLC=ETSOL*SOLC1
      ES=SOLC/COS(BETAS)
      ZENL=SKYC/P 
      TERM=SOLC/(SOLC+SKYC) 
C IF TOO LITTLE DIFFUSE LIGHT, FAULTY INPUT DATA HAVE BEEN GIVEN
      IF(TERM.LT..9) GOTO 33
      WRITE(7,8000) 
      STOP
C ROUTINE TO CALCULATE ZENL AS FUNCTION OF TURBIDITY, VAPOR & SOL.ALT.
   67 CALL ZENLCD(P,W(MONTH),BW(MONTH),HSUN,ZENL) 
C ROUTINE TO CALCULATE DIRECT SUN FROM SIMILAR DATA 
      IF(ICLD.EQ.2)SOLC=DIR(MONTH,IDAY,BW(MONTH),W(MONTH),ALT,HSUN) 
      IF(ICLD.EQ.4)SOLC=0.
C USE SKY INTEGRAL AND ZENL TO GET DIFFUSE SKY ILLUMINATION 
      SKYC=P*ZENL 
   33 II=0
      ES=SOLC/COS(BETAS)
      GRNDL=(SKYC+SOLC)*RHOGR/ZENL
      IF(NSS.LE.0) GOTO 70
C INITIALIZATION OF SUNNY-NODES IDENTIFIERS FOR OUTSIDE SURFACES
      DO 38 K=2,NSP1
      NK=NW(K)
      DO 38 J=1,NK
      II=II+1 
      KSUN(K,J)=0 
      DSOL(K,J)=0.
      DO 37 M=1,11
      DO 37 N=1,11
      KD(II,M,N)=0. 
   37 CONTINUE
   38 CONTINUE
C 
C DIRECTION COSINES FOR VECTOR TOWARDS SUN
C 
70    LS(1)=SIN(BETAS)*COS(PSIS)
      LS(2)=SIN(BETAS)*SIN(PSIS)
      LS(3)=COS(BETAS)
C 
C CALCULATION OF DIRECT SUNSHINE ILLUMINATION INSIDE ROOM 
C 
      WRITE(7,6067) 
      DO 1270 IW=1,N1 
      ISUN(IW)=0
      CBW=LS(1)*LL(1,IW,7)+LS(2)*LL(1,IW,8)+LS(3)*LL(1,IW,9)
C DOES SUN HIT WINDOW?
      IF(CBW.GE.-1.E-4) GOTO 1270 
C ANY OBSTRUCTION BETWEEN WINDOW AND SUN? 
      KK=IDS(IW)
      IF(KK.LE.0) GOTO 1240 
      NWK=NW(KK)
      DO 1250 L=1,3 
 1250 X(L)=X0(1,IW,L) 
      DO 1230 J=1,NWK 
      IF(IOB(KK,1,J).LT.0) GOTO 1230
      CALL SECT(X,LS,KK,J,KUT)
      IF(KUT.EQ.1) GOTO 1270
 1230 CONTINUE
 1240 ISUN(IW)=1
      WRITE(7,6068) IW,ALPHA(IW)
 1270 CONTINUE
C 
      IF(N0.EQ.0) GOTO 1300 
      DO 1200 I=1,NTS 
C COULD SUN POSSIBLY SHINE ON THIS SURFACE? 
      CBS=LS(1)*LL(1,I,7)+LS(2)*LL(1,I,8)+LS(3)*LL(1,I,9) 
      IF(CBS.LE.1.E-4) GOTO 1200
      KMAX=NUM(I) 
      DO 1210 IW=1,N0 
      IF(ISUN(IW).EQ.0) GOTO 1210 
      IF(IOB(1,I,IW).LT.0) GOTO 1210
      CBW=LS(1)*LL(1,IW,7)+LS(2)*LL(1,IW,8)+LS(3)*LL(1,IW,9)
C ATTENUATED AMOUNT OF DIRECT SUN LIGHT TO BE ADDED TO LUMINANCE
      WS=RHO(1,I)*ES*CBS*ALPHA(IW)*TAUB(IW,-CBW)
      DO 1260 K=1,KMAX
      DO 1220 L=1,3 
 1220 X(L)=XN(I,K,L)
      CALL SECT(X,LS,1,IW,KUT)
C DOES RAY FROM NODE TO SUN PASS THROUGH WINDOW?
      IF(KUT.NE.1) GOTO 1260
C BACK SIDE OF WINDOW ALSO? 
      CALL SECT(X,LS,1,IW+NTS,KUT)
      IF(KUT.NE.1) GOTO 1260
C HOW ABOUT INTERNAL OBSTRUCTIONS?
      IF(IOB(1,I,IW).EQ.0) GOTO 1255
      IIOB=IOB(1,I,IW)
      DO 1245 IO=1,IIOB 
      JO=JOB(I,IW,IO) 
      CALL SECT(X,LS,1,JO,KUT)
      IF(KUT.EQ.1) GOTO 1260
 1245 CONTINUE
C ANY OVERHANG IN THE WAY?
 1255 IF(IOH(IW).EQ.0) GOTO 1257
      II1=ISTART(IW)
      II2=ISTART(IW+1)-1
      DO 1256 II=II1,II2
      CALL SECT(X,LS,NOH,II,KUT)
      IF(KUT.NE.1) GOTO 1256
C IF SO, IS OVERHANG SEMI-TRANSPARENT? (IN THAT CASE ATTENUATE MORE)
      IF(KOH(IW).NE.1) GOTO 1260
      WS=WS*OHTAU(IW) 
      GOTO 1257 
 1256 CONTINUE
C ANY EXTERNAL OBSTRUCTION IN THE WAY?
1257  KK=IDS(IW)
      IF (KK.LE.0) GOTO 1275
      NWK=NW(KK)
      DO 1258 J=1,NWK 
      IF (KOB(KK,I,J).LT.0) GOTO 1258 
      CALL SECT(X,LS,KK,J,KUT)
      IF (KUT.EQ.1) GOTO 1260 
 1258 CONTINUE
C DIRECT SUN IS INITIAL VALUE OF LUMINANCE
 1275 LUMI(I,K)=WS
 1260 CONTINUE
 1210 CONTINUE
 1200 CONTINUE
C STORE THESE VALUES ALSO IN WSOL FOR PLOTTING ROUTINE
      DO 1280 I=NT1,NTS 
      KMAX=NUM(I) 
      DO 1280 K=1,KMAX
1280  WSOL(I-NT,K)=LUMI(I,K)
C 
C CALCULATION OF DIRECT SUNSHINE ON CURTAINED AND 
C DIFFUSE WINDOWS 
C 
      IF(N1.EQ.NC) GOTO 1350
 1300 DO 1310 I=N01,N1
C DOES WINDOW GET SUNSHINE? 
      IF(ISUN(I).EQ.0) GOTO 1310
      CBW=LS(1)*LL(1,I,7)+LS(2)*LL(1,I,8)+LS(3)*LL(1,I,9) 
C CALCULATE ATTENUATED CONTRIBUTION TO LUMINANCE (NOTE CBW<0) 
      WS=-ES*CBW*(1.-ALPHA(I))*TAUB(I,-CBW) 
      KMAX=NUM(I) 
C ANY OVERHANG IN THE WAY FOR THIS NODE?
      DO 1320 K=1,KMAX
      TT=1. 
      IF(IOH(I).EQ.0) GOTO 1320 
      DO 1330 L=1,3 
 1330 X(L)=XN(I,K,L)
      II1=ISTART(I) 
      II2=ISTART(I+1)-1 
      DO 1340 II=II1,II2
      CALL SECT(X,LS,NOH,II,KUT)
      IF(KUT.NE.1) GOTO 1340
      TT=0. 
C IS OVERHANG SEMI-TRANSPARENT? 
      IF(KOH(I).EQ.1) TT=OHTAU(I) 
      GOTO 1320 
 1340 CONTINUE
 1320 LUMI(I,K)=WS*TT 
 1310 CONTINUE
C 
C CALCULATION OF DIRECT  SUNSHINE ILLUMINATION FOR OUTSIDE ENCLOSURE SURFACES 
C 
 1350 IF(NSS.LE.0)GO TO 210 
      II=0
      DO 90 K=2,NSP1
      NK=NW(K)
      DO 90 J=1,NK
      II=II+1 
      FR(J)=0.
C ARE SUNNY AREAS POSSIBLE ON THIS SURFACE? 
      CBJ=LS(1)*LL(K,J,7)+LS(2)*LL(K,J,8)+LS(3)*LL(K,J,9) 
      IF(CBJ.LE.1.E-4) GO TO 90 
C 
C BREAK UP SURFACE INTO 11 BY 11 NODES AND CHECK FOR SUN (SEND RAYS TO SUN) 
      DO 95 M=1,11
      ETAM=(M-1.)/10. 
      DO 95 N=1,11
      ETAN=(N-1.)/10. 
      XPP=X0(K,J,4)-(X0(K,J,4)-X0(K,J,6)+X0(K,J,5))*ETAM
      XP=X0(K,J,5)*ETAM+XPP*ETAN
      YP=X0(K,J,7)*ETAM 
      DO 94 L=1,3 
   94 X(L)=X0(K,J,L)+LL(K,J,L)*XP+LL(K,J,L+3)*YP
C ANYTHING BLOCKING THE SUN?
      DO 93 JJ=1,NK 
      IF(J.EQ.JJ) GO TO 93
      CBL=LS(1)*LL(K,JJ,7)+LS(2)*LL(K,JJ,8)+LS(3)*LL(K,JJ,9)
      IF(CBL.GE.-1.E-4)GO TO 93 
      CALL SECT (X,LS,K,JJ,KUT) 
      IF(KUT.EQ.1) GO TO 95 
   93 CONTINUE
      KD(II,M,N)=1
C SUNNY FRACTION OF SURFACE 
      FR(J)=FR(J)+C(M,11)*C(N,11)*XPP 
   95 CONTINUE
C 
      IF(FR(J).LE.1.E-4)GO TO 90
      FR(J)=FR(J)*2/(X0(K,J,4)+X0(K,J,6)-X0(K,J,5)) 
      FS(II)=(1.-FR(J))/FR(J) 
      KSUN(K,J)=1 
C CONTRIBUTION OF DIRECT SUN TO OUTSIDE SURFACE LUMINANCE 
      DSOL(K,J)=FR(J)*CBJ*ES*RHO(K,J) 
   90 CONTINUE
C 
C CALCULATION OF F I-J IN OUTSIDE ENCLOSURES
C 
  100 ZENL0=ZENL
      IF(NSS.LE.0) GOTO 210 
      DO 400 K=2,NSP1 
      NMW=NW(K) 
      CALL CONFAC(K,NMW,F)
  400 CONTINUE
C 
C LUMINANCES FOR EXTERNAL REFLECTIONS (SOLUTION OF SIMULTANEOUS EQS.
C            BY MATRIX INVERSION) 
C 
      DO 220 K=2,NSP1 
      NKT=NW(K) 
      DO 190 I=1,NKT
      DO 180 J=1,NKT
      DIJ=(I/J)*(J/I) 
  180 COEF(I,J)=DIJ/RHO(K,J)-F(K,I,J) 
  190 COFB(I,1)=ZENL*F(K,I,NKT+1)+DSOL(K,I)/RHO(K,I)
      IF (NKT.EQ.1) GOTO 230
      CALL MATINV(COEF,NKT,COFB,1,DETERM) 
      DO 200 I=1,NKT
  200 LUM(K,I)=COFB(I,1)
      GOTO 220
230   LUM(K,1)=COFB(1,1)/COEF(1,1)
220   CONTINUE
C 
C CALCULATION OF OVERHANG LUMINANCES
C (USING MONTE CARLO, ASSUMING NO INTERFERENCE BETWEEN OVERHANGS) 
C 
  210 NBOH=NBUNDL/5 
      DO 1495 IW=1,N1 
      IF(IOH(IW).EQ.0) GOTO 1495
      K=IDS(IW) 
C NUMBER OF OVERHANG SURFACES FOR THIS WINDOW 
      II1=ISTART(IW)
      II2=ISTART(IW+1)-1
      DO 1450 II=II1,II2
      SUM=0.
C SEND OUT BUNDLES
      DO 1455 NBUN=1,NBOH 
C ************************ RANDOM NUMBER GENERATOR ********************** 
      DO 112,INDEX=1,4
112   RS(INDEX)=RANDOM(ISEED)
C CHOOSE DIRECTION
      SINB=SQRT(RS(1))
      COSB=SQRT(1.-RS(1)) 
      THET=2.*PI*RS(2)
      COST=COS(THET)
      SINT=SIN(THET)
C CHOOSE EMISSION POINT 
      DO 1460 L=1,3 
      X(L)=X0(NOH,II,L)+LL(NOH,II,L)*X0(NOH,II,4)*RS(3) 
     1                 +LL(NOH,II,L+3)*X0(NOH,II,7)*RS(4) 
      Y(L)=(LL(NOH,II,L)*COST+LL(NOH,II,L+3)*SINT)*SINB 
     1                       +LL(NOH,II,L+6)*COSB 
 1460 YM(L)=-Y(L) 
      IF(K.EQ.0) GOTO 1466
C DOES IT HIT OUTSIDE SURFACE OF WINDOW?
      NWK=NW(K) 
      CALL SECT(X,YM,1,IW+NTS,KUT)
      IF(KUT.NE.1) GOTO 1462
      SUM=SUM+RHO(1,IW)/RHO(K,1)*LUM(K,1) 
      GOTO 1455 
C ANY OUTSIDE SURFACE?
 1462 IF (NSS.LT.1) GOTO 1466 
      DO 1465 IR=1,NWK
      CALL SECT(X,Y,K,IR,KUT) 
      IF(KUT.NE.1) GOTO 1465
      SUM=SUM+LUM(K,IR) 
      GOTO 1455 
 1465 CONTINUE
C HITS THE SKY, (OR GROUND IF Y POINTS BELOW HORIZ) 
 1466 SUM=SUM+RLPZ(ICLD,Y)*ZENL 
 1455 CONTINUE
      LUM(NOH,II)=OHREF(IW)*SUM/NBOH
      IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 1470
C DOES DIRECT SUN HIT INSIDE OVERHANG SURFACES? 
      IF(ISUN(IW).EQ.0) GOTO 1470 
      CBS=LS(1)*LL(NOH,II,7)+LS(2)*LL(NOH,II,8)+LS(3)*LL(NOH,II,9)
      IF(CBS.LT.1.E-4) GOTO 1470
      IF(K.EQ.0) GOTO 1476
      DO 1475 IR=1,NWK
      CALL SECT(X,LS,K,IR,KUT)
      IF(KUT.EQ.1) GOTO 1470
 1475 CONTINUE
 1476 LUM(NOH,II)=LUM(NOH,II)+OHREF(IW)*ES*CBS
 1470 IF(KOH(IW).NE.2) GOTO 1450
C FOR TRANSLUCENT OVERHANG ALSO NEED LUMINANCE CONTRIBUTION FROM OUTSIDE
C   OF OVERHANG SURFACE; ROUTINE SAME AS THE INSIDE (SEE ABOVE) 
      SUM=0.
      DO 1480 NBUN=1,NBOH 
C ************************ RANDOM NUMBER GENERATOR ********************** 
      DO 113 INDEX=1,4
113   RS(INDEX)=RANDOM(ISEED)
      SINB=SQRT(RS(1))
      COSB=SQRT(1.-RS(1)) 
      THET=2.*PI*RS(2)
      COST=COS(THET)
      SINT=SIN(THET)
      DO 1485 L=1,3 
      X(L)=X0(NOH,II,L)+LL(NOH,II,L)*X0(NOH,II,4)*RS(3) 
     1                 +LL(NOH,II,L+3)*X0(NOH,II,7)*RS(4) 
 1485 Y(L)=-(LL(NOH,II,L)*COST+LL(NOH,II,L+3)*SINT)*SINB
     1                       -LL(NOH,II,L+6)*COSB 
      IF(K.EQ.0) GOTO 1491
      DO 1490 IR=1,NWK
      CALL SECT(X,Y,K,IR,KUT) 
      IF(KUT.NE.1) GOTO 1490
      SUM=SUM+LUM(K,IR) 
      GOTO 1480 
 1490 CONTINUE
 1491 SUM=SUM+RLPZ(ICLD,Y)*ZENL 
 1480 CONTINUE
      LUM(NOH,II)=LUM(NOH,II)+OHTAU(IW)*SUM/NBOH
      IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 1450
      IF(ISUN(IW).EQ.0) GOTO 1450 
      CBS=-LS(1)*LL(NOH,II,7)-LS(2)*LL(NOH,II,8)-LS(3)*LL(NOH,II,9) 
      IF(CBS.LT.1.E-4) GOTO 1450
      IF(K.EQ.0) GOTO 1451
      DO 1505 IR=1,NWK
      CALL SECT(X,LS,K,IR,KUT)
      IF(KUT.EQ.1) GOTO 1450
 1505 CONTINUE
 1451 LUM(NOH,II)=LUM(NOH,II)+OHTAU(IW)*ES*CBS
 1450 CONTINUE
 1495 CONTINUE
C 
C OUTSIDE ILLUMINATION ON CURTAINED AND DIFFUSE WINDOWS 
C (USING MIXED METHOD.. INTEGRATION AND MONTE CARLO)
C 
      IF(NC.EQ.N1) GOTO 1600
      DO 1500 I=N01,N1
      IO=I+NTS
      XY1=(X0(1,IO,6)-X0(1,IO,5))/X0(1,IO,4)
      ID1=1 
      IF(ABS(XY1-1.).LT.1.E-3) ID1=0
      AX=XY1*ID1
      HI=0. 
      K=IDS(I)
C INTEGRATE OVER ALL DIRECTIONS 
      DO 1550 IETA=1,21 
      ETA=(IETA-1.)/20. 
      SETA=SQRT(1.-ETA*ETA) 
      TT=C(IETA,21)*ETA*TAUB(I,ETA) 
      DO 1550 IMU=1,21
      AMU=(IMU-1.)/20.
      SWY=SETA*COS(2.*PI*AMU) 
      SWZ=SETA*SIN(2.*PI*AMU) 
C EMISSION LOCATION FIXED BY RANDOM NUMBERS 
C ************************ RANDOM NUMBER GENERATOR ********************** 
      DO 114,INDEX=1,4
114   RS(INDEX)=RANDOM(ISEED)
      YY=(1.-SQRT(1.-(1.-XY1*XY1)*RS(1)))/(1.-AX)+(1-ID1)*RS(1) 
      XX=X0(1,IO,5)*YY+X0(1,IO,4)*(1.-(1.-XY1)*YY)*RS(2)
      DO 1520 L=1,3 
      X(L)=X0(1,IO,L)+LL(1,I,L)*XX+LL(1,I,L+3)*YY*X0(1,IO,7)
 1520 Y(L)=-ETA*LL(1,I,L+6)+SWY*LL(1,I,L+3)+SWZ*LL(1,I,L) 
C ANY OVERHANG IN THE WAY?
      IF(IOH(I).EQ.0) GOTO 1545 
      II1=ISTART(I) 
      II2=ISTART(I+1)-1 
      DO 1515 II=II1,II2
      CALL SECT(X,Y,NOH,II,KUT) 
      IF(KUT.NE.1) GOTO 1515
      SS=LUM(NOH,II)
      GOTO 1550 
 1515 CONTINUE
C ANY OUTSIDE SURFACE IN THE WAY? 
 1545 IF(K.EQ.0) GOTO 1530
      NWK=NW(K) 
      DO 1540 J=1,NWK 
      IF(KOB(K,I,J)) 1540,1535,1525 
 1525 JR=KKOB(K,I,J)
      SS=LUM(K,JR)
      CALL SECT(X,Y,K,JR,KUT) 
      IF(KUT.EQ.1) GOTO 1550
 1535 CALL SECT(X,Y,K,J,KUT)
      IF(KUT.NE.1) GOTO 1540
      SS=LUM(K,J) 
      GOTO 1550 
 1540 CONTINUE
C HITS THE SKY, (OR GROUND IF Y POINTS BELOW HORIZ) 
 1530 SS=ZENL*RLPZ(ICLD,Y)
 1550 HI=HI+TT*SS*C(IMU,21) 
      KMAX=NUM(I) 
      DO 1500 K=1,KMAX
 1500 LUMI(I,K)=LUMI(I,K)+2.*(1.-ALPHA(I)*(1/ISD(I)))*HI
C 
C OUTSIDE ILLUMINATION ON INSIDE NODES
C (INCLUDES INSIDE LUMINANCE OF WINDOWS, IN PARTICULAR DIFFUSING WINDOWS) 
C 
 1600 DO 1750 I=1,NTS 
      KMAXI=NUM(I)
      FRAC=RHO(1,I)/PI
      DO 1750 J=1,N1
C CAN I SEE J?
      IF(IOB(1,I,J).LT.0) GOTO 1750 
      KMAXJ=NUM(J)
      K=IDS(J)
      IF(K.GT.0.AND.J.LE.NC) GOTO 1610
      IF(K.GT.0.AND.ISD(J).EQ.1) GOTO 1610
C AUXILIARY DATA FOR COMPUTER EFFICIENCY, FOR CASE OF NO OUTSIDE ILLUMINATION 
C   POSSIBLE ON THAT NODE 
      K=4 
      LUM(4,1)=0. 
      K1=30 
      K2=30 
      KOB(4,I,30)=-1
      GOTO 1620 
 1610 K1=1
      K2=NW(K)
 1620 IF(IOB(1,I,J).GT.0) IIOB=IOB(1,I,J) 
C THE FOLLOWING ALGORITHMS ARE FOR CASES WITH & WITHOUT INTERNAL OBSTRUCTORS
      DO 1660 KI=1,KMAXI
C START CALCULATIONS FOR NODE KI ON SURFACE I 
      DO 1630 L=1,3 
 1630 X(L)=XN(I,KI,L) 
      HIT=0.
      DO 1642 KJ=1,KMAXJ
C DIRECTION, COSINES*S, AND DISTANCE-SQUARED TO NODE KJ ON J (WINDOW) 
      SCB1=0. 
      SCB2=0. 
      SSQ=0.
      DO 1650 L=1,3 
      Y(L)=XN(J,KJ,L)-X(L)
      SCB1=SCB1+Y(L)*LL(1,I,L+6)
      SCB2=SCB2-Y(L)*LL(1,J,L+6)
 1650 SSQ=SSQ+Y(L)*Y(L) 
      IF(IOB(1,I,J).EQ.0) GOTO 1680 
C ANY INTERNAL OBSTRUCTORS IN THE WAY?
      DO 1710 IO=1,IIOB 
      JO=JOB(I,J,IO)
      CALL SECT(X,Y,1,JO,KUT) 
      IF(KUT.EQ.1) GOTO 1642
 1710 CONTINUE
 1680 S=SQRT(SSQ) 
      DO 1675 L=1,3 
 1675 Y(L)=Y(L)/S 
      ICOS1=1.+SCB1/S 
      ICOS2=1.+SCB2/S 
C CONFIGURATION FACTOR (CORRECTED BY FUDGE FACTOR FOR CLOSE NODES)
      FIJ=ICOS1*ICOS2*SCB1*SCB2/(SSQ*SSQ)*AIK(J,KJ) 
      FIJ=FIJ/(1.+.6*FIJ*FIJ) 
C ASSUME OUTSIDE-WINDOW-WALL LUMINANCE FOR WINDOW FRAME 
      HI=LUM(K,1) 
      CB2=SCB2/S
C DOES RAY GET TO THE OUTSIDE?
      CALL SECT(X,Y,1,J+NTS,KUT)
      SDM=1.
      IF(IOB(1,I,J).LE.0.AND.KUT.NE.1) GOTO 1640
      IF(KUT.NE.1) GOTO 1635
C DOES RAY THROUGH WINDOW HIT ANY OUTSIDE SURFACE?
      DO 1670 JKL=K1,K2 
      JK=JKL
      IF(KOB(K,I,JK)) 1670,1625,1615
 1615 JKK=KKOB(K,I,JK)
      CALL SECT(X,Y,K,JKK,KUT)
      IF(KUT.NE.1) GOTO 1625
      JK=JKK
      GOTO 1605 
 1625 CALL SECT(X,Y,K,JK,KUT) 
      IF(KUT.NE.1) GOTO 1670
 1605 HI=LUM(K,JK)
      IF(ICLD.NE.2.OR.KSUN(K,JK).NE.1) GOTO 1635
C IS POINT ON OUTSIDE SURFACE IN SUN OR IN SHADE? 
C INTERPOLATE OUTSIDE LUMINANCES IF CLOSE TO SUN/SHADE DIVIDE 
      II=JK 
      KM1=K-1 
      DO 1655 K3=1,KM1
 1655 II=II+NW(K3)
      TL=LUM(K,JK)-DSOL(K,JK) 
      TH=LUM(K,JK)+FS(II)*DSOL(K,JK)
      ETAM=YDP/X0(K,JK,7) 
      XPP=X0(K,JK,4)-(X0(K,JK,4)-X0(K,JK,6)+X0(K,JK,5))*ETAM
      ETAN=(XDP-X0(K,JK,5)*ETAM)/XPP
      TY=10*ETAM+1
      LY=TY 
      TX=10*ETAN+1
      LX=TX 
      F1=KD(II,LY,LX)+(TY-LY)*(KD(II,LY+1,LX)-KD(II,LY,LX)) 
      F2=KD(II,LY,LX+1)+(TY-LY)*(KD(II,LY+1,LX+1)-KD(II,LY,LX+1)) 
      PART=(F1+(TX-LX)*(F2-F1)) 
      HI=TL+PART*(TH-TL)
      GOTO 1635 
 1670 CONTINUE
C HITS SKY, (OR GROUND IF Y POINTS BELOW HORIZ) 
      HI=ZENL*RLPZ(ICLD,Y)
C ANY OVERHANG IN THE WAY?
 1635 IF(IOH(J).EQ.0) GOTO 1640 
      II1=ISTART(J) 
      II2=ISTART(J+1)-1 
      DO 1645 II=II1,II2
      CALL SECT(X,Y,NOH,II,KUT) 
      IF(KUT.NE.1) GOTO 1645
      HI=LUM(NOH,II)+OHTAU(J)*HI*(2-KOH(J)) 
      GOTO 1640 
 1645 CONTINUE
C LUMINANCE CONTRIBUTION FROM OUTSIDE (HI) AND WINDOW ITSELF
 1640 HIT=HIT+(HI*TAUB(J,CB2)*ALPHA(J)*(1/ISD(J))*SDM+LUMI(J,KJ))*FIJ 
 1642 CONTINUE
 1660 LUMI(I,KI)=LUMI(I,KI)+FRAC*HIT
1750  CONTINUE
C 
C LUMINANCES FOR INTERNAL REFLECTIONS 
C 
C KEEP DISTINCTION BETWEEN TOTAL (LUMF) AND EXTERNAL CONTRIBUTION (LUMI)
      DO 1760 I=1,NT
      KMAX=NUM(I) 
      DO 1760 K=1,KMAX
 1760 LUMF(I,K)=LUMI(I,K) 
C 
C GO THROUGH DESIRED NUMBER OF ITERATIONS 
      DO 2200 ITER=1,NUMIT
      DO 1900 I=1,NT
C IF WINDOW-REFLECTANCE IS SMALL, ITS CONTRIBUTION IS NEGLECTED 
C    (FOR COMPUTER EFFICIENCY)
C THE ALGORITHMS BELOW ARE BASICALLY THE SAME AS FOR OUTSIDE/WINDOW LUMINANCE 
C    EFFECTS COMMENTED ON ABOVE 
      IF(RHO(1,I).LE..15) GOTO 1900 
      FRAC=RHO(1,I)/PI
      KMAXI=NUM(I)
      DO 1805 KI=1,KMAXI
 1805 DELF(KI)=0. 
      DO 1890 J=N11,NT
      KMAXJ=NUM(J)
      IF(IOB(1,I,J)) 1890,1855,1850 
C 
 1850 IIOB=IOB(1,I,J) 
 1855 DO 1880 KI=1,KMAXI
      DO 1860 L=1,3 
 1860 X(L)=XN(I,KI,L) 
      DO 1880 KJ=1,KMAXJ
      SCB1=0. 
      SCB2=0. 
      SSQ=0.
      DO 1870 L=1,3 
      Y(L)=XN(J,KJ,L)-X(L)
      SCB1=SCB1+Y(L)*LL(1,I,L+6)
      SCB2=SCB2-Y(L)*LL(1,J,L+6)
 1870 SSQ=SSQ+Y(L)*Y(L) 
      IF(IOB(1,I,J).EQ.0) GOTO 1877 
      DO 1875 IO=1,IIOB 
      CALL SECT(X,Y,1,JOB(I,J,IO),KUT)
      IF(KUT.EQ.1) GOTO 1880
 1875 CONTINUE
 1877 ICOS1=1.+SCB1/(1.+SSQ)
      ICOS2=1.+SCB2/(1.+SSQ)
      FIJ=SCB1*SCB2/(SSQ*SSQ)*AIK(J,KJ)*ICOS1*ICOS2 
      FIJ=FIJ/(1.+.6*FIJ*FIJ) 
      DELF(KI)=DELF(KI)+FRAC*FIJ*LUMF(J,KJ) 
 1880 CONTINUE
 1890 CONTINUE
      DO 1895 KI=1,KMAXI
C ITERATION FINISHED FOR THIS SURFACE, IMPROVE VALUES FOR LUMF
 1895 LUMF(I,KI)=LUMI(I,KI)+DELF(KI)
 1900 CONTINUE
 2200 CONTINUE
C 
C CALCULATION OF ILLUMINATION AND DAYLIGHT FACTORS FOR WORKING SURFACE NODES
C 
C ALGORITHMS ARE THE SAME AS FOR INTERNAL REFLECTIONS, BUT NO ITERATION,
C BECAUSE OF THEIR IMPORTANCE, CONFIGURATION FACTORS FOR CLOSE NODES
C ARE CALCULATED MORE ACCURATELY
      DO 1995 I=NT1,NTS 
      KMAXI=NUM(I)
      DO 1905 KI=1,KMAXI
 1905 DELF(KI)=0. 
      DO 1990 J=N11,NT
      KMAXJ=NUM(J)
      IF(IOB(1,I,J)) 1990,1910,1950 
 1950 IIOB=IOB(1,I,J) 
 1910 DO 1980 KI=1,KMAXI
      DO 1960 L=1,3 
 1960 X(L)=XN(I,KI,L) 
      DO 1980 KJ=1,KMAXJ
      SCB1=0. 
      SCB2=0. 
      SSQ=0.
      DO 1970 L=1,3 
      Y(L)=XN(J,KJ,L)-X(L)
      SCB1=SCB1+Y(L)*LL(1,I,L+6)
      SCB2=SCB2-Y(L)*LL(1,J,L+6)
 1970 SSQ=SSQ+Y(L)*Y(L) 
      IF(IOB(1,I,J).EQ.0) GOTO 1977 
      DO 1975 IO=1,IIOB 
      CALL SECT(X,Y,1,JOB(I,J,IO),KUT)
      IF(KUT.EQ.1) GOTO 1980
 1975 CONTINUE
 1977 ICOS1=SCB1/(1.+SSQ)+1.
      ICOS2=SCB2/(1.+SSQ)+1.
      FIJ=ICOS1*ICOS2*SCB1*SCB2/(SSQ*SSQ)*AIK(J,KJ) 
      IF(AIK(J,KJ)/SSQ.LT.1.) GOTO 2160 
      KKJ=KJ
 2100 L1=(KKJ-1)/NM(J)
      L2=KKJ-L1*NM(J) 
      IF(J.LE.N2) GOTO 2140 
      ETA=1.-(2.*L1+1.)/(2.*NN(J))
      XSI=(2.*L2-1.)/(2.*NM(J)) 
      CC=X0(1,J,4)-(X0(1,J,4)+X0(1,J,5)-X0(1,J,6))*ETA
      XJP=X0(1,J,5)*ETA+CC*XSI
      YJP=X0(1,J,7)*ETA 
      XYZ=0.
      DO 2150 L=1,3 
 2150 XYZ=XYZ+X0(1,J,L)+LL(1,J,L)*XJP+LL(1,J,L+3)*YJP-XN(J,KJ,L)
      IF(ABS(XYZ).LT.1E-4) GOTO 2140
      KKJ=KKJ+1 
      GOTO 2100 
 2140 FIJ=0.
      DO 2110 J1=1,3
      ETA=1.-(2.*L1-(.5-J1)/1.5)/(2.*NN(J)) 
      CC=X0(1,J,4)-(X0(1,J,4)+X0(1,J,5)-X0(1,J,6))*ETA
      AP=X0(1,J,7)*CC/(9.*NN(J)*NM(J))
      DO 2110 J2=1,3
      XSI=(2.*L2+(.5-J2)/1.5)/(2.*NM(J))
      SCB1=0. 
      SCB2=0. 
      SSQ=0.
      XJP=X0(1,J,5)*ETA+CC*XSI
      YJP=X0(1,J,7)*ETA 
      DO 2120 L=1,3 
      Y(L)=X0(1,J,L)+LL(1,J,L)*XJP+LL(1,J,L+3)*YJP-X(L) 
      SCB1=SCB1+Y(L)*LL(1,I,L+6)
      SCB2=SCB2-Y(L)*LL(1,J,L+6)
 2120 SSQ=SSQ+Y(L)*Y(L) 
      ICOS1=SCB1/(1.+SSQ)+1.
      ICOS2=SCB2/(1.+SSQ)+1.
 2110 FIJ=FIJ+ICOS1*ICOS2*SCB1*SCB2/(SSQ*SSQ)*AP
 2160 DELF(KI)=DELF(KI)+FIJ*LUMF(J,KJ)
 1980 CONTINUE
 1990 CONTINUE
C ADD THE INTERNAL-REFLECTION CONTRIBUTION TO THE WORKING SURFACE 
C   SURFACE ILLUMINATION (WHICH SO FAR CONTAINS DIRECT SUN AND
C   CONTRIBUTION FROM OUTSIDE (INCLUDING SKY, EXTERNAL REFLECTIONS, 
C   AND DIFFUSE WINDOW LIGHT, SUCH AS TRANSMISSION THROUGH CURTAINS,
C   SLATTED BLINDS, DIFFUSING WINDOWS)
      DO 1995 KI=1,KMAXI
      LUMF(I,KI)=LUMI(I,KI)+DELF(KI)/PI 
 1995 DAYF(I-NT,KI)=100.*(LUMF(I,KI)-WSOL(I-NT,KI))/SKYC
C 
C PRINT WALL NODAL DATA AND LUMINANCES
C 
 1405 WRITE(7,6250) SOLC,SKYC,ZENL
      IF(IWRITE.EQ.1) GOTO 1440 
C 
C PRINTING OF LIGHT EXCHANGE FACTORS IN OUTSIDE ENCLOSURES
C 
      DO 1410 K=2,NSP1
      NMW=NW(K) 
      NM1=NMW+1 
      K4=K-1
      WRITE(7,6560)K4 
      WRITE(7,6570)(J,J=1,NM1)
      WRITE(7,6572) 
      DO 1410 I=1,NMW 
 1410 WRITE(7,6580) I,(F(K,I,J),J=1,NM1)
C 
C PRINTING OF OUTSIDE LUMINANCES
C 
      DO 1420 K=1,NSS 
      NKT=NW(K+1) 
      WRITE(7,7010)K
      WRITE(7,7020)(I,I=1,NKT)
 1420 WRITE(7,6200)(LUM(K+1,I),I=1,NKT) 
      DO 1415 IW=1,N1 
      IF(IOH(IW).EQ.0) GOTO 1415
      K=IDS(IW) 
      II1=ISTART(IW)
      II2=ISTART(IW+1)-1
      WRITE(7,6041) IW
      IF(IOH(IW).EQ.1) WRITE(7,6042)(LUM(NOH,II),II=II1,II2)
      IF(IOH(IW).EQ.2) WRITE(7,6043)(LUM(NOH,II),II=II1,II2)
      IF(IOH(IW).EQ.3) WRITE(7,6044)(LUM(NOH,II),II=II1,II2)
      IF(IOH(IW).EQ.4) WRITE(7,6045)(LUM(NOH,II),II=II1,II2)
      IF(IOH(IW).EQ.5) WRITE(7,6046)(LUM(NOH,II),II=II1,II2)
 1415 CONTINUE
      WRITE(7,6710) 
      DO 1400 I=1,NT
      KP=NM(I)
      KS=1+(KP-1)/12
      L=NN(I) 
      DO 1400 L1=1,L,KS 
      KL=1+(L1-1)*KP
      KR=L1*KP
      IF(KL.GT.NUM(I)) GOTO 1440
      IF(KR.GT.NUM(I)) KR=NUM(I)
      WRITE(7,6720)(K,K=KL,KR,KS) 
      WRITE(7,6730)I,(XN(I,K,1),K=KL,KR,KS) 
      WRITE(7,6740)(XN(I,K,2),K=KL,KR,KS) 
      WRITE(7,6750)(XN(I,K,3),K=KL,KR,KS) 
      WRITE(7,6760)(AIK(I,K),K=KL,KR,KS)
 1400 WRITE(7,6770)(LUMF(I,K),K=KL,KR,KS) 
C 
C PRINT WORKING SURFACE NODAL DATA, ILLUMINATION, DAYLIGHT FACTOR 
C 
 1440 WRITE(7,6715) 
      WRITE(7,6716)
      IF(INID.EQ.1)WRITE(7,6717)
      IF(INID.EQ.1)NT1=NWINDO+7 
      DO 1430 I=NT1,NTS 
      KP=NM(I)
      KS=1+(KP-1)/12
      L=NN(I) 
      DO 1430 L1=1,L,KS 
      KL=1+(L1-1)*KP
      KR=L1*KP
      WRITE(7,6720)(K,K=KL,KR,KS) 
      IF(INID.EQ.1)WRITE(7,6730)I,(ABS(XN(I,K,1)),K=KL,KR,KS) 
      IF(INID.EQ.2)WRITE(7,6730)I,(XN(I,K,1),K=KL,KR,KS)
      WRITE(7,6740)(XN(I,K,2),K=KL,KR,KS) 
      WRITE(7,6750)(XN(I,K,3),K=KL,KR,KS) 
      WRITE(7,6773) (LUMI(I,K),K=KL,KR,KS)
      WRITE(7,6774) (DELF(K)/PI,K=KL,KR,KS) 
      WRITE(7,6775)(LUMF(I,K),K=KL,KR,KS) 
      IF(INID.EQ.1)WRITE(7,6780)(DAYF(1,K),K=KL,KR,KS)
 1430 IF(INID.EQ.2)WRITE(7,6780)(DAYF(I-NT,K),K=KL,KR,KS) 
      IF(IPLOT.EQ.0) GOTO 3400
C 
C SET UP DATA FOR THE PLOTTING ROUTINE; CHOOSE LOCAL ORIGIN OF FIRST
C WORKING SURFACE AS OVERALL ORIGIN FOR PLOT, AND STORE COORDINATES 
C FOR CORNERPOINTS@ (MUST BE RECTANGULAR) 
      XY(1,1,1)=0.
      XY(1,1,2)=0.
      XY(1,2,1)=X0(1,NT1,5) 
      XY(1,2,2)=X0(1,NT1,7) 
      XY(1,3,1)=X0(1,NT1,6) 
      XY(1,3,2)=X0(1,NT1,7) 
      XY(1,4,1)=X0(1,NT1,4) 
      XY(1,4,2)=0.
C 
C SUBTRACT DIRECT SUN FROM ILLUMINATION FOR PLOTTING PURPOSES 
C AND RE-ORDER NODES
      KMAX=NUM(NT1) 
      DO 3000 K=1,KMAX
      K1=(K-1)/NM(NT1)
      K2=K+(NN(NT1)-2*K1-1)*NM(NT1) 
 3000 LUMI(1,K)=LUMF(NT1,K2)-WSOL(1,K2) 
C 
C FIND ENDPOINTS FOR WINDOW PROJECTION ON FIRST WORKING SURFACE 
      DO 2900 I=1,N1
      DO 2800 L=1,3 
      PLTTMP(1,L)=X0(1,I,L) 
 2800 PLTTMP(2,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,6)+LL(1,I,L+3)*X0(1,I,7) 
      DO 2900 I4=1,2
      DO 2900 I2=1,2
      XYW(I,I4,I2)=0. 
      DO 2900 L=1,3 
      L2=L+3*(I2-1) 
 2900 XYW(I,I4,I2)=XYW(I,I4,I2)+(PLTTMP(I4,L)-X0(1,NT1,L))*LL(1,NT1,L2) 
C 
      IF(NWS.EQ.1) GOTO 3200
C DO THE SAME AS ABOVE FOR ANY ADDITIONAL WORKING SURFACES
C     WITH RESPECT TO LOCAL AXES OF FIRST WORK SFC
      DO 3100 IW=2,NWS
      DO 3010 L=1,3 
      PLTTMP(1,L)=X0(1,NT+IW,L) 
      PLTTMP(2,L)=PLTTMP(1,L)+LL(1,NT+IW,L)*X0(1,NT+IW,5)+
     1       LL(1,NT+IW,L+3)*X0(1,NT+IW,7)
      PLTTMP(3,L)=PLTTMP(1,L)+LL(1,NT+IW,L)*X0(1,NT+IW,6)+
     1       LL(1,NT+IW,L+3)*X0(1,NT+IW,7)
 3010 PLTTMP(4,L)=PLTTMP(1,L)+LL(1,NT+IW,L)*X0(1,NT+IW,4) 
C 
      DO 3020 I4=1,4
      DO 3020 I2=1,2
      XY(IW,I4,I2)=0. 
      DO 3020 L1=1,3
      LI2=L1+(I2-1)*3 
 3020 XY(IW,I4,I2)=XY(IW,I4,I2)+(PLTTMP(I4,L1)-X0(1,NT1,L1))
     1              *LL(1,NT1,LI2)
C 
      KMAX=NUM(NT+IW) 
      DO 3100 K=1,KMAX
      K1=(K-1)/NM(NT+IW)
      K2=K+(NN(NT+IW)-2*K1-1)*NM(NT+IW) 
 3100 LUMI(IW,K)=LUMF(NT+IW,K2)-WSOL(IW,K2) 
C 
C THESE PLOTTING DATA ARE WRITTEN INTO A DATA FILE PLDAT WHICH
C UPON JOB COMPLETION IS PERMANENTLY PLACED ON PSS-DISK 
C 
 3200 WRITE(8,8200) NT,NWS,ICLD,N0,N1,IPLOT 
      WRITE(8,8200)(NM(I),NN(I),I=NT1,NTS)
      WRITE(8,8100) WSOR,SOLC,SKYC
      IF(ICLD.EQ.2.OR.ICLD.EQ.4) WRITE(8,8100) HSUND,PSID 
      WRITE(8,8100) (((XYW(I,I4,I2),I2=1,2),I4=1,2),I=1,N1) 
      DO 3150 I=1,NWS 
      KMAX=NUM(NT+I)
      WRITE(8,8100)(LUMI(I,K),K=1,KMAX) 
 3150 WRITE(8,8100)((XY(I,I4,I2),I2=1,2),I4=1,4)
      IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 3400
C 
C ADDITIONAL PLOTTING DATA ARE GENERATED IF DIRECT SUN PENETRATES 
C INTO ENCLOSURE, VZ. THE IMAGE OF THE SUNNY AREA ON THE WORKING SURFACE
C 
      DO 3280 I=1,N0
      DO 3205 IP=1,4
      XSUN(I,IP,1)=0. 
 3205 XSUN(I,IP,2)=0. 
      ANG(I)=0. 
      IF(ISUN(I).EQ.0) GOTO 3280
      DO 3210 L=1,3 
      PLTTMP(1,L)=X0(1,I,L) 
      PLTTMP(2,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,5)+LL(1,I,L+3) 
     1      *X0(1,I,7)
      PLTTMP(3,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,6)+LL(1,I,L+3) 
     1      *X0(1,I,7)
 3210 PLTTMP(4,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,4) 
      IDENT=0 
      DO 3250 IP=1,4
      DO 3220 L=1,3 
      X(L)=PLTTMP(IP,L) 
 3220 Y(L)=X(L)-(DW(I)+1E-3)*LL(1,I,L+6)
      IF(IDENT.NE.0) GOTO 3240
      CALL SECT(X,LS,1,I+NTS,KUT) 
      IF(KUT.EQ.1) IDENT=IP 
 3240 CALL SECT(X,LS,1,NT1,KUT) 
      XSUN(I,IP,1)=XDP
      XSUN(I,IP,2)=YDP
      CALL SECT(Y,LS,1,NT1,KUT) 
 3250 CONTINUE
      IAP2=IDENT+2
      IF (IAP2.GT.4) IAP2=IAP2-4
      IIII=IDENT+1
      IF (IIII.GT.4) IIII=IIII-4
      JJJJ=IDENT+3
      IF (JJJJ.GT.4) JJJJ=JJJJ-4
C 
      DX=0. 
      DY=0. 
      DO 3260 L=1,3 
      DX=DX+LS(L)*LL(1,NT1,L) 
 3260 DY=DY+LS(L)*LL(1,NT1,L+3) 
      ANG(I)=-180./PI*ATAN(DY/(DX+1E-8))
      DO 3300 IP=1,ISNGRD 
      DO 3300 L=1,ISNGRD
3300  IMG(I,IP,L)=0 
      ICNR(I,1)=IDENT 
      ICNR(I,2)=IIII
      ICNR(I,3)=JJJJ
      ICNR(I,4)=IAP2
      DO 3302 L=1,3 
      DO 3302 K=1,4 
      PLTTMP(K,L)=XSUN(I,ICNR(I,K),1)*LL(1,NT1,L) 
     2 +XSUN(I,ICNR(I,K),2)*LL(1,NT1,L+3) 
     2 +X0(1,NT1,L) 
3302  CONTINUE
      DO 3390 M=1,ISNGRD
      DO 3390 N=1,ISNGRD
      DO 3304 L=1,3 
      X(L)=(PLTTMP(3,L)-PLTTMP(1,L))*(M-1)/(ISNGRD-1) 
     1     + PLTTMP(1,L)
      DUMMY=(PLTTMP(4,L)-PLTTMP(2,L))*(M-1)/(ISNGRD-1)
     1     + PLTTMP(2,L)
      X(L)=(DUMMY-X(L))*(N-1)/(ISNGRD-1) + X(L) 
3304  CONTINUE
      CALL SECT(X,LS,1,I+NTS,KUT) 
      IF (KUT.NE.1) GOTO 3390 
      DO 3320 IWS=1,NWS 
      IWSN=NT+IWS 
      IIOB=IOB(1,IWSN,I)
      IF (IIOB.LT.1) GO TO 3320 
      DO 3310 K=1,IIOB
      JO=JOB(IWSN,I,K)
      CALL SECT(X,LS,1,JO,KUT)
      IF (KUT.NE.1) GO TO 3310
      IF ((YDP+X0(1,JO,3)).LT.X0(1,IWSN,3)) GOTO 3310 
      IF (I.NE.4) GOTO 3390 
      GOTO 3390 
3310  CONTINUE
3320  CONTINUE
      IF (IOH(I).EQ.0 .OR. KOH(I).EQ.1) GO TO 3370
      IOHK1=ISTART(I) 
      IOHK2=ISTART(I+1)-1 
      DO 3350 IOHK=IOHK1,IOHK2
      CALL SECT(X,LS,NOH,IOHK,KUT)
      IF (KUT.EQ.1) GO TO 3390
3350  CONTINUE
3370  K=IDS(I)
      IF (K.LT.1) GO TO 3385
      K2=NW(K)
      DO 3380 JK=1,K2 
      CALL SECT(X,LS,K,JK,KUT)
      IF (KUT.EQ.1) GO TO 3390
3380  CONTINUE
3385  IMG(I,M,N)=1
3390  CONTINUE
 3280 CONTINUE
      WRITE(8,8100)(((XSUN(I,IP,J),J=1,2),IP=1,4),ANG(I),I=1,N0)
      WRITE(8,8200) (((IMG(I,M,N),N=1,ISNGRD),M=1,ISNGRD),
     2   (ICNR(I,MM),MM=1,4),I=1,N0) 
C 
 3400 CONTINUE
      STOP
C 
C6005 FORMAT(1H1,10X,19HINPUT DATA FOR ROOM/11X,19H*******************/)
 6005 format(1H1,10X,'INPUT DATA FOR ROOM',/,11X,'*******************',
     1/)
 6000 FORMAT(/5X,33HNUMBER OF CLEAR WINDOWS          ,I3/ 
     1        5X,33H          SHEER CURTAIN WINDOWS  ,I3/ 
     2        5X,33H          DIFFUSE WINDOWS        ,I3/ 
     3        5X,33H          WALLS WITHOUT CUT-OUTS ,I3/ 
     4        5X,33H          WALLS WITH CUT-OUTS    ,I3/ 
     5        5X,33H          WORKING SURFACES       ,I3/ 
     6        5X,33H          OUTSIDE ENCLOSURES     ,I3///)
 6010 FORMAT( 5X,7HSURFACE,I3/5X,10H----------//
     1       10X,4HX0 =,F6.2,6H  Y0 =,F6.2,6H  Z0 =,F6.2/ 
     2       10X,4HX1 =,F6.2,6H  X2 =,F6.2,6H  X3 =,F6.2,5H  Y =,F6.2/
     3       10X,7HBETAX =,F6.1,8H  PSIX =,F6.1,9H  BETAY =,F6.1, 
     4           8H  PSIY =,F6.1/ 
     5       10X,14HREFLECTANCE  =,F4.1)
 6020 FORMAT(10X,15HGRID DIMENSION=,I3,4H  BY,I3//) 
 6030 FORMAT(10X,23HWINDOW SURROUNDING NO.=,I3,16H  MAINT.FACTOR =,F4.2/
     1       10X,12HTAU-NORMAL =,F4.2,20H  WINDOW THICKNESS =,F6.2//) 
 6031 FORMAT(15X,29HWINDOW HAS AN OVERHANG ON TOP/
     1       20X,6HWIDTH=,F4.1,18H  DIST.FROM GLASS=,F4.1,15H  REFLECTAN
     2CE =,F4.2)
 6032 FORMAT(15X,31HWINDOW HAS OVERHANGS BOTH SIDES/
     1       20X,6HWIDTH=,F4.1,18H  DIST.FROM GLASS=,F4.1,15H  REFLECTAN
     2CE =,F4.2)
 6033 FORMAT(15X,34HWINDOW HAS OVERHANG ON BOTTOM ONLY/ 
     1       20X,6HWIDTH=,F4.1,18H  DIST.FROM GLASS=,F4.1,15H  REFLECTAN
     2CE =,F4.2)
 6034 FORMAT(15X,37HWINDOW HAS OVERHANGS ON TOP AND SIDES/
     1       20X,6HWIDTH=,F4.1,18H  DIST.FROM GLASS=,F4.1,15H  REFLECTAN
     2CE =,F4.2)
 6035 FORMAT(15X,38HWINDOW HAS OVERHANGS ON ALL FOUR SIDES/ 
     1       20X,6HWIDTH=,F4.1,18H  DIST.FROM GLASS=,F4.1,15H  REFLECTAN
     2CE =,F4.2)
 6064 FORMAT(10X,47HWINDOW HAS SHEER CURTAIN, SHEER CURTAIN FACTOR=,
     1       F6.2//)
 6067 FORMAT(10X,40HDIRECT-SUN TRANSMISSION FRACTION, ALPHA ) 
 6068 FORMAT(15X,6HWINDOW,I3,8H  ALPHA=,F5.3/)
 6036 FORMAT(15X,18HOVERHANG IS OPAQUE//) 
 6037 FORMAT(15X,45HOVERHANG IS SEMI-TRANSPARENT, TRANSMISSIVITY=,F4.2
     1//) 
 6038 FORMAT(15X,40HOVERHANG IS TRANSLUCENT, TRANSMISSIVITY=,F4.2//)
 6040 FORMAT(5X,15HSURROUNDING NO.,I3/5X,18H******************//
     1      10X,23HNUMBER OF OPAQUE WALLS ,I3/) 
 6041 FORMAT(//2X,37HINSIDE-OVERHANG LUMINANCES FOR WINDOW,I2,
     1            15H (FOOT-LAMBERT)/ 
     2         2X,39H***************************************, 
     3            15H***************/)
 6042 FORMAT(10X,10HTOP       ,F7.1)
 6043 FORMAT(10X,10HLEFT SIDE ,F7.1/
     1       10X,10HRIGHT SIDE,F7.1)
 6044 FORMAT(10X,10HBOTTOM    ,F7.1)
 6045 FORMAT(10X,10HTOP       ,F7.1/
     1       10X,10HLEFT SIDE ,F7.1/
     2       10X,10HRIGHT SIDE,F7.1)
 6046 FORMAT(10X,10HTOP       ,F7.1/
     1       10X,10HLEFT SIDE ,F7.1/
     2       10X,10HRIGHT SIDE,F7.1/
     3       10X,10HBOTTOM    ,F7.1)
 6050 FORMAT(1H1////5X,30HDIRECT HORIZONTAL INSOLATION  ,F6.2,1X,
     1       16HWATTS PER SQ.FT./5X,18HLUMINOUS EFFICACY ,F6.1,1X,
     2       15HLUMENS PER WATT/5X,30HDIFFUSE HORIZONTAL INSOLATION ,
     3       F6.2,1X,16HWATTS PER SQ.FT./5X,18HLUMINOUS EFFICACY ,
     4       F6.1,1X,15HLUMENS PER WATT/ 
     5       5X,22HSKY MODEL             ,I3/ 
     6       5X,20HGROUND REFLECTANCE  ,F4.2) 
 6070 FORMAT(//5X,15HWORKING SURFACE,I3/5X,18H******************//
     1       10X,4HX0 =,F6.2,6H  Y0 =,F6.2,6H  Z0 =,F6.2/ 
     2       10X,4HX1 =,F6.2,6H  X2 =,F6.2,6H  X3 =,F6.2,5H  Y =,F6.2/
     3       10X,7HBETAX =,F6.1,8H  PSIX =,F6.1,9H  BETAY =,F6.1, 
     4           8H  PSIY =,F6.1/ 
     5       10X,15HGRID DIMENSION=,I3,4H  BY,I3//) 
 6100 FORMAT(///10X,46HSKY LUMINANCE DETERMINED BY GEOGRAPHICAL DATA /
     A          15X,24HLATITUDE                ,F6.1,5H DEG./ 
     B          15X,24HLONGITUDE               ,F6.1,5H DEG./ 
     C          15X,24HTIME ZONE               ,I4/ 
     D          15X,24HALTITUDE                ,F6.1,7H METERS/ 
     E          15X,24HSKY MODEL               ,I4/ 
     F          15X,24HGROUND REFLECTANCE      ,F7.2) 
 6110 FORMAT(   15X,24HFIRST MONTH             ,I4/ 
     A          15X,24HLAST MONTH              ,I4/ 
     B          15X,24HINCREMENT BETWEEN MONTHS,I4/ 
     C          15X,24HDAY OF MONTH            ,I4/ 
     D          15X,24HFIRST TIME OF DAY       ,F6.1,5H HRS./ 
     E          15X,24HLAST TIME OF DAY        ,F6.1,5H HRS./ 
     F          15X,24HTIME INCREMENT          ,F6.1,5H HRS.) 
 6105 FORMAT(///10X,46HSKY LUMINANCE DETERMINED BY SUN POSITON  DATA /
     A          15X,24HALTITUDE ANGLE          ,F6.1,5H DEG./ 
     B          15X,24HAZIMUTH ANGLE           ,F6.1,5H DEG./ 
     D          15X,24HALTITUDE OF STATION     ,F6.1,7H METERS/ 
     E          15X,24HSKY MODEL               ,I4/ 
     F          15X,24HGROUND REFLECTANCE      ,F7.2) 
 6120 FORMAT(15X,45HMONTHLY VALUES FOR TURBIDITY AND COND.WATER../
     A       20X,20HMONTH  C.W.    TURB.) 
 6121 FORMAT(21X,I3,2F8.4)
 6130 FORMAT(//,5X,7HMONTH =,I3,10X,6HTIME =,F5.1,5H HRS./
     1           5X,36H====================================)
 6140 FORMAT(5X,20HSUN ALTITUDE ANGLE =,F6.1,/
     15X,19HSUN AZIMUTH ANGLE =,F6.1,1X,27HDEG. OFF SOUTH TOWARDS EAST) 
 6200 FORMAT(7X,15F6.0) 
 6250 FORMAT(//10X,35HILLUMINANCE ON A HORIZONTAL SURFACE/
     1         15X,16HSOLAR COMPONENT=,F8.0,12H FOOT-CANDLE/
     2         15X,16HSKY   COMPONENT=,F8.0,12H FOOT-CANDLE/
     3         10X,17HZENITH LUMINANCE=,F8.0,1X,14H(FT.-LAMBERTS))
 6300 FORMAT(1H1/15X,47H***********************************************/
     *       15X,47H*                                             */
     1       15X,47H* DAYLIGHT ILLUMINATION IN A RECTANGULAR ROOM */
     *       15X,47H*                                             */
     2       15X,47H***********************************************// 
     3       10X,33HROOM HAS A WIDTH (WINDOW WALL) OF,F5.1,6H UNITS/
     4       10X,33H                       A DEPTH OF,F5.1,6H UNITS/
     5       10X,33H                  AND A HEIGHT OF,F5.1,6H UNITS/
     6       10X,47HROOM ORIENTATION (FROM OUTWARD WINDOW NORMAL) =,
     7               F5.1,15H DEG. OFF SOUTH/)
 6302 FORMAT(10X,8HROOM HAS,I2,28H WINDOW(S) IN THE FRONT WALL) 
 6303 FORMAT(10X,8HROOM HAS,I2,37H WINDOW(S) IN THE CEILING (SKYLIGHTS))
 6310 FORMAT(/10X,10HWINDOW NO.,I2,3H IS,F6.2,14H UNITS WIDE BY,F6.2, 
     1            11H UNITS HIGH/15X,16HWINDOW CENTER IS,F6.2,
     2            30H UNITS TO RIGHT OF WALL CENTER/
     3            31X,F6.2,17H UNITS FROM FLOOR/
     4        15X,24HWINDOW DISCRETIZATION IS,I3,3H BY,I3)
 6311 FORMAT(15X,16HWINDOW TYPE IS ",I1,16H" (CLEAR WINDOW))
 6312 FORMAT(15X,16HWINDOW TYPE IS ",I1,34H" (SHEER CURTAIN, CLEAR FRACT
     1ION =,F4.2,1H)) 
 6313 FORMAT(15X,16HWINDOW TYPE IS ",I1,20H" (DIFFUSING WINDOW))
 6320 FORMAT(15X,23HMAINTENANCE FACTOR    =,F4.2/ 
     1       15X,23HWINDOW THICKNESS      =,F4.2/ 
     2       15X,23HWINDOW TRANSMISSIVITY =,F4.2/ 
     3       15X,23HWINDOW REFLECTANCE    =,F4.2) 
 6331 FORMAT(10X,32HFRONT  WALL..     DISCRETIZATION,I3,3H BY,I3, 
     1                                16H, REFLECTANCE  =,F4.2) 
 6332 FORMAT(10X,32HLEFT SIDE..       DISCRETIZATION,I3,3H BY,I3, 
     1                                16H, REFLECTANCE  =,F4.2) 
 6333 FORMAT(10X,32HREAR WALL..       DISCRETIZATION,I3,3H BY,I3, 
     1                                16H, REFLECTANCE  =,F4.2) 
 6334 FORMAT(10X,32HRIGHT SIDE..      DISCRETIZATION,I3,3H BY,I3, 
     1                                16H, REFLECTANCE  =,F4.2) 
 6335 FORMAT(10X,32HCEILING..         DISCRETIZATION,I3,3H BY,I3, 
     1                                16H, REFLECTANCE  =,F4.2) 
 6336 FORMAT(10X,32HFLOOR..           DISCRETIZATION,I3,3H BY,I3, 
     1                                16H, REFLECTANCE  =,F4.2) 
 6337 FORMAT(10X,32HWORKING SURFACE.. DISCRETIZATION,I3,3H BY,I3, 
     1                                13H, ELEVATION =,F5.2,6H UNITS) 
 6340 FORMAT(/10X,18HSUBJECT BUILDING../
     1        15X,32HBUILDING WIDTH (FRONT  WALL)   =,F6.1,6H UNITS/
     2        15X,32HBUILDING HEIGHT                =,F6.1,6H UNITS/
     3        15X,32HOFFSET TO RIGHT OF ROOM CENTER =,F6.1,6H UNITS/
     4        15X,32H                REFLECTANCE    =,F7.2) 
 6350 FORMAT(/10X,13HOBSTRUCTION../ 
     1        15X,32HOBSTRUCTION WIDTH              =,F6.1,6H UNITS/
     2        15X,32HOBSTRUCTION HEIGHT             =,F6.1,6H UNITS/
     3        15X,32HOFFSET TO RIGHT OF ROOM CENTER =,F6.1,6H UNITS/
     4        15X,32HDISTANCE FROM SUBJECT BUILDING =,F6.1,6H UNITS/
     5        15X,32H    OBSTRUCTION REFLECTANCE    =,F7.2) 
 6560 FORMAT(//2X,69HLIGHT EXCHANGE FACTORS FROM SURFACE I TO SURFACE J 
     1IN SURROUNDING NO.,I2,/2X,71H*************************************
     2**********************************/)
 6570 FORMAT(/7X,1HJ,15I8)
 6572 FORMAT(4X,1HI)
 6580 FORMAT(2X,I3,5X,15F8.4/)
 6710 FORMAT(1H1//3X,71H************************************************
     1***********************/
     2         3X,71H* DATA FOR SURFACE NODES; I=SURFACE, K=NODE,   X,Y,
     3Z=COORDINATES (FT) */ 
     4         3X,71H*                         A=AREA (SQFT), L=LUMINANC
     5E (FT-LAMBERT)     */ 
     6         3X,71H***************************************************
     7********************) 
 6715 FORMAT(1H1//3X,71H************************************************
     1***********************/
     2         3X,71H* DATA FOR WORKING SURFACE NODES; I=SURFACE, K=NODE
     3-NO.               */ 
     4         3X,71H*      X,Y,Z=COORDINATES,                          
     5                   */ ) 
6716  FORMAT(  3X,71H*      S=ILLUMINANCE FROM EXTERNAL SOURCES (FOOT-CA
     7NDLE)              */ 
     8         3X,71H*      R=INTERNAL REFLECTED COMPONENT (FOOT-CANDLE)
     9                   */ 
     A         3X,71H*      I=TOTAL ILLUMINANCE (FOOT-CANDLE), D=DAYLIGH
     BT FACTOR (PERCENT) */ 
     6         3X,71H***************************************************
     7********************) 
 6717 FORMAT(//10X,53HX = DISTANCE FROM FRONT WALL (PARALLEL TO FRONT WA
     1LL)/
     2         10X,64HY = DISTANCE FROM ROOM CENTER LINE (PERPENDICULAR 
     3TO FRONT WALL)/ 
     4         14X,37HPOSITIVE VALUES (LEFT OF CENTER LINE)/
     5         14X,38HNEGATIVE VALUES (RIGHT OF CENTER LINE)/ 
     6         10X,24HZ = HEIGHT OF WORK PLANE) 
 6720 FORMAT(//3X,2H*K,12I10) 
 6730 FORMAT(3X,2HI*/I4,3H  X,12F10.2)
 6740 FORMAT(6X,1HY,12F10.2)
 6750 FORMAT(6X,1HZ,12F10.2)
 6760 FORMAT(6X,1HA,12F10.2)
 6770 FORMAT(6X,1HL,12F10.2)
 6773 FORMAT(6X,1HS,12F10.2)
 6774 FORMAT(6X,1HR,12F10.2)
 6775 FORMAT(6X,1HI,12F10.2)
 6780 FORMAT(6X,1HD,12F10.2)
 7010 FORMAT(//5X,41HLUMINANCE OF SURFACE I IN SURROUNDING NO.,I2,1X,14H
     1(FT.-LAMBERTS)/5X,58H=============================================
     2=============)
 7020 FORMAT(/5X,1HI,15I6)
 8000 FORMAT(2X,33HTHERE IS TOO MUCH DIRECT SUNLIGHT/ 
     1       2X,35HAS COMPARED TO TOTAL HEMISPHERICAL,/ 
     2       2X,20HEXECUTION TERMINATED)
 8002 FORMAT(32HNO. OF WINDOWS EXCEEDS MAXIMUM 5) 
 8005 FORMAT(33HILLEGAL REFLECTANCE FOR SURFACE  ,I2, 
     2        14H IN ENCLOSURE ,I2) 
 8007 FORMAT(38HNO. OF WORK SURFACES EXCEEDS MAXIMUM 3) 
 8008 FORMAT(35HNO. OF ENCLOSURES EXCEEDS MAXIMUM 2)
 8010 FORMAT(41HNO. OF INDOOR SURFACES EXCEEDS MAXIMUM 30)
 8011 FORMAT(44HNO. OF NODES EXCEEDS 400 FOR INDOOR SURFACE ,I2)
 8012 FORMAT(40HILLEGAL REFLECTANCE FOR INDOOR SURFACE  ,I2)
 8014 FORMAT(44HILLEGAL SURFACE OCCUPYING CUTOUT IN SURFACE ,I2)
 8020 FORMAT(29HILLEGAL ENCLOSURE FOR WINDOW ,I2) 
 8022 FORMAT(34HILLEGAL TRANSMISSIVITY FOR WINDOW ,I2)
 8024 FORMAT(38HILLEGAL MAINTENANCE FACTOR FOR WINDOW ,I2)
 8026 FORMAT(27HILLEGAL GROUND REFLECTANCE )
 8030 FORMAT(54HTHICKNESS OF CONDENSABLE WATER OUT OF RANGE FOR MONTH , 
     1    I2) 
 8032 FORMAT(45HTURBIDITY COEFFICIENT OUT OF RANGE FOR MONTH ,I2) 
 8034 FORMAT(41HILLEGAL OVERHANG REFLECTANCE FOR WINDOW  ,I2) 
 8036 FORMAT(43HILLEGAL OVERHANG TRANSMISSIVITY FOR WINDOW ,I2) 
 8100 FORMAT(6F10.3)
 8200 FORMAT(20I4)
      END 
      FUNCTION C(K,NK)
C THIS FUNCTION CALCULATES NEWTON-COTES COEFFICIENTS FOR NUMERICAL
C QUADRATURE
      C=1.
      IF(NK.EQ.1)RETURN 
      I=(K+1)/2-K/2 
      A=1+2*(1-I)+((K+100)/102)-(K/NK)
      C=A/(3.*(NK-1.))
      RETURN
      END 
      SUBROUTINE SECT(X,Y,K,J,KUT)
C THIS ROUTINE CHECKS WHETHER A RAY GOING FROM X INTO DIRECTION Y CUTS
C    SURFACE J IN ENCLOSURE K FROM ITS FRONT SIDE (KUT=1) OR NOT (KUT=0)
      REAL X(3),Y(3),X0(4,30,7),LL(4,30,9)
      COMMON /BLK1/ X0,LL,XDP,YDP 
C 
C INTERSECTION POINT IN LOCAL COORDINATES 
      BN=0. 
      BD=0. 
      DO 1 I=1,3
      BN=BN+(X(I)-X0(K,J,I))*LL(K,J,I+6)
    1 BD=BD+ Y(I)*LL(K,J,I+6) 
      BD=BD+1E-10+BD*1E-5 
      B=BN/BD 
      XDP=0.
      YDP=0.
      DO 2 I=1,3
      CC=X(I)-X0(K,J,I)-B*Y(I)
      XDP=XDP+CC*LL(K,J,I)
    2 YDP=YDP+CC*LL(K,J,I+3)
      KUT=0 
C WITHIN CONFINES OF SURFACE? 
      ETA=YDP/X0(K,J,7) 
      IF(ETA.LT.0..OR.ETA.GT.1.) RETURN 
      XI=(XDP-ETA*X0(K,J,5))/(X0(K,J,4)-(X0(K,J,4)-X0(K,J,6)+ 
     1   X0(K,J,5))*ETA)
      IF(XI.LT.-1.E-5.OR.XI.GT.1.00001) RETURN
C IS IT HITTING TOP OF SURFACE? 
      KUT=1-(Y(1)*LL(K,J,7)+Y(2)*LL(K,J,8)+Y(3)*LL(K,J,9))*1E-4 
      RETURN
      END 
      SUBROUTINE CONFAC(J,NT,F) 
C THIS ROUTINE CALCULATES EXCHANGE FACTORS FOR ENCLOSURE J WITH NT
C SURFACES, INCLUDING EXCHANGE FACTORS TO SKY WITH LUMINANCE VARIATION
      REAL X0(4,30,7),LL(4,30,9),F(4,30,30),X(3),Y(3),
     1     A(4,30),RS(4),LPZ
      INTEGER NW(6),IOB(4,30,30),KJOB(4,30,30)
      INTEGER*4 ISEED
      COMMON /BLK1/ X0,LL,XDP,YDP 
      COMMON /BLK2/ NW,PI,IOB,ICLD,ISEED,NBUNDL,A,KJOB
      NT1=NT+1
      DO 10 I=1,NT
      DO 5 IR=1,NT1 
    5 F(J,I,IR)=0.
      XY=(X0(J,I,6)-X0(J,I,5))/X0(J,I,4)
      ID1=1 
      IF(ABS(XY-1).LT.1E-3)ID1=0
      AX=ID1*XY 
C CALCULATION OF COORDINATES AND DIRECTION COSINES OF LIGHT BUNDLE
      DO 60 N=1,NBUNDL
      TT=1. 
C ************************ RANDOM NUMBER GENERATOR ********************** 
      DO 115,INDEX=1,4
115   RS(INDEX)=RANDOM(ISEED)
      YY=(1-SQRT(1-(1-XY*XY)*RS(1)))/(1-AX)+(1-ID1)*RS(1) 
      XX=X0(J,I,5)*YY+X0(J,I,4)*(1-(1-XY)*YY)*RS(2) 
      SINB=SQRT(RS(3))
      COSB=SQRT(1.-RS(3)) 
      THET=2.*PI*RS(4)
      COST=COS(THET)
      SINT=SIN(THET)
      DO 30 L=1,3 
      X(L)=X0(J,I,L)+LL(J,I,L)*XX+LL(J,I,L+3)*YY*X0(J,I,7)
      Y(L)=(LL(J,I,L)*COST+LL(J,I,L+3)*SINT)*SINB+LL(J,I,L+6)*COSB
  30  CONTINUE
C 
C DETERMINATION OF SURFACE INTERSECTED BY LIGHT BUNDLE
C 
      DO 80 IR=1,NT 
      IF(IOB(J,I,IR)) 80,70,50
   50 IMARK=KJOB(J,I,IR)
      CALL SECT(X,Y,J,IMARK,KUT)
      IF(KUT.EQ.1) GOTO 60
   70 IMARK=IR
      CALL SECT(X,Y,J,IR,KUT) 
      IMARK=IR
      IF(KUT.EQ.1) GOTO 60
   80 CONTINUE
C 
C EVALUATION OF F I-J IF BUNDLE INTERSECTS SKY, (OR GROUND IF Y POINTS BELOW HORIZ) 
C 
      IMARK=NT1 
      TT=RLPZ(ICLD,Y) 
   60 F(J,I,IMARK)=F(J,I,IMARK)+TT
      DO 130 IR=1,NT1 
  130 F(J,I,IR)=F(J,I,IR)/NBUNDL
   10 CONTINUE
C 
C AVERAGING OF F I-J BASED UPON SURFACE AREA
C 
      DO 220 I=1,NT 
      DO 220 IR=I,NT
      WF=1-EXP(-A(J,I)/A(J,IR))+EXP(-A(J,IR)/A(J,I))
      AF=((2-WF)*A(J,I)*F(J,I,IR)+WF*A(J,IR)*F(J,IR,I))/2.
      F(J,I,IR)=AF/A(J,I) 
  220 F(J,IR,I)=AF/A(J,IR)
      RETURN
      END 
      FUNCTION TAUB(I,CB) 
C THIS ROUTINE CALCULATES THE TRANSMISSIVITY OF WINDOW I FOR DIRECTION
C CB (COSINE OF POLAR ANGLE)
      REAL AK2(10),TAU(10)
      COMMON /BLK3/ AK2,TAU,NC
      TAUB=1.018*AK2(I)*TAU(I)*CB*(1.+(1.-CB*CB)**1.5)
      IF (I.GT.NC) TAUB=TAU(I)
      RETURN
      END 
      FUNCTION RLPZ(ICLD,Y) 
C THIS ROUTINE CALCULATES SKY-LUMINANCE (NORMALIZED BY ZENITH LUMINANCE)
C   FOR ANY GIVEN DIRECTION AND SKY CONDITION 
      REAL Y(3),LS(3) 
      COMMON /BLK4/ BETAS,LS,AUX(3),GRNDL 
      IF (Y(3).GT..001) GOTO 10 
C HITS GROUND 
      RLPZ=GRNDL
      RETURN
   10 GOTO(1,2,3,2),ICLD
C OVERCAST SKY
    1 IY=.999+Y(3)
      RLPZ=AUX(1)+IY*Y(3)*AUX(2)
      RETURN
C CLEAR SKY 
    2 CD=Y(1)*LS(1)+Y(2)*LS(2)+Y(3)*LS(3) 
      RLPZ=(1.-EXP(-.32/Y(3)))*(.91+10*EXP(-3*ACOS(CD))+.45*CD*CD)
     1     /AUX(3)
      RETURN
C UNIFORM SKY 
    3 RLPZ=1. 
      RETURN
      END 
      SUBROUTINE MATINV(A,N,B,M,DETERM) 
      DIMENSION IPIVOT(30),A(30,30),B(30,1),INDEX(30,2),PIVOT(30) 
      EQUIVALENCE (IROW,JROW), (ICOLUM,JCOLUM), (AMAX, T, SWAP) 
      DETERM=1.0
      DO 20 J=1,N 
   20 IPIVOT(J)=0 
      DO 550 I=1,N
      AMAX=0.0
      DO 105 J=1,N
      IF (IPIVOT(J)-1) 60, 105, 60
   60 DO 100 K=1,N
      IF (IPIVOT(K)-1) 80, 100, 740 
   80 IF (ABS(AMAX)-ABS(A(J,K))) 85, 100, 100 
   85 IROW=J
      ICOLUM=K
      AMAX=A(J,K) 
  100 CONTINUE
  105 CONTINUE
      IF(AMAX) 110,800,110
  110 IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1 
      IF (IROW-ICOLUM) 140, 260, 140
  140 DETERM=-DETERM
      DO 200 L=1,N
      SWAP=A(IROW,L)
      A(IROW,L)=A(ICOLUM,L) 
  200 A(ICOLUM,L)=SWAP
      IF(M) 260, 260, 210 
  210 DO 250 L=1, M 
      SWAP=B(IROW,L)
      B(IROW,L)=B(ICOLUM,L) 
  250 B(ICOLUM,L)=SWAP
  260 INDEX(I,1)=IROW 
      INDEX(I,2)=ICOLUM 
      PIVOT(I)=A(ICOLUM,ICOLUM) 
      DETERM=DETERM*PIVOT(I)
      A(ICOLUM,ICOLUM)=1.0
      DO 350 L=1,N
  350 A(ICOLUM,L)=A(ICOLUM,L)/PIVOT(I)
      IF(M) 380, 380, 360 
  360 DO 370 L=1,M
  370 B(ICOLUM,L)=B(ICOLUM,L)/PIVOT(I)
  380 DO 550 L1=1,N 
      IF(L1-ICOLUM) 400, 550, 400 
  400 T=A(L1,ICOLUM)
      A(L1,ICOLUM)=0.0
      DO 450 L=1,N
  450 A(L1,L)=A(L1,L)-A(ICOLUM,L)*T 
      IF(M) 550, 550, 460 
  460 DO 500 L=1,M
  500 B(L1,L)=B(L1,L)-B(ICOLUM,L)*T 
  550 CONTINUE
      DO 710 I=1,N
      L=N+1-I 
      IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630
  630 JROW=INDEX(L,1) 
      JCOLUM=INDEX(L,2) 
      DO 705 K=1,N
      SWAP=A(K,JROW)
      A(K,JROW)=A(K,JCOLUM) 
      A(K,JCOLUM)=SWAP
  705 CONTINUE
  710 CONTINUE
  740 RETURN
  800 DETERM = 0. 
      RETURN
      END 
      SUBROUTINE ZENLCD(P,W,BETA,H,ZENL)
C  THIS SUBROUTINE CALCULATES ZENITH LUMINANCE BASED ON 
C  LIEBELT'S EQUATION WHEN SOLAR ALTITUDE IS LESS THAN
C  60.  WHEN SOLAR ALTITUDE IS GREATER THAN 60, THE ZENITH
C  LUMINANCE IS EXTRAPOLATED TO RESULT INTEGRATED HORIZONTAL
C  ILLUMINANCE FOLLOWS SINE CURVE. 4/10/83 JJK. 
      PI=4.*ATAN(1.)
      HD=H*180./PI
      T=TUR(HD,W,BETA)
      ZENL=((1.34*T-3.46)*TAN(H)+0.1*T+0.9)*291.87
      IF(HD.LT.60.)RETURN 
      ZEN60=((1.34*T-3.46)*TAN(1.0472)+0.1*T+0.9)*291.87
      P60=SKYSUM(0.5236)
      ZENL=1.1547*ZEN60*P60*SIN(H)/P
      RETURN
      END
      SUBROUTINE SOLAR (DLAT,LNG,TZN,M,TIME,ISUN,THETA,PHI) 
C  TO FIND THE POSITION OF THE SUN
C  INPUT
C    DLAT  - LATITUDE
C    LNG  - LONGITUDE 
C    TZN  - TIME ZONE(HOURS BEHIND GREENWICH MEAN TIME
C            4- ATLANTIC
C            5- EASTERN 
C            6- CENTRAL 
C            7- MOUNTAIN
C            8- PACIFIC 
C    M    - MONTH(01-12)
C    TIME - HOURS FROM MIDNIGHT 
C 
C  OUTPUT 
C    THE ALTITUDE AND AZIMUTH ANGLES OF THE SUN.
      REAL LAT, LNG 
      INTEGER TZN 
      DIMENSION DEC(12), ET(12) 
C 
      DATA (DEC(J), J=1,12)/ -0.3491, -0.1885,  0.0   ,  0.2025,
     $                        0.3491,  0.4093,  0.3595,  0.2147,
     $                        0.0   , -0.1833, -0.3456, -0.4093/
C 
      DATA (ET(J), J=1,12)/ -0.190, -0.230, -0.123,  0.020, 
     $                        0.060, -0.025, -0.103,-0.051, 
     $                        0.113,  0.255,  0.235,  0.033/
      DATA PI/3.1415926/
      DATA CNVRT/0.017453295/ 
C  WHEN SUN IS DOWN THERE IS NO RADIATION 
      LAT = DLAT*CNVRT
      SPAN = ACOS(-TAN(LAT)*TAN(DEC(M)))
      HANG =-(15.0*(TIME - 12.0 + TZN + ET(M)) - LNG)*CNVRT 
      CHANG = COS(HANG) 
      IF(ABS(HANG) .GT. ABS(SPAN)) THEN 
          ISUN = 0
C  WHEN SUN IS UP, FIND ITS ORIENTATION 
      ELSE
          ISUN = 1
          COSZ = SIN(LAT)*SIN(DEC(M)) + 
     $           COS(LAT)*COS(DEC(M))*COS(HANG) 
          COSW = COS(DEC(M))*SIN(HANG)
      TEST = 1.0 - COSZ**2 + COSW**2
      ARG = TAN(DEC(M))/TAN(LAT)
      IF (TEST .LE. 0.0) THEN 
      COSS = 0.0
      ELSEIF (CHANG .LE. ARG) THEN
      COSS = -SQRT(TEST)
      ELSE
      COSS = SQRT(TEST) 
      ENDIF 
      THETA = ASIN(COSZ)
      CHECK = COSW/COS(THETA) 
      IF (CHECK .GE. 1.0) THEN 
      ANG = PI/2.0
      ELSEIF (CHECK .LE. -1.0) THEN
      ANG = -PI/2.0 
      ELSE
      ANG = ASIN(CHECK) 
      ENDIF 
      IF (COSS .GT. 0.0) THEN
      PHI = ANG 
      ELSE
      PHI = PI - ANG
      ENDIF 
      ENDIF 
      RETURN
      END 
      FUNCTION E(DATE,IWMOD)
      REAL OMEGA
      DATA OMEGA/.01715/
      A = OMEGA*DATE
      IF(IWMOD.NE.3) GOTO 100 
      E=128.8 
      RETURN
  100 E = 128.82 + 4.248*COS(A) + 0.0825*COS(2.0*A) 
     $-0.00043*COS(3.0*A) + 0.169*SIN(A) + 0.00914*SIN(2.0*A) 
     $+0.01726*SIN(3.0*A) 
      RETURN
      END 
      FUNCTION ABAR(T,BETA) 
      DIMENSION A(3),B(3) 
      DATA (A(I),I=1,3)/0.1512,0.1656,0.2021/ 
      DATA (B(I),I=1,3)/0.0262,0.0215,0.0193/ 
      INDEX=3 
      IF (BETA .LE. .15) INDEX=2
      IF (BETA .LE. .075) INDEX=1 
      ABAR = A(INDEX) - T*B(INDEX)
      RETURN
      END 
      FUNCTION TUR(HD,W,BETA) 
      TUR=(HD + 85.0)/(39.5*EXP(-W) + 47.4) + 0.1 
     $+ (16.0+0.22*W)*BETA
      RETURN
      END 
      FUNCTION AIR(H,ALT) 
      AIR = SIN(H) + 0.15*(H*57.29578+3.885)**(-1.253)
      AIR = 1/AIR 
      AIR = AIR*(1.0-1.E-4*ALT) 
      RETURN
      END 
      FUNCTION DIR(MONTH,IDAY,BETA,W,ALT,H) 
      DATA CFAC/92.904/ 
C CFAC CONVERTS FROM KLX TO FT-CANDLES
      HD = H*57.2958
      T = TUR(HD,W,BETA)
      A = ABAR(T,BETA)
      AMAS = AIR(H,ALT) 
      DATE = FLOAT(JDATE(MONTH,IDAY)) 
      DIR = CFAC*E(DATE,IWMOD)*EXP(-A*AMAS*T)*SIN(H)
      RETURN
      END 
      FUNCTION JDATE(MONTH,IDAY)
      DIMENSION NDAY(12)
      DATA (NDAY(I),I=1,12)/31,28,31,30,31, 
     $30,31,31,30,31,30,31/ 
      JDATE = 0 
      IF (MONTH .EQ. 1) GO TO 200 
      M1=MONTH-1
      DO 100 I=1,M1 
      JDATE = JDATE + NDAY(I) 
  100 CONTINUE
  200 JDATE = JDATE + IDAY
      RETURN
      END 
      FUNCTION SKYSUM(H)
C --  THIS SUBROUTINE CALCULATES SKY LUMINANCE INTEGRATION FOR
C--   CLEAR SKY 
      PI=4.*ATAN(1.)
      P1=.27385*(.91+10.*EXP(-3*H)+.45*COS(H)**2) 
      P=0.
      CBS=COS(H)
      SBS=SIN(H)
      DO 100 K=1,21 
      ETAK=(K-1.)/20. 
      RT=SQRT(1.-ETAK**2) 
      P2=0. 
      IF(K.GT.1)P2=(1.-EXP(-.32/ETAK))*C(K,21)*ETAK 
      DO 100 L=1,21 
      ETA=(L-1.)/20.
      CDEL=CBS*ETAK+SBS*RT*COS(PI*ETA)
  100 P=P+P2*(.91+10.*EXP(-3*ACOS(CDEL))+.45*CDEL**2)*C(L,21) 
      SKYSUM=2./P1*P
      RETURN
      END 
C
      FUNCTION RANDOM(ISEED)
      INTEGER*4 ISEED
      ISEED=DMOD(16807.0D0 * ISEED,2147483647.0D0)
      RANDOM=ISEED * 4.6566128752458D-10
      RETURN
      END
      SUBROUTINE LOGO(N)
C     THIS ROUTINE PRINTS THE SUPERLITE NOTICE TO THE UNIT
C     SPECIFIED AS 'N'
C     CURRENTLY N=6 SENDS NOTICE TO SCREEN ON VAX.  FOR BATCH
C     MACHINES THIS CAN BE SET TO THE APPROPRIATE DEVICE
      WRITE(*,100)
 100  FORMAT(/////////
     *6X,'                 S  U  P  E  R  L  I  T  E    1 . 0'/
     *//
     *6X,'                           DEVELOPED BY:'/
     */
     *6X,'                    WINDOWS AND DAYLIGHTING GROUP'/
     *6X,'                    LAWRENCE BERKELEY LABORATORY'/
     *6X,'                      UNIVERSITY OF CALIFORNIA'/
     *6X,'                        BERKELEY, CA  94720'/
     *///
     *6X,'      COPYRIGHT 1985  REGENTS OF THE UNIVERSITY OF CALIFORNIA'
     *//)
      END
