      SUBROUTINE SOLR(TIME,NDYRI,IDAY,KM,STM,DNLAT,DLONG, 
     +                DST,ICC,CLEAR,
     +                ISUN,ITH,IDHC,SAZM,SALT,
     +                IDC,IDFC,DIRECT,DIFFUS,SFST,SFTT)
C 
C  CALCULATE SOLAR FACTORS FOR TIME=TIME AND DAY = NDYRI
C 
C  INPUTS FROM FORMAL PARAMETERS IN CALL STMT 
C   TIME     TIME OF DAY,(0-24 HRS) 
C   NDYRI    DAY OF YEAR
C   IDAY     NO. OF DAYS RUN AT THIS TIME
C             (IDAY IS ONLY USED TO IDENTIFY DEBUG OUTPUT)
C   KM       MONTH OF YEAR (1 TO 12)
C  STM       TIME ZONE (SEE ETIME SUBPROGRAM) 
C  ICC     =1  CLOUD COVER DATA USED TO COMPUTE SOLAR 
C          =0  TOTAL SOLAR RADIATION SUPPLIED IN  CC  ARRAY 
C DNLAT      DEG N. LATTITUDE 
C DLONG      DEG LONGITUDE
C   DST    =0 FOR STANDARD TIME 
C          =1 FOR DAYLIGHT SAVINGS TIME 
C  CLEAR   CLEARNESS NO. FOR USE WITH CLOUD COVER DATA
C 
C  INPUTS FROM COMMON BLOCK ROOF - - - -
C   IDEXPR  INTEGER 1-4 COMPASS DIRECTION THAT ROOF FACES 
C     PAA   PITCH ANGLE BETWEEN VERTICAL AND ROOF LINE 
C 
C INPUTS FROM COMMON BLOCK NCR - - - - -
C   IFXDS  INTEGER 1-4 COMPASS DIRECTIONS OF STRUCTURE THAT ARE ALWAYS
C          SHADED.
C   WROTA  ANGLE THAT STRUCTURE IS ROTATED FROM THE CARDINAL COMPASS
C          DIRECTIONS (+CLOCKWISE, -COUNTERCLOCKWISE) 
C 
C INPUTS FROM COMMON BLOCK SIMA - - - - - - 
C REFLG    GROUND REFLECTANCE FOR EA. MONTH OF YR.
C
C SOLARH   CLOUD COVER(0-10), IF ICC> 0 
C            OR TOTAL HORIZ RADIATION, IF ICC=0, BTU/HR-SQ FT
C  OUTPUTS - - - - - - - -
C   ISUN  =1,NO SUN;  =0, SUN SUFFICIENTLY ABOVE HORIZON
C   ITH  TOTAL HORIZONTAL SOLAR RADIATION, BTU/HR-SQ FT 
C        IF ICC=1, ITH IS CLOUDLESS VALUE 
C        IF ICC=0, ITH IS MEASURED VALUE FROM CC ARRAY
C  IDHC  DIRECT SOLAR RADIATION ON HORIZONTAL SURFACE WITH
C        CLOUDY SKY, BTU/HR-SQ FT 
C  SAZM  SOLAR AZIMUTH ANGLE, RADIANS 
C  SALT  SOALR ALTITUDE, RADIANS
C 
C  VALUES FOR EACH COMPASS DIRECTION ARE GIVEN BELOW  - - - - - - - 
C  IDC  DIRECT SOLAR RADIATION ONTO VERTICAL SURFACE FOR CLOUDY SKY,
C       BTU/HR-SQ FT
C  IDFC  DIFFUSE SOLAR RADIATION ONTO VERTICAL SURFACE FOR
C        CLOUDY SKY, BTU/HR-SQ FT 
C  DIRECT  DIRECT SOLAR TRANSMITTED THROUGH WINDOWS / RATIO OF
C          UNSHADED(BY OVERHANG AND REVEAL) WINDOW AREA TO TOTAL
C          WINDOW AREA ,BTU/HR-SQ FT
C  DIFFUS  DIFFUSE SOLAR RADIATION TRANSMITTED THROUGH WINDOWS
C          (INCLUDES GROUND REFLECTANCE), BTU/HR-SQ FT
C 
C  SFST  SHADING-FROM-SIDE FACTOR (SHADING FROM REVEAL) 
C  SFTT  SHADING FROM TOP FACTOR (SHADING FROM OVERHANG AND REVEAL) 
C 
C 
CMDK NRF
CMDK ROOF
CMDK NCR
CMDK SIMA
C 
      REAL IDN,IDH,IDFH,ITH,IDHC,IDC,IDFHC,IDFC 
     + ,IZERO,IDHBS 
      CHARACTER *4 STM
      DIMENSION AJ(6),DIFFUS(4),DIRECT(4),IDC(4),IDFC(4),NDY21(12),
     + PS(12), QS(12),RS(12),SFST(4),SFTT(4),TJ(6)
