C March 1, 2008 C*********************************************************************** C** ** C** PROGRAM AtticSim II ** C** ** C*********************************************************************** C AtticSim II adds several options ---- C C C ATICSIM11 adds the option for input of hourly values for indoor air C temperature TI instead of average only. Input KFLAG(18)=1 to use C option. If KFLAG(18)=0 is entered in input file put correct average C value TIAVG in duct input. Hourly values will be set to it. C ATICSIM10 adds the option for input of hourly values for duct system C fraction ONTIME instead of average only. Input KFLAG(17)=1 to use C option. If KFLAG(17)=0 is entered in input file put correct average C value ONTIMEAVG in duct input. Hourly values will be set to it. C ATICSIM9 adds the option for input of ventilation air temperatures C different from outside air temperature if input KFLAG(16)=1. C If KFLAG(16)=0, ventilation air temperature defaults to outside air C temperature. C ATICSIM8 adds more output to the default output file and a header C line to identify current outputs; the limit that duct sections be C longer than 1 ft is changed to longer than 0 ft to allow short C dummy duct sections if needed C C THIS PROGRAM CALCULATES THE HEAT TRANSFER THROUGH THE CEILING C BELOW A GABLED ATTIC HAVING A FIVE-SIDED CROSS SECTION C ATICSIM5 HAS OPTIONS FOR TRUSSES AND ATTIC DUCTS C STEADY-STATE AND TRANSIENT DUCT OPTIONS ARE AVAILABLE C ATICSIM5 USES UNDERRELAXATION TO CALCULATE NEW TEMPERTURES FOR C TRUSSES C USED WITH WEATHER DATA PROCESSED FROM WEATHER TAPES C BY THE DOE-2 PREPROCESSOR C THIS VERSION LIMITS THE MOISTURE CONTENT OF WOOD TO 30 PERCENT C ASSUMES THE AIR WITHIN THE ATTIC IS NOT WELL-MIXED C INCLUDES APPROXIMATE CALCULATIONS FOR LATENT HEAT EFFECTS C DUE TO MOISTURE ADSORPTION AND DESORPTION C SURFACE TEMPERATURES ON INSIDE OR OUTSIDE, OR VENTILATION RATE C MAY BE SPECIFIED THROUGH KFLAG(I) = 1 C NOMENCLATURE: C N(I) = NUMBER OF CONDUCTION TRANSER FUNCTIONS FOR SURFACE I C U(I) = SURFACE-TO-SURFACE U-VALUE FOR SURFACE I C CR(I) = COMMON RATIO OF RESPONSE FACTORS FOR SURFACE I C BETA(I) = TEMPERATURE DEPENDENT COEFFICIENT FOR CONDUCTION C TRANSFER FUNCTIONS FOR SURFACE I C X(I,J) C Y(I,J) = X, Y, & Z CONDUCTION TRANSFER FUNCTIONS FOR SURFACE I C Z(I,J) C C SURFACES ARE NUMBERED AS (FOR NORTH-SOUTH RIDGE): C 1 = CEILING C 2 = EAST ROOF or EAST Roof Deck C 3 = WEST ROOF or West Roof Deck C 4 = SOUTH GABLE C 5 = NORTH GABLE C 6 = EAST EAVE WALL C 7 = WEST EAVE WALL C 8 = EAST ROOF COVER with ASV C 9 = WEST ROOF COVER with ASV C ALF(I) = SOLAR ABSORPTANCE FOR SURFACE I C EO(I) = INFRARED EMITTANCE FOR OUTSIDE OF SURFACE I C EI(I) = INFRARED EMITTANCE FOR INSIDE OF SURFACE I C TIS(I,J) = INSIDE SURFACE TEMPERATURE FOR SURFACE I,TIME J=1 or 2 C TOS(I,J) = OUTSIDE SURFACE TEMPERATURE FOR SURFACE I,TIME J=1 or 2 C TA = AVERAGE TEMPERATURE OF AIR IN ATTIC C QI(I,J) = INSIDE HEAT FLUX FOR SURFACE I,TIME J=1 or 2 C QO(I,J) = OUTSIDE HEAT FLUX FOR SURFACE I,TIME J=1 or 2 C HCI(I) = CONVECTION COEFFICIENT FOR INSIDE OF SURFACE I C HR(I,K) = RADIATION INTERCHANGE COEFFICIENT FROM SURFACE I C TO SURFACE K C AA(I,J) = MATRIX OF COEFFICIENTS FOR SIMULTANEOUS EQUATIONS C BB(I) = KNOWN VECTOR FOR SIMULTANEOUS EQUATIONS C XXX(I) = SOLUTION VECTOR FOR SIMULTANEOUS EQUATIONS C QSOL(I) = SOLAR RADIATION INCIDENT ON SURFACE I C A(I) = AREA OF SURFACE I C G(I,K) = OVERALL RADIATION VIEW FACTOR FROM SURFACE I TO C SURFACE K C HCO(I) = CONVECTION COEFFICIENT ON OUTSIDE OF SURFACE I C HRO(I) = RADIATION COEFFICIENT ON OUTSIDE OF SURFACE I C R(I) = SURFACE-TO-SURFACE RESISTANCE FOR SURFACE I C AMC(I) = MOISTURE CONTENT (WEIGHT FRACTION) FOR SURFACE I C PERM(I) = MOISTURE PERMEANCE OF SURFACE I C AWRAT(I) = RATIO OF EXPOSED SURFACE AREA TO PROJECTED SURFACE C AREA FOR SURFACE I C AMASS(I) = MASS PER UNIT AREA OF SURFACE I THAT PARTICIPATES IN C MOISTURE EXCHANGE C AMW(I) = MOISTURE FLUX TOWARD SURFACE I C BMW(I) = TERM IN TAYLOR SERIES EXPANSION FOR MOISTURE FLUX C TS(I) = TEMPERATURE OF INSIDE OF SURFACE I C LFLAG(1) = FLAG FOR ABBOVE SHEATHING VENTILATION (ASV) C = 0 implies conventional roof C = 1 implies ASV WITH 2 ADDITIONAL OUTSIDE ROOF SURFACES C = 2 implies ASV WITH 4 ADDITIONAL OUTSIDE ROOF SURFACES C LFLAG(2) = FLAG FOR ROUTINE TO CALCULATE AIRFLOW FOR ASV C 0 implies KKW Vent Subroutine used also for attic C 1 implies modified Ra Number correlation (see MODRA) C LFLAG(3) = FLAG Signaling a closed natural convection cavity C KFLAG = FLAG FOR SPECIFIED SURFACE TEMPERATURES C 0 = PARAMETER NOT SPECIFIED C 1 = PARAMETER SPECIFIED C I = 1, 7 FOR TIS(I,1) C I = 8,14 FOR TOS(I,1) C I = 15 FOR VENTILATION RATE C I = 16 FOR VENTILATION AIR TEMPERATURE C I = 17 FOR DUCT SYSTEM ONTIME C I = 18 FOR INDOOR AIR TEMPERATURE C I = 19,20 FOR TIS(I,1) ASV Roof C I = 21,22 FOR TOS(I,1) ASV Roof C AL = LENGTH OF ATTIC, FEET C W = WIDTH OF ATTIC, FEET C PITCH1 = PITCH OF EAST ROOF, DEGREES C PITCH2 = PITCH OF WEST ROOF, DEGREES C ORIENT = ORIENTATION OF HOUSE C H1 = HEIGHT OF VERTICAL WALLS AT EAVES, FEET C GAP = HEIGHT OF ASV CHANNEL NORMAL TO ROOF DECK, FEET C AI = AREA OF INLET VENT, SQUARE FEET C AO = AREA OF OUTLET VENT, SQUARE FEET C ITYPE = TYPE VENTS C 1 = SOFFIT AND RIDGE VENTS C 2 = SOFFIT AND GABLE VENTS C 3 = SOFFIT VENTS ONLY C QLAT = HEAT OF VAPORIZATION OF WATER, BTU/LBM C EXFIL = RATE OF FLOW OF AIR FROM HOUSE TO ATTIC, LB/HR C AL1 ... AL7 = CHARACTERISTIC LENGTHS OF SURFACES 1 ... 7, FEET C TO = OUTDOOR TEMPERATURE, F C QSOLH = SOLAR RADIATION INCIDENT ON HORIZONTAL SURFACE C BTU/(HR.-SQ.FT.) C WS = WIND SPEED, MPH C DIR = WIND DIRECTION C TI = INDOOR TEMPERATURE, F C HUMI = INDOOR RELATIVE HUMIDITY, PERCENT C PATM = AMOSPHERIC PRESSURE (Psia) C VDOT = VENTILATION RATE, CUBIC FEET PER HOUR C IVFLAG = FLAG FOR CALCULATING FLOW VELOCITY C 0 = FLOW THROUGH BOTH SIDE OF ATTIC C 1 = FLOW THROUGH ONE SIDE OF ATTIC C FLUX = MEASURED CEILING HEAT FLUX, BTU/(HR.-SQ.FT.) C TAIR = MEASURED ATTIC AIR TEMPERATURE, F C TEXIT = MEASURED EXIT AIR TEMPERATURE, F C H = HEIGHT, FEET C AMDOT = VENTILATION RATE, POUNDS PER HOUR C AMCDOT = VENTILATION RATE TIMES SPECIFIC HEAT C V = CRUDE ESTIMATE OF AIR SPEED THROUGH ATTIC, FEET PER HOUR C ACH = AIR CHANGE RATE FOR ATTIC, AIR CHANGES PER HOUR C NEXT LINE ADDED FOR VERSION 1A C IMPLICIT REAL*8 (A-H, O-Z) PARAMETER (Nsrf = 9, IMx = 100, Max = 900, NsrfTm = 18, & Nequ = 21, Neqmx = 441, Nduct = 100) REAL*8 MDOTIN,MDOTLK,MLEAKS,NUS(Nsrf) INTEGER*4 ICLDTY C C COMMON /CLOUD/CLDAMT,CC,ICLDTY COMMON /SUND/GUNDOG,HORANG,TDECLN,EQTIME,SOLCON, 1 ATMEXT,SKYDFF,RAYCOS(3),RDN, 2 BSUN,DECLN,CD,SD,FSUNUP,ISUNUP COMMON /SUN2/RDNCC,BSCC COMMON /HUMID/WO,PATM C Character*2 Blank Character*7 System Character*12 City C Integer CLZone C DIMENSION LFLAG(3),HM(10) DIMENSION N(Nsrf),U(Nsrf),CR(Nsrf),BETA(Nsrf) DIMENSION X(Nsrf,IMx),Y(Nsrf,IMx),Z(Nsrf,IMx) DIMENSION ALF(Nsrf),EO(Nsrf),EI(Nsrf) DIMENSION AMC(Nsrf),PERM(Nsrf),AWRAT(Nsrf),AMASS(Nsrf), & AMW(Nsrf),BMW(Nsrf) DIMENSION HCI(Nsrf),HR(Nsrf,Nsrf),HCO(Nsrf),HRO(Nsrf) DIMENSION TIS(Nsrf,IMx),TOS(Nsrf,IMx),TS(Nsrf),TA(2),TASV(2,2) DIMENSION QSOL(Nsrf),QI(Nsrf,2),QO(Nsrf,2) DIMENSION XXX(Nequ),AA(Nequ,Nequ),BB(Nequ) C DIMENSION A(Nsrf),F(7,7),G(32,32),R(Nsrf) C Next line added for model with ducts COMMON /DUCT1/ALD(25),xMDOTin(25,2),xMDOTlk(25,2),EODUCT(25), & ZONE(25),PID(25),PO(25),RD(25),xTDUCTin(2),RHOC(25),RDO(25), & RDI(25),NSUPLY,NRETRN,ND(25),NINLET(25),NEXIT(25) COMMON /DUCT2/DH(25),DCHAR(25),AD(25),ISHAPE(25) COMMON /DUCT3/TI,V,ONTIME,DuctTin,DuctMin(25), & DuctLkg(25),KFLAG(19),NTOT COMMON /DUCT4/QCON,QRAD(8),QRADtot,MLEAKS,QLEAKS,QDUCTS,QSTOR, & MDOTin(25),MDOTlk(25) COMMON /DUCT5/DX(25),Rduct,NSEC(25) C Next line added for model with trusses COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS COMMON /TRUS2/AMASST,AMCT,HCTRUS,AMWT,BMWT,QLAT DIMENSION T(27) DIMENSION TAD(25,Nduct),TWWD(25,Nduct,2),TWOD(25,Nduct) DATA TAD/2500*75.0D0/,TWWD/5000*75.0D0/,TWOD/2500*75.0D0/ C Next line modified 2-12-96 to initialize TTRNEW DATA TTRUS/2*75./,TTRNEW/75./ DATA QRAD/8*0.0D0/,QRADT/7*0.0D0/ C DATA BLANK / ' ' / DATA KHEADER / 0 / ! PRINT FLAG for OUTPUT FILE HEADER C DATA AA/ Neqmx*0.0D0/, BB/Nequ*0.0D0/ DATA TIS/ Max*75.0D0/, TOS/Max*75.0D0/ DATA TA/ 2*75.0D0/, TASV/4*75.0D0/, TR/75.0D0/ DATA QI/NsrfTm*0.0D0/, QO/NsrfTm*0.0D0/, BMW/Nsrf*0.0D0/ C C OPEN DATA FILES OPEN(UNIT= 5,FILE= & 'IN') OPEN(UNIT=10,FILE='WEA') OPEN(UNIT=11,FILE='OUT') OPEN(UNIT=15,FILE='COEF') C C READ IN LFLAGS READ (5,*) (LFLAG(I),I=1,3),IVFLAG,ISURF,Relax C IVFLAG = FLAG FOR CALCULATING FLOW VELOCITY C 0 = FLOW THROUGH BOTH SIDES OF ATTIC C 1 = FLOW THROUGH ONLY ONE SIDE OF ATTIC IF(LFLAG(1) .EQ. 0) then Nsurf = Nsrf - 2 ! No ASV reduces surfaces to 7 Narray = Nequ - 6 ! No ASV reduces solution to 15 Equations System = 'No ASV' ELSE IF (LFLAG(1) .EQ. 1 .and. IVFLAG .eq. 1) then Nsurf = Nsrf - 1 Narray = Nequ - 1 System = 'Yes ASV' else Nsurf = Nsrf Narray = Nequ System = 'Yes ASV' END IF C C READ IN KFLAGS READ (5,*) (KFLAG(I),I=1,7) READ (5,*) (KFLAG(I),I=8,14) READ (5,*) (KFLAG(I),I=15,19) JFLAG = 0 DO 1 I = 1,18 IF(KFLAG(I).EQ.0) GO TO 1 JFLAG = 1 1 CONTINUE C********** THIS INCLUDES LEAKS IN THE DUCT IF(JFLAG.NE.0) OPEN(UNIT=12, & FILE='BCD') C C READ IN CONDUCTION TRANSFER FUNCTIONS DO 100 I = 1,Nsurf READ (5,*) N(I),U(I),CR(I),BETA(I) M = N(I) DO 100 J = 1,M READ (5,*) X(I,J),Y(I,J),Z(I,J) 100 CONTINUE DO 101 I = 1,Nsurf 101 R(I) = 1.0D0/U(I) C CONVENTION FOR CTF'S IS SAME AS USED IN TARP. C CONVENTION FOR ATTIC MODEL HAS THE INSIDE SURFACES FACING THE ATTIC SPACE C C CTF'S FOR THE CEILING ARE CALCULATED FOR USE IN A WHOLE HOUSE MODEL, C WHERE THE OUTSIDE SURFACE FACES THE ATTIC. THEREFORE, FOR USE INTERNAL C TO THE ATTIC MODEL, THE CODE INTERCHANGES THE X'S AND THE Z'S M = N(1) DO 110 J = 1,M XX = X(1,J) X(1,J) = Z(1,J) Z(1,J) = XX 110 CONTINUE C READ IN SOLAR ABSORPTANCES OF OUTSIDE SURFACES AND C EMITTANCES OF BOTH OUTSIDE AND INSIDE SURFACES READ(5,*) (ALF(I), I = 1,Nsurf) READ(5,*) (EO(I), I = 1,Nsurf) READ(5,*) (EI(I), I = 1,Nsurf) C READ IN ATTIC GEOMETRY: LENGTH, WIDTH, PITCHES OF ROOF C SURFACES(DEGREES), ORIENTATION ANGLE OF HOUSE C HEIGHT OF SIDE WALLS AND GAP FOR ASV READ (5,*) AL,W,PITCH1,PITCH2,ORIENT,H1,GAP C READ VENT INLET AND OUTLET AREAS, AND VENT TYPE C TYPE 1 = SOFFIT AND RIDGE VENTS C TYPE 2 = SOFFIT AND GABLE VENTS C TYPE 3 = SOFFIT VENTS ONLY READ(5,*) AI,AO,ITYPE C READ WATER VAPOR PERMEANCES OF ATTIC SURFACES C PERM = GRAINS PER (HR-FT2-INCH HG) READ (5,*) (PERM(I),I=1,Nsurf) C CONVERT PERM VALUES TO POUNDS PER (HR-FT2-PSIA) DO 120 I = 1,Nsurf 120 PERM(I) = PERM(I)*(29.921/14.696)/7000. C READ RATIO OF TOTAL AREA OF EXPOSED WOOD TO GEOMETRICAL C PROJECTION OF SURFACE AREAS READ (5,*) (AWRAT(I),I=1,Nsurf) C READ WEIGHT OF WOOD PER UNIT PROJECTED AREA FOR EACH SURFACE READ (5,*) (AMASS(I),I=1,Nsurf) C READ INITIAL MOISTURE CONTENTS OF WOOD, WEIGHT FRACTION READ (5,*) (AMC(I),I=1,Nsurf) C READ LATENT HEAT OF VAPORIZATION, BTU/LB READ (5,*) QLAT C READ RATE OF FLOW OF HOUSE AIR INTO ATTIC, POUNDS PER HOUR READ (5,*) EXFIL C NEXT SECTION ADDED FOR MODEL WITH DUCTS C Read flag (IDUCT) for ducts; 0 for no duct, 1 for ducts C read flag (IDMOD) for duct model: 0 for steady-state, 1 for transient READ(5,*) IDUCT,IDMOD IF(IDUCT.NE.0) CALL DUCTREAD CONTINUE C NEXT SECTION ADDED FOR MODEL WITH TRUSSES C Read flag (ITRUS) for trusses: 0 for no trusses, 1 for trusses READ(5,*) ITRUS IF(ITRUS.NE.0) CALL TRUSREAD C NEXT SECTION FOR INDOOR AMBIENT, ATTIC VENTILATION AND COMPRESSOR ON-TIME READ(5,*) TIAVG,HUMI,ONTIMEAVG C C*********************************************************************** C CALCULATE CHARACTERISTIC LENGTHS AND AREAS OF ATTIC SURFACES C CHARACTERISTIC LENGTHS ARE: C CEILING - AVERAGE OF LENGTH AND WIDTH C ROOF - DISTANCE FROM EAVE TO RIDGE C GABLE - AVERAGE HEIGHT C EAVE SIDE WALLS - HEIGHT AL1 = (AL+W)/2. P1 = PITCH1*3.14159265/180.D0 P2 = PITCH2*3.14159265/180.D0 P3 = 3.14159265 - P1 - P2 IP1 = PITCH1 IP2 = PITCH2 IP3 = 180.0D0 - PITCH1 - PITCH2 IF(IP3.EQ.180.0D0) GO TO 5 AL2 = W*DSIN(P2)/DSIN(P3) ! Typically P1 = P2 (If shed type roof P2 = 90) AL3 = W*DSIN(P1)/DSIN(P3) AL4 = 0.5*(AL2*DSIN(P1)) + H1 GO TO 6 5 AL2 = W/2.0D0 AL3 = W/2.0D0 AL4 = H1 6 AL5 = AL4 AL6 = H1 AL7 = H1 if(LFLAG(1) .eq. 1) then AL8 = AL2 AL9 = AL3 end if A(1) = AL*W ! AREAS OF RESPECTIVE ATTIC SURFACES A(2) = AL*AL2 A(3) = AL*AL3 A(4) = W*AL4 A(5) = A(4) A(6) = AL*H1 A(7) = AL*H1 if(LFLAG(1) .eq. 1) then A(8) = AL*AL8 A(9) = AL*AL9 end if C NEXT 7 LINES ADDED FOR MODEL WITH DUCTS AND TRUSSES IF(IDUCT.EQ.0.AND.ITRUS.EQ.0) THEN CALL VIEW2(AL,W,PITCH1,PITCH2,H1,EI,G,7,F) ELSE IF(IDUCT.NE.0.AND.ITRUS.EQ.0) THEN CALL VIEW2D(AL,W,PITCH1,PITCH2,H1,EI,G,7,A) ELSE IF(ITRUS.NE.0) THEN CALL VIEW2T(AL,W,PITCH1,PITCH2,H1,EI,IDUCT,G) END IF C*********************************************************************** C READ LATITUDE, LONGITUDE, TIME ZONE, CLEARNESS NUMBER, C GROUND REFLECTANCE, AND FLAG FOR AMOUNT OF SOLAR DATA C TIME ZONE = 5 FOR EASTERN, = 6 FOR CENTRAL, ETC. C ISUN = 1, MEASURED TOTAL HORIZONTAL AND DIRECT AVAILABLE C ISUN = 2, MEASURED TOTAL HORIZONTAL ONLY AVAILABLE C ISUN = 3, NO MEASURED SOLAR, BUT MEASURED CLOUD COVER C NOTE, THIS VERSION EXPECTS ISUN = 1 C IN THIS VERSION CLRNES WILL BE READ IN EACH HOUR READ(10,*) STALAT,STALON,ITIMZ,RHOG,ISUN STALAT = STALAT*3.14159265/180. SSTALA = DSIN(STALAT) CSTALA = DCOS(STALAT) STALON = STALON*3.14159265/180. BAZIM = ORIENT*3.14159265/180. SBAZIM = DSIN(BAZIM) CBAZIM = DCOS(BAZIM) C IF (LFLAG(1) .eq. 0) then WRITE( 6,8) System,1.0D0-ALF(2),EO(2),EI(2),R(1),Rduct WRITE(11,9) System,1.0D0-ALF(2),EO(2),EI(2),R(1),Rduct WRITE(15,9) System,1.0D0-ALF(2),EO(2),EI(2),R(1),Rduct END IF IF (LFLAG(1) .eq. 1) then WRITE( 6,8) System,1.0D0-ALF(8),EO(8),EI(2),R(1),Rduct WRITE(11,9) System,1.0D0-ALF(8),EO(8),EI(2),R(1),Rduct WRITE(15,9) System,1.0D0-ALF(8),EO(8),EI(2),R(1),Rduct END IF 8 FORMAT(1x,'System,',A7,/, & 1x,'SR',F10.3,5x,'TE',f10.2,5x,'RB',f10.2,/, & 5x,'R-Value Ceiling = ',f10.3,/, & 5x,'R-value Duct = ',f10.1) 9 FORMAT(1x,'System,',A7,/, & 'SR,',F10.3,/,'TE,',f10.2,',',f10.2,/, & 'R-Ceiling,',f10.1,/,'R-Duct,',f10.1) C*********************************************************************** C READ LINES OF HOURLY WEATHER DATA C*********************************************************************** QIR = 0.0D0 10 READ(10,*,END=999) IDOY,IHR,TO,RHod,PATM,DIR,WS,QIR,QSOLH, & DIRSOL,DIFSOL,CLDAMT,CLRNES,ICLDTY C c write(6,1129) IDOY,IHR,TO,RHod,PATM,DIR,WS, c & QIR,QSOLH,DIRSOL,DIFSOL c1129 format(' DOY : ',2i3,/, c & ' TO : ',5f7.2,/, c & ' QIR : ',4f7.2) c write(6,1131) CLDAMT,CLRNES,ICLDTY c1131 format(' Cldamt ',1pD10.2,' Clrnew ',f10.2,' Icldty ',i4) c pause Call xMoist(TO,RHod,PATM,TodDwPt,WO,Pwod) if(WO .lt. 0.0001) WO = 0.00005 C Next section added for cases where boundary data are supplied IF(JFLAG.EQ.0) GO TO 11 READ (12,*,END = 999) (TIS(I,1),I=1,7),(TOS(I,1),I=1,7), & Tattic,TsoffS,TsoffN,Tridge,Ontime, & TI,RHi,S_OSB_RH,xN_OSB_RH Tvent = (TsoffS + TsoffN)/2.0 Call xMoist(TIS(2,1), S_OSB_RH,PATM,TsRfDwPt,W1,PwsRf) Call xMoist(TIS(3,1),xN_OSB_RH,PATM,TnRfDwPt,W2,PwnRf) c write(6,1129) (TIS(I,1),I=1,7),(TOS(I,1),I=1,7) c1129 format(' TIS : ',7F7.2,' TOS ',7F7.2) c pause c write(6,1131) Tattic,TsoffS,TsoffN,Tridge,Ontime, c & TI,RHi,S_OSB_RH,xN_OSB_RH c1131 format(' Tattic:',4f7.2,' Ontime',f7.2,' TI ',f7.2,' RHi ',3f7.2) c pause C RESTORE ESTIMATES FOR TEMPERATURES NOT SPECIFIED DO 15 I = 1,7 IF(KFLAG(I).EQ.0) TIS(I,1) = TIS(I,2) 15 IF(KFLAG(I+7).EQ.0) TOS(I,1) = TOS(I,2) 11 CONTINUE IF (KFLAG(16).EQ.0) TVENT = TO IF (KFLAG(17).EQ.0) ONTIME=ONTIMEAVG IF (KFLAG(18).EQ.0) TI=TIAVG C CONVERT WIND DIRECTION TO BUILDING FRAME OF REFERENCE DIR = ORIENT - DIR C*********************************************************************** IF(IHR.EQ.1) CALL SUN1(STALAT,IDOY) CALL SKY(TO,TodDwPt,IHR,QIR,CLDAMT,TSKY) TS2 = (((1.0D0 - DCOS(P1))/2.0D0*(TO+459.67)**4) & + ((1.0D0 + DCOS(P1))/2.0D0*(TSKY+459.67)**4))**0.25 - 459.67 TS3 = (((1.0D0 - DCOS(P2))/2.0D0*(TO+459.67)**4) & + ((1.0D0 + DCOS(P2))/2.0D0*(TSKY+459.67)**4))**0.25 - 459.67 TSV = (0.5*(TO+459.67)**4 + 0.5*(TSKY+459.67)**4)**0.25 - 459.67 IF(ISUN.NE.1) DIRSOL = 0.0D0 CALL WDTSUN(IHR,ITIMZ,STALON,SSTALA,CSTALA,SBAZIM,CBAZIM, & CLRNES,QSOLH,DIRSOL,ISUN) C CALCULATE SOLAR LOAD ON EACH SURFACE QSOL(1) = 0.0D0 ! Ceiling if (LFLAG(1) .eq. 0) then QSOL(2) = SUN3OR( 90.0D0,PITCH1,RHOG) ! East Facing Roof for E-W Ridge QSOL(3) = SUN3OR(270.0D0,PITCH2,RHOG) ! North Facing Roof for E-W Ridge else QSOL(2) = 0.0 QSOL(3) = 0.0 QSOL(8) = SUN3OR( 90.0D0,PITCH1,RHOG) QSOL(9) = SUN3OR(270.0D0,PITCH2,RHOG) end if QSOL(4) = SUN3OR(180.0D0,90.0D0,RHOG) ! West gable QSOL(5) = SUN3OR( 0.0D0,90.0D0,RHOG) ! East gable QSOL(6) = SUN3OR( 90.0D0,90.0D0,RHOG) ! South eave QSOL(7) = SUN3OR(270.0D0,90.0D0,RHOG) ! North eave C*********************************************************************** NIT = 1 ! Iterative loop per time step (limit = 100) C*********************************************************************** C CALCULATE ATTIC VENTILATION RATE AND FLOW VELOCITY H = AL4 + AL2*DSIN(P1)/2.0D0 20 CONTINUE IF(KFLAG(15).EQ.0) THEN CALL VENT(TVENT,TA(1),WS,DIR,H,AI,AO,ITYPE,AMCDOT,AMDOT,VDOT) ELSE C ASSUME VENTILATION RATES ARE BASED ON VENTILATION AIR TEMPERATURE AMDOT = VDOT*22.0493/((TVENT+459.67)/1.8) AMCDOT = AMDOT*(3.4763 + 1.066D-4*((TVENT+459.67)/1.8))*0.068559 END IF C ESTIMATE CRUDE FLOW VELOCITY IF(IVFLAG.EQ.0) V = VDOT/(AL4*AL) IF(IVFLAG.EQ.1) V = 2.0D0*VDOT/(AL4*AL) ACH = VDOT/A(4)/AL C IF(LFLAG(1).eq. 0) GO TO 30 ! Logic for Conventional Attic C IF(LFLAG(2).eq. 1) then CALL MODRA(TASV(1,1),TOS(2,1),TIS(8,1),TO,P1,AL,GAP,TMDOT, & TCMDOT,TV,RA,RE) AMCDOT2 = TCMDOT AMDOT2 = TMDOT VDOT2 = TV Ray2 = RA Rey2 = RE V2 = VDOT2/(AL*GAP) ACH2 = VDOT2/(A(8)*GAP) ELSE IF(LFLAG(2) .eq. 0) then ! Option for future algorithm (not presently used) CALL VENT(TO,TASV(1,1),WS,DIR,H,AI,AO,1,AMCDOT2,AMDOT2,VDOT2) V2 = VDOT2/(AL*GAP) ACH2 = VDOT2/(A(8)*GAP) ELSE GO TO 998 END IF C IF(IVFLAG.EQ.1) GO TO 30 ! FLAG FOR SHED TYPE ROOF C IF(LFLAG(2).eq. 1) then CALL MODRA(TASV(2,1),TOS(3,1),TIS(9,1),TO,P1,AL,GAP,TMDOT, & TCMDOT,TV,RA,RE) AMCDOT3 = TCMDOT AMDOT3 = TMDOT VDOT3 = TV Ray3 = RA Rey3 = RE V3 = VDOT3/(AL*GAP) ACH3 = VDOT3/(A(9)*GAP) ELSE IF(LFLAG(2) .eq. 0) then ! Option for future algorithm (not presently used) CALL VENT(TO,TASV(2,1),WS,DIR,H,AI,AO,1,AMCDOT3,AMDOT3,VDOT3) V3 = VDOT3/(AL*GAP) ACH3 = VDOT3/(A(9)*GAP) ELSE GO TO 998 END IF 30 CONTINUE C*********************************************************************** C CALCULATE CONVECTION COEFFICIENTS AT INSIDE SURFACES CALL HCON(TIS(1,1),TA(1), 0.0D0,AL1,1,1,V,HCI(1)) CALL HCON(TIS(2,1),TA(1),PITCH1,AL2,2,1,V,HCI(2)) CALL HCON(TIS(3,1),TA(1),PITCH2,AL3,2,1,V,HCI(3)) CALL HCON(TIS(4,1),TA(1), 90.0D0,AL4,1,1,V,HCI(4)) CALL HCON(TIS(5,1),TA(1), 90.0D0,AL5,1,1,V,HCI(5)) CALL HCON(TIS(6,1),TA(1), 90.0D0,AL6,1,1,V,HCI(6)) CALL HCON(TIS(7,1),TA(1), 90.0D0,AL7,1,1,V,HCI(7)) C if (LFLAG(1) .eq. 1 .and. LFLAG(3) .eq. 1) then CALL HGapOpn(TIS(8,1),TASV(1,1),TOS(2,1),TO,EI(8),PITCH1,AL,AL8, & GAP,2,0.0D0,HCI(8),AMDOT2,NUS(8)) if (IVFLAG .eq. 0) &CALL HGapOpn(TIS(9,1),TASV(2,1),TOS(3,1),TO,EI(9),PITCH2,AL,AL9, & GAP,2,0.0D0,HCI(9),AMDOT3,NUS(9)) else if (LFLAG(1) .eq. 1 .and. LFLAG(3) .eq. 0) then CALL HGapCld(TIS(8,1),TASV(1,1),TOS(2,1),TO,PITCH1,AL,AL8,GAP,2, & 0.0D0,HCI(8),AMDOT2,NUS(8)) if (IVFLAG .eq. 0) &CALL HGapCld(TIS(9,1),TASV(2,1),TOS(3,1),TO,PITCH2,AL,AL9,GAP,2, & 0.0D0,HCI(9),AMDOT3,NUS(9)) end if C*********************************************************************** C CALCULATE CONVECTION COEFFICIENTS AT OUTSIDE SURFACES WS1 = WS*5280. ! WIND VELOCITY FT PER HR (POINT FOR TMY2 ADJUSTMENT) C CALL HCON(TOS(1,1),TI,0.0D0,AL1,2,1,0.0D0,HCO(1)) IF (LFLAG(1) .eq. 0) then CALL HCON(TOS(2,1),TO,PITCH1,AL2,1,ISURF,WS1,HCO(2)) if (IVFLAG .eq. 0) CALL HCON(TOS(3,1),TO,PITCH2,AL3,1,ISURF, & WS1,HCO(3)) go to 150 END IF IF (LFLAG(1) .eq. 1 .and. LFLAG(3) .eq. 1) then CALL HGapOpn(TIS(8,1),TASV(1,1),TOS(2,1),TO,EO(2),PITCH1,AL,AL2, & GAP,1,0.0D0,HCO(2),AMDOT2,NUS(2)) CALL HCON(TOS(8,1),TO,PITCH1,AL2,1,ISURF,WS1,HCO(8)) if (IVFLAG .eq. 0) &CALL HGapOpn(TIS(9,1),TASV(2,1),TOS(3,1),TO,EO(3),PITCH2,AL,AL3, & GAP,1,0.0D0,HCO(3),AMDOT3,NUS(3)) if (IVFLAG .eq. 0) CALL HCON(TOS(9,1),TO,PITCH2,AL3,1,ISURF, & WS1,HCO(9)) go to 150 END IF IF (LFLAG(1) .eq. 1 .and. LFLAG(3) .eq. 0) then CALL HGapCld(TIS(8,1),TASV(1,1),TOS(2,1),TO,PITCH1,AL,AL2,GAP,1, & 0.0D0,HCO(2),AMDOT2,NUS(2)) CALL HCON(TOS(8,1),TO,PITCH1,AL2,1,ISURF,WS1,HCO(8)) if (IVFLAG .eq. 0) &CALL HGapCld(TIS(9,1),TASV(2,1),TOS(3,1),TO,PITCH2,AL,AL3,GAP,1, & 0.0D0,HCO(3),AMDOT3,NUS(3)) if (IVFLAG .eq. 0) CALL HCON(TOS(9,1),TO,PITCH2,AL3,1,ISURF, & WS1,HCO(9)) END IF 150 CALL HCON(TOS(4,1),TO,90.0D0,AL4,1,1,WS1,HCO(4)) CALL HCON(TOS(5,1),TO,90.0D0,AL5,1,1,WS1,HCO(5)) CALL HCON(TOS(6,1),TO,90.0D0,AL6,1,1,WS1,HCO(6)) CALL HCON(TOS(7,1),TO,90.0D0,AL7,1,1,WS1,HCO(7)) C CALCULATE RADIATION COEFFICIENTS AT INSIDE SURFACES DO 200 I = 1,7 DO 200 J = 1,7 HR(I,J) = HRAD(G(I,J),TIS(I,1),TIS(J,1)) 200 CONTINUE C CALCULATE RADIATION COEFFICIENTS AT OUTSIDE SURFACES HRO(1) = HRAD(EO(1),TOS(1,1),TI) IF (LFLAG(1) .eq. 0) then HRO(2) = HRAD(EO(2),TOS(2,1),TS2) HRO(3) = HRAD(EO(3),TOS(3,1),TS3) ELSE E_2 = (1.0D0-EI(8))/EI(8) + (1.0D0-EO(2))/EO(2) + 1.0D0/FP(AL,AL2,GAP) E_3 = (1.0D0-EI(9))/EI(9) + (1.0D0-EO(3))/EO(3) + 1.0D0/FP(AL,AL3,GAP) HRO(2) = HRAD(1.0D0,TOS(2,1),TIS(8,1))/E_2 HRO(3) = HRAD(1.0D0,TOS(3,1),TIS(9,1))/E_3 HRO(8) = HRAD(EO(8),TOS(8,1),TS2) HRO(9) = HRAD(EO(9),TOS(9,1),TS3) END IF DO 250 I = 4,7 HRO(I) = HRAD(EO(I),TOS(I,1),TSV) 250 CONTINUE C C CALCULATE MOISTURE BALANCES, LATENT HEATS IF(QLAT.LE.1.) THEN DO 259 I = 1,7 AMW(I) = 0.0D0 BMW(I) = 0.0D0 259 CONTINUE IF(ITRUS.NE.0) AMWT = 0.0 ELSE DO 260 I = 1,7 260 TS(I) = TIS(I,1) CALL MOIST(TS,TI,TA(1),HCI,A,AWRAT,AMC,HUMI,PERM, & AMDOT,EXFIL,AMW,BMW,HM,7) DO 270 I = 1,7 IF(AMC(I).LE.1.D-6.AND.AMW(I).LT.0.0) THEN AMW(I) = 0.0D0 BMW(I) = 0.0D0 ELSE IF(AMW(I).LT.0.0.AND.ABS(AMW(I)).GT.(AMC(I)*AMASS(I))) THEN FACT = ABS(AMC(I)*AMASS(I)/AMW(I)) AMW(I) = AMW(I)*FACT BMW(I) = BMW(I)*FACT ELSE CONTINUE END IF 270 CONTINUE IF(ITRUS.NE.0) THEN IF(AMCT.LE.1.D-6.AND.AMWT.LT.0.0D0) THEN AMWT = 0.0D0 BMWT = 0.0D0 ELSE IF(AMWT.LT.0.0D0 .AND. DABS(AMWT).GT.(AMCT*AMASST*ATRUS)) 7 THEN FACT = DABS(AMCT*AMASST*ATRUS/AMWT) AMWT = AMWT*FACT BMWT = BMWT*FACT END IF END IF END IF C next section added for model with attic trusses IF(ITRUS.NE.0) THEN CALL TRUSMOD(G,TIS,TA,V,IDUCT,TTRNEW) END IF C next line added for duct model IF(IDUCT.NE.0) THEN IF(IDMOD.EQ.0) CALL DUCTSS(G,TIS,TA,T) END IF C C*********************************************************************** C SET UP MATRIX OF COEFFICIENTS IN HEAT BALANCE EQUATIONS DO 300 I = 1,7 DO 300 J = 1,7 AA(I,J) = -HR(I,J) 300 CONTINUE IF (LFLAG(1) .eq. 1) then E8_2 = (1.0D0-EI(8))/EI(8) + (1.0D0-EO(2))/EO(2) & + 1.0D0/FP(AL,AL2,GAP) E9_3 = (1.0D0-EI(9))/EI(9) + (1.0D0-EO(3))/EO(3) & + 1.0D0/FP(AL,AL3,GAP) HRasv2 = HRAD(1.0D0,TIS(8,1),TOS(2,1))/E8_2 HRasv3 = HRAD(1.0D0,TIS(9,1),TOS(3,1))/E9_3 END IF C IF (LFLAG(1) .EQ. 1 .AND. LFLAG(3) .EQ. 1) THEN AA(16,9) = -HRasv2 ! Used in Interior Heat Balance for Surface 8 AA(9,16) = -HRasv2 ! Used in Exterior Heat Balance for Surface 2 IF(IVFLAG.EQ.0) THEN AA(17,10) = -HRasv3 ! Used in Interior Heat Balance for Surface 9 AA(10,17) = -HRasv3 ! Used in Exterior Heat Balance for Surface 3 END IF ELSE IF (LFLAG(1) .EQ. 1 .AND. LFLAG(3) .EQ. 0) THEN AA(16,9) = -(HRasv2 + HCI(8)) ! Used in Interior Heat Balance for Surface 8 AA(9,16) = -(HRasv2 + HCO(2)) ! Used in Exterior Heat Balance for Surface 2 IF(IVFLAG.EQ.0) THEN AA(17,10) = -(HRasv3 + HCI(9)) ! Used in Interior Heat Balance for Surface 9 AA(10,17) = -(HRasv3 + HCO(3)) ! Used in Exterior Heat Balance for Surface 3 END IF END IF C DO 310 I = 1,7 DO 305 J = 1,7 305 AA(I,I) = AA(I,I) + HR(I,J) 310 AA(I,I) = AA(I,I) + Z(I,1) + HCI(I) c C ADD IN LATENT HEAT TERMS DO 315 I = 1,7 315 AA(I,I) = AA(I,I) + BMW(I)*QLAT C*********************************************************************** DO 320 I = 1,7 AA(I,I+7) = - Y(I,1) AA(I,15) = -HCI(I) AA(I+7,I) = -Y(I,1) AA(I+7,I+7) = X(I,1) + HCO(I) + HRO(I) 320 AA(15,I) = A(I)*HCI(I) C IF (LFLAG(1) .eq. 1) then AA(16,16) = Z(8,1) + HRasv2 + HCI(8) ! Used in Interior Heat Balance for Surface 8 AA(16,18) = -Y(8,1) ! Used in Interior Heat Balance for Surface 8 AA(18,16) = -Y(8,1) ! Used in Exterior Heat Balance for Surface 8 AA(18,18) = X(8,1) + HRO(8) + HCO(8) ! Used in Exterior Heat Balance for Surface 8 IF(IVFLAG.EQ.0) THEN AA(17,17) = Z(9,1) + HRasv3 + HCI(9) ! Used in Interior Heat Balance for Surface 9 AA(17,19) = -Y(9,1) ! Used in Interior Heat Balance for Surface 9 AA(19,17) = -Y(9,1) ! Used in Exterior Heat Balance for Surface 9 AA(19,19) = X(9,1) + HCO(9) + HRO(9) ! Used in Exterior Heat Balance for Surface 9 END IF END IF C IF (LFLAG(1) .eq. 1 .and. LFLAG(3) .eq. 1) then AA(9,20) = -HCO(2) ! Used in Exterior Heat Balance for Surface 2 AA(16,20) = -HCI(8) ! Used in Interior Heat Balance for Surface 8 IF(IVFLAG.EQ.0) then AA(10,21) = -HCO(3) ! Used in Exterior Heat Balance for Surface 3 AA(17,21) = -HCI(9) ! Used in Interior Heat Balance for Surface 9 END IF ELSE IF (LFLAG(1) .eq. 1 .and. LFLAG(3) .eq. 0) then AA(9,20) = 0.0D0 ! Used in Exterior Heat Balance for Surface 2 AA(16,20) = 0.0D0 ! Used in Interior Heat Balance for Surface 8 IF(IVFLAG.EQ.0) then AA(10,21) = 0.0D0 ! Used in Exterior Heat Balance for Surface 3 AA(17,21) = 0.0D0 ! Used in Interior Heat Balance for Surface 9 END IF END IF C C*********************************************************************** C1 = 0.0D0 DO 330 I = 1,7 330 C1 = C1 + A(I)*HCI(I) IF (LFLAG(1) .eq. 1 .and. Lflag(3) .eq. 1) then AA(20,9) = A(2)*HCO(2) ! Used in ASV Heat Balance between Surfaces 2 & 8 AA(20,16) = A(8)*HCI(8) ! Used in ASV Heat Balance between Surfaces 2 & 8 IF(IVFLAG.EQ.0) THEN AA(21,10) = A(3)*HCO(3) ! Used in ASV Heat Balance between Surfaces 3 & 9 AA(21,17) = A(9)*HCI(9) ! Used in ASV Heat Balance between Surfaces 3 & 9 END IF ELSE IF (LFLAG(1) .eq. 1 .and. Lflag(3) .eq. 0) then AA(20,9) = -0.5 ! Used in ASV Heat Balance between Surfaces 2 & 8 closed airspace AA(20,16) = -0.5 ! Used in ASV Heat Balance between Surfaces 2 & 8 IF(IVFLAG.EQ.0) THEN AA(21,10) = -0.5 ! Used in ASV Heat Balance between Surfaces 3 & 9 AA(21,17) = -0.5 ! Used in ASV Heat Balance between Surfaces 3 & 9 closed airspace END IF END IF C C Next line added by K. Wilkes on 2-19-97 if(c1.le.1.D-10) c1 = 1.D-6 C2I = (AMCDOT + EXFIL*0.24)/C1 IF(C2I.GT.0.02D0) C3 = DEXP(-1./C2I) - 1.0 IF(C2I.LE.0.02D0) C3 = -1.0 AA(15,15) = -C1/(1. + C2I*C3) C IF (LFLAG(1) .eq. 1) then CV21 = 0.0 CV21 = A(8)*HCI(8) + A(2)*HCO(2) if(CV21.le.1.D-10) CV21 = 1.D-6 CV22I = AMCDOT2/CV21 IF(CV22I .GT. 0.02D0) CV23 = DEXP(-1./CV22I) - 1.0 IF(CV22I .LE. 0.02D0) CV23 = -1.0D0 AA(20,20) = -CV21/(1. + CV22I*CV23) ! Used in ASV Heat Balance between Surfaces 2 & 8 IF (LFLAG(3) .EQ. 0) AA(20,20) = 1.0D0 ! Used in ASV Heat Balance between Surfaces 2 & 8, closed airspace END IF C IF (LFLAG(1) .EQ. 1 .AND. IVFLAG.EQ.0) THEN CV31 = 0.0 CV31 = A(9)*HCI(9) + A(3)*HCO(3) if(CV31.le.1.D-10) CV31 = 1.D-6 CV32I = AMCDOT3/CV31 IF(CV32I .GT. 0.02D0) CV33 = DEXP(-1./CV32I) - 1.0 IF(CV32I .LE. 0.02D0) CV33 = -1.0D0 AA(21,21) = -CV31/(1. + CV32I*CV33) ! Used in ASV Heat Balance between Surfaces 3 & 9 IF (LFLAG(3) .EQ. 0) AA(21,21) = 1.0D0 ! Used in ASV Heat Balance between Surfaces 3 & 9, closed airspace END IF C*********************************************************************** C SET UP VECTOR FOR RIGHT HAND SIDE OF HEAT BALANCE EQUATIONS DO 400 I = 1,7 BB(I) = TR*(Z(I,1)-Y(I,1)) - CR(I)*QI(I,2) BB(I) = BB(I) -BETA(I)/2.*(Z(I,1)*(TIS(I,1)-TR)**2 & -Y(I,1)*(TOS(I,1)-TR)**2) DO 400 J = 2,N(I) 400 BB(I)=BB(I)-Z(I,J)*(TIS(I,J)-TR)+Y(I,J)*(TOS(I,J)-TR) & -BETA(I)/2.*(Z(I,J)*(TIS(I,J)-TR)**2 & -Y(I,J)*(TOS(I,J)-TR)**2) C ADD IN LATENT HEAT EFFECTS DO 405 I = 1,7 405 BB(I) = BB(I) + AMW(I)*QLAT + BMW(I)*QLAT*TIS(I,1) DO 410 I = 1,7 BB(I+7) = TR*(X(I,1)-Y(I,1)) + CR(I)*QO(I,2) BB(I+7) = BB(I+7) -BETA(I)/2.*(X(I,1)*(TOS(I,1)-TR)**2 & -Y(I,1)*(TIS(I,1)-TR)**2) DO 410 J = 2,N(I) 410 BB(I+7) = BB(I+7)-X(I,J)*(TOS(I,J)-TR)+Y(I,J)*(TIS(I,J)-TR) & -BETA(I)/2.*(X(I,J)*(TOS(I,J)-TR)**2 & -Y(I,J)*(TIS(I,J)-TR)**2) C BB(8) = BB(8) + HCO(1)*TI + HRO(1)*TI ! Source for inside house C IF (LFLAG(1) .eq. 0) then BB(9) = BB(9) + HCO(2)*TO + HRO(2)*TS2 + ALF(2)*QSOL(2) ! East Roof NO ASV BB(10) = BB(10) + HCO(3)*TO + HRO(3)*TS3 + ALF(3)*QSOL(3) ! West Roof NO ASV END IF C DO 420 I = 4,7 420 BB(I+7) = BB(I+7) + HCO(I)*TO + HRO(I)*TSV + ALF(I)*QSOL(I) BB(15) = (AMCDOT*TVENT + EXFIL*0.24*TI)*C3/(1.0 + C2I*C3) C C SET UP VECTOR FOR RIGHT HAND SIDE OF ASV HEAT BALANCE EQUATIONS (Surfaces 8 and 9) IF (LFLAG(1) .eq. 1) then DO 430 I = 8,9 ! Inside Surface 8 and 9 BB(I+8) = TR*(Z(I,1)-Y(I,1)) - CR(I)*QI(I,2) BB(I+8) = BB(I+8) -BETA(I)/2.*(Z(I,1)*(TIS(I,1)-TR)**2 & -Y(I,1)*(TOS(I,1)-TR)**2) DO 430 J = 2,N(I) 430 BB(I+8)=BB(I+8)-Z(I,J)*(TIS(I,J)-TR)+Y(I,J)*(TOS(I,J)-TR) & -BETA(I)/2.*(Z(I,J)*(TIS(I,J)-TR)**2 & -Y(I,J)*(TOS(I,J)-TR)**2) DO 440 I = 8,9 ! Outside Surface 8 and 9 BB(I+10) = TR*(X(I,1)-Y(I,1)) + CR(I)*QO(I,2) BB(I+10) = BB(I+10) -BETA(I)/2.*(X(I,1)*(TOS(I,1)-TR)**2 & -Y(I,1)*(TIS(I,1)-TR)**2) DO 440 J = 2,N(I) 440 BB(I+10) = BB(I+10)-X(I,J)*(TOS(I,J)-TR)+Y(I,J)*(TIS(I,J)-TR) & -BETA(I)/2.*(X(I,J)*(TOS(I,J)-TR)**2 & -Y(I,J)*(TIS(I,J)-TR)**2) C BB(18) = BB(18) + HCO(8)*TO + HRO(8)*TS2 + ALF(8)*QSOL(8) IF(IVFLAG.EQ.0) & BB(19) = BB(19) + HCO(9)*TO + HRO(9)*TS3 + ALF(9)*QSOL(9) C BB(20) = (AMCDOT2*TO)*CV23/(1.0 + CV22I*CV23) IF(LFLAG(3) .EQ. 0) BB(20) = 0.0D0 IF(IVFLAG.EQ.0) & BB(21) = (AMCDOT3*TO)*CV33/(1.0 + CV32I*CV33) IF (LFLAG(3) .EQ. 0) BB(21) = 0.0D0 END IF C C NEXT SECTION MODIFIES THE BB MATRIX TO ACCOUNT FOR DUCTS IF(IDUCT.NE.0) THEN DO 445 I = 1,7 BB(I) = BB(I) + QRAD(I)/A(I) 445 CONTINUE C Next line modified 2-14-96 to correct sign convention on QCON and QLEAKS BB(15) = BB(15) - QCON - QLEAKS END IF C Next section modifies the BB matrix to account for trusses IF(ITRUS.NE.0) THEN DO 450 I = 1,7 BB(I) = BB(I) + QRADT(I)/A(I) 450 CONTINUE C Next line modified 2-14-96 to correct sign convention on QCONT BB(15) = BB(15) - QCONT END IF C*********************************************************************** SKBIG = 1.E20 DO 455 I = 1,7 IF(KFLAG(I).EQ.0) GO TO 460 AA(I,I) = AA(I,I)*SKBIG BB(I) = TIS(I,1)*AA(I,I) 460 IF(KFLAG(I+7).EQ.0) GO TO 455 AA(I+7,I+7) = AA(I+7,I+7)*SKBIG BB(I+7) = TOS(I,1)*AA(I+7,I+7) 455 CONTINUE C C*********************************************************************** C SOLVE SYSTEM OF EQUATIONS c If(ihr .eq. 2) then c do 1232 i = 1,21 c 1232 write(15,1233) IDOY,IHR,(AA(i,j),j=1,21),BB(i) c1233 format(2(i7,','),22(1pE12.5,',')) c end if CALL SOLVP(Narray,AA,BB,XXX,Nequ) c write(6,1244) Narray,Nequ,(i,XXX(i),i=1,Narray) c1244 format(' Narray : ',i3,' Nequ : ',i3,/,1(i3,1pE12.5)) c pause C*********************************************************************** C TEST FOR CONVERGENCE OF SOLUTION EPS = 1.D-3 IFLAG = 0 DO 500 I = 1,7 500 IF(ABS(XXX(I)-TIS(I,1)).GT.EPS) GO TO 600 DO 505 I = 1,7 505 IF(ABS(XXX(I+7)-TOS(I,1)).GT.EPS) GO TO 600 IF(ABS(XXX(15)-TA(1)).GT.EPS) GO TO 600 if(LFLAG(1) .eq. 0) go to 520 DO 510 I = 8,9 510 IF(ABS(XXX(I+8)-TIS(I,1)).GT.EPS) GO TO 600 DO 515 I = 8,9 515 IF(ABS(XXX(I+10)-TOS(I,1)).GT.EPS) GO TO 600 IF(ABS(XXX(20)-TASV(1,1)).GT.EPS) GO TO 600 IF(ABS(XXX(21)-TASV(2,1)).GT.EPS) GO TO 600 520 CONTINUE C Next line added for model with trusses IF(ITRUS.NE.0) THEN IF(ABS(TTRUS(1) - TTRNEW).GT.EPS) GO TO 600 END IF C*********************************************************************** GO TO 700 600 IFLAG = 1 700 CONTINUE C*********************************************************************** DO 800 I = 1,7 TIS(I,1) = TIS(I,1) + Relax*(XXX(I) - TIS(I,1)) 800 TOS(I,1) = TOS(I,1) + Relax*(XXX(I+7) - TOS(I,1)) TA(1) = TA(1) + Relax*(XXX(15) - TA(1)) IF(LFLAG(1) .EQ.0) GO TO 810 DO 805 I = 8,9 TIS(I,1) = TIS(I,1) + Relax*(XXX(I+8) - TIS(I,1)) 805 TOS(I,1) = TOS(I,1) + Relax*(XXX(I+10) - TOS(I,1)) TASV(1,1) = TASV(1,1) + Relax*(XXX(20) - TASV(1,1)) IF(IVFLAG.EQ.0) &TASV(2,1) = TASV(2,1) + Relax*(XXX(21) - TASV(2,1)) 810 CONTINUE C Next line added for model with trusses IF(ITRUS.NE.0) TTRUS(1) = TTRNEW C*********************************************************************** S1 = 0.0 DO 850 I = 1,7 850 S1 = S1 + A(I)*HCI(I)*TIS(I,1) TE = AMCDOT*TVENT + EXFIL*0.24 + (S1 - C1*TA(1)) TE = TE/(AMCDOT + EXFIL) NIT = NIT + 1 IF(IFLAG.EQ.1.AND.NIT.LE.100) GO TO 20 C CALCULATE HEAT FLUXES DO 900 I = 1,Nsurf QI(I,1) = CR(I)*QI(I,2) 900 QO(I,1) = CR(I)*QO(I,2) DO 950 I = 1,Nsurf DO 950 J = 1,N(I) c if(idoy .eq. 253 .and. ihr .eq. 14 .and. i .eq. 8) c &write(15,949) ihr,i,j,QI(i,1),QO(i,1),TR,Beta(i),X(i,j),Y(i,j), c & Z(i,j),TIS(i,j),TOS(i,j) c 949 format(3(i3,","),9(1PE12.3,",")) QI(I,1) = QI(I,1) + Z(I,J)*(TIS(I,J)-TR)-Y(I,J)*(TOS(I,J)-TR) & +BETA(I)/2.*(Z(I,J)*(TIS(I,J)-TR)**2-Y(I,J)*(TOS(I,J)-TR)**2) 950 QO(I,1) = QO(I,1) + Y(I,J)*(TIS(I,J)-TR)-X(I,J)*(TOS(I,J)-TR) & +BETA(I)/2.*(Y(I,J)*(TIS(I,J)-TR)**2-X(I,J)*(TOS(I,J)-TR)**2) C CALCULATE TOTAL CEILING HEAT FLOW QCEIL = QO(1,1)*A(1) C CALCULATE Above-Sheating VENTILATION HEAT FLOW PER HR [Btu/hr] IF (LFLAG(1) .EQ. 0) GO TO 955 Qasv2 = AMCDOT2*(0.5*(TIS(8,1) + TOS(2,1)) - TO) Qasv_ch2 = -A(2)*(QI(8,1) - QO(2,1)) IF (IVFLAG .EQ. 1) &Qasv3 = (0.5*(TIS(9,1) + TOS(3,1)) - TO) Qasv_ch3 = -A(3)*(QI(9,1) - QO(3,1)) C C Determination of Conduction, Convection and Radiation Effects Across Air Space TK = (TASV(1,1)+459.67)/1.8 C THERMAL COND. IN BTU/(HR-FT-F) AirK = 0.6325D-5*DSQRT(TK)/(1.+(245.4*10.**(-12./TK))/TK)*241.77 Rcond = Airk/Gap Rrad = HRasv2 Rconv = 1.0/(1.0/HCO(2) + 1.0/HCI(8)) Rtotal = 1.0/(Rcond + Rrad + Rconv) Qrdgap = Rrad*(TIS(8,1)-TOS(2,1)) Qcvtop = HCI(8)*(TIS(8,1)-Tasv(1,1)) Qcvosb = HCO(2)*(Tasv(1,1)-TOS(2,1)) Qvtgap = Qasv2/A(2) c write (15,1189) IDOY,IHR,TO,Tsky,TA(1),ACH,WS,Qsolr,HCI(8),NUS(8), c & HCO(2),NUS(2),TIS(8,1),TOS(2,1),(TIS(8,1)-TOS(2,1)),Tasv(1,1), c & -QI(8,1),Qrdgap,Qcvtop,-QO(2,1),Qcvosb,Qvtgap,QO(1,1) c1189 FORMAT(1X,2(I4,','),3(F10.1,','),2(f10.2,','),1(f10.1,','), c & 4(f10.2,','),4(f10.1,','),7(f10.2,',')) 955 continue C CALCULATE NEW MOISTURE CONTENTS DO 960 I = 1,7 AMC(I) = AMC(I) + AMW(I)/AMASS(I) C NEXT LINE ADDED 10-1-89 IF(AMC(I).GT.0.30) AMC(I) = 0.3 960 IF(AMC(I).LE.0.00) AMC(I) = 0.0 IF(ITRUS.NE.0) THEN AMCT = AMCT + AMWT/(AMASST*ATRUS) IF(AMCT.GT.0.30) AMCT = 0.3 IF(AMCT.LE.0.0) AMCT = 0.0 END IF C Next lines added for duct model IF(IDUCT.NE.0) THEN IF(IDMOD.NE.0.) & CALL DUCTTR(G,TIS,TO,TA,T,TAD,TWWD,TWOD) END IF C C Account for Heat Transfer across Leaky Attic Floor C QDUCTS Computed as +Heat (Tin - Tex) in Heating Mode; Attic Heat Loss C QDUCTS Computed as -Heat (Tin - Tex) in Cooling Mode; Attic Heat Gain c Qfloor = 0.24*Exfil*(TI - TA(1)) c QDUCTx = QDUCTS + Qfloor C C WRITE OUT RESULTS IN COMMA DELIMITED FILE (Include duct if IDUCT=1) IF(KHEADER.EQ.0) THEN WRITE(15,2007) (Blank,I,I,I,I,I, I = 1,7) 2007 FORMAT(1X,'Day,','HOUR,',7(a2,'HCO(',i2,'),','HCI(',i2,'),', & 'HRO(',i2,'),','HR(',i2,' 1),','HM(',i2,'),')) IF (IDUCT.EQ.0) WRITE(11,1005) IF (IDUCT.EQ.1) WRITE(11,1007) (Blank,I,I=1,NSUPLY+NRETRN+2) KHEADER=1 END IF IF (IDUCT.EQ.0) THEN WRITE(11,1006) IDOY,IHR,TI,RHi,TO,TodDwPt,Patm,Tsky, & TA(1), ACH, WS,DIR,QIR,QSOLH,DIRSOL,DIFSOL, & QSOL(2), QSOL(3), QSOL(8), QSOL(9), & QO(1,1), QI(1,1),TOS(1,1),TIS(1,1), & -1*QO(2,1),-1*QI(2,1),TOS(2,1),TIS(2,1),TsRfDwPt, & -1*QO(3,1),-1*QI(3,1),TOS(3,1),TIS(3,1),TnRfDwPt, & (-1*QO(I,1),-1*QI(I,1),TOS(I,1),TIS(I,1),I=4,9) ELSE WRITE(11,1008) IDOY,IHR,TO,TDwPt,Patm,QSOLH,DIRSOL,DIFSOL, & Tsky, TA(1), DIR,WS,HCO(1),HCI(1), & QSOL(2), QSOL(3), QSOL(8), QSOL(9), & QO(1,1), QI(1,1),TOS(1,1),TIS(1,1), & (-1*QO(I,1),-1*QI(I,1),TOS(I,1),TIS(I,1),I=2,9), & (T(I),I=1,12),QRADtot,Qcon,Qleaks, & Qstor,QDUCTx,Qfloor,TI END IF 1005 FORMAT(1X,'Day,','HOUR,', & 'Tair indr,',' RH indr,','Tair outd,','TairDwPt,', & ' Patm,',' Tsky,',' T attic,',' ACH,', & 'Wind(mph),','Wind dir,',' QIR,',' QSOLH,', & ' DIRSOL,',' DIFSOL,', & ' Qsol(2),',' Qsol(3),',' Qsol(8),',' Qsol(9),', & 'Qext ceil,','Qint ceil,','Text ceil,','Tint ceil,', & 'Qext sRf,','Qint sRf,','Text sRf,','Tint sRf,', & 'TDwPt sRf,', & 'Qext nRf,','Qint nRf,','Text nRf,','Tint nRf,', & 'TDwPt nRf,', & 'Qext eGab,','Qint eGab,','Text eGab,','Tint eGab,', & 'Qext wGab,','Qint wGab,','Text wGab,','Tint wGab,', & 'Qext sEv,','Qint sEv,','Text sEv,','Tint sEv,', & 'Qext nEv,','Qint nEv,','Text nEv,','Tint nEv,', & 'Qext sASV,','Qint sASV,','Text sASV,','Tint sASV,', & 'Qext nASV,','Qint nASV,','Text nASV,','Tint nASV,') C C 1006 FORMAT(1X,2(I4,','),4(F10.1,','),F10.2,',',2(F10.1,','), & 2(F10.3,','),9(F10.1,','), & 1(F10.3,',',F10.3,',',F9.2,',',F9.2,','), & 2(F10.3,',',F10.3,',',F9.2,',',F9.2,',',F9.2,','), & 6(F10.3,',',F10.3,',',F9.2,',',F9.2,',')) C 2006 FORMAT(1X,7(A2,F10.3,',')) 1007 FORMAT(1X,'Day,','HOUR,', & 'Tair outd,','TairDwPt,',' Patm,', & ' QSOLH,',' DIRSOL,',' DIFSOL,', & ' Tsky,','T attic,',' Wind dir,', & 'Wind(mph),',' HCO,',' HCI,', & ' Qsol(2),',' Qsol(3),',' Qsol(8),',' Qsol(9),', & 'Qext ceil,','Qint ceil,','Text ceil,','Tint ceil,', & 'Qext eOSB,','Qint eOSB,','Text eOSB,','Tint eOSB,', & 'Qext wOSB,','Qint wOSB,','Text wOSB,','Tint wOSB,', & 'Qext nOSB,','Qint nOSB,','Text nOSB,','Tint nOSB,', & 'Qext sOSB,','Qint sOSB,','Text sOSB,','Tint sOSB,', & 'Qext eEav,','Qint eEAV,','Text eEAV,','Tint eEAV,', & 'Qext wEav,','Qint wEAV,','Text wEAV,','Tint wEAV,', & 12(a2,'Tduct(',i2,'),'), & 'QradDucts,','QconDucts,','QleakDucts,', & 'QstorDucts,','QDUCTS,','QFloor,','Tair indr') 1008 FORMAT(1X,2(I4,','),8(F10.1,','),F10.1,',',7(F10.1,','), & 7(F10.3,',',F10.3,',',F9.2,',',F9.2,','), & 12(F9.2,','),6(f10.2,','),f10.1) WRITE(15,2008) IDOY,IHR,(blank,HCO(I),HCI(I),HRO(I), & HR(I,1),HM(I),I=1,7) 2008 FORMAT(1X,2(I4,','),7(a2,F10.3,',',F10.3,',',F10.3,',', & F10.3,',',F10.3,',')) c WRITE(11,1008) IDOY,IHR,TO,TDwPt,Patm,QSOLH,DIRSOL,DIFSOL, c & Tsky, TA(1), DIR,WS,Vdot,ACH, c & QSOL(2), QSOL(3), QSOL(8), QSOL(9), c & QO(1,1), QI(1,1),TOS(1,1),TIS(1,1), c & (-1*QO(I,1),-1*QI(I,1),TOS(I,1),TIS(I,1),I=2,3), c & (-1*QO(I,1),-1*QI(I,1),TOS(I,1),TIS(I,1),I=8,9), c & Qasv_ch2,Qasv2,AMDOT2,V2,ACH2,Ray2,Rey2, c & Qasv_ch3,Qasv3,AMDOT3,V3,ACH3,Ray3,Rey3,ONTIME, c & (T(I),I=1,12),QRADtot,Qcon,Qleaks, c & Qstor,QDUCTS,Qfloor,TI C IF(IHR .EQ. 1) WRITE(*,*) IDOY,IHR C UPDATE TEMPERATURES AND HEAT FLUXES FOR NEXT HOUR DO 1100 I = 1,Nsurf QI(I,2) = QI(I,1) QO(I,2) = QO(I,1) DO 1100 J = 1,N(I)-1 TIS(I,N(I)-J+1) = TIS(I,N(I)-J) 1100 TOS(I,N(I)-J+1) = TOS(I,N(I)-J) TA(2) = TA(1) TASV(1,2) = TASV(1,1) IF(IVFLAG.EQ.0) TASV(2,2) = TASV(2,1) C Next line added for model with trusses IF(ITRUS.NE.0) TTRUS(2) = TTRUS(1) GO TO 10 998 Write(*,1000) (blank,i,LFLAG(i),i=1,2) 999 STOP 1000 format(' Problem Occurred in Specifing ASV Options',/, & a2,'LFLAG(',i,') = ',i2,/, & a2,'LFLAG(',i,') = ',i2) END C*********************************************************************** C** ** C** FUNCTION FMN ** C** ** C*********************************************************************** FUNCTION FMN(A,B,C,PHI) C THIS FUNCTION CALCULATES THE VIEW FACTOR BETWEEN TWO FINITE C RECTANGULAR PLATES SHARING AN EDGE C SEE APPENDIX A OF SPARROW AND CESS, CONFIGURATION 2 C SEE ALSO A. FEINGOLD, PROC. ROY. SOC., A, VOL. 292, PP. 51-60 (1965) C VIEW FACTOR FROM PLATE 1 TO PLATE 2 C A = WIDTH OF PLATE 2, C = WIDTH OF PLATE 1 C B = LENGTH OF BOTH PLATES C PHI = ANGLE BETWEEN PLATES (RADIANS) C NUMERICAL INTEGRATION PERFORMED BY SIMPSON'S RULE IMPLICIT REAL*8 (A-H, O-Z) DIMENSION Q(65) X = A/B Y = C/B SP = DSIN(PHI) S2P = DSIN(2.*PHI) CP = DCOS(PHI) C2P = DCOS(2.*PHI) Z = X*X + Y*Y - 2.*X*Y*CP DXI = Y/64. DO 1 N = 1,65 ANM1 = N - 1 XI = DXI*ANM1 E = DSQRT(1. + XI*XI*SP*SP) F = DATAN((X-XI*CP)/E) + DATAN(XI*CP/E) 1 Q(N) = E*F AINT = Q(1) + 4.*Q(64) + Q(65) DO 2 N = 1,31 K = 2*N 2 AINT = AINT + 4.*Q(K) + 2.*Q(K+1) FMN = CP*AINT*DXI/3. E = DSQRT(1.+X*X*SP*SP) F = DATAN(X*CP/E) + DATAN((Y-X*CP)/E) FMN = FMN + SP*S2P/2.*X*E*F FMN = FMN + X*DATAN(1./X) + Y*DATAN(1./Y) & - DSQRT(Z)*DATAN(1./DSQRT(Z)) E = 1. + X*X F = 1. + Y*Y FMN = FMN + (0.5 - SP*SP/4.)*DLOG(E*F/(1.+Z)) FMN = FMN + SP*SP/4.*(Y*Y*DLOG(Y*Y*(1.+Z)/F/Z) & + X*X*(DLOG(X*X/Z) + C2P*DLOG(E/(1.+Z)))) E = Y*Y*DATAN((X-Y*CP)/Y/SP) + X*X*DATAN((Y-X*CP)/X/SP) FMN = FMN - S2P/4.*(X*Y*SP + (3.14159265/2. - PHI) & * (X*X + Y*Y) + E) FMN = FMN/3.14159265/Y RETURN END C*********************************************************************** C** ** C** FUNCTION FP ** C** ** C*********************************************************************** FUNCTION FP(A,B,C) C THIS FUNCTION CALCULATES THE VIEW FACTOR BETWEEN TWO FINITE EQUAL C PARALLEL PLATES C SEE APPENDIX A OF SPARROW AND CESS, CONFIGURATION 1 C A AND B = DIMENSIONS OF PLATES, C = SEPARATION OF PLATES IMPLICIT REAL*8 (A-H, O-Z) X = A/C Y = B/C D = DSQRT(1.+X*X) E = DSQRT(1.+Y*Y) F = DSQRT(1.+X*X+Y*Y) FP = DLOG(D*E/F) + Y*D*DATAN(Y/D) + X*E*DATAN(X/E) & - Y*DATAN(Y) - X*DATAN(X) FP = FP*2./(3.14159265*X*Y) RETURN END C*********************************************************************** C** ** C** SUBROUTINE HGap "Open Cavity" ** C** ** C*********************************************************************** SUBROUTINE HGapOpn(TS1,TA,TS2,TO,EO,PHI,Width,Length,GAP,IFLAG,V, & HC,AMDOT,NUSx) C THIS SUBROUTINE CALCULATES NATURAL CONVECTION HEAT TRANSFER C COEFFICIENTS BASED ON ISOLATED ISOTHERMAL FLAT PLATES C CALCULATES FORCED CONVECTION COEFFICIENTS AND TOTAL C CONVECTION HEAT TRANSFER COEFFICIENTS C SEE HOLMAN, 5TH EDITION, PP. 272-286 C TS = SURFACE TEMPERATURE, F C TA = AIR TEMPERATURE, F C PHI = TILT ANGLE, DEGREES, 0 FOR HORIZONTAL, 90 FOR VERTICAL C AW = CHARACTERISTIC WIDTH OF AIRGAP C AL = CHARACTERISTIC LENGTH OF AIRGAP C GAP = HEIGHT OF AIRGAP C IFLAG = 1 FOR SURFACE FACING UPWARD C IFLAG = 2 FOR SURFACE FACING DOWNWARD C V = AIR SPEED, FEET PER HOUR C HCF = FORCED CONVECTION COEFFICIENT C HCN = NATURAL CONVECTION COEFFICIENT C HC = TOTAL CONVECTION COEFFICIENT C IMPLICIT REAL*8 (A-H, O-Z) REAL*8 NUS,K,MU,NU, Length,NUSx Data PI /3.14159265/ C IF(IFLAG.EQ.1) DT = TA - TS2 IF(IFLAG.EQ.2) DT = TS1 - TA C C CALCULATE FILM TEMPERATURE IF (IFLAG .EQ. 1) TF = (TS2 + TA)/2.0 IF (IFLAG .EQ. 2) TF = (TS1 + TA)/2.0 TK = (TF+459.67)/1.8 C THERMAL CONDUCTIVITY OF AIR FROM EQUATION IN NBS CIRC. 564 C THERMAL COND. IN BTU/(HR-FT-F) K = 0.6325D-5*DSQRT(TK)/(1.+(245.4*10.**(-12./TK))/TK)*241.77 C DYNAMIC VISCOSITY OF AIR FROM EQUATION IN NBS CIRC. 564 C LB/(HR-FT) MU = (145.8*TK*DSQRT(TK)/(TK+110.4))*241.90D-7 C PRANDTL NUMBER, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564 PR = 0.7880 - 2.631D-4*TK C VOLUME EXPANSION COEFFICIENT OF AIR, 1/TABS, PERFECT GAS BETA = 1./(TF+459.67) C DENSITY OF AIR, PERFECT GAS, NBS CIRC. 564, LB/CF RHO = 22.0493/TK C KINEMATIC VISCOSITY, FT2/HR NU = MU/RHO C SPECIFIC HEAT OF AIR, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564, BTU/(LB-F) CP = (3.4763 + 1.066D-4*TK)*0.068559 C RAYLEIGH NUMBER, LEADING COEFFICIENT IS C 32.174*3600*3600, FT/HR2 RAf =(4.16975E8)*BETA*RHO*CP*(GAP**3)/NU/K C IF(ABS(DT) .LE. 2.0) THEN NUS = 1.0 HCN = NUS*K/GAP NUSx = NUS GO TO 20 END IF C IF(DT.GE. 0.0) GO TO 10 C Holland's Nusselt equation for heating mode RAh = RAf*ABS(DT)*DCOS(PI*PHI/180.) X = 0.0 if(Rah .gt. 0.0) X = 1.0 - 1708.0/RAh IF(X .LE. 0.E0) X = 0.0 Y = ((RAh/5830.)**(1./3.) - 1.0) IF(Y .LE. 0.E0) Y = 0.0 if(Rah .gt. 0.0) then Z = (1.0-1708.0*(DSIN(1.8*PI*PHI/180.))**(1.6)/RAh) else Z = 0.0 end if NUS = 1.0 + 1.44*X*Z + Y HCN = NUS*K/GAP NUSx = NUS GO TO 20 10 CONTINUE C AZEVEDO & SPARROW Nusselt equation for open channel RAaz = RAf*ABS(DT)*DSIN(PI*PHI/180.) NUS = 0.673*((Gap/Length)*RAaz)**0.30 !0.25 HCN = NUS*K/GAP NUSx = NUS C C Brinkworth's Nusselt correlation for cooling mode c WP = 2.0*(GAP+Width) c DH = 4.0*(GAP*Width)/WP c RE = 4.0*AMDOT/(MU*WP) c NUS = 5.385 + 0.071*(Length/(DH*RE)+0.002) c HCN = NUS*K/DH c NUSx = NUS c IF(DT.GE. 0.0) GO TO 10 c C ELENBAAS CORRELATION C RAeb = (Gap/Length)*RA*SIN(PI*PHI/180.) C NUS = (1.0/24.0)*RAeb*(1 - EXP(-35.0/RAeb))**0.75 C HCN = NUS*K/GAP c C CALCULATE FORCED CONVECTION COEFFICIENT C SEE HOLMAN, 5TH EDITION, PG. 191, EQN. 5-44, 5-46 C AND PG. 202, EQN. 5-85 20 RE = V*Length/NU WP = 2.0*(GAP+Width) DH = 4.0*(GAP*Width)/WP RE = V*DH/NU IF(RE.LT.2.2E3) NUS = 0.023*(PR**(0.4))*RE**0.8 IF(RE.GT.2.2E3) NUS = 0.023*(PR**(0.4))*RE**0.8 HCF = NUS*K/Length C COMBINE NATURAL AND FORCED CONVECTION COEFFICIENTS C USING CHURCHILL'S CORRELATION C SEE J. HEAT TRANS., VOL. 108, PP. 835-840 (1986) C ASSUME ASSISTING FLOW IN ALL CASES HC = (HCF**3 + HCN**3)**(1./3.) RETURN END C*********************************************************************** C** ** C** SUBROUTINE HGap "Closed Cavity" ** C** ** C*********************************************************************** SUBROUTINE HGapCld(TS1,TA,TS2,TO,PHI,Width,Length,GAP,IFLAG,V, & HCN,AMDOT,NUSx) C THIS SUBROUTINE CALCULATES NATURAL CONVECTION HEAT TRANSFER C COEFFICIENTS BASED ON ISOLATED ISOTHERMAL FLAT PLATES C CALCULATES FORCED CONVECTION COEFFICIENTS AND TOTAL C CONVECTION HEAT TRANSFER COEFFICIENTS C SEE HOLMAN, 5TH EDITION, PP. 272-286 C TS = SURFACE TEMPERATURE, F C TA = AIR TEMPERATURE, F C PHI = TILT ANGLE, DEGREES, 0 FOR HORIZONTAL, 90 FOR VERTICAL C AW = CHARACTERISTIC WIDTH OF AIRGAP C AL = CHARACTERISTIC LENGTH OF AIRGAP C GAP = HEIGHT OF AIRGAP C IFLAG = 1 FOR SURFACE FACING UPWARD C IFLAG = 2 FOR SURFACE FACING DOWNWARD C V = AIR SPEED, FEET PER HOUR C HCF = FORCED CONVECTION COEFFICIENT C HCN = NATURAL CONVECTION COEFFICIENT C HC = TOTAL CONVECTION COEFFICIENT C IMPLICIT REAL*8 (A-H, O-Z) REAL*8 NUS,K,MU,NU, Length,NUSx Data PI /3.14159265/ C C C CALCULATE FILM TEMPERATURE TF = (TS1 + TS2)/2.0 TK = (TF+459.67)/1.8 C THERMAL CONDUCTIVITY OF AIR FROM EQUATION IN NBS CIRC. 564 C THERMAL COND. IN BTU/(HR-FT-F) K = 0.6325D-5*DSQRT(TK)/(1.+(245.4*10.**(-12./TK))/TK)*241.77 C DYNAMIC VISCOSITY OF AIR FROM EQUATION IN NBS CIRC. 564 C LB/(HR-FT) MU = (145.8*TK*DSQRT(TK)/(TK+110.4))*241.90D-7 C PRANDTL NUMBER, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564 PR = 0.7880 - 2.631D-4*TK C VOLUME EXPANSION COEFFICIENT OF AIR, 1/TABS, PERFECT GAS BETA = 1./(TF+459.67) C DENSITY OF AIR, PERFECT GAS, NBS CIRC. 564, LB/CF RHO = 22.0493/TK C KINEMATIC VISCOSITY, FT2/HR NU = MU/RHO C SPECIFIC HEAT OF AIR, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564, BTU/(LB-F) CP = (3.4763 + 1.066D-4*TK)*0.068559 C RAYLEIGH NUMBER, LEADING COEFFICIENT IS C 32.174*3600*3600, FT/HR2 DT = TS1 - TS2 C C Holland's Nusselt equation RA =(4.16975E8)*BETA*RHO*CP*ABS(DT)*DCOS(PI*PHI/180.) & *(GAP**3)/NU/K X = 0.0 if(Ra .gt. 0.0) X = 1.0 - 1708.0/RA IF(X .LE. 0.E0) X = 0.0 Y = ((RA/5830.)**(1./3.) - 1.0) IF(Y .LE. 0.E0) Y = 0.0 if(Ra .gt. 0.0) then Z = (1.0-1708.0*(DSIN(1.8*PI*PHI/180.))**(1.6)/RA) else Z = 0.0 end if NUS = 1.0 + 1.44*X*Z + Y c write(6,202) x,y,z,RA,NUS c202 format(' X = ',f10.3,/, c & ' Y = ',f10.3,/, c & ' Z = ',f10.3,/, c & ' Ra = ',f10.1,/, c & 'NUS = ',f10.1) c pause HCN = NUS*K/GAP NUSx = NUS RETURN END C*********************************************************************** C** ** C** SUBROUTINE HCON ** C** ** C*********************************************************************** SUBROUTINE HCON(TS,TA,PHI,AL,IFLAG,Isurf,V,HC) C THIS SUBROUTINE CALCULATES NATURAL CONVECTION HEAT TRANSFER C COEFFICIENTS BASED ON ISOLATED ISOTHERMAL FLAT PLATES C CALCULATES FORCED CONVECTION COEFFICIENTS AND TOTAL C CONVECTION HEAT TRANSFER COEFFICIENTS C SEE HOLMAN, 5TH EDITION, PP. 272-286 C TS = SURFACE TEMPERATURE, F C TA = AIR TEMPERATURE, F C PHI = TILT ANGLE, DEGREES, 0 FOR HORIZONTAL, 90 FOR VERTICAL C AL = CHARACTERISTIC LENGTH OF SURFACE C IFLAG = 1 FOR SURFACE FACING UPWARD C IFLAG = 2 FOR SURFACE FACING DOWNWARD C V = AIR SPEED, FEET PER HOUR C HCF = FORCED CONVECTION COEFFICIENT C HCN = NATURAL CONVECTION COEFFICIENT C HC = TOTAL CONVECTION COEFFICIENT C RExc = Critical Reynolds number for transition to turbulence C RExc(1) = 500,000 for smooth plate; Isurf = 1 C RExc(2) = 300,000 for course surface:Isurf = 2 C RExc(3) = 100,000 for rough surface; Isurf = 3 C RExc(4) = 50,000 for rough surface; Isurf = 4 C RExc(5) = 0 for blunt surface; Isurf = 5 C IMPLICIT REAL*8 (A-H, O-Z) REAL*8 NUS,K,MU,NU DIMENSION RExc(5) DATA RExc(1), RExc(2), RExc(3), RExc(4), RExc(5) & / 500000., 300000., 100000., 50000., 0.0 / C DT = TS - TA IF (IFLAG.EQ.2) DT = -DT C CALCULATE FILM TEMPERATURE TF = (TS+TA)/2.0 TF1 = TF IF(ABS(PHI).GT.1.D-3.AND.ABS(PHI-90.).GT.1.D-3) & TF = TS - 0.25*(TS-TA) IF(ABS(PHI).GT.1.D-3.AND.ABS(PHI-90.).GT.1.D-3) & TF1 = TA + 0.25*(TS-TA) TK = (TF+459.67)/1.8 C THERMAL CONDUCTIVITY OF AIR FROM EQUATION IN NBS CIRC. 564 C THERMAL COND. IN BTU/(HR-FT-F) K = 0.6325D-5*DSQRT(TK)/(1.+(245.4*10.**(-12./TK))/TK)*241.77 C DYNAMIC VISCOSITY OF AIR FROM EQUATION IN NBS CIRC. 564 C LB/(HR-FT) MU = (145.8*TK*DSQRT(TK)/(TK+110.4))*241.90D-7 C PRANDTL NUMBER, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564 PR = 0.7880 - 2.631D-4*TK C VOLUME EXPANSION COEFFICIENT OF AIR, 1/TABS, PERFECT GAS BETA = 1./(TF1+459.67) C DENSITY OF AIR, PERFECT GAS, NBS CIRC. 564, LB/CF RHO = 22.0493/TK C KINEMATIC VISCOSITY, FT2/HR NU = MU/RHO C SPECIFIC HEAT OF AIR, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564, BTU/(LB-F) CP = (3.4763 + 1.066D-4*TK)*0.068559 C RAYLEIGH NUMBER, LEADING COEFFICIENT IS C 32.174*3600*3600, FT/HR2 RA = (4.16975E8)*BETA*RHO*CP*ABS(DT)*(AL**3)/NU/K C BRANCH TO DIFFERENT CORRELATIONS DEPENDING UPON C SURFACE ORIENTATION IF(ABS(PHI).LE.1.D-3) GO TO 100 IF(ABS(PHI-90.).LE.1.D-3) GO TO 200 IF(ABS(PHI).GT.1.D-3.AND.ABS(PHI).LT.2.) GO TO 300 IF(ABS(PHI).GT.2.0.AND.ABS(PHI-90.).GT.1.D-3) GO TO 400 C FOR HORIZONTAL SURFACES 100 IF(DT.LT.0.0) GO TO 150 NUS = 0.15*RA**(1./3.) IF(RA.LT.8.E6) NUS = 0.54*RA**0.25 GO TO 1000 150 NUS = 0.58*RA**0.2 GO TO 1000 C FOR VERTICAL SURFACES 200 NUS = 0.10*RA**(1./3.) IF(RA.LT.1.E9) NUS = 0.59*RA**0.25 GO TO 1000 C FOR NEARLY HORIZONTAL(UP TO 2 DEGREES TILT) SURFACES 300 IF(DT.GT.0.0) GO TO 450 NUS = 0.58*RA**0.2 GO TO 1000 C FOR TILTED SURFACES 400 IF(DT.GT.0.0) GO TO 450 NUS = 0.56*(RA*DCOS((90.-PHI)*3.14159265/180.))**0.25 GO TO 1000 450 GRC = 10.0**(PHI/(1.1870+0.0870*PHI)) IF(ABS(PHI).LT.15.) GRC = 1.E6 IF(ABS(PHI).GT.75.) GRC = 5.E9 GR = RA/PR IF(GR.LE.GRC) NUS=0.56*(RA*DCOS((90.-PHI)*3.14159265/180.))**0.25 IF(GR.GT.GRC) NUS = 0.14*(RA**(1./3.) - (GRC*PR)**(1./3.)) & +0.56*(GRC*PR*DCOS((90.-PHI)*3.14159265/180.))**0.25 GO TO 1000 1000 IF (AL .LE. 0.042) then ! CHECK FOR EAVE HEIGHT (if < 0.5-in) set NCN = 0.0 HCN = 0.0 else HCN = NUS*K/AL END IF c write(6,10) DT,RA,TS,TA,NUS,HCN c10 format(' DT = ',1pE12.3,/, c & ' ra = ', e12.3,/, c & 'TS = ', e12.3,/, c & 'TA = ', e12.3,/, c & 'NU = ', e12.3,/, c & 'hcv = ', e12.3) c pause C CALCULATE FORCED CONVECTION COEFFICIENT C SEE HOLMAN, 5TH EDITION, PG. 191, EQN. 5-44, 5-46 C AND PG. 202, EQN. 5-85 RE = V*AL/NU IF(RE.LT.RExc(Isurf)) then NUSf = 0.664*(PR**(1./3.))*DSQRT(RE) GO TO 2000 END IF Afactor = 0.037*RExc(Isurf)**0.8 - 0.664*DSQRT(RExc(Isurf)) NUSf = (PR**(1./3.))*(0.037*(RE**0.8) - Afactor) 2000 IF (AL .LE. 0.042) THEN ! CHECK FOR EAVE HEIGHT (if < 0.5-in) set HDF = 0.0 HCF = 0.0 ELSE HCF = NUSf*K/AL END IF C COMBINE NATURAL AND FORCED CONVECTION COEFFICIENTS C USING CHURCHILL'S CORRELATION C SEE J. HEAT TRANS., VOL. 108, PP. 835-840 (1986) C ASSUME ASSISTING FLOW IN ALL CASES HC = (HCF**3 + HCN**3)**(1./3.) c write(6,10) RA,NUS,NUSf c10 format(' ra = ',1pe12.3,/, c & ' NUS = ', e12.3,/, c & ' NUSf= ', e12.3,/) RETURN END C*********************************************************************** C** ** C** SUBROUTINE HMASS ** C** ** C*********************************************************************** SUBROUTINE HMASS(TS,TA,HC,HM) C THIS SUBROUTINE CALCULATES THE MASS TRANSFER COEFFICIENT FOR C MOISTURE TRANSFER BETWEEN A SURFACE AND MOIST AIR C USING THE LEWIS RELATIONSHIP, AND EVALUATING ALL THERMOPHYSICAL C PROPERTIES AS THOSE OF DRY AIR C SEE HOLMAN, 5TH EDITION, PP. 494 C TS = SURFACE TEMPERATURE, F C TA = AIR TEMPERATURE, F C HC = CONVECTION HEAT TRANSFER COEFFICIENT C HM = MASS TRANSFER COEFFICIENT, LB/(HR-FT2) C IMPLICIT REAL*8 (A-H, O-Z) REAL*8 K DIMENSION RExc(5) C CALCULATE FILM TEMPERATURE TF = (TS+TA)/2.0 TK = (TF+459.67)/1.8 C THERMAL CONDUCTIVITY OF AIR FROM EQUATION IN NBS CIRC. 564 C THERMAL COND. IN BTU/(HR-FT-F) K = 0.6325D-5*DSQRT(TK)/(1.+(245.4*10.**(-12./TK))/TK)*241.77 C DENSITY OF AIR, PERFECT GAS, NBS CIRC. 564, LB/CF RHO = 22.0493/TK C SPECIFIC HEAT OF AIR, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564, BTU/(LB-F) CP = (3.4763 + 1.066D-4*TK)*0.068559 C CALCULATE THERMAL DIFFUSIVITY, FT2/HR ALF = K/CP/RHO C CALCULATE DIFFUSION COEFFICIENT FOR WATER VAPOR THROUGH AIR C D = FT2/HR, SEE 1985 ASHRAE HANDBOOK, PG. 5.2 D = 0.035883/101.32*(TK**2.5)/(TK+245.0) HM = HC/CP/((ALF/D)**(2.0/3.0)) RETURN END C*********************************************************************** C** ** C** FUNCTION HRAD ** C** ** C*********************************************************************** FUNCTION HRAD(G,T1,T2) C*********************************************************************** C* FUNCTION: HRAD C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculates radiative coffficients for heat balance C* equations. C*********************************************************************** C* INPUT VARIABLES C* G Radiation interchange factor, C* including effects of surface C* emittances and view factors (dimensionless) C* T1 Temperature of surface 1 (°F) C* T2 Temperature of surface 2 (°F) C* C* OUTPUT VARIABLES C* HRAD Radiation heat transfer coefficient (Btu/h·ft²·°F) C*********************************************************************** C* MAJOR RESTRICTION: None C* C* DEVELOPER: Kenneth E. Wilkes, PhD, PE C* C* INCLUDE FILES: None C* C* SUBROUTINES CALLED: None C* C* FUNCTIONS REQUIRED: None C* C* REVISION HISTORY: None C* C* REFERENCE: None C* C*********************************************************************** C* INTERNAL VARIABLES C* T1R Absolute temperature of surface 1, divided by 100 C* T2R Absolute temperature of surface 2, divided by 100 C*********************************************************************** IMPLICIT REAL*8 (A-H, O-Z) T1R = (T1 + 459.67)/100. T2R = (T2 + 459.67)/100. HRAD = 1.712D-3*(T1R*T1R+T2R*T2R)*(T1R+T2R)*G RETURN END C*********************************************************************** C** ** C** SUBROUTINE MOIST ** C** ** C*********************************************************************** SUBROUTINE MOIST(TS,TI,TA,HCI,A,AWRAT,AMC,HUMI,PERM, & AMDOT,EXFIL,AMW,BMW,HM,K) C NEXT LINE ADDED FOR VERSION 1B IMPLICIT REAL*8 (A-H, O-Z) COMMON /HUMID/WO,PATM COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS COMMON /TRUS2/AMASST,AMCT,HCTRUS,AMWT,BMWT,QLAT CC THIS SUBROUTINE PERFORMS A MOISTURE BALANCE ON THE ATTIC SPACE C TS = SURFACE TEMPERATURE, DEGREES F C TO = OUTSIDE AIR TEMPERATURE C TI = INSIDE AIR TEMPERATURE C TA = ATTIC AIR TEMPERATURE C HCI = CONVECTION HEAT TRANSFER COEFFICIENTS C A = PROJECTED SURFACE AREAS C AWRAT = RATIO OF EXPOSED WOOD SURFACE AREAS TO PROJECTED AREAS C AMC = WOOD MOISTURE CONTENTS, FRACTION OF DRY WEIGHT C HUMI = RELATIVE HUMIDITY OF INSIDE AIR C PERM = WATER VAPOR PERMEANCES C AMDOT = MASS FLOW OF VENTILATION AIR C EXFIL = MASS FLOW OF AIR FROM HOUSE INTO ATTIC C AMW = MOISTURE ADSORBED PER UNIT PROJECTED AREA C BMW = COEFFICIENT OF LINEAR TERM IN TAYLOR SERIES C EXPANSION OF WS(I) C K = NUMBER OF SURFACES DIMENSION TS(K),HCI(K),A(K),AWRAT(K),AMC(K),PERM(K),AMW(K) DIMENSION WS(10),P(10),HM(10),BMW(K) EPS = 1.D-3 C CALCULATE HUMIDITY RATIOS AT WOOD SURFACES DO 10 I = 1,K 10 WS(I) = WDHUM(AMC(I),TS(I)) IF(ITRUS.NE.0) WSTRUS = WDHUM(AMCT,TTRUS(1)) C CALCULATE HUMIDITY RATIOS OF INSIDE AND OUTSIDE AIR CALL PSY(TI,HUMI,PATM,PI,WI) XX = WO/0.62198 PO = PATM*XX/(1.0 + XX) C CALCULATE MASS TRANSFER COEFFICIENTS DO 20 I = 1,K 20 CALL HMASS(TS(I),TA,HCI(I),HM(I)) IF(ITRUS.NE.0) THEN CALL HMASS(TTRUS(1),TA,HCTRUS,HMTRUS) END IF C SET UP WATER VAPOR PARTIAL PRESSURES P(1) = PI DO 30 I = 2,K 30 P(I) = PO C CALCULATE HUMIDITY RATIO OF ATTIC AIR FROM MOISTURE BALANCE PA = PO X = 0.0 DO 40 I = 1,K X = X + A(I)*PERM(I)*P(I) 40 X = X + A(I)*AWRAT(I)*HM(I)*WS(I) X = X + AMDOT*WO + EXFIL*WI IF(ITRUS.NE.0) X = X + ATRUS*HMTRUS*WSTRUS NIT = 0 100 Y = 0.0 DO 50 I = 1,K 50 Y = Y + A(I)*AWRAT(I)*HM(I) Y = Y + AMDOT + EXFIL IF(ITRUS.NE.0) Y = Y + ATRUS*HMTRUS C IN FOLLOWING, PATM HAS BEEN SUBSTITUTED FOR 14.696 IN VERSION 1B DO 60 I = 1,K C 60 Y = Y + A(I)*PERM(I)*(14.696 - PA)/0.62198 60 Y = Y + A(I)*PERM(I)*(PATM - PA)/0.62198 WA = X/Y XX = WA/0.62198 C PANEW = 14.696*XX/(1.0 + XX) PANEW = PATM*XX/(1.0 + XX) NIT = NIT + 1 IF (NIT.GT.10) GO TO 200 IF(ABS(1.-PANEW/PA).LT.EPS) GO TO 200 PA = PANEW GO TO 100 200 CONTINUE C CALCULATE WATER ADSORBED BY EACH SURFACE DO 300 I = 1,K 300 AMW(I) = HM(I)*AWRAT(I)*(WA-WS(I)) IF(ITRUS.NE.0) AMWT = HMTRUS*ATRUS*(WA - WSTRUS) C CALCULATE CORRECTION TERM FOR TAYLOR SERIES EXPANSION OF C ADSORBED WATER DO 310 I = 1,K 310 BMW(I) = HM(I)*AWRAT(I)*WS(I)/28.6 IF(ITRUS.NE.0) BMWT = HMTRUS*ATRUS*WSTRUS/28.6 RETURN END C*********************************************************************** C** ** C** SUBROUTINE PSY ** C** ** C*********************************************************************** SUBROUTINE PSY(T,HUM,PATM,P,W) C THIS SUBROUTINE CALCULATES THE WATER VAPOR PARTIAL PRESSURE AND C THE HUMIDITY RATIO, GIVEN THE TEMPERATURE AND RELATIVE HUMIDITY C USING EQUATION FROM 1985 ASHRAE FUNDAMENTALS P.6.6 C AND ERRATA FROM 1988 ASHRAE HANDBOOK C T = DRYBULB TEMPERATURE, DEGREES F C HUM = RELATIVE HUMIDITY, PERCENT C P = WATER VAPOR PARTIAL PRESSURE, PSIA C W = HUMIDITY RATIO, POUNDS OF WATER PER POUNDS OF DRY AIR IMPLICIT REAL*8 (A-H, O-Z) C1 = -10214.16 C2 = -4.8932631 C3 = -0.53769056E-2 C4 = 0.19202377E-6 C5 = 0.35575832E-9 C6 = 0.090344688E-12 C7 = 4.1635019 C8 = -10440.397 C9 = -11.2946496 C10 = -0.027022355 C11 = 0.12890360E-4 C12 = -0.2478068E-8 C13 = 6.5459673 TR = T + 459.67 IF(T.GE.32.0) GO TO 10 X = C1/TR + C2 + C3*TR + C4*TR*TR +C5*TR**3+C6*TR**4+C7*DLOG(TR) GO TO 20 10 X = C8/TR + C9 + C10*TR + C11*TR*TR + C12*TR**3 + C13*DLOG(TR) 20 PSAT = EXP(X) P = (HUM/100.0)*PSAT W = 0.62198*P/(PATM - P) RETURN END C*********************************************************************** C** ** C** SUBROUTINE Humidity Ratio ** C** ** C*********************************************************************** SUBROUTINE xMoist(T,RH,PATM,DewPt,HUM,Pw) C THIS SUBROUTINE CALCULATES THE WATER VAPOR PARTIAL PRESSURE AND C THE DEW POINT AND HUMIDITY RATIO GIVEN THE TEMPERATURE AND RELATIVE HUMIDITY C USING EQUATION FROM 1985 ASHRAE FUNDAMENTALS P.6.6 C AND ERRATA FROM 1988 ASHRAE HANDBOOK C T = DRYBULB TEMPERATURE, DEGREES F C RH = RELATIVE HUMIDITY, PERCENT C P = WATER VAPOR PARTIAL PRESSURE, PSIA C HUM = HUMIDITY RATIO, POUNDS OF WATER PER POUNDS OF DRY AIR IMPLICIT REAL*8 (A-H, O-Z) C1 = -10214.16462D0 ! Pws Over Ice C2 = -4.89350301D0 C3 = -0.00537657944D0 C4 = 0.000000192023769D0 C5 = 3.55758316D-10 C6 = -9.03446883D-14 C7 = 4.1635019D0 C8 = -10440.397D0 C9 = -11.2946496D0 C10 = -0.027022355D0 C11 = 0.12890360D-4 C12 = -0.2478068D-8 C13 = 6.5459673D0 C13 = 6.5459673D0 C14 = 90.12D0 ! Tdp over Ice C15 = 26.412D0 C16 = 0.8927D0 C17 = 100.45D0 ! Tdp over Liquid water C18 = 33.193D0 C19 = 2.319D0 C20 = 0.17074D0 C21 = 1.2063D0 C TR = T + 459.67 P = 0.491154*Patm C IF(T.GE.32.0) GO TO 10 X = C1/TR + C2 + C3*TR + C4*TR*TR +C5*TR**3 + C6*TR**4 & + C7*DLOG(TR) Pws = Exp(X) Pw = Pws * RH / 100.0D0 SpHUM = 0.62198 * (Pw / (P - Pw)) C alpha = DLOG(Pw) DewPt = C14 + C15 * alpha + C16 * alpha ** 2 C GO TO 20 10 X = C8/TR + C9 + C10*TR + C11*TR*TR + C12*TR**3 + C13*DLOG(TR) Pws = Exp(X) Pw = Pws * RH / 100.0D0 SpHUM = 0.62198 * (Pw / (P - Pw)) alpha = DLOG(Pw) DewPt = C17 + C18 * alpha + C19 * alpha ** 2 + C20 * alpha ** 3 & + C21 * (Pw) ** 0.1984D0 20 CONTINUE RETURN END C*********************************************************************** SUBROUTINE SOLVP (I,C,D,X,NM) C MODIFIED FOR DUCT VERSION TO ADD DIMENSION NM C SUBROUTINE TAKEN FROM NBSLD C THIS SUBROUTINE SOLVES A SET OF SIMULTANEOUS LINEAR EQUATIONS C USING GAUSSIAN ELIMINATION C I = NUMBER OF EQUATIONS TO BE SOLVED C C = MATRIX OF COEFFICIENTS OF SET OF EQUATIONS C D = VECTOR ON RIGHT HAND SIDE OF EQUATIONS C X = SOLUTION VECTOR, CX = D IMPLICIT REAL*8 (A-H, O-Z) DIMENSION A(100,101),C(NM,NM),D(NM),X(NM) M = I N = M + 1 DO 10 IX = 1,M DO 10 IY = 1,M 10 A(IX,IY) = C(IX,IY) DO 20 IZ = 1,M 20 A(IZ,N) = D(IZ) L = 1 30 AA = A(L,L) if(AA .eq. 0.0) go to 65 DO 40 K = L,N A(L,K) = A(L,K)/AA c write(6,100) L,K,A(l,k) c100 format('A(',i2,',',i2,') = ',1(1Pe12.3)) c pause 40 continue DO 60 K = 1,M IF (K.EQ.L) GO TO 60 AA = -A(K,L) DO 50 IA = L,N 50 A(K,IA) = A(K,IA) + AA*A(L,IA) 60 CONTINUE 65 CONTINUE L = L + 1 IF (L.LE.M) GO TO 30 DO 70 IP = 1,M 70 X(IP) = A(IP,N) RETURN END C*********************************************************************** C** ** C** SUBROUTINE VENT ** C** ** C*********************************************************************** SUBROUTINE VENT(TVENT,TA,WS,DIR,H,AI,AO,ITYPE,AMCDOT,AMDOT,V) C THIS SUBROUTINE CALCULATES THE VENTILATION RATE FOR A GABLED ATTIC C TVENT = VENTILATION AIR TEMPERATURE, F C TVENT = OUTDOOR AIR TEMPERATURE UNLESS INPUT OTHERWISE C TA = AVERAGE ATTIC AIR TEMPERATURE, F C WS = WIND SPEED, MPH C DIR = WIND DIRECTION C H = THE ATTIC INLET TO OUTLET HEIGHT C AI = AREA OF VENT INLET, SQUARE FEET C AO = AREA OF VENT OUTLET, SQUARE FEET C ITYPE = 1, SOFFIT AND RIDGE VENTS C ITYPE = 2, SOFFIT AND GABLE VENTS C ITYPE = 3, SOFFIT VENTS C AMCDOT = AIR MASS FLOW RATE TIMES SPECIFIC HEAT, BTU/HR-F C V = VOLUME FLOW RATE, CUBIC FEET PER HOUR C CALCULATE ABSOLUTE TEMPERATURES, R IMPLICIT REAL*8 (A-H, O-Z) TVR = TVENT + 459.67 TAR = TA + 459.67 C CALCULATE DENSITIES AND SPECIFIC HEATS OF AIR RHOO = 22.0493/(TVR/1.8) RHOA = 22.0493/(TAR/1.8) CP = (3.4763 + 1.066D-4*(TAR/1.8))*0.068559 c write(12,100) Tvent,TVR,RHOO,TA,TAR,RHOA c100 format(//,' Tvent = ',2f10.2,1x,f10.4,/, c & ' TA = ',2f10.2,1x,f10.4,/) C DETERMINE MINIMUM VENT AREA AND MAXIMUM VENT AREAS AMIN = DMIN1(AI,AO) AMAX = DMAX1(AI,AO) C CALCULATE THE EFFECTS OF WIND DIRECTION ON WIND AT VENT Arg = DABS(DSIN((DIR + 90.0D0)*3.14159265/180.0D0)) if (Arg .eq. 0.0) then WSi = 0.087*WS else WSi = 0.45*WS*(0.087 + & + 0.130*Arg**2.5) end if C CALCULATE HEIGHT OF NEUTRAL PRESSURE LEVEL FROM LOWER OPENING C SEE PG. 22.3 OF 1985 ASHRAE HANDBOOK OF FUNDAMENTALS HNPL = H/(1.+AI*AI/AO/AO*TAR/TVR) IF(TAR.LT.TVR) HNPL = H/(1.+AI*AI/AO/AO*TVR/TAR) c write(12,110) AMIN,AMAX,H,HNPL c110 format(' AMIN, AMAX = ',2f10.2,/, c & ' NEUTRAL Pr = ',2f10.3,/) C CALCULATE FLOW DUE TO STACK EFFECT, SEE ASHRAE PG. 22.7 IF(TAR.GE.TVR) QSTACK = 0.65*AMIN*DSQRT(2.0*32.174*3600.*3600. & *HNPL*(TAR-TVR)/TAR) IF(TAR.GE.TVR) AMSTACK = QSTACK*RHOO IF(TAR.LT.TVR) QSTACK = 0.65*AMIN*DSQRT(2.0*32.174*3600.*3600. & *HNPL*(TVR-TAR)/TVR) IF(TAR.LT.TVR) AMSTACK = QSTACK*RHOA IF(ITYPE.EQ.3) QSTACK = 0.0 IF(ITYPE.EQ.3) AMSTACK = 0.0 C CALCULATE FLOW DUE TO WIND PRESSURE, SEE ASHRAE PG. 22.6 QWIND = 5280.0*0.6*AMIN*WSi c write(12,120) Qstack,AMstack,Qwind c120 format(' Qstack AMstack Qwind',/,1x,3f10.2) C MULTIPLY BY EFFECTIVENESS COEFFICIENT DERIVED FROM C BURCH AND TREADO'S DATA FOR ATTIC VENTILATION C TYPE 1 = SOFFIT AND RIDGE VENTS C TYPE 2 = SOFFIT AND ROOF VENTS, ASSUMED TO BE SAME AS SOFFIT AND C GABLE VENTS C TYPE 3 = SOFFIT VENTS ONLY IF(ITYPE.EQ.1) CF = 0.38 IF(ITYPE.EQ.2) CF = 0.54 IF(ITYPE.EQ.3) THEN IF(DIR.EQ.0.OR.DIR.EQ.180.) CF = 0.089 IF(DIR.NE.0.0.AND.DIR.NE.180.0) & CF = 0.089 + 0.132*(DABS(DSIN(DIR*3.14159265/180.D0)))**2.5 END IF QWIND = QWIND/0.6*CF AMWIND = QWIND*RHOO c write(12,130) WS,WSi,Qwind,AMWIND c130 format(' WS, WSi, Qwind, AMWIND = ',4f10.3,/) C COMBINE STACK AND WIND FLOWS, SEE ASHRAE PG. 22.5 AMDOT = DSQRT(AMSTACK*AMSTACK + AMWIND*AMWIND) C CORRECT FOR UNEQUAL INLET AND OUTLET AREAS (ASHRAE PG. 22.7) ARAT = AMIN/AMAX AMDOT = AMDOT * (1.0 + 0.4077*(1.0-ARAT**1.5)) AMCDOT = CP*AMDOT C CALCULATE VOLUME FLOW RATE, CUBIC FEET PER HOUR V = AMDOT/RHOA c write(12,140) amstack,amwind,amdot,v c140 format(' AMstack AMwind AMdot V',/,1x,4f10.2) c pause RETURN END C*********************************************************************** C** ** C** SUBROUTINE VIEW2 ** C** ** C*********************************************************************** SUBROUTINE VIEW2(AL,W,PITCH1,PITCH2,H1,EI,G,K,F) C THIS SUBROUTINE CALCULATES THE VIEW FACTORS AMONG THE SURFACES C IN A GABLED ATTIC, HAVING A FIVE-SIDED CROSS-SECTION, AND C HAVING ARBITRARILY PITCHED ROOF SURFACES C AL = LENGTH OF ATTIC C W = WIDTH OF ATTIC C PITCH1 = PITCH OF SOUTH ROOF, DEGREES C PITCH2 = PITCH OF NORTH ROOF, DEGREES C H1 = HEIGHT OF SIDE WALL C SURFACE 1 = CEILING C SURFACE 2 = SOUTH ROOF (FOR E-W RIDGE) C SURFACE 3 = NORTH ROOF C SURFACE 4 = EAST GABLE C SURFACE 5 = WEST GABLE C SURFACE 6 = SOUTH SIDE WALL C SURFACE 7 = NORTH SIDE WALL C EI = EMITTANCE OF SURFACE FACING ATTIC SPACE C F(I,J) = RADIATION VIEW FACTOR FROM SURFACE I TO C SURFACE J C G(I,J) = RADIATION FACTOR IMPLICIT REAL*8 (A-H, O-Z) DIMENSION F(7,7),CHI(7,7),PSI(7,7),E(7),B(7),EI(K),G(32,32) PI = 3.14159265 PIO2 = PI/2.0 P1 = PITCH1*PI/180. P2 = PITCH2*PI/180. F(1,1) = 0.0 F(2,2) = 0.0 F(3,3) = 0.0 F(4,4) = 0.0 F(5,5) = 0.0 F(6,6) = 0.0 F(7,7) = 0.0 P3 = (180.-PITCH1-PITCH2)*PI/180. IP1 = PITCH1 IP2 = PITCH2 IP3 = 180.0 - PITCH1 - PITCH2 IF(IP1.EQ.90.AND.IP2.EQ.90) STOP IF(IP1.EQ.0.AND.IP2.NE.0) STOP IF(IP1.NE.0.AND.IP2.EQ.0) STOP IF(IP1.EQ.0) GO TO 20 A2 = W*DSIN(P2)/DSIN(P3) A3 = W*DSIN(P1)/DSIN(P3) A2EXT = H1/DSIN(P1) A3EXT = H1/DSIN(P2) IF(IP1.NE.90) W2EXT = H1/DTAN(P1) IF(IP1.EQ.90) W2EXT = 0.0 IF(IP2.NE.90) W3EXT = H1/DTAN(P2) IF(IP2.EQ.90) W3EXT = 0.0 A2P = A2 + A2EXT A3P = A3 + A3EXT W2P = W + W2EXT W3P = W + W3EXT IF(IP1.NE.90) F(1,2) = (W2P*FMN(A2P,AL,W2P,P1) & +W2EXT*FMN(A2EXT,AL,W2EXT,P1) & -W2P*FMN(A2EXT,AL,W2P,P1)-W2EXT*FMN(A2P,AL,W2EXT,P1))/W IF(IP1.EQ.90) F(1,2) = FMN(A2P,AL,W,P1)-FMN(H1,AL,W,P1) IF(IP2.NE.90) F(1,3) = (W3P*FMN(A3P,AL,W3P,P2) & +W3EXT*FMN(A3EXT,AL,W3EXT,P2) & -W3P*FMN(A3EXT,AL,W3P,P2)-W3EXT*FMN(A3P,AL,W3EXT,P2))/W IF(IP2.EQ.90) F(1,3) = FMN(A3P,AL,W,P2)-FMN(H1,AL,W,P2) F(1,6) = FMN(H1,AL,W,PIO2) F(1,7) = FMN(H1,AL,W,PIO2) F(1,4) = (1.0 - F(1,2) - F(1,3) - F(1,6) - F(1,7))/2.0 F(1,5) = F(1,4) F(2,1) = W*AL*F(1,2)/A2/AL F(2,3) = FMN(A3,AL,A2,P3) ANG = P1 + PIO2 IF(IP1.NE.90) F(2,6) = FMN(H1,AL,A2,ANG) IF(IP1.EQ.90) F(2,6) = 0.0 ANG = PIO2 - P1 IF(IP1.EQ.90) GO TO 5 A2P = W/DSIN(ANG) A2EXT = A2P - A2 H1EXT = W/DTAN(ANG) H1P = H1 + H1EXT C NEXT LINE ADDED 4-6-88 IF(IP2.NE.90) THEN F(2,7) = (A2P*FMN(H1P,AL,A2P,ANG)+A2EXT*FMN(H1EXT,AL,A2EXT,ANG) & -A2P*FMN(H1EXT,AL,A2P,ANG)-A2EXT*FMN(H1P,AL,A2EXT,ANG))/A2 C NEXT 3 LINES ADDED 4-6-88 ELSE F(2,7) = FMN(H1P,AL,A2,ANG) - FMN(H1EXT,AL,A2,ANG) END IF GO TO 6 5 F(2,7) = ((A2+H1)*FP((A2+H1),AL,W)-H1*FP(H1,AL,W)-A2*FP(A2,AL,W))/ & 2.0/A2 6 F(2,4) = (1.0 - F(2,1) - F(2,3) - F(2,6) - F(2,7))/2.0 F(2,5) = F(2,4) F(3,1) = W*AL*F(1,3)/A3/AL F(3,2) = A2*AL*F(2,3)/A3/AL ANG = P2 + PIO2 IF(IP2.NE.90) F(3,7) = FMN(H1,AL,A3,ANG) IF(IP2.EQ.90) F(3,7) = 0.0 ANG = PIO2 - P2 IF(IP2.EQ.90) GO TO 7 A3P = W/DSIN(ANG) A3EXT = A3P - A3 H1EXT = W/DTAN(ANG) H1P = H1 + H1EXT C NEXT LINE ADDED 4-6-88 IF(IP1.NE.90) THEN F(3,6) = (A3P*FMN(H1P,AL,A3P,ANG)+A3EXT*FMN(H1EXT,AL,A3EXT,ANG) & -A3P*FMN(H1EXT,AL,A3P,ANG)-A3EXT*FMN(H1P,AL,A3EXT,ANG))/A3 C NEXT 3 LINES ADDED 4-6-88 ELSE F(3,6) = FMN(H1P,AL,A3,ANG) - FMN(H1EXT,AL,A3,ANG) END IF GO TO 8 7 F(3,6) = ((A3+H1)*FP((A3+H1),AL,W)-H1*FP(H1,AL,W)-A3*FP(A3,AL,W))/ & 2.0/A3 8 F(3,4) = (1.0D0 - F(3,1) - F(3,2) - F(3,6) - F(3,7))/2.0 F(3,5) = F(3,4) F(6,1) = W*AL*F(1,6)/H1/AL F(6,2) = A2*AL*F(2,6)/H1/AL F(6,3) = A3*AL*F(3,6)/H1/AL F(6,7) = FP(H1,AL,W) F(6,4) = (1.0D0 - F(6,1) - F(6,2) - F(6,3) - F(6,7))/2.0 F(6,5) = F(6,4) F(7,1) = W*AL*F(1,7)/H1/AL F(7,2) = A2*AL*F(2,7)/H1/AL F(7,3) = A3*AL*F(3,7)/H1/AL F(7,6) = F(6,7) F(7,4) = F(6,4) F(7,5) = F(7,4) F(4,1) = W*AL*F(1,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,2) = A2*AL*F(2,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,3) = A3*AL*F(3,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,6) = H1*AL*F(6,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,7) = H1*AL*F(7,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,5) = 1.0 - F(4,1) - F(4,2) - F(4,3) - F(4,6) - F(4,7) F(5,1) = F(4,1) F(5,2) = F(4,2) F(5,3) = F(4,3) F(5,4) = F(4,5) F(5,6) = F(4,6) F(5,7) = F(4,7) GO TO 50 20 F(1,2) = 0.5*FP(W,AL,H1) F(1,3) = F(1,2) F(1,4) = FMN(H1,W,AL,PIO2) F(1,5) = F(1,4) F(1,6) = FMN(H1,AL,W,PIO2) F(1,7) = F(1,6) F(2,1) = 2.0*F(1,2) F(2,3) = 0.0 X = W/2.0 F(2,6) = FMN(H1,AL,X,PIO2) X = W/2.0 F(2,7) = H1*(FMN(W,AL,H1,PIO2)-FMN(X,AL,H1,PIO2))/X F(2,4) = (1.0D0 - F(2,1) - F(2,3) - F(2,6) - F(2,7))/2.0 F(2,5) = F(2,4) F(3,1) = F(2,1) F(3,2) = 0.0D0 F(3,4) = F(2,4) F(3,5) = F(2,5) F(3,6) = F(2,7) F(3,7) = F(2,6) F(4,1) = AL*F(1,4)/H1 F(4,2) = (AL*W/2.0)*F(2,4)/(H1*W) F(4,3) = F(4,2) F(4,5) = FP(H1,W,AL) F(4,6) = FMN(AL,H1,W,PIO2) F(4,7) = F(4,6) F(5,1) = F(4,1) F(5,2) = F(4,2) F(5,3) = F(4,3) F(5,4) = F(4,5) F(5,6) = F(4,6) F(5,7) = F(4,7) F(6,1) = W*F(1,6)/H1 F(6,2) = (AL*W/2.0)*F(2,6)/(AL*H1) F(6,3) = (AL*W/2.0)*F(3,6)/(AL*H1) F(6,4) = W*F(4,6)/AL F(6,5) = F(6,4) F(6,7) = FP(H1,AL,W) F(7,1) = F(6,1) F(7,2) = F(6,3) F(7,3) = F(6,2) F(7,4) = F(6,4) F(7,5) = F(6,5) F(7,6) = F(6,7) C 50 CONTINUE DO 100 I = 1,7 DO 100 J = 1,7 100 CHI(I,J) = -(1. - EI(I))*F(I,J)/EI(I) DO 200 I = 1,7 200 CHI(I,I) = CHI(I,I) + 1./EI(I) C INVERT CHI MATRIX TO OBTAIN PSI MATRIX DO 500 L = 1,7 DO 300 I = 1,7 E(I) = 0.0D0 300 IF(I.EQ.L) E(I) = 1.0D0 CALL SOLVP(7,CHI,E,B,7) DO 400 I = 1,7 400 PSI(I,L) = B(I) 500 CONTINUE DO 600 I = 1,7 DO 600 J = 1,7 600 G(I,J) = PSI(I,J)*EI(I)/(1. - EI(I)) C C DO 1000 I = 1,7 C1000 WRITE(11,2000) (F(I,J),J=1,7) C WRITE(11,2000) C DO 1001 I = 1,7 C1001 WRITE(11,2000) (CHI(I,J),J=1,7) C WRITE(11,2000) C DO 1002 I = 1,7 C1002 WRITE(11,2000) (PSI(I,J),J=1,7) C WRITE(11,2000) C DO 1003 I = 1,7 C1003 WRITE(11,2000) (G(I,J),J=1,7) C WRITE(11,2000) C2000 FORMAT(1X,7F10.4) C Pause RETURN END C*********************************************************************** C** ** C** FUNCTION WDHUM ** C** ** C*********************************************************************** FUNCTION WDHUM(AMC,T) C THIS FUNCTION CALCULATES THE HUMIDITY RATIO OF A WOOD SURFACE C USING THE EQUATION FROM P. G. CLEARY (1985). C WDHUM = HUMIDITY RATIO, POUNDS OF WATER PER POUND OF DRY AIR C AMC = WOOD MOISTURE CONTENT, FRACTION OF DRY WEIGHT OF WOOD C T = WOOD SURFACE TEMPERATURE, DEGREES F IMPLICIT REAL*8 (A-H, O-Z) A = 28.6 B = -0.00049 C = 0.0172 D = -0.060 E = 0.076 WDHUM = B + C*AMC + D*AMC*AMC + E*AMC**3 C NEXT LINE ADDED 11-2-88 IF(WDHUM.LT.0.0) WDHUM = 0.0 WDHUM = WDHUM * EXP(T/A) RETURN END C*********************************************************************** C** ** C** SUBROUTINE SUN1 ** C** ** C*********************************************************************** SUBROUTINE SUN1(STALAT,IDOY) C KW SUBROUTINE SUN1 C CALCULATES DAILY DATA ON SOLAR RADIATION CKW SUBROUTINE BORROWED FROM DOE 2.1C WITH MODIFICATIONS C CKW COMMON /LOCALD/STALAT,STALON,ITIMZ,BAZIM,BALTIT, CKW 1 SSTALA,CSTALA,TSTALA,SBAZIM,CBAZIM, CKW 2 BXORG,BYORG,SHCOEF,TP1,TP2,WSTP1,WSTP2, CKW 3 WSHGT CKW DIMENSION LOCALD(18) CKW EQUIVALENCE (LOCALD(1), STALAT) IMPLICIT REAL*8 (A-H, O-Z) COMMON /SUND/GUNDOG,HORANG,TDECLN,EQTIME,SOLCON, 1 ATMEXT,SKYDFF,RAYCOS(3),RDN, 2 BSUN,DECLN,CD,SD,FSUNUP,ISUNUP CKW COMMON /TIME/IDOY,IDOW,IDSTF,ISCHR,ISCDAY, CKW 1 IMO,IDAY,IYR,IHR,CLOCK(10),IHLFLG,IDSFLG, CKW 2 MONDSC(12),MONLEN(12),MONSDA(12),IEODMR C GET SIN, COS OF DAY OF YEAR/365 C1 = DCOS(0.01721*DFLOTJ(IDOY)) S1 = DSIN(0.01721*DFLOTJ(IDOY)) S2 = 2.*S1*C1 C2 = C1*C1 - S1*S1 C3 = C1*C2 - S1*S2 S3 = C1*S2 + S1*C2 C CALC TANGENT OF DECLINATION ANGLE TDECLN = 0.00527 - 0.4001*C1 - 0.003996*C2 - 0.004240*C3 & + 0.0672*S1 C CALC EQUATION OF TIME EQTIME = 0.696D-4 + 0.706D-2*C1 - 0.0533*C2 - 0.157D-2*C3-0.122*S1 & -0.156*S2 - 0.556D-2*S3 C CALCULATE SOLAR CONSTANT SOLCON = 368.44 + 24.52*C1 - 1.14*C2 - 1.09*C3 + 0.58*S1 - 0.18*S2 & +0.28*S3 C IF SOUTHERN HEMISPHERE, SHIFT AIRMASS AND ATMOSPHERIC C EXTINCTION CURVES BY 6 MONTHS IF(STALAT.GE.0.) GO TO 10 C1 = -C1 S1 = -S1 C3 = -C3 S3 = -S3 10 CONTINUE C CALCULATE ATMOSPHERIC EXTINCTION COEFFICIENT ATMEXT = 0.1717 - 0.0344*C1 + 0.0032*C2 + 0.0024*C3 - 0.0043*S1 & -0.0008*S3 C CALCULATE SKY DIFFUSIVITY SKYDFF = 0.0905 - 0.0410*C1 + 0.0073*C2 + 0.0015*C3 - 0.0034*S1 & +0.0004*S2 - 0.0006*S3 C CALCULATE HOUR ANGLE OF SUNRISE C NEXT LINE ADDED BY K. WILKES TSTALA = DTAN(STALAT) GUNDOG = DACOS(-TSTALA*TDECLN) DECLN = DATAN(TDECLN) CD = DCOS(DECLN) SD = DSIN(DECLN) RETURN END C*********************************************************************** C** ** C** SUBROUTINE WDTSUN ** C** ** C*********************************************************************** SUBROUTINE WDTSUN(IHR,ITIMZ,STALON,SSTALA,CSTALA,SBAZIM, & CBAZIM,CLRNES,SOLRAD,DIRSOL,ISUN) CKW SUBROUTINE WDTSUN CKW CALCULATE WEATHER DATA AND SUN HOURLY DATA CKW SUBROUTINE BORROWED FROM DOE 2.1C, WITH MODIFICATIONS CKW IHR = HOUR OF DAY CKW ITIMZ = TIME ZONE INDICATOR, = 5 FOR EASTERN, = 6 FOR CENTRAL,ETC. CKW STALON = LONGITUDE CKW STALAT = LATITUDE CKW SSTALA = SIN OF LATITUDE CKW CSTALA = COSINE OF LATITUDE CKW SBAZIM = SIN OF BUILDING AZIMUTH CKW CBAZIM = COSINE OF BUILDING AZIMUTH CKW CLRNES = ATMOSPHERIC CLEARNESS NUMBER CKW QSOLH = SOLAR RADIATION ON HORIZONTAL SURFACE CKW CLDAMT = CLOUD AMOUNT CKW ICLDTY = CLOUD TYPE CKW ISUN = FLAG FOR TYPE OF SOLAR DATA AVAILABLE CKW ISUN = 1, MEASURED TOTAL HORIZONTAL AND DIRECT AVAILABLE CKW ISUN = 2, MEASURED TOTAL HORIZONTAL ONLY AVAILABLE CKW ISUN = 3, NO MEASURED SOLAR, BUT MEASURED CLOUD AMOUNT CKW RDNCC = DIRECT NORMAL SOLAR, WITH CLOUD COVER CKW BSCC = DIFFUSE WITH CLOUD COVER CKW COMMON /CLCFLG/IFSTHR,IPRDFL,IGOLGE,NDD,IDDFLG,LDSTYP CKW COMMON /INFPAR/PTWV,DUMG(8),PSE CKW COMMON /LOCALD/STALAT,STALON,ITIMZ,BAZIM,BALTIT, CKW 1 SSTALA,CSTALA,TSTALA,SBAZIM,CBAZIM, CKW 2 BXORG,BYORG,SHCOEF,TP1,TP2,WSTP1,WSTP2, CKW 3 WSHGT CKW DIMENSION LOCALD(18) CKW EQUIVALENCE (LOCALD(1), STALAT) IMPLICIT REAL*8 (A-H, O-Z) COMMON /SUND/GUNDOG,HORANG,TDECLN,EQTIME,SOLCON, 1 ATMEXT,SKYDFF,RAYCOS(3),RDN, 2 BSUN,DECLN,CD,SD,FSUNUP,ISUNUP CKW COMMON /TIME/IDOY,IDOW,IDSTF,ISCHR,ISCDAY, CKW 1 IMO,IDAY,IYR,IHR,CLOCK(10),IHLFLG,IDSFLG, CKW 2 MONDSC(12),MONLEN(12),MONSDA(12),IEODMR CKW COMMON /WEATH/IWDID(5),LRECX,WLAT,WLONG,LTIMZ,IFX,IWSIZ, CKW 1 CLRNES,TGNDR,WBT,DBT,PATM,CLDAMT,ISNOW,IRAIN, CKW 2 IWNDDR,HUMRAT,DENSTY,ENTHAL,DIFSOL,DIRSOL,SOLRAD,ICLDTY, CKW 3 WNDSPD,IDUMMY,DPT,WNDDRR,CLDCOV,RDNCC,BSCC,SKYA, CKW 4 DBTR,GTEMP(12),CLR(12) CKW NEXT TWO COMMONS ADDED BY K. WILKES COMMON /SUN2/RDNCC,BSCC COMMON /CLOUD/CLDAMT,CC,ICLDTY CKW IF(NDD.GT.0) GO TO 50 C IF NOT DESIGN DAY, GET WEATHER DATA FROM THE HOURLY WEATHER FILE CKW CALL WEATHI CKW GO TO 60 C OTHERWISE, GET DESIGN DAY WEATHER CKW50 CALL DESWTH CKW60 CONTINUE C CONVERT WIND DIRECTION FROM 16THS OF CIRCLE TO RADIANS CKW WNDDRR = FLOAT(IWNDDR)*0.39269908 C CALCULATE PRESSURE DUE TO WIND VELOCITY CKW PTWV = 0.000638*WNDSPD*WNDSPD C SET RADIATION LOSS TO SKY FOR TILT = 0 CKW SKYA = (10. - CLDAMT)*2. C CALCULATE HOUR ANGLE HORANG = 0.2618*(DFLOTJ(IHR-12+ITIMZ) + EQTIME - 0.5) - STALON C SET TEST TO BE THE HOUR ANGLE OF THE BIN EDGE NEAREST NOON TEST = HORANG + 0.1309 IF(IHR.GE.12) TEST = HORANG - 0.1309 C IF SUN IS DOWN, SKIP IF(ABS(TEST).GT.ABS(GUNDOG)) GO TO 1800 C SUN IS UP. SET FLAG ISUNUP = 1 C TEST TO SEE IF THIS HOUR BIN CONTAINS SUNRISE OR SUNSET FSUNUP = 1 DIFF = ABS(GUNDOG) - ABS(TEST) IF((DIFF.LT.0.) .OR. (DIFF.GE.0.2618)) GO TO 200 C RESET THE HOUR ANGLE HALF WAY BETWEEN SUNRISE OR C SUNSET AND THE BIN EDGE NEAREST NOON IF(IHR.LT.12) HORANG = HORANG + 0.5*(0.2618 - DIFF) IF(IHR.GE.12) HORANG = HORANG - 0.5*(0.2618 - DIFF) C SET FSUNUP TO BE THE FRACTION OF THE HOUR THE SUN WAS UP FSUNUP = 3.8197*DIFF 200 CONTINUE CHCD = DCOS(HORANG)*CD SHCD = DSIN(HORANG)*CD CLSD = SD*CSTALA CHCDSC = CHCD*SSTALA - CLSD C CALCULATE SOLAR DIRECTION COSINES RAYCOS(1) = CHCDSC*SBAZIM - SHCD*CBAZIM RAYCOS(2) = -CHCDSC*CBAZIM - SHCD*SBAZIM RAYCOS(3) = SSTALA*SD + CSTALA*CHCD IF(RAYCOS(3).LT.0.001) GO TO 1800 CKW IF((IFX.GE.5) .AND. (NDD.LE.0)) GO TO 2000 CKW NEXT LINE ADDED BY K. WILKES IF(ISUN.EQ.1) GO TO 2000 C CALCULATE DIRECT NORMAL AND DIFFUSE RADIATION FOR SKY WITH C NO CLOUDS C IBM MACHINE PROBLEM WITH EXPONENT UNDERFLOW FOR EXP( ). C PROBABLE CAUSE < RAYCOS(3) VERY SMALL --> ARGUMENT VERY LARGE. C FIX < CHECK ARGUMENT FOR LARGE NEGATIVE NUMBER. RDN = 0 IF SO. ARGUE = -ATMEXT/RAYCOS(3) IF(ARGUE.LT.-50.0) GO TO 1 RDN = SOLCON*CLRNES*EXP(ARGUE) GO TO 2 C1 WRITE(6,3) C2 FORMAT('>>>>NOTE 0 MEANS NO DIRECT SOLAR 240 IF(ETA.GT.0.) GO TO 260 RDIR = 0. GO TO 270 C CALCULATE DIRECT RADIATION ON WALL 260 RDIR = RDNCC*ETA 270 FFG = (1.0 - GAMMA)*0.5 FFS = 1.0 - FFG RDIF = FFS*BSCC + FFG*BG SUN3OR = RDIF + RDIR RETURN END C*********************************************************************** C** ** C** SUBROUTINE CCMINV ** C** ** C*********************************************************************** SUBROUTINE CCMINV(CLDCOV,AL,CLDAMT) C THIS SUBROUTINE CALCULATES THE CLOUD AMOUNT FROM THE CLOUD COVER, C ASSUMING TYPE 2 CLOUDS IMPLICIT REAL*8 (A-H, O-Z) IF(AL.GT.0.707) GO TO 100 C LOW SUN IF(CLDCOV.GE.1.0) THEN CLDAMT = 0.0 ELSE IF(CLDCOV.GE.0.71919) THEN CLDAMT = (1.0 - CLDCOV)/(1.0 - 0.71919) ELSE IF(CLDCOV.LE.0.3825) THEN CLDAMT = 10.0 ELSE Y0 = CLDCOV K = 1 XK = 9.995 10 XKP1 = XK - (0.724-0.00625*XK+0.00191*XK*XK-0.00047*XK*XK*XK & -Y0)/(-0.00625+2.*0.00191*XK-3.*0.00047*XK*XK) IF(ABS(XKP1-XK).GE.0.01) THEN XK = XKP1 K = K + 1 IF(K.GT.25) STOP 'ERROR IN CCMINV' GO TO 10 ELSE CLDAMT = XKP1 END IF END IF GO TO 200 C HIGH SUN 100 IF(CLDCOV.GE.1.0) THEN CLDAMT = 0.0 ELSE IF(CLDCOV.GE.0.94292) THEN CLDAMT = (1.0 - CLDCOV)/(1.0 - 0.94292) ELSE IF(CLDCOV.LE.0.6056) THEN CLDAMT = 10.0 ELSE Y0 = CLDCOV K = 1 XK = 9.995 20 XKP1 = XK - (0.959-0.02304*XK+0.00787*XK*XK-0.00091*XK*XK*XK & -Y0)/(-0.02304+2.*0.00787*XK-3.*0.00091*XK*XK) IF(ABS(XKP1-XK).GE.0.01) THEN XK = XKP1 K = K + 1 IF(K.GT.25) STOP 'ERROR IN CCMINV' GO TO 20 ELSE CLDAMT = XKP1 END IF END IF 200 CONTINUE RETURN END C*********************************************************************** C** ** C** SUBROUTINE SKY ** C** ** C*********************************************************************** SUBROUTINE SKY(TO,TDP,IHR,QIR,CLDAMT,TS) IMPLICIT REAL*8 (A-H, O-Z) c COMMON /HUMID/WO,PATM IF (QIR .le. 0.0) go to 10 TS = (QIR/1.714D-09)**0.25 - 459.67 go to 20 10 continue c XX = WO/0.62198 c P = PATM*XX/(1. + XX) c ALF = DLOG(P*29.921/14.696) c IF(P.GT.0.08865) TDP = 79.047 + 30.5790*ALF + 1.8893*ALF*ALF c IF(P.LE.0.08865) TDP = 71.98 + 24.873*ALF + 0.8927*ALF*ALF TD = (TDP+459.67)/1.8 - 273.15 EPS = 0.711 + 0.56*TD/100. + 0.73*(TD/100.)*(TD/100.) + 0.013* & DCOS(2.0*3.14159265/24.*DFLOTJ(IHR)) EPS = EPS + (1.0 - EPS)*CLDAMT/10.*0.784 TS = (TO+459.67)*DSQRT(DSQRT(EPS)) - 459.67 20 RETURN END C*********************************************************************** C** ** C** SUBROUTINE DUCTREAD ** C** ** C*********************************************************************** C This subroutine reads input for the duct model, and also calculates C inside and outside perimeters and overall thermal resistance, C (actually, resistance divided by perimeter). C Next line added for model with ducts SUBROUTINE DUCTREAD C Nomenclature: C NSUPLY = number of supply ducts C NRETRN = number of return ducts C NTOT = total number of ducts C xTDUCTIN(I,J) = temperature of air entering supply duct J=1 Heating Mode C TI = temperature of air inside house J=2 Cooling Mode C = temperature of air entering return duct C ND(I) = number of duct I, only used for convenience in input C NINLET(I) = node number at inlet of duct I C NEXIT(I) = node number at exit of duct I C ISHAPE(I) = flag for shape of duct: 0 for round, 1 for rectangular C AKI(I) = thermal conductivity of duct wall, inside layer of duct I C AKO(I) = thermal conductivity of duct wall, outside layer of duct I C ALD(I) = length of duct I C MDOTIN(I) = mass flow rate entering duct I C MDOTLK(I) = mass leakage per unit length for duct I C EODUCT(I) = emittance of outside of duct I C ZONE(I) = zone of house that duct I supplies C DI(I) = inside diameter of round duct C DW(I) = inside diameter of outer layer of round duct wall C DO(I) = outside diameter of round duct C XI(I), Y(I) = inside width and height of rectangular duct C XW(I), YW(I) = inside width and height of outer layer of C rectangular duct C XO(I), YO(I) = outside width and height of rectangular duct C PID(I) = inside perimeter of duct C PW(I) = inside perimeter of outer layer of duct C PO(I) = outside perimeter of duct C RD(I) = thermal resistance of duct wall (usual resistance divided C by perimeter C DH(I) = hydraulic diameter of duct C DCHAR(I) = characteristic dimension of duct for convection: C outside diameter for round duct, C height for rectangular duct C RHOC(I) = heat capacity per unit length of duct, Btu/F-ft C DX(I) = length of duct section for transient model C NSEC(I) = number of duct sections for duct I IMPLICIT REAL*8 (A-H, O-Z) COMMON /DUCT1/ALD(25),xMDOTin(25,2),xMDOTlk(25,2),EODUCT(25), & ZONE(25),PID(25),PO(25),RD(25),xTDUCTin(2),RHOC(25),RDO(25), & RDI(25),NSUPLY,NRETRN,ND(25),NINLET(25),NEXIT(25) COMMON /DUCT2/DH(25),DCHAR(25),AD(25),ISHAPE(25) COMMON /DUCT3/TI,V,ONTIME,DuctTin,DuctMin(25), & DuctLkg(25),KFLAG(19),NTOT COMMON /DUCT5/DX(25),Rduct,NSEC(25) DIMENSION AKI(25),AKO(25),DI(25),DW(25),DO(25),PW(25) DIMENSION XI(25),YI(25),XW(25),YW(25),XO(25),YO(25) PI = 3.14159265 READ(5,*) NSUPLY,NRETRN NTOT = NSUPLY + NRETRN READ(5,*) (xTDUCTin(i), i=1,2) DO 10 I = 1,NTOT READ(5,*) ND(I),NINLET(I),NEXIT(I),ISHAPE(I),AKI(I),AKO(I), & RHOC(I),ALD(I),(xMDOTin(I,J),xMDOTlk(I,J), J=1,2), &EODUCT(I),ZONE(I) IF(ISHAPE(I).EQ.0) THEN READ(5,*) DI(I),DW(I),DO(I) ELSE READ(5,*) XI(I),YI(I),XW(I),YW(I),XO(I),YO(I) END IF 10 CONTINUE C Calculate perimeters DO 20 I = 1,NTOT IF(ISHAPE(I).EQ.0) THEN PID(I) = PI*DI(I) PW(I) = PI*DW(I) PO(I) = PI*DO(I) Rduct = 0.5*DO(1)*(DLOG(DW(1)/DI(1))/AKI(1) & + DLOG(DO(1)/DI(1))/AKO(1)) ! Duct R-Value used for Excel Identifier ELSE PID(I) = 2.*(XI(I) + YI(I)) PO(I) = 2.*(XO(I) + YO(I)) END IF 20 CONTINUE C Calculate wall resistance DO 30 I = 1,NTOT IF(ISHAPE(I).EQ.0) THEN RDO(I) = DLOG(PO(I)/PW(I))/2./PI/AKO(I) RDI(I) = DLOG(PW(I)/PID(I))/2./PI/AKI(I) RD(I) = RDO(I) + RDI(I) ELSE RHOR1 = (XO(I) - XW(I))/(YO(I)+YW(I))*2./AKO(I) RHOR2 = (XW(I) - XI(I))/(YW(I)+YI(I))*2./AKI(I) RVER1 = (YO(I) - YW(I))/(XO(I)+XW(I))*2./AKO(I) RVER2 = (YW(I) - YI(I))/(XW(I)+XI(I))*2./AKI(I) RDO(I) = 1./(2./(RHOR1) + 2./(RVER1)) RDI(I) = 1./(2./(RHOR2) + 2./(RVER2)) RD(I) = 1./(2./(RHOR1+RHOR2) + 2./(RVER1 + RVER2)) END IF 30 CONTINUE C Set characteristic dimensions for convection and hydraulic C diameter on inside DO 40 I = 1,NTOT IF(ISHAPE(I).EQ.0) THEN DH(I) = DI(I) DCHAR(I) = DO(I) ELSE DH(I) = 4.*XI(I)*YI(I)/(2.*(XI(I)+YI(I))) DCHAR(I) = YO(I) END IF 40 CONTINUE C Calculate flow areas DO 50 I = 1,NTOT IF(ISHAPE(I).EQ.0) THEN AD(I) = PI/4.*DI(I)*DI(I) ELSE AD(I) = XI(I)*YI(I) END IF 50 CONTINUE C Calculate number of sections and section length for transient model DO 60 I = 1,NTOT C 12-5-96 Increase NSEC by 1 to allow ALD less than 1 ft NSEC(I) = IDINT(ALD(I))+1 ANSEC = NSEC(I) DX(I) = ALD(I)/ANSEC 60 CONTINUE RETURN END C*********************************************************************** C** ** C** SUBROUTINE VIEW2D ** C** ** C*********************************************************************** SUBROUTINE VIEW2D(AL,W,PITCH1,PITCH2,H1,EI,G,K,A) C THIS SUBROUTINE CALCULATES THE VIEW FACTORS AMONG THE SURFACES C IN A GABLED ATTIC, HAVING A FIVE-SIDED CROSS-SECTION, AND C HAVING ARBITRARILY PITCHED ROOF SURFACES C VIEW2D INCLUDES ATTIC DUCTS C ASSUMES THAT THE DUCTS DO NOT SEE EACH OTHER C ASSUMES THAT THE VIEW FACTOR FROM A DUCT TO THE ATTIC FLOOR IS 0.5 C ASSUMES THAT THE VIEW FACTORS AMONG THE ATTIC SURFACES OTHER THAN C THE ATTIC FLOOR ARE NOT INFLUENCED BY THE PRESENCE OF THE DUCTS C ASSUMES THAT THE VIEW FACTORS FROM THE ATTIC FLOOR TO THE OTHER C ATTIC SURFACES IN THE PRESENCE OF DUCTS ARE PROPORTIONAL TO THOSE C IN THE ABSENCE OF DUCTS C AL = LENGTH OF ATTIC C W = WIDTH OF ATTIC C PITCH1 = PITCH OF SOUTH ROOF, DEGREES C PITCH2 = PITCH OF NORTH ROOF, DEGREES C H1 = HEIGHT OF SIDE WALL C SURFACE 1 = CEILING C SURFACE 2 = SOUTH ROOF (FOR E-W RIDGE) C SURFACE 3 = NORTH ROOF C SURFACE 4 = EAST GABLE C SURFACE 5 = WEST GABLE C SURFACE 6 = SOUTH SIDE WALL C SURFACE 7 = NORTH SIDE WALL C EI = EMITTANCE OF SURFACE FACING ATTIC SPACE C A(J) = AREAS OF ATTIC SURFACES C F(I,J) = RADIATION VIEW FACTOR FROM SURFACE I TO C SURFACE J C G(I,J) = RADIATION FACTOR C NEXT LINE DELETED AND FOLLOWING LINE ADDED FOR MODEL WITH DUCTS C DIMENSION F(7,7),CHI(7,7),PSI(7,7),E(7),B(7),EI(K),G(7,7) IMPLICIT REAL*8 (A-H, O-Z) DIMENSION F(32,32),CHI(32,32),PSI(32,32),E(32) DIMENSION B(32),EI(7),G(32,32),A(K),EMIT(32) REAL MDOTIN,MDOTLK COMMON /DUCT1/ALD(25),xMDOTin(25,2),xMDOTlk(25,2),EODUCT(25), & ZONE(25),PID(25),PO(25),RD(25),xTDUCTin(2),RHOC(25),RDO(25), & RDI(25),NSUPLY,NRETRN,ND(25),NINLET(25),NEXIT(25) C PI = 3.14159265 PIO2 = PI/2.0 P1 = PITCH1*PI/180. P2 = PITCH2*PI/180. F(1,1) = 0.0 F(2,2) = 0.0 F(3,3) = 0.0 F(4,4) = 0.0 F(5,5) = 0.0 F(6,6) = 0.0 F(7,7) = 0.0 P3 = (180.-PITCH1-PITCH2)*PI/180. IP1 = PITCH1 IP2 = PITCH2 IP3 = 180.0 - PITCH1 - PITCH2 IF(IP1.EQ.90.AND.IP2.EQ.90) STOP IF(IP1.EQ.0.AND.IP2.NE.0) STOP IF(IP1.NE.0.AND.IP2.EQ.0) STOP IF(IP1.EQ.0) GO TO 20 A2 = W*DSIN(P2)/DSIN(P3) A3 = W*DSIN(P1)/DSIN(P3) A2EXT = H1/DSIN(P1) A3EXT = H1/DSIN(P2) IF(IP1.NE.90) W2EXT = H1/DTAN(P1) IF(IP1.EQ.90) W2EXT = 0.0 IF(IP2.NE.90) W3EXT = H1/DTAN(P2) IF(IP2.EQ.90) W3EXT = 0.0 A2P = A2 + A2EXT A3P = A3 + A3EXT W2P = W + W2EXT W3P = W + W3EXT IF(IP1.NE.90) F(1,2) = (W2P*FMN(A2P,AL,W2P,P1) & +W2EXT*FMN(A2EXT,AL,W2EXT,P1) & -W2P*FMN(A2EXT,AL,W2P,P1)-W2EXT*FMN(A2P,AL,W2EXT,P1))/W IF(IP1.EQ.90) F(1,2) = FMN(A2P,AL,W,P1)-FMN(H1,AL,W,P1) IF(IP2.NE.90) F(1,3) = (W3P*FMN(A3P,AL,W3P,P2) & +W3EXT*FMN(A3EXT,AL,W3EXT,P2) & -W3P*FMN(A3EXT,AL,W3P,P2)-W3EXT*FMN(A3P,AL,W3EXT,P2))/W IF(IP2.EQ.90) F(1,3) = FMN(A3P,AL,W,P2)-FMN(H1,AL,W,P2) F(1,6) = FMN(H1,AL,W,PIO2) F(1,7) = FMN(H1,AL,W,PIO2) F(1,4) = (1.0 - F(1,2) - F(1,3) - F(1,6) - F(1,7))/2.0 F(1,5) = F(1,4) F(2,1) = W*AL*F(1,2)/A2/AL F(2,3) = FMN(A3,AL,A2,P3) ANG = P1 + PIO2 IF(IP1.NE.90) F(2,6) = FMN(H1,AL,A2,ANG) IF(IP1.EQ.90) F(2,6) = 0.0 ANG = PIO2 - P1 IF(IP1.EQ.90) GO TO 5 A2P = W/DSIN(ANG) A2EXT = A2P - A2 H1EXT = W/DTAN(ANG) H1P = H1 + H1EXT C NEXT LINE ADDED 4-6-88 IF(IP2.NE.90) THEN F(2,7) = (A2P*FMN(H1P,AL,A2P,ANG)+A2EXT*FMN(H1EXT,AL,A2EXT,ANG) & -A2P*FMN(H1EXT,AL,A2P,ANG)-A2EXT*FMN(H1P,AL,A2EXT,ANG))/A2 C NEXT 3 LINES ADDED 4-6-88 ELSE F(2,7) = FMN(H1P,AL,A2,ANG) - FMN(H1EXT,AL,A2,ANG) END IF GO TO 6 5 F(2,7) = ((A2+H1)*FP((A2+H1),AL,W)-H1*FP(H1,AL,W)-A2*FP(A2,AL,W))/ & 2.0/A2 6 F(2,4) = (1.0 - F(2,1) - F(2,3) - F(2,6) - F(2,7))/2.0 F(2,5) = F(2,4) F(3,1) = W*AL*F(1,3)/A3/AL F(3,2) = A2*AL*F(2,3)/A3/AL ANG = P2 + PIO2 IF(IP2.NE.90) F(3,7) = FMN(H1,AL,A3,ANG) IF(IP2.EQ.90) F(3,7) = 0.0 ANG = PIO2 - P2 IF(IP2.EQ.90) GO TO 7 A3P = W/DSIN(ANG) A3EXT = A3P - A3 H1EXT = W/DTAN(ANG) H1P = H1 + H1EXT C NEXT LINE ADDED 4-6-88 IF(IP1.NE.90) THEN F(3,6) = (A3P*FMN(H1P,AL,A3P,ANG)+A3EXT*FMN(H1EXT,AL,A3EXT,ANG) & -A3P*FMN(H1EXT,AL,A3P,ANG)-A3EXT*FMN(H1P,AL,A3EXT,ANG))/A3 C NEXT 3 LINES ADDED 4-6-88 ELSE F(3,6) = FMN(H1P,AL,A3,ANG) - FMN(H1EXT,AL,A3,ANG) END IF GO TO 8 7 F(3,6) = ((A3+H1)*FP((A3+H1),AL,W)-H1*FP(H1,AL,W)-A3*FP(A3,AL,W))/ & 2.0/A3 8 F(3,4) = (1.0 - F(3,1) - F(3,2) - F(3,6) - F(3,7))/2.0 F(3,5) = F(3,4) F(6,1) = W*AL*F(1,6)/H1/AL F(6,2) = A2*AL*F(2,6)/H1/AL F(6,3) = A3*AL*F(3,6)/H1/AL F(6,7) = FP(H1,AL,W) F(6,4) = ( 1.0 - F(6,1) - F(6,2) - F(6,3) - F(6,7))/2.0 F(6,5) = F(6,4) F(7,1) = W*AL*F(1,7)/H1/AL F(7,2) = A2*AL*F(2,7)/H1/AL F(7,3) = A3*AL*F(3,7)/H1/AL F(7,6) = F(6,7) F(7,4) = F(6,4) F(7,5) = F(7,4) F(4,1) = W*AL*F(1,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,2) = A2*AL*F(2,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,3) = A3*AL*F(3,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,6) = H1*AL*F(6,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,7) = H1*AL*F(7,4)/(0.5*A2*DSIN(P1)*W + W*H1) F(4,5) = 1.0 - F(4,1) - F(4,2) - F(4,3) - F(4,6) - F(4,7) F(5,1) = F(4,1) F(5,2) = F(4,2) F(5,3) = F(4,3) F(5,4) = F(4,5) F(5,6) = F(4,6) F(5,7) = F(4,7) GO TO 50 20 F(1,2) = 0.5*FP(W,AL,H1) F(1,3) = F(1,2) F(1,4) = FMN(H1,W,AL,PIO2) F(1,5) = F(1,4) F(1,6) = FMN(H1,AL,W,PIO2) F(1,7) = F(1,6) F(2,1) = 2.0*F(1,2) F(2,3) = 0.0 X = W/2.0 F(2,6) = FMN(H1,AL,X,PIO2) X = W/2.0 F(2,7) = H1*(FMN(W,AL,H1,PIO2)-FMN(X,AL,H1,PIO2))/X F(2,4) = (1.0 - F(2,1) - F(2,3) - F(2,6) - F(2,7))/2.0 F(2,5) = F(2,4) F(3,1) = F(2,1) F(3,2) = 0.0 F(3,4) = F(2,4) F(3,5) = F(2,5) F(3,6) = F(2,7) F(3,7) = F(2,6) F(4,1) = AL*F(1,4)/H1 F(4,2) = (AL*W/2.0)*F(2,4)/(H1*W) F(4,3) = F(4,2) F(4,5) = FP(H1,W,AL) F(4,6) = FMN(AL,H1,W,PIO2) F(4,7) = F(4,6) F(5,1) = F(4,1) F(5,2) = F(4,2) F(5,3) = F(4,3) F(5,4) = F(4,5) F(5,6) = F(4,6) F(5,7) = F(4,7) F(6,1) = W*F(1,6)/H1 F(6,2) = (AL*W/2.0)*F(2,6)/(AL*H1) F(6,3) = (AL*W/2.0)*F(3,6)/(AL*H1) F(6,4) = W*F(4,6)/AL F(6,5) = F(6,4) F(6,7) = FP(H1,AL,W) F(7,1) = F(6,1) F(7,2) = F(6,3) F(7,3) = F(6,2) F(7,4) = F(6,4) F(7,5) = F(6,5) F(7,6) = F(6,7) C 50 CONTINUE C NEXT SECTION ADDED FOR MODEL WITH DUCTS NTOT = NSUPLY + NRETRN NTOT1 = 7 + NTOT DO 51 I = 1,7 51 EMIT(I) = EI(I) DO 52 I = 1,NTOT 52 EMIT(I+7) = EODUCT(I) DO 60 I = 1,NTOT F(I+7,1) = 0.5 60 F(1,I+7) = PO(I)*ALD(I)/A(1)*F(I+7,1) SUM = 0.0 DO 70 I = 1,NTOT 70 SUM = SUM + F(1,I+7) DO 80 J = 2,7 F(1,J) = (1. - SUM)*F(1,J) 80 F(J,1) = (1. - SUM)*F(J,1) SUMA = 0.0 DO 85 I = 1,NTOT 85 SUMA = SUMA + PO(I)*ALD(I) DO 95 I = 1,NTOT DO 95 J = 2,7 F(I+7,J) = A(J)/SUMA*F(J,1)*(1./(1.-SUM) - 1.) 95 F(J,I+7) = F(I+7,J)*PO(I)*ALD(I)/A(J) DO 97 I = 1,NTOT DO 97 J = 1,NTOT 97 F(I+7,J+7) = 0.0 C ABOVE SECTION ADDED FOR MODEL WITH DUCTS DO 100 I = 1,NTOT1 DO 100 J = 1,NTOT1 100 CHI(I,J) = -(1. - EMIT(I))*F(I,J)/EMIT(I) DO 200 I = 1,NTOT1 200 CHI(I,I) = CHI(I,I) + 1./EMIT(I) C INVERT CHI MATRIX TO OBTAIN PSI MATRIX DO 500 L = 1,NTOT1 DO 300 I = 1,NTOT1 E(I) = 0.0 300 IF(I.EQ.L) E(I) = 1.0 CALL SOLVP(NTOT1,CHI,E,B,32) DO 400 I = 1,NTOT1 400 PSI(I,L) = B(I) 500 CONTINUE DO 600 I = 1,NTOT1 DO 600 J = 1,NTOT1 600 G(I,J) = PSI(I,J)*EMIT(I)/(1. - EMIT(I)) c DO 1000 I = 1,NTOT1 c1000 WRITE(11,2000) (F(I,J),J=1,NTOT1) C Pause c WRITE(11,2000) c DO 1001 I = 1,NTOT1 c1001 WRITE(11,2000) (CHI(I,J),J=1,NTOT1) c WRITE(11,2000) c DO 1002 I = 1,NTOT1 c1002 WRITE(11,2000) (PSI(I,J),J=1,NTOT1) c WRITE(11,2000) c DO 1003 I = 1,NTOT1 c1003 WRITE(11,2000) (G(I,J),J=1,NTOT1) c WRITE(11,2000) c2000 FORMAT(1X,10F10.4) RETURN END C*********************************************************************** C** ** C** SUBROUTINE DUCTSS ** C** ** C*********************************************************************** C This subroutine contains algorithms for steady-state duct model, C with or without leakage in any or all of the duct segments. SUBROUTINE DUCTSS(G,TIS,TA,T) C Duct and node numbering scheme: C Supply ducts are numbered 1 through NSUPLY C Return ducts are numbered NSUPLY+1 through NSUPLY+NRETRN C The inlet from the supply source is numbered node 1 C Outlets of supply nodes are numbered 2 through NSUPLY+1 C The inlet to the return duct(s) is node NSUPLY+2 C The outlet(s) from return duct(s) are nodes NSUPLY+2 C through NSUPLY+NRETRN+2 C Nomenclature: C NSUPLY = number of supply ducts C NRETRN = number of return ducts C NTOT = total number of ducts C xTDUCTIN(i,j) = temperature of air entering supply duct J=1 Heating Mode C TI = temperature of air inside house J=2 Cooling Mode C = temperature entering return duct C ND(I) = number of duct I, only used for convenience in input C NINLET(I) = node number at inlet of duct I C NEXIT(I) = node number at exit of duct I C ISHAPE(I) = flag for shape of duct: 0 for round, 1 for rectangular C ALD(I) = length of duct I C MDOTIN(I,J) = mass flow rate entering duct I, J=1 Heating Mode C MDOTLK(I,J) = mass leakage per unit length for duct I J=2 Cooling Mode C EODUCT(I) = emittance of outside of duct I C ZONE(I) = zone of house that duct I supplies C DI(I) = inside diameter of round duct C DW(I) = inside diameter of outer layer of round duct wall C DO(I) = outside diameter of round duct C XI(I), Y(I) = inside width and height of rectangular duct C XW(I), YW(I) = inside width and height of outer layer of C rectangular duct C XO(I), YO(I) = outside width and height of rectangular duct C PID(I) = inside perimeter of duct C PW(I) = inside perimeter of outer layer of duct C PO(I) = outside perimeter of duct C RD(I) = thermal resistance of duct wall (usual resistance divided C by perimeter C DH(I) = hydraulic diameter of duct C DCHAR(I) = characteristic dimension of duct for convection: C outside diameter for round duct, C height for rectangular duct C TAIR(I), TAIRNEW(I) = average air temperature inside duct C TWO(I), TWONEW(I) = average exterior temperature of duct wall C TIN(I) = temperature at inlet to duct I C TEXIT(I) = temperature at exit of duct I C TSURR(I) = effective temperature of surroundings of duct I C TEXTER(I) = effective temperature of surroundings of return duct I C T(I) = air temperature at node I C UPRIME = effective U-value between inside duct air and TSURR(I), C including leakage C HI(I) = convection coefficient inside duct C HO(I) = convection coefficient on outside of duct C HR(I,J) = radiation coefficient between duct I and attic surface J C QCON = total heat convected from all ducts to attic air C QRAD(J) = total heat radiated from all ducts to attic surface J C MLEAKS = total mass of air leaked from ducts C QLEAKS = total heat flow due to leakage from ducts to attic air C ONTIME = fraction that HVAC equipment runs during hour time step IMPLICIT REAL*8 (A-H, O-Z) REAL*8 MDOTin,MDOTlk,MLEAKS,MDOTAVG COMMON /DUCT1/ALD(25),xMDOTin(25,2),xMDOTlk(25,2),EODUCT(25), & ZONE(25),PID(25),PO(25),RD(25),xTDUCTin(2),RHOC(25),RDO(25), & RDI(25),NSUPLY,NRETRN,ND(25),NINLET(25),NEXIT(25) COMMON /DUCT2/DH(25),DCHAR(25),AD(25),ISHAPE(25) COMMON /DUCT3/TI,V,ONTIME,DuctTin,DuctMin(25), & DuctLkg(25),KFLAG(19),NTOT COMMON /DUCT4/QCON,QRAD(8),QRADtot,MLEAKS,QLEAKS,QDUCTS,QSTOR, & MDOTin(25),MDOTlk(25) DIMENSION TIS(7,100),TA(2) DIMENSION G(32,32),TAIR(25),TWO(25),HI(25),HO(25),HR(25,8) DIMENSION UPRIME(25),TSURR(25),TEXTER(25),TIN(25),TEXIT(25),T(27) DIMENSION TAIRNEW(25),TWONEW(25) C Next line added for model with trusses COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS C Initialize duct average air temperature, average exterior temperature C NTOT = NSUPLY + NRETRN IF(TA(1) .le. 65.0) then ! Heating Mode Jmode = 1 TDUCTin = xTDUCTin(Jmode) ELSE IF(TA(1) .gt. 65.0) then ! Cooling Mode Jmode = 2 TDUCTin = xTDUCTin(Jmode) End if Do 3 I = 1,NTOT MDOTin(i) = xMDOTin(I,Jmode) MDOTlk(i) = xMDOTlk(I,Jmode) 3 CONTINUE DO 10 I = 1,NTOT IF(I.LE.NSUPLY) THEN TAIR(I) = TDUCTIN ELSE TAIR(I) = TI END IF TWO(I) = TA(1) 10 CONTINUE NIT = 0 15 CONTINUE C Calculate interior film coefficient C See SP43 reports for equation DO 20 I = 1,NTOT IF (I.LE.NSUPLY) THEN MDOTAVG = MDOTIN(I) - 0.5*MDOTLK(I)*ALD(I) ELSE MDOTAVG = MDOTIN(I) + 0.5*MDOTLK(I)*ALD(I) END IF HI(I) = (0.00368+1.5D-6*(TAIR(I)-80.))* & ((MDOTAVG/AD(I))**0.8)/(DH(I)**0.2) 20 CONTINUE C Calculate exterior convection coefficient DO 30 I = 1,NTOT CALL HCOND(TWO(I),TA(1),DCHAR(I),V,ISHAPE(I),HO(I)) 30 CONTINUE C Calculate radiation coefficients between duct and attic surfaces IF(ITRUS.EQ.0) THEN NSURF = 7 ELSE NSURF = 8 END IF DO 40 I = 1,NTOT DO 40 J = 1,7 HR(I,J) = HRAD(G(I+NSURF,J),TWO(I),TIS(J,1)) 40 CONTINUE IF(ITRUS.EQ.1) THEN DO 41 I = 1,NTOT C Next line modified 2-12-96 to use index 8 in HRAD instead of J HR(I,8) = HRAD(G(I+NSURF,8),TWO(I),TTRUS(1)) 41 CONTINUE END IF C Calculate summations of radiation coefficients C Calculate Uprime, temperature of surrounding surfaces, effective C temperature of surroundings DO 50 I = 1,NTOT SUMHR = 0. SUMHRT = 0. DO 45 J = 1,7 SUMHR = SUMHR + HR(I,J) SUMHRT = SUMHRT + HR(I,J)*TIS(J,1) 45 CONTINUE IF(ITRUS.NE.0) THEN SUMHR = SUMHR + HR(I,8) SUMHRT = SUMHRT + HR(I,8)*TTRUS(1) END IF UPRIME(I) = 1./(1./HI(I)/PID(I) + RD(I) + 1./(HO(I)+SUMHR)/PO(I)) TSURR(I) = (HO(I)*TA(1) + SUMHRT)/(HO(I) + SUMHR) TEXTER(I) = (MDOTLK(I)*0.24*TA(1) + UPRIME(I)*TSURR(I))/ & (MDOTLK(I)*0.24 + UPRIME(I)) 50 CONTINUE C Calculate duct exit and inlet temperatures N = 0 DO 55 I = 1,NTOT+2 55 T(I) = -460. DO 60 I = 1,NTOT IF(NINLET(I).EQ.1) THEN TIN(I) = TDUCTIN T(1) = TDUCTIN IF(MDOTLK(I).LE.1.D-6) THEN TEXIT(I) = TSURR(I) + (TIN(I) - TSURR(I))* & EXP(-UPRIME(I)*ALD(I)/(MDOTIN(I)*0.24)) ELSE TEXIT(I) = TSURR(I) + (TIN(I) - TSURR(I))* & ((MDOTIN(I) - MDOTLK(I)*ALD(I))/MDOTIN(I))** & (UPRIME(I)/(0.24*MDOTLK(I))) END IF T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 60 CONTINUE IF(N.GE.NSUPLY) GO TO 75 65 CONTINUE DO 70 I = 1,NTOT IF(T(NINLET(I)).GT.-460..AND.T(NEXIT(I)).LT.-459.) THEN TIN(I) = T(NINLET(I)) IF(MDOTLK(I).LE.1.D-6) THEN TEXIT(I) = TSURR(I) + (TIN(I) - TSURR(I))* & EXP(-UPRIME(I)*ALD(I)/(MDOTIN(I)*0.24)) ELSE TEXIT(I) = TSURR(I) + (TIN(I) - TSURR(I))* & ((MDOTIN(I) - MDOTLK(I)*ALD(I))/MDOTIN(I))** & (UPRIME(I)/(0.24*MDOTLK(I))) END IF T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 70 CONTINUE IF(N.LT.NSUPLY) GO TO 65 75 CONTINUE DO 80 I = 1,NTOT IF(NINLET(I).EQ.(NSUPLY+2)) THEN TIN(I) = TI T(NSUPLY+2) = TI IF(ABS(MDOTLK(I)).LE.1.D-6) THEN TEXIT(I) = TSURR(I) + (TIN(I) - TSURR(I))* & EXP(-UPRIME(I)*ALD(I)/(MDOTIN(I)*0.24)) ELSE TEXIT(I) = TEXTER(I) + (TIN(I) - TEXTER(I))* & ((MDOTIN(I) + MDOTLK(I)*ALD(I))/MDOTIN(I))** & ((UPRIME(I)+0.24*MDOTLK(I))/(0.24*MDOTLK(I))) END IF T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 80 CONTINUE IF(N.EQ.NTOT) THEN GO TO 95 ELSE WRITE(11,1000) 1000 FORMAT(1X,'ERROR IN DUCTSS, COUNT OF DUCTS IS WRONG',/) STOP END IF 95 CONTINUE C Calculate average duct air temperatures and average duct exterior C temperatures. DO 100 I = 1,NTOT IF(ABS(MDOTLK(I)).LE.1.D-6) THEN TAIRNEW(I) = TSURR(I) + (TIN(I) - TSURR(I))/ALD(I)/UPRIME(I)* & MDOTIN(I)*0.24*(1.-EXP(-UPRIME(I)*ALD(I)/(MDOTIN(I)*0.24))) ELSE IF(I.LE.NSUPLY) THEN C1 = UPRIME(I)/(0.24*MDOTLK(I)) + 1. C2 = MDOTIN(I)/(ALD(I)*MDOTLK(I)) TAIRNEW(I) = TSURR(I) + (TIN(I) - TSURR(I))*C2* & (1.- (1.-1./C2)**C1)/C1 ELSE C1 = UPRIME(I)/(0.24*MDOTLK(I)) + 2. C2 = MDOTIN(I)/(ALD(I)*MDOTLK(I)) TAIRNEW(I) = TEXTER(I) + (TIN(I) - TEXTER(I))*C2* & (1.- (1.+1./C2)**C1)/C1 END IF END IF TWONEW(I) = TAIRNEW(I) - (TAIRNEW(I) - TSURR(I))*UPRIME(I)* & (1./HI(I)/PID(I) + RD(I)) 100 CONTINUE C Test for convergence of duct air and exterior surface temperatures IFLAGIT = 0 DO 110 I = 1,NTOT IF(ABS(TAIRNEW(I) - TAIR(I)).GT.1.D-3) GO TO 120 IF(ABS(TWONEW(I) - TWO(I)).GT.1.D-3) GO TO 120 110 CONTINUE GO TO 130 120 IFLAGIT = 1 130 CONTINUE DO 200 I = 1,NTOT TAIR(I) = TAIRNEW(I) TWO(I) = TWONEW(I) 200 CONTINUE NIT = NIT + 1 IF(IFLAGIT.EQ.1.AND.NIT.LE.15) GO TO 15 C Calculate total heat losses from ducts QCON = 0. DO 300 I = 1,NTOT QCON = QCON + HO(I)*PO(I)*ALD(I)*(TWO(I) - TA(1))*ONTIME 300 CONTINUE DO 350 J = 1,7 QRAD(J) = 0. DO 350 I = 1,NTOT QRAD(J) = QRAD(J) & + HR(I,J)*PO(I)*ALD(I)*(TWO(I) - TIS(J,1))*ONTIME 350 CONTINUE IF(ITRUS.NE.0) THEN QRAD(8) = 0. DO 360 I = 1,NTOT QRAD(8) = QRAD(8) & + HR(I,8)*PO(I)*ALD(I)*(TWO(I) - TTRUS(1))*ONTIME 360 CONTINUE END IF QRADTOT = 0. DO 375 J = 1,NSURF QRADTOT = QRADTOT + QRAD(J) 375 CONTINUE MLEAKS = 0. QLEAKS = - QCON - QRADTOT DO 400 I = 1,NTOT MLEAKS = MLEAKS + MDOTLK(I)*ALD(I)*ONTIME IF(I.LE.NSUPLY) THEN QLEAKS = QLEAKS + MDOTIN(I)*0.24*TIN(I)*ONTIME & -(MDOTIN(I) - MDOTLK(I)*ALD(I))*0.24*TEXIT(I)*ONTIME ELSE QLEAKS = QLEAKS + MDOTIN(I)*0.24*TIN(I)*ONTIME & -(MDOTIN(I) + MDOTLK(I)*ALD(I))*0.24*TEXIT(I)*ONTIME END IF 400 CONTINUE C Calulate average temperature of air returning to equipment SUM1 = 0. SUM2 = 0. DO 500 I = NSUPLY+1,NTOT SUM1 = SUM1 + (MDOTIN(I) + MDOTLK(I)*ALD(I))*TEXIT(I) SUM2 = SUM2 + MDOTIN(I) + MDOTLK(I)*ALD(I) 500 CONTINUE TRETRN = SUM1/SUM2 RETURN END C*********************************************************************** C** ** C** SUBROUTINE HCOND ** C** ** C*********************************************************************** C This subroutine calculates the convection coefficient on the exterior C of a duct. C See Holman, pp. 243,244,247,248,275 SUBROUTINE HCOND(TS,TA,AL,V,IFLAG,HC) IMPLICIT REAL*8 (A-H, O-Z) REAL*8 NUS,K,MU,NU DT = TS - TA TF = (TS + TA)/2. TF1 = TF TK = (TF + 459.67)/1.8 C THERMAL CONDUCTIVITY OF AIR FROM EQUATION IN NBS CIRC. 564 C THERMAL COND. IN BTU/(HR-FT-F) K = 0.6325D-5*DSQRT(TK)/(1.+(245.4*10.**(-12./TK))/TK)*241.77 C DYNAMIC VISCOSITY OF AIR FROM EQUATION IN NBS CIRC. 564 C LB/(HR-FT) MU = (145.8*TK*DSQRT(TK)/(TK+110.4))*241.90D-7 C PRANDTL NUMBER, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564 PR = 0.7880 - 2.631D-4*TK C VOLUME EXPANSION COEFFICIENT OF AIR, 1/TABS, PERFECT GAS BETA = 1./(TF1+459.67) C DENSITY OF AIR, PERFECT GAS, NBS CIRC. 564, LB/CF RHO = 22.0493/TK C KINEMATIC VISCOSITY, FT2/HR NU = MU/RHO C SPECIFIC HEAT OF AIR, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564, BTU/(LB-F) CP = (3.4763 + 1.066D-4*TK)*0.068559 C RAYLEIGH NUMBER, LEADING COEFFICIENT IS C 32.174*3600*3600, FT/HR2 RA = (4.16975E8)*BETA*RHO*CP*ABS(DT)*(AL**3)/NU/K IF(RA.LE.1.E9) THEN NUS = 0.53*RA**(0.25) ELSE NUS = 0.13*RA**(1./3.) END IF HCN = NUS*K/AL RE = V*AL/NU IF(IFLAG.NE.0) THEN C = 0.102 AN = 0.675 ELSE IF(RE.LE.4.) THEN C = 0.989 AN = 0.330 ELSE IF(RE.LE.40.) THEN C = 0.911 AN = 0.385 ELSE IF(RE.LE.4.E3) THEN C = 0.683 AN = 0.466 ELSE IF(RE.LE.4.E4) THEN C = 0.193 AN = 0.618 ELSE C = 0.0266 AN = 0.805 END IF END IF NUS = C*(RE**AN)*(PR**(1./3.)) HCF = NUS*K/AL HC = (HCF**3 + HCN**3)**(1./3.) RETURN END C*********************************************************************** C** ** C** SUBROUTINE DUCTTR ** C** ** C*********************************************************************** C This subroutine contains algorithms for transient duct model, C with or without leakage in any or all of the duct segments. SUBROUTINE DUCTTR(G,TIS,TO,TA,T,TAD,TWWD,TWOD) C Duct and node numbering scheme: C Supply ducts are numbered 1 through NSUPLY C Return ducts are numbered NSUPLY+1 through NSUPLY+NRETRN C The inlet from the supply source is numbered node 1 C Outlets of supply nodes are numbered 2 through NSUPLY+1 C The inlet to the return duct(s) is node NSUPLY+2 C The outlet(s) from return duct(s) are nodes NSUPLY+2 C through NSUPLY+NRETRN+2 C Nomenclature: C NSUPLY = number of supply ducts C NRETRN = number of return ducts C NTOT = total number of ducts C xTDUCTIN(I,J) = temperature of air entering supply duct J=1 Heating Mode C TI = temperature of air inside house J=2 Cooling Mode C = temperature entering return duct C ND(I) = number of duct I, only used for convenience in input C NINLET(I) = node number at inlet of duct I C NEXIT(I) = node number at exit of duct I C ISHAPE(I) = flag for shape of duct: 0 for round, 1 for rectangular C ALD(I) = length of duct I C MDOTIN(I) = mass flow rate entering duct I J=1 Heating Mode C MDOTLK(I) = mass leakage per unit length for duct I J=2 Cooling Mode C EODUCT(I) = emittance of outside of duct I C ZONE(I) = zone of house that duct I supplies C DI(I) = inside diameter of round duct C DW(I) = inside diameter of outer layer of round duct wall C DO(I) = outside diameter of round duct C XI(I), Y(I) = inside width and height of rectangular duct C XW(I), YW(I) = inside width and height of outer layer of C rectangular duct C XO(I), YO(I) = outside width and height of rectangular duct C PID(I) = inside perimeter of duct C PW(I) = inside perimeter of outer layer of duct C PO(I) = outside perimeter of duct C RD(I) = thermal resistance of duct wall (usual resistance divided C by perimeter C DH(I) = hydraulic diameter of duct C DCHAR(I) = characteristic dimension of duct for convection: C outside diameter for round duct, C height for rectangular duct C TAIR(I), TAIRNEW(I) = average air temperature inside duct C TWO(I), TWONEW(I) = average exterior temperature of duct wall C TIN(I) = temperature at inlet to duct I C TEXIT(I) = temperature at exit of duct I C TSURR(I) = effective temperature of surroundings of duct I C TEXTER(I) = effective temperature of surroundings of return duct I C T(I) = air temperature at node I, averaged over ontime C UPRIME = effective U-value between inside duct air and TSURR(I), C including leakage C HI(I) = convection coefficient inside duct C HO(I) = convection coefficient on outside of duct C HR(I,J) = radiation coefficient between duct I and attic surface J C QCON = total heat convected from all ducts to attic air C QRAD(J) = total heat radiated from all ducts to attic surface J C MLEAKS = total mass of air leaked from ducts C QLEAKS = total heat flow due to leakage from ducts to attic air C ONTIME = fraction that HVAC equipment runs during hour time step IMPLICIT REAL*8 (A-H, O-Z) REAL*8 MDOTin,MDOTlk,MLEAKS,MDOTAVG Character*2 Blank COMMON /DUCT1/ALD(25),xMDOTin(25,2),xMDOTlk(25,2),EODUCT(25), & ZONE(25),PID(25),PO(25),RD(25),xTDUCTin(2),RHOC(25),RDO(25), & RDI(25),NSUPLY,NRETRN,ND(25),NINLET(25),NEXIT(25) COMMON /DUCT2/DH(25),DCHAR(25),AD(25),ISHAPE(25) COMMON /DUCT3/TI,V,ONTIME,DuctTin,DuctMin(25), & DuctLkg(25),KFLAG(19),NTOT COMMON /DUCT4/QCON,QRAD(8),QRADtot,MLEAKS,QLEAKS,QDUCTS,QSTOR, & MDOTin(25),MDOTlk(25) COMMON /DUCT5/DX(25),Rduct,NSEC(25) DIMENSION TIS(7,100),TA(2) DIMENSION G(32,32),TAIR(25),TWO(25),HI(25),HO(25),HR(25,8) DIMENSION TSURR(25),TIN(25),TEXIT(25),T(27) DIMENSION SUMT(27),TADSUM(27) Data Blank / ' '/ C TAD(I,J) = air temperature in duct I, section J C TWWD(I,J,K) = mid duct wall temperature in duct I, section J, at time K C TWOD(I,J) = outer duct wall temperature in duct I, section J DIMENSION TAD(25,100),TWWD(25,100,2),TWOD(25,100),HPRIME(25) C Next line added for model with trusses COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS C Initialize variables C NTOT = NSUPLY + NRETRN DT = 1./60./3. TWWDSUM = 0.0 TIME = 0.0 TIMEON = 0.0 QCON = 0.0 DO 1 J = 1,7 QRAD(J) = 0. 1 CONTINUE IF(ITRUS.NE.0) THEN QRAD(8) = 0. END IF QRADTOT = 0.0 MLEAKS = 0.0 QLEAKS = 0.0 QSTOR = 0.0 QDUCTS = 0.0 DO 2 I = 1,NTOT+2 TADSUM(I) = 0.0 SUMT(I) = 0.0 2 CONTINUE IF (KFLAG(19) .EQ. 1) go to 4 IF(TO .le. 65.0) then ! Heating Mode Jmode = 1 TDUCTin = xTDUCTin(Jmode) ELSE IF(TO .gt. 65.0) then ! Cooling Mode Jmode = 2 TDUCTin = xTDUCTin(Jmode) End if Do 3 I = 1,NTOT MDOTin(i) = xMDOTin(I,Jmode) MDOTlk(i) = xMDOTlk(I,Jmode) 3 CONTINUE GO TO 6 C 4 CONTINUE TDUCTin = DuctTin ! Input from Auxiliary file Do 5 I = 1,NTOT MDOTin(i) = DuctMin(I) MDOTlk(i) = DuctLkg(I) 5 CONTINUE 6 CONTINUE C C CALCULATE CYCLE TIME ICYCL = 1 ANCYCL = 4.*3.*ONTIME*(1.-ONTIME) IF(ONTIME.GE.1.) THEN DTON = 1.0 DTOFF = 0.0 ELSE IF (ONTIME .LE. 0.0) THEN ICYCL = 0 DTON = 0.0 DTOFF = 1.0 ELSE DTON = ONTIME/ANCYCL DTOFF = (1.-ONTIME)/ANCYCL END IF C CALCULATE FLOW RATES BASED ON ON-OFF CYCLES IF(TIME.GT.DTON.AND.TIME.LE.(DTON+DTOFF)) THEN ICYCL = 0 END IF IF(TIME.GT.(2.*DTON+DTOFF).AND.TIME.LE.(2.*(DTON+DTOFF)))THEN ICYCL = 0 END IF IF(TIME.GT.(3.*DTON+2.*DTOFF).AND.TIME.LE.(3.*(DTON+DTOFF))) & THEN ICYCL = 0 END IF C CALCULATE AVERAGE TAIR AND TWO DO 10 I = 1,NTOT SUM1 = 0.0 NP1 = NSEC(I) + 1 IF(NP1 .GT. 100) THEN WRITE(6,16) NP1 16 FORMAT(' 1-FT DUCT SECTIONS > 100! NP1 = ',I3) STOP END IF DO 12 J = 1,NP1 SUM1 = SUM1 + TAD(I,J) 12 CONTINUE ANP1 = NP1 TAIR(I) = SUM1/ANP1 SUM2 = 0.0 DO 14 J = 1,NSEC(I) SUM2 = SUM2 + TWOD(I,J) 14 CONTINUE ANSEC = NSEC(I) TWO(I) = SUM2/ANSEC 10 CONTINUE C Calculate interior film coefficient C See SP43 reports for equation DO 20 I = 1,NTOT IF (I.LE.NSUPLY) THEN MDOTAVG = MDOTIN(I) - 0.5*MDOTLK(I)*ALD(I) ELSE MDOTAVG = MDOTIN(I) + 0.5*MDOTLK(I)*ALD(I) END IF HI(I) = (0.00368+1.5D-6*(TAIR(I)-80.))* & ((MDOTAVG/AD(I))**0.8)/(DH(I)**0.2) 20 CONTINUE C Calculate exterior convection coefficient DO 30 I = 1,NTOT CALL HCOND(TWO(I),TA(1),DCHAR(I),V,ISHAPE(I),HO(I)) 30 CONTINUE C Calculate radiation coefficients between duct and attic surfaces IF(ITRUS.EQ.0) THEN NSURF = 7 ELSE NSURF = 8 END IF DO 40 I = 1,NTOT DO 40 J = 1,7 HR(I,J) = HRAD(G(I+NSURF,J),TWO(I),TIS(J,1)) 40 CONTINUE IF(ITRUS.EQ.1) THEN DO 41 I = 1,NTOT C Next line modified 2-12-96 to use index 8 in HRAD instead of J HR(I,8) = HRAD(G(I+NSURF,8),TWO(I),TTRUS(1)) 41 CONTINUE END IF C Calculate summations of radiation coefficients C Calculate Hprime, temperature of surrounding surfaces DO 50 I = 1,NTOT SUMHR = 0. SUMHRT = 0. DO 45 J = 1,7 SUMHR = SUMHR + HR(I,J) SUMHRT = SUMHRT + HR(I,J)*TIS(J,1) 45 CONTINUE IF(ITRUS.NE.0) THEN SUMHR = SUMHR + HR(I,8) SUMHRT = SUMHRT + HR(I,8)*TTRUS(1) END IF HPRIME(I) = 1./(RDO(I) + 1./(HO(I)+SUMHR)/PO(I)) TSURR(I) = (HO(I)*TA(1) + SUMHRT)/(HO(I) + SUMHR) 50 CONTINUE C Calculate duct exit and inlet temperatures IF(ICYCL.EQ.1) THEN N = 0 DO 55 I = 1,NTOT+2 55 T(I) = -460. DO 60 I = 1,NTOT IF(NINLET(I).EQ.1) THEN TIN(I) = TDUCTIN T(1) = TDUCTIN TAD(I,1) = T(1) NP1 = NSEC(I) + 1 DO 61 J = 2,NP1 AJ = J TAD(I,J) = TAD(I,J-1) + (TWWD(I,J-1,1)-TAD(I,J-1))* & 1./(1./HI(I)/PID(I)+RDI(I))*DX(I)/0.24/ & (MDOTIN(I) - (AJ-1.5)*DX(I)*MDOTLK(I)) 61 CONTINUE TEXIT(I) = TAD(I,NP1) T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 60 CONTINUE IF(N.GE.NSUPLY) GO TO 75 65 CONTINUE DO 70 I = 1,NTOT IF(T(NINLET(I)).GT.-460..AND.T(NEXIT(I)).LT.-459.) THEN TIN(I) = T(NINLET(I)) TAD(I,1) = TIN(I) NP1 = NSEC(I) + 1 DO 71 J = 2,NP1 AJ = J TAD(I,J) = TAD(I,J-1) + (TWWD(I,J-1,1)-TAD(I,J-1))* & 1./(1./HI(I)/PID(I)+RDI(I))*DX(I)/0.24/ & (MDOTIN(I) - (AJ-1.5)*DX(I)*MDOTLK(I)) 71 CONTINUE TEXIT(I) = TAD(I,NP1) T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 70 CONTINUE IF(N.LT.NSUPLY) GO TO 65 75 CONTINUE DO 80 I = 1,NTOT IF(NINLET(I).EQ.(NSUPLY+2)) THEN TIN(I) = TI T(NSUPLY+2) = TI TAD(I,1) = TIN(I) NP1 = NSEC(I) + 1 DO 81 J = 2,NP1 AJ = J TAD(I,J) = TAD(I,J-1) + & ((TWWD(I,J-1,1)-TAD(I,J-1))*(1./(1./HI(I)/PID(I)+RDI(I))) & +(TA(1) - TAD(I,J-1))*MDOTLK(I)*0.24)*DX(I)/0.24/ & (MDOTIN(I) + (AJ-1.5)*DX(I)*MDOTLK(I)) 81 CONTINUE TEXIT(I) = TAD(I,NP1) T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 80 CONTINUE IF(N.EQ.NTOT) THEN GO TO 95 ELSE WRITE(11,1000) 1000 FORMAT(1X,'ERROR IN DUCTTR, COUNT OF DUCTS IS WRONG',/) STOP END IF ELSE N = 0 DO 155 I = 1,NTOT+2 155 T(I) = -460. DO 160 I = 1,NTOT IF(NINLET(I).EQ.1) THEN TIN(I) = TWWD(I,1,1) T(1) = TWWD(I,1,1) TAD(I,1) = T(1) NP1 = NSEC(I) + 1 DO 161 J = 2,NP1 AJ = J TAD(I,J) = TWWD(I,J-1,1) 161 CONTINUE TEXIT(I) = TAD(I,NP1) T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 160 CONTINUE IF(N.GE.NSUPLY) GO TO 175 165 CONTINUE DO 170 I = 1,NTOT IF(T(NINLET(I)).GT.-460..AND.T(NEXIT(I)).LT.-459.) THEN TIN(I) = T(NINLET(I)) TAD(I,1) = TIN(I) NP1 = NSEC(I) + 1 DO 171 J = 2,NP1 AJ = J TAD(I,J) = TWWD(I,J-1,1) 171 CONTINUE TEXIT(I) = TAD(I,NP1) T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 170 CONTINUE IF(N.LT.NSUPLY) GO TO 165 175 CONTINUE DO 180 I = 1,NTOT IF(NINLET(I).EQ.(NSUPLY+2)) THEN TIN(I) = TWWD(I,1,1) T(NSUPLY+2) = TWWD(I,1,1) TAD(I,1) = TIN(I) NP1 = NSEC(I) + 1 DO 181 J = 2,NP1 AJ = J TAD(I,J) = TWWD(I,J-1,1) 181 CONTINUE TEXIT(I) = TAD(I,NP1) T(NEXIT(I)) = TEXIT(I) N = N + 1 END IF 180 CONTINUE IF(N.EQ.NTOT) THEN GO TO 95 ELSE WRITE(11,1000) STOP END IF END IF 95 CONTINUE C Calculate duct wall temperatures IF(ICYCL.EQ.1) THEN DO 2000 I = 1,NTOT DO 2000 J = 1,NSEC(I) AJ = J IF(I.LE.NSUPLY) THEN TWWD(I,J,2) = TWWD(I,J,1) + DT/RHOC(I)/DX(I)* & ((MDOTIN(I) - (AJ-0.5)*DX(I)*MDOTLK(I))*0.24 & *(TAD(I,J)-TAD(I,J+1)) & + DX(I)*HPRIME(I)*(TSURR(I) - TWWD(I,J,1))) ELSE TWWD(I,J,2) = TWWD(I,J,1) + DT/RHOC(I)/DX(I)* & ((MDOTIN(I) + (AJ-0.5)*DX(I)*MDOTLK(I))*0.24 & *(TAD(I,J)-TAD(I,J+1)) & + DX(I)*HPRIME(I)*(TSURR(I) - TWWD(I,J,1)) & + MDOTLK(I)*DX(I)*0.24*(TA(1) - (TAD(I,J)+TAD(I,J+1))/2.)) END IF 2000 CONTINUE ELSE DO 2001 I = 1,NTOT DO 2001 J = 1,NSEC(I) AJ = J TWWD(I,J,2) = TWWD(I,J,1) + DT/RHOC(I)* & HPRIME(I)*(TSURR(I) - TWWD(I,J,1)) 2001 CONTINUE END IF C Calculate duct exterior temperature DO 2100 I = 1,NTOT DO 2100 J = 1,NSEC(I) TWOD(I,J) = TWWD(I,J,2) & + RDO(I)*HPRIME(I)*(TSURR(I) - TWWD(I,J,2)) 2100 CONTINUE C Calculate total heat losses from ducts DO 300 I = 1,NTOT DO 300 J = 1,NSEC(I) QCON = QCON + HO(I)*PO(I)*DX(I)*(TWOD(I,J) - TA(1))*DT 300 CONTINUE DO 350 J = 1,7 DO 350 I = 1,NTOT DO 350 JJ = 1,NSEC(I) QRAD(J) = QRAD(J) + & HR(I,J)*PO(I)*DX(I)*(TWOD(I,JJ) - TIS(J,1))*DT 350 CONTINUE IF(ITRUS.NE.0) THEN DO 360 I = 1,NTOT DO 360 JJ = 1,NSEC(I) QRAD(8) = QRAD(8) & + HR(I,8)*PO(I)*DX(I)*(TWOD(I,JJ) - TTRUS(1))*DT 360 CONTINUE END IF IF(ICYCL.EQ.1) THEN DO 400 I = 1,NTOT MLEAKS = MLEAKS + MDOTLK(I)*ALD(I)*DT IF(I.LE.NSUPLY) THEN QDUCTS = QDUCTS + MDOTIN(I)*0.24*TIN(I)*DT & -(MDOTIN(I) - MDOTLK(I)*ALD(I))*0.24*TEXIT(I)*DT ELSE QDUCTS = QDUCTS + MDOTIN(I)*0.24*TIN(I)*DT & -(MDOTIN(I) + MDOTLK(I)*ALD(I))*0.24*TEXIT(I)*DT END IF 400 CONTINUE END IF DO 450 I = 1,NTOT DO 450 J = 1,NSEC(I) QSTOR = QSTOR + RHOC(I)*DX(I)*(TWWD(I,J,2) - TWWD(I,J,1)) 450 CONTINUE C Calulate average temperature of air returning to equipment IF(ICYCL.EQ.1) THEN SUM1 = 0. SUM2 = 0. DO 500 I = NSUPLY+1,NTOT SUM1 = SUM1 + (MDOTIN(I) + MDOTLK(I)*ALD(I))*TEXIT(I) SUM2 = SUM2 + MDOTIN(I) + MDOTLK(I)*ALD(I) 500 CONTINUE TRETRN = SUM1/SUM2 END IF C Calculate average temperature of air at nodes IF(ICYCL.EQ.1) THEN DO 510 I = 1,NTOT+2 SUMT(I) = SUMT(I) + T(I)*DT 510 CONTINUE TIMEON = TIMEON + DT END IF C Shift index on TWWD for next time step DO 550 I = 1,NTOT DO 550 J = 1,NSEC(I) TWWD(I,J,1) = TWWD(I,J,2) 550 CONTINUE TIME = TIME + DT IF(TIME.LT.1.) GO TO 6 DO 600 J = 1,7 QRADTOT = QRADTOT + QRAD(J) 600 CONTINUE IF(ITRUS.NE.0) THEN QRADTOT = QRADTOT + QRAD(8) END IF QLEAKS = QDUCTS - (QCON + QRADTOT + QSTOR) IF (ONTIME .LE. 0.0) THEN QDUCTS = 0.0 QLEAKS = 0.0 END IF C Calculate average temperature of air at nodes during on cycles DO 650 I = 1,NTOT+2 IF (ONTIME .GT. 0.0) T(I) = SUMT(I)/TIMEON 650 CONTINUE IF (ONTIME .LE. 0.0) THEN DO 4002 K = 1,NTOT DO 4001 J = 1,NSEC(K) 4001 TADSUM(K) = TADSUM(K) + TAD(K,J) 4002 T(K) = TADSUM(K)/DFLOTJ(NSEC(K)) T(NTOT+1) = TI T(NTOT+2) = T(NTOT) END IF RETURN END C*********************************************************************** C** ** C** SUBROUTINE TRUSREAD ** C** ** C*********************************************************************** C This subroutine reads input for the trusses in the attic SUBROUTINE TRUSREAD C Nomenclature: C AOPEN = open area of trusses between framing C ATRUS = exposed surface area of truss framing C VTRUS = volume of truss framing C ALTRUS = characteristic length for convection from truss framing C ETRUS = emittance of truss framing C CPTRUS = specific heat of truss framing C RHOTRUS = density of truss framing C AMASST = mass per unit area of truss that participates in moisture C exchange C AMCT = moisture content (weight fraction) of trusses IMPLICIT REAL*8 (A-H, O-Z) COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS COMMON /TRUS2/AMASST,AMCT,HCTRUS,AMWT,BMWT,QLAT READ(5,*) AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS READ(5,*) AMASST,AMCT RETURN END C********************************************************************** C** ** C** SUBROUTINE VIEW2T ** C** ** C********************************************************************** SUBROUTINE VIEW2T(AL,W,PITCH1,PITCH2,H1,EI,IDUCT,G) C F(I,J) = VIEW FACTOR FROM SURFACE I TO SURFACE J C FF(I,J,K) = VIEW FACTOR FROM SURFACE I TO SURFACE J OF ATTIC WITH C K JOIST SPACES C FFF(I,J,K) = VIEW FACTOR FROM STRIP OF SURFACE I TO STRIP OF SURFACE C J. K = 1 INDICATES SURFACES IN SAME JOIST SPACE, K = 2 SURFACES C IN ADJACENT JOIST SPACES, K = 3 SURFACES WITH 1 JOIST SPACE IN C BETWEEN, ETC. IMPLICIT REAL*8 (A-H, O-Z) REAL*8 MDOTIN,MDOTLK DIMENSION EI(7),A(8),G(32,32) DIMENSION F(7,7),FF(7,7,30),FFF(7,7,30),X(7,7) DIMENSION H(30),T(30),S(7,7),ST(7,7) DIMENSION FNEW(33,33) DIMENSION CHI(33,33),PSI(33,33),E(33) DIMENSION B(33),EMIT(33) COMMON /DUCT1/ALD(25),xMDOTin(25,2),xMDOTlk(25,2),EODUCT(25), & ZONE(25),PID(25),PO(25),RD(25),xTDUCTin(2),RHOC(25),RDO(25), & RDI(25),NSUPLY,NRETRN,ND(25),NINLET(25),NEXIT(25) COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS C NSPAC = IDINT((AL + 1.)/2.) ANSPAC = NSPAC K = 1 C LOOP FOR CALCULATING VIEW FACTORS FOR ATTICS OF WITH 1, 2, C 3, 4 .... NSPAC JOIST SPACES 10 CONTINUE AKK = K ALL = AL*AKK/ANSPAC C 10 WRITE(11,1000) AL C1000 FORMAT(1X,'AL = ',F10.0) C CALCULATE CHARACTERISTIC LENGTHS AND AREAS OF SURFACES C CHAR. LENGTHS ARE TAKEN TO BE: C CEILING - AVERAGE OF LENGTH AND WIDTH C ROOF - DISTANCE FROM EAVE TO RIDGE C GABLE - AVERAGE HEIGHT C SIDE WALLS - HEIGHT AL1 = (ALL+W)/2. P1 = PITCH1*3.14159265/180. P2 = PITCH2*3.14159265/180. P3 = 3.14159265 - P1 - P2 IP1 = PITCH1 IP2 = PITCH2 IP3 = 180.0 - PITCH1 - PITCH2 IF(IP3.EQ.180) GO TO 5 AL2 = W*DSIN(P2)/DSIN(P3) AL3 = W*DSIN(P1)/DSIN(P3) AL4 = 0.5*(AL2*DSIN(P1)) + H1 GO TO 6 5 AL2 = W/2.0 AL3 = W/2.0 AL4 = H1 6 AL5 = AL4 AL6 = H1 AL7 = H1 A(1) = ALL*W A(2) = ALL*AL2 A(3) = ALL*AL3 A(4) = W*AL4 A(5) = A(4) A(6) = ALL*H1 A(7) = ALL*H1 CALL VIEW2(ALL,W,PITCH1,PITCH2,H1,EI,G,7,F) DO 20 I = 1,7 DO 20 J = 1,7 FF(I,J,K) = F(I,J) 20 CONTINUE C WRITE(11,2000) C WRITE(11,*) K C DO 30 I = 1,7 C WRITE(11,2000) (FF(I,J,K),J=1,7) C 30 CONTINUE C2000 FORMAT(1X,7F10.8) K = K + 1 IF(K.LE.NSPAC) GO TO 10 C LOOP FOR CALCULATING VIEW FACTORS FROM ONE JOIST SPACE TO C ANOTHER JOIST SPACE, FFF(I,J,K) DO 100 M = 1,7 DO 100 N = 1,7 IF(M.EQ.N) GO TO 100 FFF(M,N,1) = FF(M,N,1) FFF(M,N,2) = FF(M,N,2) - FFF(M,N,1) KK = 3 40 X(M,N) = 0.0 DO 50 K = 2,KK-1 X(M,N) = X(M,N) + (KK - K + 1)*FFF(M,N,K) 50 CONTINUE FFF(M,N,KK) = KK*(FF(M,N,KK) - FF(M,N,1))/2. - X(M,N) KK = KK + 1 IF(KK.LE.NSPAC) GO TO 40 100 CONTINUE DO 150 M = 1,7 IF(M.EQ.4.OR.M.EQ.5) GO TO 150 FFF(M,4,1) = FF(M,4,1) FFF(M,5,1) = FF(M,5,1) KK = 2 60 X(M,4) = 0.0 X(M,5) = 0.0 DO 70 K = 1,KK-1 X(M,4) = X(M,4) + FFF(M,4,K) X(M,5) = X(M,5) + FFF(M,5,K) 70 CONTINUE FFF(M,4,KK) = KK*FF(M,4,KK) - X(M,4) FFF(M,5,KK) = KK*FF(M,5,KK) - X(M,5) KK = KK + 1 IF(KK.LE.NSPAC) GO TO 60 150 CONTINUE DO 101 K = 1,NSPAC DO 101 M = 1,7 FFF(M,M,K) = 0.0 101 CONTINUE C DO 170 M = 1,7 C DO 170 K = 1,NSPAC C WRITE(11,*) (FFF(M,N,K), N = 1,7) C 170 CONTINUE C CALCULATE INTEGRATED VIEW FACTORS WITHOUT AND WITH TRUSS WEBS C LOOP 200 IS FOR ALL RECTANGULAR SURFACES DO 200 M = 1,7 DO 200 N = 1,7 IF(M.EQ.N) GO TO 200 IF(M.EQ.4.OR.M.EQ.5.OR.N.EQ.4.OR.N.EQ.5) GO TO 200 DO 220 I = 1,NSPAC H(I) = 0.0 T(I) = 0.0 DO 230 J = 1,NSPAC+1-I H(I) = H(I) + FFF(M,N,J) T(I) = T(I) + FFF(M,N,J)*AOPEN**(J-1) 230 CONTINUE IF(I.NE.1) THEN DO 240 J = 2,I H(I) = H(I) + FFF(M,N,J) T(I) = T(I) + FFF(M,N,J)*AOPEN**(J-1) 240 CONTINUE END IF 220 CONTINUE S(M,N) = 0.0 ST(M,N) = 0.0 DO 250 I = 1,NSPAC S(M,N) = S(M,N) + H(I) ST(M,N) = ST(M,N) + T(I) 250 CONTINUE S(M,N) = S(M,N)/NSPAC ST(M,N) = ST(M,N)/NSPAC 200 CONTINUE DO 201 M = 1,7 S(M,M) = 0.0 ST(M,M) = 0.0 201 CONTINUE C LOOP 300 IS FOR RECTANGULAR SURFACES TO THE GABLES DO 300 M = 1,7 IF(M.EQ.4.OR.M.EQ.5) GO TO 300 S(M,4) = 0.0 ST(M,4) = 0.0 DO 310 J = 1,NSPAC S(M,4) = S(M,4) + FFF(M,4,J) ST(M,4) = ST(M,4) + FFF(M,4,J)*AOPEN**(J - 1) 310 CONTINUE S(M,4) = S(M,4)/NSPAC ST(M,4) = ST(M,4)/NSPAC S(M,5) = S(M,4) ST(M,5) = ST(M,4) 300 CONTINUE C LOOP 400 IS FOR GABLES TO RECTANGULAR SURFACES C NOTE THAT ATTIC SURFACE AREAS CALCULATED IN THE LAST LOOP THROUGH C LOOP 10 ARE THE VALUES FOR THE ENTIRE ATTIC DO 400 M = 1,7 IF(M.EQ.4.OR.M.EQ.5) GO TO 400 S(4,M) = S(M,4)*A(M)/A(4) S(5,M) = S(4,M) ST(4,M) = ST(M,4)*A(M)/A(4) ST(5,M) = ST(4,M) 400 CONTINUE C CALCULATE VIEW FACTOR BETWEEN TWO GABLES S(4,5) = FF(4,5,NSPAC) S(5,4) = S(4,5) ST(4,5) = FF(4,5,NSPAC)*AOPEN**(NSPAC - 1) ST(5,4) = ST(4,5) C LOOP 500 IS FOR CALCULATION OF NEW SET OF VIEW FACTORS C AMONG ATTIC SURFACES AND TRUSSES DO 500 M = 1,7 SUM = 0.0 DO 501 N = 1,7 FNEW(M,N) = ST(M,N) SUM = SUM + FNEW(M,N) 501 CONTINUE FNEW(M,8) = 1. - SUM 500 CONTINUE SUM = 0.0 DO 510 M = 1,7 FNEW(8,M) = FNEW(M,8)*A(M)/ATRUS SUM = SUM + FNEW(8,M) 510 CONTINUE FNEW(8,8) = 1. - SUM C DO 600 I = 1,7 C WRITE (11,*) (S(I,J), J = 1,7) C 600 CONTINUE C WRITE(11,*) C DO 601 I = 1,7 C WRITE(11,*) (ST(I,J),J = 1,7) C 601 CONTINUE C WRITE(11,*) C DO 602 I = 1,8 C WRITE(11,*) (FNEW(I,J),J = 1,8) C 602 CONTINUE A(8) = ATRUS DO 701 I = 1,7 701 EMIT(I) = EI(I) EMIT(8) = ETRUS C NEXT SECTION ADDED FOR MODEL WITH DUCTS IF(IDUCT.EQ.0) GO TO 3000 NTOT = NSUPLY + NRETRN NTOT1 = 8 + NTOT DO 702 I = 1,NTOT 702 EMIT(I+8) = EODUCT(I) DO 760 I = 1,NTOT FNEW(I+8,1) = 0.5 760 FNEW(1,I+8) = PO(I)*ALD(I)/A(1)*FNEW(I+8,1) SUM = 0.0 DO 770 I = 1,NTOT 770 SUM = SUM + FNEW(1,I+8) DO 780 J = 2,8 FNEW(1,J) = (1. - SUM)*FNEW(1,J) 780 FNEW(J,1) = (1. - SUM)*FNEW(J,1) SUMA = 0.0 DO 785 I = 1,NTOT 785 SUMA = SUMA + PO(I)*ALD(I) DO 795 I = 1,NTOT DO 795 J = 2,8 FNEW(I+8,J) = A(J)/SUMA*FNEW(J,1)*(1./(1.-SUM) - 1.) 795 FNEW(J,I+8) = FNEW(I+8,J)*PO(I)*ALD(I)/A(J) DO 797 I = 1,NTOT DO 797 J = 1,NTOT 797 FNEW(I+8,J+8) = 0.0 GO TO 3001 C ABOVE SECTION ADDED FOR MODEL WITH DUCTS 3000 CONTINUE NTOT1 = 8 3001 CONTINUE DO 800 I = 1,NTOT1 DO 800 J = 1,NTOT1 800 CHI(I,J) = -(1. - EMIT(I))*FNEW(I,J)/EMIT(I) DO 900 I = 1,NTOT1 900 CHI(I,I) = CHI(I,I) + 1./EMIT(I) C INVERT CHI MATRIX TO OBTAIN PSI MATRIX DO 1200 L = 1,NTOT1 DO 1000 I = 1,NTOT1 E(I) = 0.0 1000 IF(I.EQ.L) E(I) = 1.0 CALL SOLVP(NTOT1,CHI,E,B,33) DO 1100 I = 1,NTOT1 1100 PSI(I,L) = B(I) 1200 CONTINUE DO 1300 I = 1,NTOT1 DO 1300 J = 1,NTOT1 1300 G(I,J) = PSI(I,J)*EMIT(I)/(1. - EMIT(I)) C DO 1400 I = 1,NTOT1 C1400 WRITE(11,2000) (FNEW(I,J),J=1,NTOT1) C WRITE(11,2000) C DO 1500 I = 1,NTOT1 C1500 WRITE(11,2000) (CHI(I,J),J=1,NTOT1) C WRITE(11,2000) C DO 1600 I = 1,NTOT1 C1600 WRITE(11,2000) (PSI(I,J),J=1,NTOT1) C WRITE(11,2000) C DO 1700 I = 1,NTOT1 C1700 WRITE(11,2000) (G(I,J),J=1,NTOT1) C WRITE(11,2000) C2000 FORMAT(1X,7F10.4) RETURN END C********************************************************************** C** ** C** SUBROUTINE TRUSMOD ** C** ** C********************************************************************** SUBROUTINE TRUSMOD(G,TIS,TA,V,IDUCT,TTRNEW) C This subroutine performs the energy balance on the truss webs C Next line added 12-27-94 IMPLICIT REAL*8 (A-H, O-Z) REAL*8 MDOTIN,MDOTLK,MLEAKS COMMON /TRUS1/AOPEN,ATRUS,VTRUS,ALTRUS,ETRUS,CPTRUS,RHOTRUS, & TTRUS(2),QRADT(7),QCONT,AMTRUS,ITRUS COMMON /TRUS2/AMASST,AMCT,HCTRUS,AMWT,BMWT,QLAT COMMON /DUCT4/QCON,QRAD(8),QRADtot,MLEAKS,QLEAKS,QDUCTS,QSTOR, & MDOTin(25),MDOTlk(25) DIMENSION G(32,32),TIS(7,100),TA(2),HRT(7) AMTRUS = RHOTRUS*VTRUS CALL HCON(TTRUS(1),TA(1),90.0D0,ALTRUS,1,1,V,HCTRUS) DO 100 J = 1,7 HRT(J) = HRAD(G(8,J),TTRUS(1),TIS(J,1)) 100 CONTINUE QINC = ATRUS*HCTRUS*(TA(1) - TTRUS(1)) DO 150 J = 1,7 QINC = QINC + ATRUS*HRT(J)*(TIS(J,1) - TTRUS(1)) 150 CONTINUE IF(IDUCT.NE.0) QINC = QINC + QRAD(8) QINC = QINC + AMWT*QLAT + BMWT*QLAT*TTRUS(1) BB = QINC + AMTRUS*CPTRUS*TTRUS(2) AA = AMTRUS*CPTRUS + BMWT*QLAT TTROLD = BB/AA TTRNEW = TTRNEW + 0.15*(TTROLD-TTRNEW) C write(11,*) bb,aa,ttrnew DO 200 J = 1,7 HRT(J) = HRAD(G(8,J),TTRNEW,TIS(J,1)) QRADT(J) = ATRUS*HRT(J)*(TTRNEW - TIS(J,1)) 200 CONTINUE CALL HCON(TTRNEW,TA(1),90.0D0,ALTRUS,1,1,V,HCTRUS) QCONT = ATRUS*HCTRUS*(TTRNEW - TA(1)) RETURN END C C************************************************************************* C** ** C** Subroutine ModRa ** C** ** C************************************************************************* SUBROUTINE MODRA(TASV,TS1,TS2,TA,P,AL,GAP,TMDOT,TCMDOT,TV,RA,RE) C THIS SUBROUTINE CALCULATES THE MODIFIED RAYLEIGH NUMBER SIMILAR TO SUBROUTINE C HCON. THE DT VARIABLE IS MODIFIED SPECIFICALLY FOR THE AIR VENT BETWEEN THE C RADIANT BARRIER AND THE ROOF UNDERSIDE. THE REYNOLDS NUMBER IS CALCULATED FROM C THE RAYLEIGH NUMBER ACCORDING TO A CORRILATION FROM ADAM YOUNGQUIST AND BILL MILLER. C TS1 = SURFACE TEMP OF DECK FACING ROOF F C TS2 = SURFACE TEMP OF ROOF UNDERSIDE F C TA = AIR TEMP ENTERING VENT C P = SLOPE OF AIR CHANNEL RADIANS C AL = CROSS-SECTIONAL LENGTH OF ATTIC (ft) C GAP = HEIGHT OF AIR CHANNEL (ft) IMPLICIT REAL*8 (A-H, O-Z) REAL*8 K,K1,K2,MU1,MU2,MU,NU C CALCULATE FILM TEMPERATURE TF1 = (TS1+TASV)/2.0 TF2 = (TS2+TASV)/2.0 C CALCULATE MODIFIED DT DT = TASV - TA C THERMAL CONDUCTIVITY OF AIR FROM EQUATION IN NBS CIRC. 564 C THERMAL COND. IN BTU/(HR-FT-F) TK1 = (TF1+459.67)/1.8 TK2 = (TF2+459.67)/1.8 K1 = 0.6325D-5*DSQRT(TK1)/(1.+(245.4*10.**(-12./TK1))/TK1)*241.77 K2 = 0.6325D-5*DSQRT(TK2)/(1.+(245.4*10.**(-12./TK2))/TK2)*241.77 K = (K1+K2)/2.0 C DYNAMIC VISCOSITY OF AIR FROM EQUATION IN NBS CIRC. 564 C LB/(HR-FT) MU1 = (145.8*TK1*DSQRT(TK1)/(TK1+110.4))*241.90D-7 MU2 = (145.8*TK2*DSQRT(TK2)/(TK2+110.4))*241.90D-7 MU = (MU1+MU2)/2.0 C VOLUME EXPANSION COEFFICIENT OF AIR, 1/TABS, PERFECT GAS BETA1 = 1./(TF1+459.67) BETA2 = 1./(TF2+459.67) BETA = (BETA1+BETA2)/2.0 C DENSITY OF AIR, PERFECT GAS, NBS CIRC. 564, LB/CF RHO1 = 22.0493/TK1 RHO2 = 22.0493/TK2 RHO = (RHO1+RHO2)/2.0 C KINEMATIC VISCOSITY, FT2/HR NU = MU/RHO C SPECIFIC HEAT OF AIR, LINEAR REGRESSION FROM 220 K TO 360 K C FROM NBS CIRC. 564, BTU/(LB-F) CP1 = (3.4763 + 1.066D-4*TK1)*0.068559 CP2 = (3.4763 + 1.066D-4*TK2)*0.068559 CP = (CP1+CP2)/2.0 C RAYLEIGH NUMBER, LEADING COEFFICIENT IS C 32.174*3600*3600, FT/HR2 RA = (4.16975E8)*BETA*RHO*CP*ABS(DT)*(GAP**3)/NU/K C REYNOLDS NUMBER FROM YOUNGQUIST CORRELATION IF(RA.GT.1.0E4) RE = 3.0668*(RA*DSIN(P))**0.4282 IF(RA.LE.1.0E4) RE = 0.7114*(RA*DSIN(P))**0.5839 C CALCULATE MASS FLOW RATE FROM REYNOLDS NUMBER TMDOT = RE*MU*AL TCMDOT = TMDOT*CP TV = TMDOT/RHO RETURN END