LOGICAL FUNCTION ASHWAT_Thermal(FS, TIN, TOUT, HCIN, HCOUT, &
TRMOUT, TRMIN, ISOL, SOURCE, TOL, &
QOCF, QOCFRoom, T, Q, JF, JB, HC, &
UCG, SHGC, HCInFlag)
!
! SUBROUTINE INFORMATION:
! AUTHOR JOHN L. WRIGHT (University of Waterloo, Mechanical Engineering)
! Chip Barnaby (WrightSoft)
!
! DATE WRITTEN LATEST MODIFICATIONS, February 2008
!
! MODIFIED Bereket Nigusse, June 2013
! added standard 155099 inside convection
! coefficient calculation for U-Factor
!
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
! Subroutine to calculate the glazing temperatures of the
! various elements of a window/shade array while solving an energy
! balance which accounts for absorbed solar radiation, indoor-
! outdoor temperature difference, any combination of hemispherical
! IR optical properties of the various glazings/shading layers.
! Mean radiant temperatures can differ from air temperature on
! both the indoor and outdoor sides.
! It is also possible to allow air-flow between the two layers
! adjacent to the indoor side and/or the two layers adjacent the
! outdoor side. U-factor and SHGC calculations are also included (optional)
! METHODOLOGY EMPLOYED:
! Uses the net radiation method developed for ASHWAT fenestration
! model by John Wright, the University of WaterLoo
! REFERENCES:
! ASHRAE RP-1311
! USE STATEMENTS:
! na
IMPLICIT NONE
! FUNCTION ARGUMENT DEFINITIONS:
TYPE( CFSTY), INTENT( IN) :: FS ! fenestration system
! FS%NL determines # of layers modelled
REAL(r64), INTENT( IN) :: TIN, TOUT ! indoor / outdoor air temperature, K
REAL(r64), INTENT( IN) :: HCIN, HCOUT ! indoor / outdoor convective heat transfer
! coefficient, W/m2K
REAL(r64), INTENT( IN) :: TRMIN, TRMOUT ! indoor / outdoor mean radiant temp, K
REAL(r64), INTENT( IN) :: ISOL ! total incident solar, W/m2 (values used for SOURCE derivation)
! = outside direct + outside diffuse + inside diffuse
REAL(r64), INTENT( IN) :: SOURCE( FS%NL+1) ! absorbed solar by layer, W/m2
REAL(r64), INTENT( IN) :: TOL ! convergence tolerance, usually
! 0.001 (good) or 0.0001 (tight)
REAL(r64), INTENT( OUT) :: QOCF( FS%NL) ! returned: heat flux to layer i from gaps i-1 and i
! due to open channel flow, W/m2
REAL(r64), INTENT( OUT) :: QOCFRoom ! returned: open channel heat gain to room, W/m2
REAL(r64), INTENT( OUT) :: T( FS%NL) ! returned: layer temperatures, 1=outside-most layer, K
REAL(r64), INTENT( OUT) :: Q(0:) ! returned: heat flux at ith gap (betw layers i and i+1), W/m2
! + = heat flow indoor to outdoor
REAL(r64), INTENT( OUT) :: JF( FS%NL+1) ! returned: front (outside facing) radiosity of surfaces, W/m2
! JF( NL+1) = room radiosity
REAL(r64), INTENT( OUT) :: JB(0:FS%NL) ! returned: back (inside facing) radiosity, W/m2
! JB( 0) = outside environment radiosity
REAL(r64), INTENT( OUT) :: HC(0:FS%NL) ! returned: gap convective heat transfer coefficient, W/m2K
! 0=outside, 1=betw layer 1-2, ..., NL=inside
REAL(r64), INTENT( OUT) :: UCG ! returned: center-glass U-factor, W/m2-K
REAL(r64), INTENT( OUT) :: SHGC ! returned: center-glass SHGC (Solar Heat Gain Coefficient)
LOGICAL, OPTIONAL,INTENT( IN) :: HCInFlag ! If true uses ISO Std 150099 routine for HCIn calc
! FUNCTION PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: Height = 1.0d0 ! Window height (m) for standard ratings calculation
CHARACTER(len=MaxNameLength), PARAMETER:: RoutineName= 'ASHWAT_Thermal: '
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! FUNCTION LOCAL VARIABLE DECLARATIONS:
REAL(r64) :: ALPHA, HCOCFout
REAL(r64) :: A( 3*FS%NL+2, 3*FS%NL+4), XSOL( 3*FS%NL+2)
REAL(r64) :: MAXERR
REAL(r64) :: TNEW( FS%NL) ! latest estimate of layer temperatures, K
REAL(r64) :: EB( 0:FS%NL+1) ! black emissive power by layer, W/m2
! EB( 0) = outdoor environment, EB( NL+1) = indoor environment
REAL(r64) :: HHAT(0:FS%NL) ! convective heat transfer coefficient (W/m2.K4)
! based on EB, NOT temperature difference
REAL(r64) :: RHOF_ROOM, TAU_ROOM, EPSF_ROOM ! effective longwave room-side properties
REAL(r64) :: RHOB_OUT, TAU_OUT, EPSB_OUT ! effective longwave outdoor environment properties
REAL(r64) :: QNET( FS%NL) ! checksum - net heat flux to a layer - should be zero - not needed
INTEGER :: ADIM ! dimension of the A matrix
INTEGER :: CONVRG
INTEGER :: NL, I, J, L, ITRY
INTEGER :: hin_scheme ! flags different schemes for indoor convection coefficients
INTEGER ISDL(0:FS%NL+1) ! Flag to mark diathermanous layers, 0=opaque
INTEGER NDLIAR ! Number of Diathermanous Layers In A Row (i.e., consecutive)
INTEGER IB, IE ! Counter begin and end limits
INTEGER IDV ! Integer dummy variable, general utility
INTEGER IM_ON ! Turns on calculation of Indices of Merit if IM_ON=1
REAL(r64) :: QOCF_F(FS%NL) ! heat flux to outdoor-facing surface of layer i, from gap i-1,
! due to open channel flow, W/m2
REAL(r64) :: QOCF_B(FS%NL) ! heat flux to indoor-facing surface of layer i, from gap i,
! due to open channel flow, W/m2
REAL(r64) :: Rvalue ! R-value in IP units [hr.ft2.F/BTU]
REAL(r64) :: TAE_IN, TAE_OUT ! Indoor and outdoor effective ambient temperatures [K]
REAL(r64) :: HR(0:FS%NL) ! Radiant heat transfer coefficient [W/m2K]
REAL(r64) :: HJR(FS%NL),HJC(FS%NL) ! radiative and convective jump heat transfer coefficients
REAL(r64) :: FHR_OUT, FHR_IN ! hre/(hre+hce) fraction radiant h, outdoor or indoor, used for TAE
REAL(r64) :: Q_IN ! net gain to the room [W/m2], including transmitted solar
REAL(r64) :: RHOF(0:FS%NL+1) ! longwave reflectance, front ! these variables help simplify
REAL(r64) :: RHOB(0:FS%NL+1) ! longwave reflectance, back ! the code because it is useful to
REAL(r64) :: EPSF(0:FS%NL+1) ! longwave emisivity, front ! increase the scope of the arrays
REAL(r64) :: EPSB(0:FS%NL+1) ! longwave emisivity, back ! to include indoor and outdoor
REAL(r64) :: TAU(0:FS%NL+1) ! longwave transmittance ! nodes - more general
REAL(r64) :: RTOT ! total resistance from TAE_OUT to TAE_IN [m2K/W]
REAL(r64) :: HC2D(6, 6) ! convective heat transfer coefficients between layers i and j
REAL(r64) :: HR2D(6, 6) ! radiant heat transfer coefficients between layers i and j
REAL(r64) :: HCIout(6), HRIout(6) ! convective and radiant heat transfer coefficients between
! layer i and outdoor air or mean radiant temperature, resp.
REAL(r64) :: HCIin(6), HRIin(6) ! convective and radiant heat transfer coefficients between
! layer i and indoor air or mean radiant temperature, resp.
REAL(r64) :: HCinout, HRinout ! convective and radiant heat transfer coefficients between
! indoor and outdoor air or mean radiant temperatures
! (almost always zero)
! Indoor side convection coefficients - used for Open Channel Flow on indoor side
REAL(r64) :: HFS ! nominal height of fen system (assumed 1 m)
REAL(r64) :: TOC_EFF ! effective thickness of open channel, m
Real(r64) :: ConvF ! convection factor: accounts for enhanced convection
! for e.g. VB adjacent to open channel
Real(r64) :: HC_GA ! convection - glass to air
Real(r64) :: HC_SA ! convection - shade (both sides) to air
Real(r64) :: HC_GS ! convection - glass to shade (one side)
REAL(r64) :: TINdv, TOUTdv ! dummy variables used
REAL(r64) :: TRMINdv, TRMOUTdv ! for boundary conditions in calculating
REAL(r64) :: SOURCEdv( FS%NL+1) ! indices of merit
REAL(r64) :: SUMERR ! error summation used to check validity of code/model
REAL(r64) :: QGAIN ! total gain to conditioned space [[W/m2]
REAL(r64) :: SaveHCNLm ! place to save HC(NL-1) - two resistance networks differ
REAL(r64) :: SaveHCNL ! place to save HC(NL) - two resistance networks differ
! in their definitions of these heat transfer coefficients
LOGICAL DoPrint ! set true to print debugging info
! Flow
ASHWAT_Thermal = .FALSE. ! init to failure
NL = FS%NL ! working copy
IF (NL < 1) RETURN
HCOCFout = HCOUT ! outdoor side
IM_ON = 1
HHAT=0.0d0
HC = 0.0d0
HR = 0.0d0
T=0.0d0
TNEW=0.0d0
EB=0.0d0
JF=0.0d0
JB=0.0d0
Q=0.0d0
QOCF_F=0.0d0
QOCF_B=0.0d0
QOCF = 0.0d0
QOCFRoom = 0.0d0
QNET=0.0d0
QGAIN=0.0d0
TAU = 0.0d0
RHOF = 0.0d0
RHOB = 0.0d0
EPSF = 0.0d0
EPSB = 0.0d0
HC_GA = 0.0d0
HC_SA = 0.0d0
HC_GS = 0.0d0
ITRY=0
EB( 0) = StefanBoltzmann * TOUT**4
EB( NL+1) = StefanBoltzmann * TIN**4
ADIM=3*NL + 2 ! DIMENSION OF A-MATRIX
! organize longwave radiant properties - book-keeping
TAU_ROOM = 0.0d0 ! must always be zero
RHOF_ROOM = 0.0d0 ! almost always zero
EPSF_ROOM = 1.d0- TAU_ROOM - RHOF_ROOM ! almost always unity
RHOF(NL+1) = RHOF_ROOM
EPSF(NL+1) = EPSF_ROOM
TAU (NL+1) = TAU_ROOM
DO I=1,NL
EPSF(I) = FS%L(I)%LWP_EL%EPSLF
EPSB(I) = FS%L(I)%LWP_EL%EPSLB
TAU (I) = FS%L(I)%LWP_EL%TAUL
RHOF(I) = 1.d0 - FS%L(I)%LWP_EL%EPSLF - FS%L(I)%LWP_EL%TAUL
RHOB(I) = 1.d0 - FS%L(I)%LWP_EL%EPSLB - FS%L(I)%LWP_EL%TAUL
END DO
TAU_OUT = 0.0d0 ! must always be zero
RHOB_OUT = 0.0d0 ! DON'T CHANGE
EPSB_OUT = 1.d0 - TAU_OUT - RHOB_OUT ! should always be unity
TAU (0) = TAU_OUT
EPSB(0) = EPSB_OUT
RHOB(0) = RHOB_OUT
! Later could add RHOF_ROOM to the parameter list
! Relaxation needed to keep solver stable if OCF is present
ALPHA=1.0d0
IF (NL.GE.2) THEN
IF (FS%G(NL-1)%GTYPE .EQ. gtyOPENin) ALPHA=0.5d0
IF (FS%G( 1) %GTYPE .EQ. gtyOPENout) ALPHA=0.1d0
ENDIF
! FIRST ESTIMATE OF GLAZING TEMPERATURES AND BLACK EMISSIVE POWERS
DO I=1,NL
T(I) = TOUT + (REAL(I,r64))/(REAL(NL+1,r64)) * (TIN-TOUT)
EB(I)= StefanBoltzmann * T(I)**4
END DO
CONVRG=0
! Start the solver
! ITERATION RE-ENTRY
DO ITRY=1,100
! CALCULATE GAS LAYER CONVECTIVE HEAT TRANSFER COEFFICIENTS
hin_scheme = 3 ! different schemes for calculating convection
! coefficients glass-to-air and shade-to-air
! if open channel air flow is allowed
! see the corresponding subroutines for detail
! = 1 gives dependence of height, spacing, delta-T
! = 2 gives dependence of spacing, delta-T but
! returns unrealistic values for large spacing
! = 3 glass-shade spacing dependence only on HCIN
! = negative, applies HCIN without adjusting for
! temperature, height, spacing, slat angle
! Recommended -> hin_scheme=3 for use with HBX,
! simplicity, right trends wrt spacing
! start by assuming no open channel flow on indoor side
HC(NL)=HCIN ! default - HC(NL) supplied by calling routine
! use this for HBX
! or
! trigger calculation of HC(NL) using ASHRAE correlation
! HC(NL) = HIC_ASHRAE(1.0d0, T(NL), TIN) ! h - flat plate
!
! Add by BAN June 2013 for standard ratings U-value and SHGC calc only
IF (PRESENT(HCInFlag)) THEN
IF ( HCInFlag) HC(NL) = HCInWindowStandardRatings(Height, T(NL), TIN)
ENDIF
HC(0) = HCOUT ! HC(0) supplied by calling routine as HCOUT
! Check for open channels - only possible with at least two layers
IF (NL >= 2) THEN
DO I=1,NL-1 ! Scan gaps between layers
! DEAL WITH INDOOR OPEN CHANNEL FLOW HERE
IF ((I .EQ. NL-1).AND.(FS%G( I)%GTYPE .EQ. gtyOPENin)) THEN
!TOC_EFF = FS%G( I)%TAS_EFF / 1000. ! effective thickness of OC gap, m
TOC_EFF = FS%G( I)%TAS_EFF ! effective thickness of OC gap, m Modified by BAN May 9, 2013
HFS = 1. ! nominal height of system (m)
! convection - glass to air
CALL GLtoAMB( TOC_EFF, HFS, T( NL-1), TIN, HCIN, HC_GA, hin_scheme)
! CALL GLtoAMB( 1.0, HFS, T( NL-1), TIN, HCIN, HC_GA, hin_scheme)
! ^ VERY WIDE GAP
! convection - shade (both sides) to air
ConvF = ConvectionFactor( FS%L( I+1))
HC_SA = ConvF * SLtoAMB( TOC_EFF, HFS, T(NL), TIN, HCIN, hin_scheme)
! HC_SA = ConvF * SLtoAMB( 1.0, HFS, T(NL), TIN, HCIN, hin_scheme)
! ^ VERY WIDE GAP
! convection - glass to shade (one side)
CALL SLtoGL( TOC_EFF, T(NL), T(NL-1), HC_GS, 1)
! CALL SLtoGL( 1.0, T(NL), T(NL-1), HC_GS, 2) ! REMOVE LATER
! ^ VERY WIDE GAP, should return near zero
! Don't use hin_scheme as last parameter - set manually
! 1 = Conduction, 2 = slight Ra penalty
! Set negative for default HC_GS=0
! Recommended: 2
HC(NL-1) = HC_GS
HC(NL ) = HCIN * ConvF
QOCF_B(NL-1) = (TIN - T(NL-1))*HC_GA
QOCF_F(NL ) = (TIN - T(NL ))*(HC_SA - HC( NL))
QOCFRoom = - QOCF_B( NL-1) - QOCF_F( NL)
! end of gap open to indoor side
ELSE IF ((I .EQ. 1).AND.(FS%G( I)%GTYPE .EQ. gtyOPENout)) THEN
! outdoor open channel
QOCF_B(1) = (TOUT - T(1))*HCOCFout
QOCF_F(2) = (TOUT - T(2))*HCOCFout
HC(1)=0.0d0
HC(0)=HCOCFout
ELSE
! normal gap
HC(I) = HConvGap( FS%G(I), T(I), T(I+1))
ENDIF
END DO ! end scan through gaps
! total OCF gain to each layer
QOCF = QOCF_F + QOCF_B
ENDIF ! end IF (NL .GE. 2)
! CONVERT TEMPERATURE POTENTIAL CONVECTIVE COEFFICIENTS to
! BLACK EMISSIVE POWER POTENTIAL CONVECTIVE COEFFICIENTS
HHAT(0) = HC(0)*(1.d0/StefanBoltzmann)/( ((TOUT**2+T(1)**2)) * ((TOUT+T(1))) )
DO I=1,NL-1 ! Scan the cavities
HHAT(I) = HC(I)*(1.d0/StefanBoltzmann)/( ((T(I)**2+T(I+1)**2)) * ((T(I)+T(I+1))) )
END DO
HHAT(NL) = HC(NL)*(1.d0/StefanBoltzmann)/( ((T(NL)**2+TIN**2)) * ((T(NL)+TIN)) )
! SET UP MATRIX
XSOL = 0.0d0
A = 0.0d0
L=1
A(L,1)=1.d0
A(L,2)= -1.d0 * RHOB(0) ! -1.0 * RHOB_OUT
A(L,ADIM+1)= EPSB_OUT * StefanBoltzmann * TRMOUT**4
DO I=1,NL
L=3*I-1
A(L,3*I-2) = RHOF( I)
A(L,3*I-1) = -1.0d0
A(L,3*I ) = EPSF(I) ! LWP( I)%EPSLF
A(L,3*I+2) = TAU (I) ! LWP( I)%TAUL
A(L,ADIM+1) = 0.0d0
L=3*I
IF (NL == 1) THEN
A(L,1)= 1.0d0 ! Single layer
A(L,2)=-1.0d0
A(L,3)=-1.0d0*(HHAT(0)+HHAT(1))
A(L,4)=-1.0d0
A(L,5)= 1.0d0
A(L,ADIM+1)= -1.0d0*(HHAT(0)*EB( 0) + HHAT(1)*EB(2) + &
SOURCE(1)+QOCF(1) )
ELSE IF (I == 1 ) THEN
A(L,1)= 1.0d0 ! Outdoor layer
A(L,2)= -1.0d0
A(L,3)= -1.0d0*(HHAT(0) + HHAT(1))
A(L,4)= -1.0d0
A(L,5)= 1.0d0
A(L,6)= HHAT(1)
A(L,ADIM+1)=-1.0d0*( HHAT(0)*EB( 0) + SOURCE(1) + QOCF(1) )
ELSE IF (I == NL) THEN
A(L,3*NL-3)= HHAT(NL-1) ! Indoor layer
A(L,3*NL-2)= 1.0d0
A(L,3*NL-1)=-1.0d0
A(L,3*NL )=-1.0d0*(HHAT(NL)+HHAT(NL-1))
A(L,3*NL+1)=-1.0d0
A(L,3*NL+2)= 1.0d0
A(L,ADIM+1)=-1.0d0*( HHAT(NL)*EB( NL+1) + SOURCE(NL)+ QOCF(NL) )
ELSE
A(L,3*I-3) = HHAT(I-1)
A(L,3*I-2) = 1.0d0
A(L,3*I-1) =-1.0d0
A(L,3*I ) =-1.0d0*(HHAT(I)+HHAT(I-1))
A(L,3*I+1) =-1.0d0
A(L,3*I+2) = 1.0d0
A(L,3*I+3) = HHAT(I)
A(L,ADIM+1) =-1.0d0*(SOURCE(I)+QOCF(I))
END IF
L=3*I+1
A(L,3*I-2) = TAU(I) ! LWP( I)%TAUL
A(L,3*I ) = EPSB(I) ! LWP( I)%EPSLB
A(L,3*I+1) = -1.0d0
A(L,3*I+2) = RHOB(I)
A(L,ADIM+1) = 0.0d0
END DO
L=3*NL+2
A(L,3*NL+1)= - 1.d0 * RHOF(NL+1) ! - 1.0 * RHOF_ROOM
A(L,3*NL+2)= 1.0d0
A(L,ADIM+1)= EPSF_ROOM * StefanBoltzmann * TRMIN**4
! SOLVE MATRIX
! Call SOLMATS for single precision matrix solution
CALL SOLMATS( ADIM, A, XSOL)
! UNPACK SOLUTION VECTOR AND RECORD LARGEST TEMPERATURE CHANGE
JB(0)=XSOL(1)
MAXERR=0.0d0
DO I=1,NL
J=3*I-1
JF(I)=XSOL(J)
J=J+1
EB(I)=MAX( 1.d0, XSOL(J)) ! prevent impossible temps
TNEW(I)=((EB(I)/StefanBoltzmann)**0.25d0)
J=J+1
JB(I)=XSOL(J)
MAXERR = MAX( MAXERR, ABS( TNEW(I)-T(I))/TNEW(I))
END DO
JF(NL+1)=XSOL(ADIM)
! CALCULATE HEAT FLUX AT EACH GAP, Q
DO I=0,NL ! Loop gaps (including inside and outside
Q(I)=JF(I+1)-JB(I) + HHAT(I)*(EB(I+1)-EB(I))
END DO
! A CHECK - NET HEAT FLUX INTO ANY LAYER, AT STEADY-STATE,
! SHOULD BE ZERO
DO I=1,NL
QNET(I)=SOURCE(I) + QOCF(I) + Q(I) - Q(I-1)
END DO
! UPDATE GLAZING TEMPERATURES AND BLACK EMISSIVE POWERS
DO I=1,NL
T(I) = T(I) + ALPHA * (TNEW(I)-T(I))
EB(I)= StefanBoltzmann * T(I)**4
END DO
! CHECK FOR CONVERGENCE
IF (CONVRG > 0) EXIT
IF (MAXERR < TOL) CONVRG=CONVRG+1
END DO ! main iteration
IF (CONVRG == 0) THEN
CALL ShowSevereError(RoutineName//'Net radiation analysis did not converge for '//TRIM( FS%Name))
CALL ShowContinueError('...Maximum error is = '//TRIM(TrimSigDigits(MAXERR,6)))
CALL ShowContinueError('...Convergence tolerance is = '//TRIM(TrimSigDigits(TOL,6)))
END IF
! NOTE: HC_SA, HC_GA and HC_SG are only available if there is
! an open channel on the indoor side and the calculation of
! these coefficients was triggered earlier
QGAIN = SOURCE(NL+1) + HC(NL)*(T(NL)-TIN) + JB(NL) - JF(NL+1)
! Modified by BAN May 3, 2013 to avoid zero layer index
IF ( NL .GE. 2 ) THEN
IF (FS%G(NL-1)%GTYPE .EQ. gtyOPENin) THEN
QGAIN = SOURCE(NL+1) + (HC_SA/2.0d0)*(T(NL)-TIN) + JB(NL) - JF(NL+1)
QGAIN = QGAIN + HC_GA * (T(NL-1) - TIN) &
+ (HC_SA/2.0d0) * (T(NL ) - TIN)
ENDIF
END IF
ASHWAT_Thermal = .TRUE.
! New code follows from here - for calculating Ucg and SHGC
! NOTE: This code can be bypassed if
! indices of merit are not needed
IF (IM_ON .NE. 1) RETURN
! Initialize various things
HR = 0.0d0
HJC = 0.0d0
HJR = 0.0d0
TAE_OUT = 0.0d0
TAE_IN = 0.0d0
FHR_OUT = 0.0d0
FHR_IN = 0.0d0
Q_IN=0.0d0
RTOT=0.0d0
UCG=0.0d0
SHGC=0.0d0
Rvalue = 0.0d0
HC2D=0.0d0
HR2D=0.0d0
HCIout=0.0d0
HRIout=0.0d0
HCIin=0.0d0
HRIin=0.0d0
HCinout=0.0d0
HRinout=0.0d0
TNEW=0.0d0
TINdv=0.0d0
TOUTdv=0.0d0
TRMINdv=0.0d0
TRMOUTdv=0.0d0
SOURCEdv=0.0d0
! Identify the diathermanous layers
ISDL = 0
DO I=1,NL
IF(FS%L( I)%LWP_EL%TAUL .GT. 0.001d0) ISDL(I)=1 ! layer is diathermanous
! of tau_lw > 0.001 (ie 0.1%)
! Note: ISDL(0) and ISDL(NL+1)
! must both be zero
END DO ! end loop to calculate ISDL(i)
! determine the largest number of consecutive diathermanous layers, NDLIAR
! i.e., the number of diathermanous layers in a row
NDLIAR = 0
IDV = 0
DO I=1,NL
IF(ISDL(I) .EQ. 1) THEN
IDV=IDV+1
ELSE
IDV=0
ENDIF
IF(IDV.GT. NDLIAR) NDLIAR=IDV
END DO ! end loop to calculate NDLIAR
IF(NDLIAR .GT. 1)RETURN ! cannot handle two (or more) consecutive
! diathermanous layers, U/SHGC calculation
! will be skipped entirely
! CHANGE TO (NDLIAR .GT. 2) ONCE
! SUBROUTINE DL2_RES IS AVAILABLE
! calculate radiant heat transfer coefficents between adjacent opaque
! layers
DO I=0,NL ! scan through all gaps - including indoor/outdoor
IF ((ISDL(I) .EQ. 0) .AND. (ISDL(I+1) .EQ. 0)) THEN
IF (I .EQ.0) THEN ! outdoor side
HR(I) = HRadPar( T(1), TRMOUT, EPSF(1), EPSB(0))
ELSE IF (I .EQ. NL) THEN ! indoor side
HR(I) = HRadPar( T(NL), TRMIN,EPSF(NL+1),EPSB(NL))
ELSE ! cavities
HR(I) = HRadPar( T(I), T(I+1), EPSF(I+1), EPSB(I))
ENDIF
ENDIF
END DO ! end loop through gaps
! calculate radiant heat transfer coefficents at single diathermanous
! layers,three coefficients in each case
DO I=0,NL-1 ! scan through all layers - look for single DL
! layers between two opaque layers
IF((ISDL(I) .EQ. 0) .AND. (ISDL(I+1) .EQ. 1) &
.AND. (ISDL(I+2) .EQ. 0)) THEN
IF (I .EQ. 0) THEN ! outdoor layer is diathermanous
IF(NL .EQ. 1) THEN
CALL DL_RES_r2 (TRMOUT, T(1), TRMIN, &
RHOB(0),RHOF(1),RHOB(1),TAU(1),RHOF(2),&
HJR(1), HR(0), HR(1))
ELSE
CALL DL_RES_r2 (TRMOUT, T(1), T(2), &
RHOB(0),RHOF(1),RHOB(1),TAU(1),RHOF(2),&
HJR(1), HR(0), HR(1))
ENDIF
ELSE ! with IF (I .EQ. 0) i.e., i != 0
IF(I .EQ. NL-1) THEN ! indoor layer is diathermanous
CALL DL_RES_r2 (T(NL-1), T(NL), TRMIN, &
RHOB(NL-1),RHOF(NL),RHOB(NL),TAU(NL),RHOF(NL+1),&
HJR(NL), HR(NL-1), HR(NL))
ELSE ! some intermediate layer is diathermanous
CALL DL_RES_r2 (T(I), T(I+1), T(I+2), &
RHOB(I),RHOF(I+1),RHOB(I+1),TAU(I+1),RHOF(I+2),&
HJR(I+1), HR(I), HR(I+1))
ENDIF ! end of IF/ELSE (I .EQ. NL-1)
ENDIF ! end of IF/ELSE (I .EQ. 0)
ENDIF ! end of IF(ISDL(I) .EQ. 0) .AND. .....
END DO ! end of scan through all layers
! calculate radiant heat transfer coefficents at double diathermanous
! layers,six coefficients in each case
! THIS SECTION NOT ACTIVE YET
IF (NL .GE. 2) THEN
DO I=0,NL-2 ! scan through all layers - look for double DL
! layers between two opaque layers
IF((ISDL(I) .EQ. 0) .AND. (ISDL(I+1) .EQ. 1) &
.AND. (ISDL(I+2) .EQ. 1) &
.AND. (ISDL(I+3) .EQ. 0)) THEN
IF (I .EQ. 0) THEN
! CALL DL2_RES (TRMOUT, T(1), T(2), T(3) etc)
ELSE
IF(I .EQ. NL-2) THEN
! CALL DL2_RES (T(NL-2), T(NL-1), T(NL), TRMIN, etc)
ELSE
! CALL DL2_RES (T(I), T(I+1), T(I+2), T(I+3) etc)
ENDIF ! end of IF/ELSE (I .EQ. NL-1)
ENDIF ! end of IF/ELSE (I .EQ. 0)
ENDIF ! end of IF(ISDL(I) .EQ. 0) .AND. .....
END DO ! end of scan through all layers
ENDIF
! calculate convective OCF/jump heat transfer coefficients
IF(NL .GE. 2) THEN ! no OCF unless at least two layers exist
! It is not possible for both of the following cases to be
! true for the same gap (i.e., for NL=2)
IF (FS%G(NL-1)%GTYPE .EQ. gtyOPENin) THEN
SaveHCNLm = HC(NL-1)
SaveHCNL = HC(NL)
HC(NL-1) = HC_GS
HC(NL) = HC_SA
HJC(NL) = HC_GA
ENDIF
HC(0) = HCOUT
IF (FS%G(1)%GTYPE .EQ. gtyOPENout) THEN
HC(0) = HCOUT + HCOCFout
HJC(1) = HCOCFout
ENDIF
ENDIF
! copy convective heat transfer coefficients to 2D arrays
! adjacent layers
IB=1
IE=NL-1
IF (IB .LE. IE) THEN
DO I=IB,IE
HC2D(I,I+1) = HC(I)
HC2D(I+1,I) = HC2D(I,I+1)
END DO
ENDIF
! jumpers
IB=2
IE=NL-1
IF (IB .LE. IE) THEN
DO I=IB,IE
HC2D(I-1,I+1) = HJC(I)
HC2D(I+1,I-1) = HC2D(I-1,I+1)
END DO
ENDIF
! double jumpers - NOT ACTIVE YET
IB=2
IE=NL-2
IF (IB .LE. IE) THEN
DO I=IB,IE
! HC2D(I-1,I+2) = H2JC(I)
! HC2D(I+2,I-1) = HC2D(I-1,I+2)
END DO
ENDIF
! outdoor side
HCIout(1)=HC(0)
IF(NL .GE. 2) HCIout(2)=HJC(1)
! indoor side
HCIin(NL) = HC(NL)
IF(NL .GE. 2) HCIin(NL-1)=HJC(NL)
! special case - indoor-to-outdoor convection (?)
HCinout=0.0d0
! copy radiative heat transfer coefficients to 2D arrays
! adjacent layers
IB=1
IE=NL-1
IF (IB .LE. IE) THEN
DO I=IB,IE
HR2D(I,I+1) = HR(I)
HR2D(I+1,I) = HR2D(I,I+1)
END DO
ENDIF
! jumpers
IB=2
IE=NL-1
IF (IB .LE. IE) THEN
DO I=IB,IE
HR2D(I-1,I+1) = HJR(I)
HR2D(I+1,I-1) = HR2D(I-1,I+1)
END DO
ENDIF
! double jumpers
IB=2
IE=NL-2
IF (IB .LE. IE) THEN
DO I=IB,IE
! HR2D(I-1,I+2) = H2JR(I)
! HR2D(I+2,I-1) = HR2D(I-1,I+2)
END DO
ENDIF
! outdoor side
HRIout(1)=HR(0)
IF(NL .GE. 2) HRIout(2)=HJR(1)
! indoor side
HRIin(NL)=HR(NL)
IF(NL .GE. 2) HRIin(NL-1)=HJR(NL)
! special case - indoor-to-outdoor radiation
IF(NL .EQ. 1) HRinout=HJR(1)
! IF(NL .EQ. 2) HRinout=H2JR(1)
! CONFIRM VALIDITY OF CODE
IF (1 .EQ. 0) THEN ! was used for debugging - successfully
! and can now be bypassed
! - code in this section generates the
! same solution of temperatures (TNEW(i))
! that was found by the net radiation
! solver above (T(i))
ADIM = NL
A=0.0d0
XSOL=0.0d0
TOUTdv =TOUT ! solution for TNEW should
TRMOUTdv =TRMOUT ! match existing solution
TINdv =TIN ! for T
TRMINdv =TRMIN
SOURCEdv =SOURCE
DO I=1,NL
A(I,ADIM+1) = HCIout(I)*TOUTdv + HRIout(I)*TRMOUTdv + &
HCIin(I) *TINdv + HRIin(I) *TRMINdv + &
SOURCEdv(I)
A(I,I) = HCIout(I) + HRIout(I) + HCIin(I) + HRIin(I)
DO J=1,NL
IF(J .NE. I) THEN
A(I,I) = A(I,I) + HC2D(I,J) + HR2D(I,J)
A(I,J) = -1.0d0 * (HC2D(I,J) + HR2D(I,J))
END IF
END DO
END DO
! SOLVE MATRIX
! Call SOLMATS for single precision matrix solution
CALL SOLMATS( ADIM, A, XSOL)
! UNPACK SOLUTION VECTOR
SUMERR=0.0d0
DO I=1,NL
TNEW(I)=XSOL(I)
SUMERR=SUMERR+ABS(TNEW(I)-T(I))
END DO
END IF ! end (1 .EQ. 0) code disabled
! calculate U-factor
ADIM = NL
A=0.0d0
XSOL=0.0d0
TNEW =0.0d0
TOUTdv =1.0d0
TRMOUTdv =1.0d0
TINdv =0.0d0
TRMINdv =0.0d0
SOURCEdv =0.0d0
DO I=1,NL
A(I,ADIM+1) = HCIout(I)*TOUTdv + HRIout(I)*TRMOUTdv + &
HCIin(I) *TINdv + HRIin(I) *TRMINdv + &
SOURCEdv(I)
A(I,I) = HCIout(I) + HRIout(I) + HCIin(I) + HRIin(I)
DO J=1,NL
IF(J .NE. I) THEN
A(I,I) = A(I,I) + HC2D(I,J) + HR2D(I,J)
A(I,J) = -1.d0 * (HC2D(I,J) + HR2D(I,J))
END IF
END DO
END DO
! SOLVE MATRIX
! Call SOLMATS for single precision matrix solution
CALL SOLMATS( ADIM, A, XSOL)
! UNPACK SOLUTION VECTOR
DO I=1,NL
TNEW(I)=XSOL(I)
END DO
Q_IN = HCinout* (TOUTdv-TINdv) + HRinout * &
(TRMOUTdv-TRMINdv)
DO I=1,NL
Q_IN = Q_IN + HCIin(I)*(TNEW(I)-TINdv ) &
+ HRIin(I)*(TNEW(I)-TRMINdv)
END DO
Q_IN = Q_IN + SOURCEdv(NL+1) ! this line not needed
UCG = Q_IN
Rvalue = 5.678d0/UCG
! calculate SHGC
SHGC = 0.0d0
IF(ABS(ISOL) .GT. 0.01d0) THEN
ADIM = NL
A=0.0d0
XSOL =0.0d0
TNEW =0.0d0
TOUTdv =0.0d0
TRMOUTdv =0.0d0
TINdv =0.0d0
TRMINdv =0.0d0
DO I=1,NL+1
SOURCEdv(I) = SOURCE(I)
END DO
DO I=1,NL
A(I,ADIM+1) = HCIout(I)*TOUTdv + HRIout(I)*TRMOUTdv + &
HCIin(I) *TINdv + HRIin(I) *TRMINdv + &
SOURCEdv(I)
A(I,I) = HCIout(I) + HRIout(I) + HCIin(I) + HRIin(I)
DO J=1,NL
IF(J .NE. I) THEN
A(I,I) = A(I,I) + HC2D(I,J) + HR2D(I,J)
A(I,J) = -1.0d0 * (HC2D(I,J) + HR2D(I,J))
END IF
END DO
END DO
! SOLVE MATRIX
! Call SOLMATS for single precision matrix solution
CALL SOLMATS( ADIM, A, XSOL)
! UNPACK SOLUTION VECTOR
DO I=1,NL
TNEW(I)=XSOL(I)
END DO
Q_IN = HCinout* (TOUTdv-TINdv) + HRinout * &
(TRMOUTdv-TRMINdv)
DO I=1,NL
Q_IN = Q_IN + HCIin(I)*(TNEW(I)-TINdv ) &
+ HRIin(I)*(TNEW(I)-TRMINdv)
END DO
Q_IN = Q_IN + SOURCEdv(NL+1)
SHGC = Q_IN/ISOL ! only executed if ISOL > 0.01 [W/m2]
ENDIF ! end if (ABS(ISOL) .GT. 0.01)
! calculate FHR_OUT
ADIM = NL
A=0.0d0
XSOL =0.0d0
TNEW =0.0d0
TOUTdv =1.0d0
TRMOUTdv =0.0d0
TINdv =0.0d0
TRMINdv =0.0d0
SOURCEdv =0.0d0
DO I=1,NL
A(I,ADIM+1) = HCIout(I)*TOUTdv + HRIout(I)*TRMOUTdv + &
HCIin(I) *TINdv + HRIin(I) *TRMINdv + &
SOURCEdv(I)
A(I,I) = HCIout(I) + HRIout(I) + HCIin(I) + HRIin(I)
DO J=1,NL
IF(J .NE. I) THEN
A(I,I) = A(I,I) + HC2D(I,J) + HR2D(I,J)
A(I,J) = -1.d0 * (HC2D(I,J) + HR2D(I,J))
END IF
END DO
END DO
! SOLVE MATRIX
! Call SOLMATS for single precision matrix solution
CALL SOLMATS( ADIM, A, XSOL)
! UNPACK SOLUTION VECTOR
DO I=1,NL
TNEW(I)=XSOL(I)
END DO
Q_IN = HCinout* (TOUTdv-TINdv) + HRinout * &
(TRMOUTdv-TRMINdv)
DO I=1,NL
Q_IN = Q_IN + HCIin(I)*(TNEW(I)-TINdv ) &
+ HRIin(I)*(TNEW(I)-TRMINdv)
END DO
Q_IN = Q_IN + SOURCEdv(NL+1)
FHR_OUT = 1.0d0 - (Q_IN/UCG)
TAE_OUT = FHR_OUT*TRMOUT + (1.d0 - FHR_OUT)*TOUT
! calculate FHR_IN
ADIM = NL
A=0.0d0
XSOL =0.0d0
TNEW =0.0d0
TOUTdv =0.0d0
TRMOUTdv =0.0d0
TINdv =1.0d0
TRMINdv =0.0d0
SOURCEdv =0.0d0
DO I=1,NL
A(I,ADIM+1) = HCIout(I)*TOUTdv + HRIout(I)*TRMOUTdv + &
HCIin(I) *TINdv + HRIin(I) *TRMINdv + &
SOURCEdv(I)
A(I,I) = HCIout(I) + HRIout(I) + HCIin(I) + HRIin(I)
DO J=1,NL
IF(J .NE. I) THEN
A(I,I) = A(I,I) + HC2D(I,J) + HR2D(I,J)
A(I,J) = -1.d0 * (HC2D(I,J) + HR2D(I,J))
END IF
END DO
END DO
! SOLVE MATRIX
! Call SOLMATS for single precision matrix solution
CALL SOLMATS( ADIM, A, XSOL)
! UNPACK SOLUTION VECTOR
DO I=1,NL
TNEW(I)=XSOL(I)
END DO
Q_IN = HCinout* (TOUTdv-TINdv) + HRinout * &
(TRMOUTdv-TRMINdv)
DO I=1,NL
Q_IN = Q_IN + HCIin(I)*(TNEW(I)-TINdv ) &
+ HRIin(I)*(TNEW(I)-TRMINdv)
END DO
Q_IN = Q_IN + SOURCEdv(NL+1)
FHR_IN = 1.0d0 + (Q_IN/UCG)
TAE_IN = FHR_IN*TRMIN + (1.0d0 - FHR_IN)*TIN
! double check heat gain to room
! Q_IN calculated this way should be equal to QGAIN calculated
! above with raw results from the net radiation solution
! The difference between the two is printed below
! Both include the directly transmitted solar gain
Q_IN = UCG*(TAE_OUT - TAE_IN) + SHGC*ISOL
! End of new code - for calculating Ucg and SHGC
! restore convective heat transfer coefficients if alterred earlier
! for more general resistor network - otherwise mainline will
! receive faulty data
IF(NL .GE. 2) THEN ! no OCF unless at least two layers exist
IF (FS%G(NL-1)%GTYPE .EQ. gtyOPENin) THEN
HC(NL-1) = SaveHCNLm
HC(NL) = SaveHCNL
ENDIF
ENDIF
RETURN
END FUNCTION ASHWAT_Thermal