C 
      DATA  RAD/57.29577/ 
      DATA  TJ/-0.00885, 2.71235,-0.62062,-7.07329, 9.75995,-3.89922/ 
      DATA  AJ/ 0.01154, 0.77674,-3.94657, 8.57881,-8.38135, 3.01188/ 
      DATA  PS/1.14,1.06,1.06,1.06,.96,.96,.96,.95,.95,.95,1.14,1.14/ 
      DATA  QS/.003,.012,.012,.012,.033,.033,.033,.030,.030,
     $         .030,.003,.003/
      DATA  RS/-.0082,-.0084,-.0084,-.0084,-.0106,-.0106,-.0106,
     $         -.0108,-.0108,-.0108,-.0082,-.0082/
      DATA  PI/3.14159/ 
C ***** SOLAR CONSTANT, BTU/FT2-HR
      DATA SC/428./ 
C  NDY21 IS DAY OF YR AT 21ST OF EACH MONTH 
      DATA NDY21/21,52,80,111,141,172,202,233,264,
     +294,325,355/
C 
      IDBG=0
      ITM=INT(TIME+1.001) 
C ***** TIME = TIME OF DAY (0 - 24),HR
      IF (TIME.GT.23.98) GO TO 9996 
C  SET SUN INDICATOR FOR SUNSHINE 
      ISUN=0
C  ZERO OUTPUTS FOR NO SUN
      ITH=0.
      IDHC=0. 
      SAZM=0. 
      SALT=0. 
      DO 10 I=1,4 
      IDC(I)=0. 
      IDFC(I)=0.
      DIRECT(I)=0.
      DIFFUS(I)=0.
      SFST(I)=0.
      SFTT(I)=0.
   10 CONTINUE
      DO 15 I=1,NRF 
      IDR(I)=0. 
      IDIFR(I)=0. 
   15 CONTINUE
C 
      NDYR=NDYRI
C  IF NDYRI <= 0, USE DAY OF YEAR AT 21ST DAY OF KM MONTH 
      IF(NDYRI.EQ.-1)NDYR=NDY21(KM)
      CALL ETIME (NDYR,STM,ET,TZN)
C 
C **** CALCULATE SOLAR FACTORS
C ***** D = DECLINATION  (TRNSYS) 
C ***** SEE "ALGORITHMS FOR BUILDING HEAT TRANSFER SUBROUTINES" 
C *****   P.17 ASHRAE 1976
      D=23.45*SIN(360.*(284.+NDYR)/(365.*RAD))/RAD
      XIM=1.0 
      DO 44 IXM=1,12
      IF (NDYR.LE.NDY21(IXM))GO TO 46 
   44 CONTINUE
   46 XIM=FLOAT(KM)
      ABC=COS(15.*(2.*XIM-1.)/RAD)
      A=367.5+22.5*ABC
      B=0.174-.032*ABC
      C=.0965-.0385*ABC 
C 
      DNL=DNLAT/RAD 
