!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! MAIN ITERATION LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | nlayer | |||
| integer, | intent(in) | :: | iwd | |||
| real(kind=r64), | intent(inout) | :: | tout | |||
| real(kind=r64), | intent(inout) | :: | tind | |||
| real(kind=r64), | intent(in) | :: | wso | |||
| real(kind=r64), | intent(in) | :: | wsi | |||
| real(kind=r64), | intent(in) | :: | VacuumPressure | |||
| real(kind=r64), | intent(in) | :: | VacuumMaxGapThickness | |||
| real(kind=r64), | intent(in) | :: | dir | |||
| real(kind=r64), | intent(inout) | :: | Ebsky | |||
| real(kind=r64), | intent(in) | :: | Gout | |||
| real(kind=r64), | intent(in) | :: | Trmout | |||
| real(kind=r64), | intent(in) | :: | Trmin | |||
| real(kind=r64), | intent(inout) | :: | Ebroom | |||
| real(kind=r64), | intent(in) | :: | Gin | |||
| real(kind=r64), | intent(in), | dimension(maxlay2) | :: | tir | ||
| real(kind=r64), | intent(in), | dimension(maxlay2) | :: | rir | ||
| real(kind=r64), | intent(in), | dimension(maxlay2) | :: | emis | ||
| real(kind=r64), | intent(in), | dimension(MaxGap) | :: | gap | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | thick | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | scon | ||
| real(kind=r64), | intent(in) | :: | tilt | |||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | asol | ||
| real(kind=r64), | intent(in) | :: | height | |||
| real(kind=r64), | intent(in) | :: | heightt | |||
| real(kind=r64) | :: | width | ||||
| integer, | intent(in), | dimension(maxlay1, maxgas) | :: | iprop | ||
| real(kind=r64), | intent(in), | dimension(maxlay1, maxgas) | :: | frct | ||
| real(kind=r64), | intent(in), | dimension(maxlay1) | :: | presure | ||
| integer, | intent(in), | dimension(maxlay1) | :: | nmix | ||
| real(kind=r64), | intent(in), | dimension(maxgas) | :: | wght | ||
| real(kind=r64), | intent(in), | dimension(maxgas, 3) | :: | gcon | ||
| real(kind=r64), | intent(in), | dimension(maxgas, 3) | :: | gvis | ||
| real(kind=r64), | intent(in), | dimension(maxgas, 3) | :: | gcp | ||
| real(kind=r64), | intent(in), | dimension(maxgas) | :: | gama | ||
| integer, | intent(in), | dimension(maxlay) | :: | SupportPillar | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | PillarSpacing | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | PillarRadius | ||
| real(kind=r64), | intent(inout), | dimension(maxlay2) | :: | theta | ||
| real(kind=r64), | intent(out), | dimension(maxlay3) | :: | q | ||
| real(kind=r64), | intent(out), | dimension(maxlay1) | :: | qv | ||
| real(kind=r64), | intent(out) | :: | flux | |||
| real(kind=r64), | intent(out) | :: | hcin | |||
| real(kind=r64), | intent(inout) | :: | hrin | |||
| real(kind=r64), | intent(inout) | :: | hcout | |||
| real(kind=r64), | intent(inout) | :: | hrout | |||
| real(kind=r64), | intent(inout) | :: | hin | |||
| real(kind=r64), | intent(inout) | :: | hout | |||
| real(kind=r64), | intent(out), | dimension(maxlay1) | :: | hcgas | ||
| real(kind=r64), | intent(out), | dimension(maxlay1) | :: | hrgas | ||
| real(kind=r64), | intent(out) | :: | ufactor | |||
| integer, | intent(inout) | :: | nperr | |||
| character(len=*) | :: | ErrorMessage | ||||
| real(kind=r64), | intent(inout) | :: | Tamb | |||
| real(kind=r64), | intent(inout) | :: | Troom | |||
| integer, | intent(in), | dimension(2) | :: | ibc | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | Atop | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | Abot | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | Al | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | Ar | ||
| real(kind=r64), | intent(in), | dimension(maxlay) | :: | Ah | ||
| real(kind=r64), | intent(in), | dimension(maxlay1) | :: | vvent | ||
| real(kind=r64), | intent(in), | dimension(maxlay1) | :: | tvent | ||
| integer, | intent(in), | dimension(maxlay) | :: | LayerType | ||
| real(kind=r64), | intent(out), | dimension(maxlay) | :: | Ra | ||
| real(kind=r64), | intent(out), | dimension(maxlay) | :: | Nu | ||
| real(kind=r64), | intent(out), | dimension(maxlay1) | :: | vfreevent | ||
| real(kind=r64), | intent(out), | dimension(maxlay1) | :: | qcgas | ||
| real(kind=r64), | intent(out), | dimension(maxlay1) | :: | qrgas | ||
| real(kind=r64), | intent(inout), | dimension(maxlay) | :: | Ebf | ||
| real(kind=r64), | intent(inout), | dimension(maxlay) | :: | Ebb | ||
| real(kind=r64), | intent(inout), | dimension(maxlay) | :: | Rf | ||
| real(kind=r64), | intent(inout), | dimension(maxlay) | :: | Rb | ||
| real(kind=r64), | intent(out) | :: | ShadeEmisRatioOut | |||
| real(kind=r64), | intent(out) | :: | ShadeEmisRatioIn | |||
| real(kind=r64), | intent(out) | :: | ShadeHcModifiedOut | |||
| real(kind=r64), | intent(out) | :: | ShadeHcModifiedIn | |||
| integer, | intent(in) | :: | ThermalMod | |||
| integer, | intent(in) | :: | Debug_mode | |||
| real(kind=r64), | intent(out) | :: | AchievedErrorTolerance | |||
| integer, | intent(out) | :: | TotalIndex |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
subroutine therm1d(nlayer, iwd, tout, tind, wso, wsi, VacuumPressure, VacuumMaxGapThickness, &
dir, Ebsky, Gout, Trmout, Trmin, Ebroom, Gin, tir, rir, emis, &
gap, thick, scon, tilt, asol, height, heightt, width, &
iprop, frct, presure, nmix, wght, gcon, gvis, gcp, gama, &
SupportPillar, PillarSpacing, PillarRadius, &
theta, q, qv, flux, hcin, hrin, hcout, hrout, hin, hout, hcgas, hrgas, ufactor, nperr, ErrorMessage, &
tamb, troom, ibc, Atop, Abot, Al, Ar, Ah, vvent, tvent, LayerType, &
Ra, Nu, vfreevent, qcgas, qrgas, Ebf, Ebb, Rf, Rb, &
ShadeEmisRatioOut, ShadeEmisRatioIn, ShadeHcModifiedOut, ShadeHcModifiedIn, &
ThermalMod, Debug_mode, AchievedErrorTolerance, TotalIndex)
!********************************************************************************
! Main subroutine for calculation of 1-D heat transfer in the center of glazing.
!********************************************************************************
! Inputs
!
! nlayer number of solid layers
! iwd wind direction
! tout outside temp in k
! tind inside temp in k
! wso wind speed in m/s
! wsi inside forced air speed m/s
! Ebsky ir flux from outside
! Gout back facing radiosity from outside
! Trmout Mean outdoor radiant temperature
! Trmin Mean indoor radiant temperature
! Ebroom ir flux from room
! Gin front facing radiosity from room
! tir ir transmittance of each layer
! rir ir reflectance of each surface
! emis ir emittances of each surface
! gap array of gap widths in meters
! thick thickness of glazing layers (m)
! scon Vector of conductivities of 'glazing' layers
! tilt Window tilt (deg). vert: tilt=90, hor out up: tilt=0, hor out down: tilt=180
! sol absorbed solar energy for each layer in w/m2
! height glazing cavity height
! heightt
! iprop
! frct
! presure
! nmix vector of number of gasses in a mixture for each gap
! hin convective indoor film coefficient (if non-zero hin input)
! hout convective outdoor film coeff. (if non-zero hout input)
!
! outputs
!
! theta temp distribution in k
! flux net heat flux between room and window
! rtot overall thermal resistance
! rs ?
! hcin convective indoor film coeff.
! hrin radiative part of indoor film coeff.
! hcout convective outdoor film coeff.
! hrout radiative part of outdoor film coeff.
! hin convective indoor film coefficient
! hout convective outdoor film coeff.
! ufactor overall u-factor
! qcgap vector of convective/conductive parts of flux in gaps
! qrgap vector of radiative parts of flux in gaps
! nperr
! *Inactives**
! wa - window azimuth (degrees, clockwise from south)
! hgas matrix of gap film coefficients
! Locals
! Ebb Vector
! Ebf Vector
! Rb Vector
! Rf Vector
! a Array
! b Array
! hhat Vector
! err iteration tolerance
! dtmax max temp dfference after iteration
! index iteration step
integer, intent(in) :: nlayer, iwd, ThermalMod
integer, intent(in) :: Debug_mode ! Switch for debug output files:
! 0 - don't create debug output files
! 1 - append results to existing debug output file
! 2 - store results in new debug output file
! 3 - save in-between results (in all iterations) to existing debug file
integer, dimension(maxlay), intent(in) :: LayerType
integer, dimension(maxlay1, maxgas), intent(in) :: iprop
integer, dimension(maxlay1), intent(in) :: nmix
integer, dimension(2), intent(in) :: ibc
real(r64), intent(in) :: dir, wso, wsi, VacuumPressure, VacuumMaxGapThickness, Gout, Trmout, Trmin, Gin
real(r64), intent(inout) :: tout, tind
integer, intent(out) :: TotalIndex
real(r64), dimension(MaxGap), intent(in) :: gap
real(r64), dimension(maxlay), intent(in) :: thick
real(r64), dimension(maxlay2), intent(in) :: rir, emis, tir
real(r64), dimension(maxlay), intent(in) :: asol, scon
real(r64), dimension(maxlay1, maxgas), intent(in) :: frct
real(r64), dimension(maxlay1), intent(in) :: presure
real(r64), dimension(maxgas), intent(in) :: wght, gama
real(r64), dimension(maxgas, 3), intent(in) :: gcon, gvis, gcp
real(r64), intent(in) :: tilt, height, heightt
real(r64), dimension(maxlay), intent(in) :: Atop, Abot, Al, Ar, Ah
real(r64), dimension(maxlay1), intent(in) :: vvent, tvent
integer, dimension(maxlay), intent(in) :: SupportPillar
real(r64), dimension(maxlay), intent(in) :: PillarSpacing, PillarRadius
real(r64), intent(inout) :: Ebroom, Ebsky
! real(r64), intent(in) :: sumsol(maxlay)
!integer, intent(in) :: nslice(maxlay)
real(r64), intent(out) :: flux
real(r64), dimension(maxlay3), intent(out) :: q
real(r64), dimension(maxlay1), intent(out) :: qv, qcgas, qrgas
real(r64), dimension(maxlay), intent(out) :: Ra, Nu
real(r64), dimension(maxlay1), intent(out) :: vfreevent
real(r64), intent(out) :: ShadeEmisRatioOut, ShadeEmisRatioIn, ShadeHcModifiedOut, ShadeHcModifiedIn
real(r64), intent(out) :: ufactor, hcin
real(r64), dimension(maxlay1), intent(out) :: hcgas, hrgas
real(r64), intent(inout) :: Tamb, Troom, hin, hout, hcout
real(r64), dimension(maxlay2), intent(inout) :: theta
real(r64), intent(inout) :: hrin, hrout
real(r64), dimension(maxlay), intent(inout) :: Ebb, Ebf, Rb, Rf
real(r64), intent(out) :: AchievedErrorTolerance
integer, intent(inout) :: nperr
character(len=*) :: ErrorMessage
real(r64) :: width, glsyswidth
!real(r64) :: Ebbold(maxlay), Ebfold(maxlay), Rbold(maxlay), Rfold(maxlay)
!real(r64) :: rs(maxlay3)
real(r64), dimension(4*nlayer, 4*nlayer) :: a
real(r64), dimension(4*nlayer) :: b
real(r64), dimension(maxlay1) :: hgas
!real(r64) :: hhatv(maxlay3),hcv(maxlay3), Ebgap(maxlay3), Tgap(maxlay1)
real(r64), dimension(maxlay1) :: Tgap
!real(r64) :: alpha
integer :: maxiter
real(r64) :: qr_gap_out, qr_gap_in
real(r64), dimension(2*nlayer) :: told
! Simon: parameters used in case of JCFN iteration method
!real(r64) :: Dvector(maxlay4) ! store diagonal matrix used in JCFN iterations
real(r64), allocatable :: FRes(:) ! store function results from current iteration
real(r64), allocatable :: FResOld(:) ! store function results from previous iteration
real(r64), allocatable :: FResDiff(:) ! save difference in results between iterations
real(r64), allocatable :: Radiation(:) ! radiation on layer surfaces. used as temporary storage during iterations
real(r64), allocatable :: x(:) ! temporary vector for storing results (theta and Radiation). used for easier handling
real(r64), allocatable :: dX(:) ! difference in results
real(r64), allocatable :: Jacobian(:,:) ! diagonal vector for jacobian comuptation-free newton method
real(r64), allocatable :: DRes(:) ! used in jacobian forward-difference approximation
! This is used to store matrix before equation solver. It is important because solver destroys
! content of matrices
real(r64), allocatable :: LeftHandSide(:,:)
real(r64), allocatable :: RightHandSide(:)
! Simon: Keep best achieved convergence
real(r64) :: prevDifference, Relaxation
real(r64), allocatable :: RadiationSave(:)
real(r64), allocatable :: thetaSave(:)
integer :: currentTry
integer :: LayerTypeSpec(maxlay)
integer :: SDLayerIndex
integer :: CSMFlag
integer :: i, j, k
real(r64) :: curDifference
integer :: index
integer :: curTempCorrection
real(r64) :: qc_gap_out, qcgapout2, hc_modified_out, qc_gap_in, hc_modified_in
integer :: CalcOutcome
logical :: iterationsFinished ! To mark whether or not iterations are finished
logical :: saveIterationResults
logical :: updateGapTemperature
!logical :: TurnOnNewton
SDLayerIndex = -1
! Allocate arrays
if (.not. allocated(FRes)) allocate(FRes(1:4*nlayer))
if (.not. allocated(FResOld)) allocate(FResOld(1:4*nlayer))
if (.not. allocated(FResDiff)) allocate(FResDiff(1:4*nlayer))
if (.not. allocated(Radiation)) allocate(Radiation(1:2*nlayer))
if (.not. allocated(RadiationSave)) allocate(RadiationSave(1:2*nlayer))
if (.not. allocated(thetaSave)) allocate(thetaSave(1:2*nlayer))
if (.not. allocated(X)) allocate(X(1:4*nlayer))
if (.not. allocated(dX)) allocate(dX(1:4*nlayer))
if (.not. allocated(Jacobian)) allocate(Jacobian(1:4*nlayer, 1:4*nlayer))
if (.not. allocated(DRes)) allocate(DRes(1:4*nlayer))
if (.not. allocated(LeftHandSide)) allocate(LeftHandSide(1:4*nlayer, 1:4*nlayer))
if (.not. allocated(RightHandSide)) allocate(RightHandSide(1:4*nlayer))
dX = 0.0d0
! Simon: This is set to zero until it is resolved what to do with modifier
ShadeHcModifiedOut = 0.0d0
!BuffIndex = 0
CSMFlag = 0
CalcOutcome = CALC_UNKNOWN
curTempCorrection = 0
AchievedErrorTolerance = 0.0d0
curDifference = 0.0d0
!TurnOnNewton = .true.
currentTry = 0
index = 0
TotalIndex = 0
iterationsFinished = .false.
qv=0.0d0
Ebb = 0.0d0
Ebf = 0.0d0
Rb = 0.0d0
Rf = 0.0d0
a = 0.0d0
b = 0.0d0
!Dvector = 0.0
FRes = 0.0d0
FResOld = 0.0d0
FResDiff = 0.0d0
Radiation = 0.0d0
Relaxation = RelaxationStart
!alpha = PicardRelaxation
maxiter = NumOfIterations
!call MarkVentilatedGaps(nlayer, isVentilated, LayerType, vvent)
if (Debug_mode == saveIntermediateResults) then
saveIterationResults = .true.
else
saveIterationResults = .false.
end if
!call guess(tout, tind, nlayer, gap, thick, glsyswidth, theta, Ebb, Ebf, Tgap)
do i=1, nlayer
k = 2*i
Radiation(k) = Ebb(i)
Radiation(k-1) = Ebf(i)
told(k-1) = 0.0d0
told(k) = 0.0d0
!told(k-1) = theta(k-1)
!told(k) = theta(k)
end do
!bi...Set LayerTypeSpec array - need to treat venetians AND woven shades as glass:
if (ThermalMod == THERM_MOD_CSM) then
do i = 1, nlayer
if (IsShadingLayer(LayerType(i))) then
LayerTypeSpec(i) = 0
SDLayerIndex = i
else
LayerTypeSpec(i) = LayerType(i)
end if
end do
end if
!first store results before iterations begin
if (saveIterationResults) then
call storeIterationResults(nlayer, index, theta, Trmout, Tamb, Trmin, Troom, Ebsky, Ebroom, &
hcin, hcout, hrin, hrout, hin, hout, Ebb, Ebf, Rb, Rf, nperr)
end if
Tgap(1) = tout
Tgap(nlayer+1) = tind
do i =2, nlayer
Tgap(i) = (theta(2*i-1) + theta(2*i-2) ) / 2
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! MAIN ITERATION LOOP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do while (.not.(iterationsFinished))
do i=1, 2*nlayer
if (theta(i).lt.0) then
theta(i) = 1.0d0 * i
end if
end do
!do i=1,nlayer+1
! if (i == 1) then
! Tgap(i) = tout
! else if (i == nlayer+1) then
! Tgap(i) = tind
! else
! Tgap(i) = (theta(2*i-2) + theta(2*i-1)) / 2.0d0
! end if
!end do
! skip updating gap temperatures for shading devices. Gap temperature in that case is not simply average
! between two layer temperatures
do i =2, nlayer
updateGapTemperature = .false.
if ((.not.(IsShadingLayer(LayerType(i-1)))).and.(.not.(IsShadingLayer(LayerType(i))))) then
updateGapTemperature = .true.
end if
if (updateGapTemperature) then
Tgap(i) = (theta(2*i-1) + theta(2*i-2) ) / 2
end if
end do
! evaluate convective/conductive components of gap
call hatter(nlayer, iwd, tout, tind, wso, wsi, VacuumPressure, VacuumMaxGapThickness, Ebsky, Tamb, Ebroom, Troom, &
gap, height, heightt, scon, tilt, theta, Tgap, Radiation, Trmout, Trmin, &
iprop, frct, presure, nmix, wght, gcon, gvis, gcp, gama, &
SupportPillar, PillarSpacing, PillarRadius, &
hgas, hcgas, hrgas, hcin, hcout, hin, hout, &
index, ibc, nperr, ErrorMessage, hrin, hrout, Ra, Nu)
! exit on error
if (.not.(GoAhead(nperr))) return
!bi...Override hhat values near SHADING DEVICE layer(s), but only for CSM thermal model:
if ((ThermalMod.eq.THERM_MOD_CSM).and.(SDLayerIndex.gt.0)) then
! adjust hhat values
!call adjusthhat(SDLayerIndex, ibc, tout, tind, nlayer, theta, wso, wsi, iwd, height, heightt, tilt, &
! & thick, gap, hout, hrout, hin, hrin, iprop, frct, presure, nmix, wght, gcon, gvis, gcp, &
! index, SDScalar, Ebf, Ebb, hgas, hhat, nperr, ErrorMessage)
!do i = 1, maxlay3
!hhatv(i) = 0.0
!Ebgap(i) = 0.0
!qv(i) = 0.0
!hcv(i) = 0.0
!end do
call matrixQBalance(nlayer, a, b, scon, thick, hcgas, hcout, hcin, asol, qv, &
Tind, Tout, Gin, Gout, theta, tir, rir, emis)
else
!bi...There are no Venetian layers, or ThermalMod is not CSM, so carry on as usual:
call shading(theta, gap, hgas, hcgas, hrgas, frct, iprop, presure, nmix, wght, gcon, gvis, gcp, nlayer, &
width, height, tilt, tout, tind, Atop, Abot, Al, Ar, Ah, vvent, tvent, LayerType, Tgap, qv, nperr, &
ErrorMessage, vfreevent)
! exit on error
if (.not.(GoAhead(nperr))) return
call matrixQBalance(nlayer, a, b, scon, thick, hcgas, hcout, hcin, asol, qv, &
Tind, Tout, Gin, Gout, theta, tir, rir, emis)
end if ! end if
FResOld = FRes
! Pack results in one array
do i = 1, nlayer
k=4*i-3
j=2*i-1
x(k) = theta(j)
x(k+1) = theta(j+1)
x(k+2) = Radiation(j)
x(k+3) = Radiation(j+1)
end do
call CalculateFuncResults(nlayer, a, b, x, FRes)
FResDiff = FRes - FResOld
!if (TurnOnNewton) then
! do i = 1, 4*nlayer
! temp = x(i)
! h = ConvergenceTolerance * abs(x(i))
! if (h == 0) then
! h = ConvergenceTolerance
! end if
! x(i) = temp + h
! h = x(i) - temp ! trick to reduce finite precision error
! call CalculateFuncResults(nlayer, a, b, x, DRes, nperr, ErrorMessage)
! do j = 1, 4*nlayer
! Jacobian(j,i) = (DRes(j) - FRes(j)) / h
! end do
! x(i) = temp
! end do
!end if
!if (TurnOnNewton) then
! LeftHandSide = Jacobian
! RightHandSide = -FRes
!else
LeftHandSide = a
RightHandSide = b
!end if
call EquationsSolver(LeftHandSide, RightHandSide, 4*nlayer, nperr, ErrorMessage)
!if (TurnOnNewton) then
! dx = RightHandSide
!end if
! Simon: This is much better, but also much slower convergence criteria. Think of how to make this flexible and allow
! user to change this from outside (through argument passing)
!curDifference = abs(FRes(1))
!do i = 2, 4*nlayer
!curDifference = max(curDifference, abs(FRes(i)))
!curDifference = curDifference + abs(FRes(i))
!end do
curDifference = abs(theta(1) - told(1))
!curDifference = abs(FRes(1))
do i = 2, 2*nlayer
!do i = 2, 4*nlayer
curDifference = max(curDifference, abs(theta(i) - told(i)))
!curDifference = MAX(abs(FRes(i)), curDifference)
end do
do i=1, nlayer
k=4*i-3
j=2*i-1
!if (TurnOnNewton) then
! theta(j) = theta(j) + Relaxation*dx(k)
! theta(j+1) = theta(j+1) + Relaxation*dx(k+1)
! Radiation(j) = Radiation(j) + Relaxation*dx(k+2)
! Radiation(j+1) = Radiation(j+1) + Relaxation*dx(k+3)
!else
! dX(k) = RightHandSide(k) - theta(j)
! dX(k+1) = RightHandSide(k + 1) - theta(j+1)
! dX(k+2) = RightHandSide(k + 2) - Radiation(j)
! dX(k+3) = RightHandSide(k + 3) - Radiation(j+1)
told(j) = theta(j)
told(j+1) = theta(j+1)
theta(j) = (1 - Relaxation) * theta(j) + Relaxation * RightHandSide(k)
theta(j+1) = (1 - Relaxation) * theta(j+1) + Relaxation * RightHandSide(k+1)
Radiation(j) = (1 - Relaxation) * Radiation(j) + Relaxation * RightHandSide(k+2)
Radiation(j+1) = (1 - Relaxation) * Radiation(j+1) + Relaxation * RightHandSide(k+3)
!end if
end do
! it is important not to update gaps around shading layers since that is already calculated by
! shading routines
do i=1, nlayer+1
updateGapTemperature = .true.
if ((i.eq.1).or.(i.eq.nlayer+1)) then
! update gap array with interior and exterior temperature
updateGapTemperature = .true.
else
! update gap temperature only if gap on both sides
updateGapTemperature = .false.
if ((.not.(IsShadingLayer(LayerType(i-1)))).and.(.not.(IsShadingLayer(LayerType(i))))) then
updateGapTemperature = .true.
end if
end if
j = 2*(i-1)
if (updateGapTemperature) then
if (i.eq.1) then
Tgap(1) = tout
else if (i.eq.(nlayer+1)) then
Tgap(i) = tind
else
Tgap(i) = (theta(j) + theta(j+1) ) / 2
end if
end if
end do
!and store results during iterations
if (saveIterationResults) then
call storeIterationResults(nlayer, index + 1, theta, Trmout, Tamb, Trmin, Troom, Ebsky, Ebroom, &
hcin, hcout, hrin, hrout, hin, hout, Ebb, Ebf, Rb, Rf, nperr)
end if
if (.not.(GoAhead(nperr))) return
prevDifference = curDifference
if ((index == 0).or.(curDifference < AchievedErrorTolerance)) then
AchievedErrorTolerance = curDifference
currentTry = 0
do i=1,2*nlayer
RadiationSave(i) = Radiation(i)
thetaSave(i) = theta(i)
end do
else
! This is case when program solution diverged
currentTry = currentTry + 1
if (currentTry >= NumOfTries) then
currentTry = 0
do i = 1, 2*nlayer
Radiation(i) = RadiationSave(i)
theta(i) = thetaSave(i)
end do
!if (.not.TurnOnNewton) then
! TurnOnNewton = .true.
!else
Relaxation = Relaxation - RelaxationDecrease
TotalIndex = TotalIndex + index
index = 0
! Start from best achieved convergence
if (Relaxation <= 0.0d0) then ! cannot continue with relaxation equal to zero
iterationsFinished = .true.
end if
! TurnOnNewton = .true.
!end if ! if (.not.TurnOnNewton) then
end if ! f (currentTry == NumOfTries) then
end if
! Chek if results were found:
if (curDifference < ConvergenceTolerance) then
CalcOutcome = CALC_OK
TotalIndex = TotalIndex + index
iterationsFinished = .true.
end if
if (index >= maxiter) then
Relaxation = Relaxation - RelaxationDecrease
TotalIndex = TotalIndex + index
index = 0
!TurnOnNewton = .true.
! Start from best achieved convergence
do i = 1, 2*nlayer
Radiation(i) = RadiationSave(i)
theta(i) = thetaSave(i)
end do
if (Relaxation <= 0.0d0) then ! cannot continue with relaxation equal to zero
iterationsFinished = .true.
end if
end if
index = index + 1
end do
! Get results from closest iteration and store it
if (CalcOutcome == CALC_OK) then
do i=1,2*nlayer
Radiation(i) = RadiationSave(i)
theta(i) = thetaSave(i)
end do
do i =2, nlayer
updateGapTemperature = .false.
if ((.not.(IsShadingLayer(LayerType(i-1)))).and.(.not.(IsShadingLayer(LayerType(i))))) then
updateGapTemperature = .true.
end if
if (updateGapTemperature) then
Tgap(i) = (theta(2*i-1) + theta(2*i-2) ) / 2
end if
end do
! Simon: It is important to recalculate coefficients from most accurate run
call hatter(nlayer, iwd, tout, tind, wso, wsi, VacuumPressure, VacuumMaxGapThickness, Ebsky, Tamb, Ebroom, Troom, &
gap, height, heightt, scon, tilt, theta, Tgap, Radiation, Trmout, Trmin, &
iprop, frct, presure, nmix, wght, gcon, gvis, gcp, gama, &
SupportPillar, PillarSpacing, PillarRadius, hgas, hcgas, hrgas, hcin, hcout, hin, hout, &
index, ibc, nperr, ErrorMessage, hrin, hrout, Ra, Nu)
call shading(theta, gap, hgas, hcgas, hrgas, frct, iprop, presure, nmix, wght, gcon, gvis, gcp, nlayer, &
width, height, tilt, tout, tind, Atop, Abot, Al, Ar, Ah, vvent, tvent, LayerType, Tgap, qv, nperr, &
ErrorMessage, vfreevent)
end if
if (CalcOutcome.eq.CALC_UNKNOWN) then
ErrorMessage = 'Tarcog failed to converge'
nperr = 2 ! error 2: failed to converge...
end if
! Get radiation results first
!if (curEquationsApproach.eq.eaQBalance) then
do i=1, nlayer
k=2*i-1
Rf(i) = Radiation(k)
Rb(i) = Radiation(k+1)
Ebf(i) = StefanBoltzmann*theta(k)**4
Ebb(i) = StefanBoltzmann*theta(k+1)**4
end do
!end if
! Finishing calcs:
call resist(nlayer, Trmout, Tout, Trmin, tind, hcgas, hrgas, Theta, q, &
& qv, LayerType, thick, scon, ufactor, flux, &
& qcgas, qrgas)
!bi... Set T6-related quantities - ratios for modified epsilon, hc for modelling external SDs:
! (using non-solar pass results)
if ((dir.eq.0.0d0).and.(nlayer.gt.1)) then
qr_gap_out = Rf(2) - Rb(1)
qr_gap_in = Rf(nlayer) - Rb(nlayer-1)
if (IsShadingLayer(LayerType(1))) then
ShadeEmisRatioOut = qr_gap_out / ( emis(3) * StefanBoltzmann * (theta(3)**4 - Trmout**4) )
!qc_gap_out = qprim(3) - qr_gap_out
!qcgapout2 = qcgas(1)
!Hc_modified_out = (qc_gap_out / (theta(3) - tout))
!ShadeHcModifiedOut = Hc_modified_out
end if
if (IsShadingLayer(LayerType(nlayer))) then
ShadeEmisRatioIn = qr_gap_in / ( emis(2*nlayer-2) * StefanBoltzmann * (Trmin**4 - theta(2*nlayer-2)**4) )
qc_gap_in = q(2*nlayer-1) - qr_gap_in
Hc_modified_in = (qc_gap_in / (Tind - theta(2*nlayer-2)))
ShadeHcModifiedIn = Hc_modified_in
end if
end if ! IF dir = 0
if (allocated(FRes)) deallocate(FRes)
if (allocated(FResOld)) deallocate(FResOld)
if (allocated(FResDiff)) deallocate(FResDiff)
if (allocated(Radiation)) deallocate(Radiation)
if (allocated(RadiationSave)) deallocate(RadiationSave)
if (allocated(thetaSave)) deallocate(thetaSave)
if (allocated(X)) deallocate(X)
if (allocated(dX)) deallocate(dX)
if (allocated(Jacobian)) deallocate(Jacobian)
if (allocated(DRes)) deallocate(DRes)
if (allocated(LeftHandSide)) deallocate(LeftHandSide)
if (allocated(RightHandSide)) deallocate(RightHandSide)
!do i=1, nlayer-1
! if (((LayerType(i).eq.VENETBLIND) &
! & .and.(ThermalMod.ne.THERM_MOD_CSM)) &
! & .or.(LayerType(i).eq.WOVSHADE)) then
! !hcgas(i+1)=hcv(2*i+1)
! else
! !hcgas(i+1)=hgas(i+1)
! end if
!end do
1111 format('Outdoor: ', F9.6,' ; alt2: ', F9.6, ' ; alt3: ', F9.6, ' ; alt4: ', F9.6)
1112 format('Indoor: ', F9.6,' ; alt2: ', F9.6, ' ; alt3: ', F9.6, ' ; alt4: ', F9.6)
!110 format(' Theta(',I1,') = ',F12.6)
!111 format(' T(',I1,')=',F15.9)
!112 format(' ',A3,' =',F15.9)
end subroutine therm1d