Simtel MSDOS 1992 June/SIMTEL_0692.cdr/msdos/progjorn/pj_6_1.arc/SLITE.ARC/SLITE3.FOR

Text File | 1987-09-07 | 77KB | 2,416 lines

      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