C ***** SUNRISE HOUR ANGLE (HPRM), DEGREES FROM SOLAR NOON
      HPRM=ACOS(-TAN(DNL)*TAN(D)) 
      HLIM=ABS(HPRM)-15./RAD
      H=15.*(TIME-12.+TZN+ET-DST)-DLONG 
      HRAD=H/RAD
      IF(IDBG.NE.0)WRITE(60,500)ITM,TIME,HPRM,HLIM,H,HRAD
      IF (ABS(HRAD).LT.ABS(HLIM)) GO TO 50
      ISUN=1
      GO TO 9996
C *** CALCULATE DIRECTION COSINES OF SOLAR BEAM 
C ***** Z = ZENITH ANGLE, ANGLE BETWEEN BEAM FROM SUN AND VERTICAL
50    COSZ=SIN(DNL)*SIN(D)+COS(DNL)*COS(D)*COS(HRAD)
      COSW=COS(D)*SIN(HRAD) 
      COSS=(1-COSZ**2-COSW**2)**.5
      COSH=COS(HRAD)
      TANDL=TAN(D)/TAN(DNL) 
      ISGN=1
      IF (COSH.LT.TANDL) ISGN=-1
      COSS=ISGN*COSS
C 
C *** CALCULATE SOLAR ALTITUDE ANGLE
      SALT=ASIN(COSZ) 
C *** CALCULATE SOLAR AZIMUTH ANGLE 
      IF (COSS.LT.0.) GO TO 52
      SAZM=ASIN(COSW/COS(SALT)) 
      GO TO 54
52    SAZM=PI-ASIN(COSW/COS(SALT))
C 
C *** CALCULATE DIRECT NORMAL,SKY DIFFUSE AND GROUND REFLECTED RADIATION
C *** IDN IS DIRECT NORMAL SOLAR RADIATION ON HORIZONTAL
C *** BS IS DIFFUSE SOLAR RADIATION ON HORIZONTAL SURFACE 
C *** BG IS PORTION OF TOTAL SOLAR RADIATION REFLECTED TO 
C ***              VERTICAL WALL
   54 IF(ICC.GT.0)GO TO 56
C ***** TRNSYS MANUAL (TYPE 16 UNIT)
      ITH=SOLARH(ITM)
      IF(ITH.EQ.1.1E10)THEN 
         WRITE(60,504) 
         STOP 'SOLR:SOLAR BAD, SEE NOTE ON OUTPUT'
         ENDIF
      IF(IDBG.NE.0)WRITE(60,503)IDAY,ITM,KM,ITH
      IZERO=SC*COSZ*(1.+.033*COS(360.*NDYR/(365.*RAD))) 
      XKT=ITH/IZERO 
      RAT=1.0045+XKT*(0.04349+XKT*((-3.5227)+2.6313*XKT)) 
      IF(XKT.LT.0.0) RAT=1.0
      IF(XKT.GT.0.8) RAT=0.165
      IDHC=ITH*(1.0-RAT)
      GO TO 60
56    CLOUDC=SOLARH(ITM)
      IDN=A*CLEAR*EXP(-B/COSZ)
      BS=C*IDN/(CLEAR)**2 
C      BG=REFLG(KM)*(BS+IDN*COSZ)
C 
C *** CALCULATE CLOUDY SKY RADIATION
      X=SIN(SALT) 
      Y1=.309-.137*X+.394*X**2
      XKSOL=X/(C+X)+(PS(KM)-1.)/(1.-Y1) 
      IDH=IDN*COSZ
      IDFH=BS 
      ITH=IDH+IDFH
      IDHC=ITH*XKSOL*(1.-CLOUDC/10.)
C 
C 
C    M= COMPASS DIRECTION NO. 
   60 DO 220 M=1,4
C 
C 
C *** CALCULATE CLOUDY SKY RADIATION FOR WALL 
      CALL SHADE (M,SAZM,SALT,WROTA,SFS,SFT,WAZ)
      SFST(M)=SFS 
      SFTT(M)=SFT 
C 
C---- WAZ IS WALL AZIMUTH ANGLE FROM SUN POSITION TO HORIZ. PLANE 
C     COSINC IS SAME AS THAT IN EQ.5B IN ASHRAE F27.3 1981
      COSINC=COS(SALT)*COS(WAZ/RAD) 
      CT=COSINC 
      VAL = 1.
      IF (CT.LE.0.) VAL=0.
      IDC(M)=VAL*IDHC*CT/COSZ 
         Y2=0.55+CT*(0.437+CT*0.313)
         IF (CT.LE. -0.2) Y2=.45
      IF (ICC.GT.0)GO TO 92 
      IDFHC=ITH-IDHC
      IDFC(M)= Y2*IDFHC+0.5*REFLG(KM)*ITH 
      GO TO 94
   92 CONTINUE
C *** IDF = DIFFUSE RADIATION FOR CLOUDLESS SKY ON VERTICAL SURFACE 
C      IDF=BS*Y2+BG/2. 
      CCF=PS(KM)+QS(KM)*CLOUDC+RS(KM)*CLOUDC**2 
      IDFHC=ITH*CCF-IDHC
      IDHBS=COSZ*CLEAR**2/C 
      IDFC(M)=(Y2+0.5*REFLG(KM)*(1.+IDHBS))*IDFHC 
C *** IDC=DIRECT
C *** IDFC=DIFFUSE
   94 CONTINUE
C 
C     CALC INSTANTANEOUS HEAT TRANSMISSION THRU WINDOWS 
      TDIR=0. 
      TDIF=0. 
      ADIR=0. 
      ADIF=0. 
      DO 219 K=1,6
      J=K-1 
      TDIR=TDIR+TJ(K)*CT**J 
      TDIF=TDIF+2.*TJ(K)/(J+2)
      ADIR=ADIR+AJ(K)*CT**J 
219   ADIF=ADIF+ 2.0*AJ(K)/(J+2)
      DIRECT(M)=IDC(M)*(TDIR+.267*ADIR) 
      IF(ICC.LE.0) THEN 
          DIFFUS(M)=0.5*(IDFHC+REFLG(KM)*ITH)*(TDIF+0.267*ADIF) 
      ELSE
          DIFFUS(M)=0.5*(1.+REFLG(KM)*(1.+IDHBS))*(TDIF+.267*ADIF)*IDFHC
      ENDIF 
  220 CONTINUE
C 
C---- CALCULATE SOLAR LOADING ON ROOF 
C 
      DO 360 I=1,NSECTA 
      M=IDEXPR(I) 
      PA   =PAA   (I) 
C 
C *** DETERMINE SOLAR AZIMUTH FOR EACH ROOF SURFACE 
      IF (M.EQ.1) RAZ=ABS(180.+SAZM*RAD-WROTA)
      IF (M.EQ.2) RAZ=ABS( 90.+SAZM*RAD-WROTA)
      IF (M.EQ.3) RAZ=ABS(     SAZM*RAD-WROTA)
      IF (M.EQ.4) RAZ=ABS( 90.+WROTA-SAZM*RAD)
      PARAD=PA/RAD
      CSINC=COS(SALT)*COS(RAZ/RAD)*COS(PARAD)+
     $       SIN(SALT)*SIN(PARAD) 
      IF (CSINC.LE.0.) GO TO 346
      IDR(I)=IDHC*CSINC/COSZ
      GO TO 347 
346   IDR(I)=0. 
347   IDIFR(I)=IDFHC*((1.+SIN(PARAD))/2.) 
      IF(IDBG.NE.0)THEN 
      WRITE(60,502)I,M,PA,RAZ,CSINC,IDR(I),IDHC,IDIFR(I),IDFHC 
          ENDIF 
360   CONTINUE
9996  CONTINUE
          IF(IDBG.NE.0)THEN 
      WRITE(60,501)ITH,IDHC,SAZM,SALT,IDC,IDFC,DIRECT,DIFFUS,SFST,SFTT 
          ENDIF 
      RETURN
  500 FORMAT(1X,'SOLR: ITM,TIME,HPRM,HLIM,H,HRAD= ',I5,5G13.5)
  501 FORMAT(1X,'SOLR OUTPUTS:ITH,IDHC,SAZM,SALT= ',4G13.5/ 
     +' IDC= ',4G13.5/' IDFC= ',4G13.5/ 
     +' DIRECT= ',4G13.5/' DIFFUS= ',4G13.5/
     +' SFST= ',4G13.5/' SFTT= ',4G13.5)
  502 FORMAT(1X,'I,M,PA,RAZ,CSINC= ',2I3,3G13.5/
     +3X,'IDR,IDHC,IDIFR,IDFHC= ',4G13.5) 
  503 FORMAT(1X,'SOLR:IDAY,ITM,KM,ITH= ',3I5,G13.5) 
  504 FORMAT(1X,'SOLR:SOLAR BAD, CHECK WETHER DATA OR VALUE OF MONTH',
     +' READ IN STRRUN')
      END 
      SUBROUTINE ETIME(NDYR,STM,ET,TZN) 
C  CALC EQUATION OF TIME
C      REAL LSTM(6)
      CHARACTER *4 TMER,STM 
      DIMENSION E(10),TMER(6),TZNS(6) 
      DATA E/-4.319780295E+00 , -6.560560432E-03 , -1.768353820E-05 , 
     +        1.880630945E-08 , -3.612378625E-12 , -1.243267535E-15 , 
     +        6.119539396E-19 , -9.877736221E-23 ,  7.200722198E-27 , 
     +       -2.009408139E-31/
C 
C *** INPUTS
C     NDYR  DAY OF YEAR(NON-LEAP YEAR), 1-365 
C     STM  TIME ZONE (EST,CST,MST,PST,YST,ANST) 
C 
C *** OUTPUTS 
C     ET  EQN OF TIME, HR (SEE SOLAR ENERGY BY DUFFIE AND BECKMAN)
C     TZN  TIME ZONE NO., 1-6 
  
C    EQN OF TIME VALUES ARE FOR THE 21ST DAY OF MONTH 
      DATA TMER/'EST','CST','MST','PST','YST','AHST'/ 
C      DATA LSTM/75.,90.,105.,120.,135.,150./
      DATA TZNS/5.,6.,7.,8.,9.,10./ 
      DO 10 I=1,6 
        IF (STM.EQ.TMER(I)) GO TO 20
10      CONTINUE
      WRITE(60,15)
15    FORMAT ('0ERROR IN INPUT FOR STANDARD TIME MERIDIAN') 
      STOP 'ETIME:ERROR IN STM FOR STD TIME ZONE' 
20    CONTINUE
C      LONGM=LSTM(I) 
C  EQUA OF TIME IS ON P. 19 OF DUFFIE AND BECKMAN 
      T=(NDYR-0.5)*24.
      ET=E(1)+T*(E(2)+T*(E(3)+T*(E(4)+T*(E(5)+T*(E(6)+
     +T*(E(7)+T*(E(8)+T*(E(9)+T*E(10))))))))) 
C  ET IS IN MINUTES SO CONVERT TO HRS 
      ET=ET/60. 
      TZN=TZNS(I) 
      RETURN
      END 
      SUBROUTINE SHADE (M,SAZM,SALT,WA,SFS,SFT,WAZ) 
C  CALCULATE WINDOW SHADING FACTORS 
      DATA RAD/57.29577/
  
C 
      SAZ=SAZM*RAD
      IF (M.EQ.1) WAZ=ABS(180.+SAZ-WA)
      IF (M.EQ.2) WAZ=ABS( 90.+SAZ-WA)
      IF (M.EQ.3) WAZ=ABS(     SAZ-WA)
      IF (M.EQ.4) WAZ=ABS( 90.+WA-SAZ)
      IF ((SALT.EQ. 0.0 .AND. (SAZM.EQ. 0.0)))GO TO 20
      IF (WAZ.GE.90.) GO TO 10
      WAZRAD=WAZ/RAD
      SFS=TAN(WAZRAD) 
      SFT=TAN(SALT)/COS(WAZRAD) 
      GO TO 20
10    SFS=0.0 
      SFT=0.0 
20    RETURN
      END 
