| blob | d9ce496f31dc3f0efc1203e58e6ebab4a10918c1 |
1 MODULE module_sf_bem
2 #if defined(mpas)
3 use mpas_atmphys_utilities, only: physics_error_fatal
4 #define FATAL_ERROR(M) call physics_error_fatal( M )
5 #else
6 use module_wrf_error
7 #define FATAL_ERROR(M) call wrf_error_fatal( M )
8 #endif
9 ! -----------------------------------------------------------------------
10 ! Variables and constants used in the BEM module
11 ! -----------------------------------------------------------------------
13 real emins !emissivity of the internal walls
14 parameter (emins=0.9)
15 real albins !albedo of the internal walls
16 !! parameter (albins=0.5)
17 parameter (albins=0.3)
19 real thickwin !thickness of the window [m]
20 parameter (thickwin=0.006)
21 real cswin !Specific heat of the windows [J/(m3.K)]
22 parameter(cswin= 2.268e+06)
24 real temp_rat !power of the A.C. heating/cooling the indoor air [K/s]
25 parameter(temp_rat=0.001)
27 real hum_rat !power of the A.C. drying/moistening the indoor air [(Kg/kg)/s]
28 parameter(hum_rat=1.e-06)
30 real,parameter :: effpv=0.19 ! Efficiency of PV panels installed at the roofs, typical values [0.10,0.15]
32 CONTAINS
34 !====6================================================================72
35 !====6================================================================72
37 subroutine BEM(nzcanm,nlev,nhourday,dt,bw,bl,dzlev, &
38 nwal,nflo,nrof,ngrd,hswalout,gswal, &
39 hswinout,hsrof,lsrof,gsrof,hspv, &
40 latent,sigma,albwal,albwin,albrof, &
41 emrof,emwal,emwin,rswal,rlwal,rair,cp, &
42 rhoout,tout,humout,press, &
43 rs,swddif,rl,dzwal,cswal,kwal,pwin,cop,beta,sw_cond, &
44 timeon,timeoff,targtemp,gaptemp,targhum,gaphum, &
45 perflo,gr_frac_roof,pv_frac_roof,gr_flag, &
46 uout,vout, &
47 hsesf,hsequip,dzflo, &
48 csflo,kflo,dzgrd,csgrd,kgrd,dzrof,csrof, &
49 krof,tlev,shumlev,twal,twin,tflo,tgrd,trof, &
50 hsout,hlout,consump,eppv,tpv,hsvent,hlvent,hfgr,&
51 tr_av,tpv_print,sfpv,sfr_indoor)
56 ! ---------------------------------------------------------------------
57 implicit none
59 ! ---------------------------------------------------------------------
60 ! TOP
61 ! ---------------------
62 ! ! ----------------- !--->roof (-) : level number
63 ! ! ! ! ! rem: the windows are given
64 ! ! !---------------! ! with respect to the
65 ! ! !---------------! ! vertical walls-->win(2)
66 ! (n)! !(1) (1)!-!(n)
67 ! ! !---------------! ! 2D vision of the building
68 ! WEST ! !-------4-------! ! EAST
69 ! I ! ! 1 ilev 2! ! II ^
70 ! ! !-------3--------! ! !
71 ! ! !---------------! !--->floor 1 !
72 ! ! ! ! ! !
73 ! ! ! ! ! !
74 ! ! ----------------- ! <--------------(n)
75 ! ------------------------>ground ------------(1)
76 ! BOTTOM
77 ! i(6)
78 ! i
79 ! +---------v-----+
80 ! /| /| 3D vision of a room
81 ! / | 4 / |
82 ! / | / |
83 ! / | / |
84 ! / | / |
85 ! +---------------+ |
86 ! | 1 | | 2 |
87 ! | +---------|-----+
88 ! dzlev | / | /
89 ! | / 3 | /
90 ! | /bw | /
91 ! | / | /
92 ! |/ |/
93 ! +---------------+
94 ! ^ bl
95 ! i
96 ! i
97 ! (5)
98 !-----------------------------------------------------------------------
101 ! Input:
102 ! -----
104 real dt !time step [s]
106 integer nzcanm !Maximum number of vertical levels in the urban grid
107 integer nlev !number of floors in the building
108 integer nwal !number of levels inside the wall
109 integer nrof !number of levels inside the roof
110 integer nflo !number of levels inside the floor
111 integer ngrd !number of levels inside the ground
112 real dzlev !vertical grid resolution [m]
113 real bl !Building length [m]
114 real bw !Building width [m]
116 real albwal !albedo of the walls
117 real albwin !albedo of the windows
118 real albrof !albedo of the roof
120 real emwal !emissivity of the walls
122 real emrof !emissivity of the roof
123 real emwin !emissivity of the windows
125 real pwin !window proportion
126 real, intent(in) :: cop !Coefficient of performance of the A/C systems
127 real, intent(in) :: beta !Thermal efficiency of the heat exchanger
128 integer, intent(in) :: sw_cond ! Air Conditioning switch
129 real, intent(in) :: timeon ! Initial local time of A/C systems
130 real, intent(in) :: timeoff ! Ending local time of A/C systems
131 real, intent(in) :: targtemp ! Target temperature of A/C systems
132 real, intent(in) :: gaptemp ! Comfort range of indoor temperature
133 real, intent(in) :: targhum ! Target humidity of A/C systems
134 real, intent(in) :: gaphum ! Comfort range of specific humidity
135 real, intent(in) :: perflo ! Peak number of occupants per unit floor area
136 real gr_frac_roof
137 real pv_frac_roof
138 integer gr_flag
139 real uout(nzcanm)
140 real vout(nzcanm)
141 real, intent(in) :: hsesf !
142 real, intent(in) :: hsequip(24) !
144 real cswal(nwal) !Specific heat of the wall [J/(m3.K)]
146 real csflo(nflo) !Specific heat of the floor [J/(m3.K)]
147 real csrof(nrof) !Specific heat of the roof [J/(m3.K)]
148 real csgrd(ngrd) !Specific heat of the ground [J/(m3.K)]
150 real kwal(nwal+1) !Thermal conductivity in each layers of the walls (face) [W/(m.K)]
151 real kflo(nflo+1) !Thermal diffusivity in each layers of the floors (face) [W/(m.K)]
152 real krof(nrof+1) !Thermal diffusivity in each layers of the roof (face) [W/(m.K)]
153 real kgrd(ngrd+1) !Thermal diffusivity in each layers of the ground (face) [W/(m.K)]
155 real dzwal(nwal) !Layer sizes of walls [m]
156 real dzflo(nflo) !Layer sizes of floors [m]
157 real dzrof(nrof) !Layer sizes of roof [m]
158 real dzgrd(ngrd) !Layer sizes of ground [m]
159 real tpv
160 real tr_av
161 real latent !latent heat of evaporation [J/Kg]
163 real swddif
164 real rs !external short wave radiation [W/m2]
165 real rl !external long wave radiation [W/m2]
166 real rswal(4,nzcanm) !short wave radiation reaching the exterior walls [W/m2]
167 real rlwal(4,nzcanm) !long wave radiation reaching the walls [W/m2]
168 real rhoout(nzcanm) !exterior air density [kg/m3]
169 real tout(nzcanm) !external temperature [K]
170 real humout(nzcanm) !absolute humidity [Kgwater/Kgair]
171 real press(nzcanm) !external air pressure [Pa]
173 real hswalout(4,nzcanm) !outside walls sensible heat flux [W/m2]
174 real hswinout(4,nzcanm) !outside window sensible heat flux [W/m2]
175 real hsrof !Sensible heat flux at the roof [W/m2]
176 real lsrof
177 real rair !ideal gas constant [J.kg-1.K-1]
178 real sigma !parameter (wall is not black body) [W/m2.K4]
179 real cp !specific heat of air [J/kg.K]
180 real hfgr !Green roof heat flux
182 !Input-Output
183 !------------
184 real tlev(nzcanm) !temperature of the floors [K]
185 real shumlev(nzcanm) !specific humidity of the floor [kg/kg]
186 real twal(4,nwal,nzcanm) !walls temperatures [K]
187 real twin(4,nzcanm) !windows temperatures [K]
188 real tflo(nflo,nzcanm-1) !floor temperatures [K]
189 real tgrd(ngrd) !ground temperature [K]
190 real trof(nrof) !roof temperature [K]
191 real hsout(nzcanm) !sensible heat emitted outside the floor [W]
192 real hlout(nzcanm) !latent heat emitted outside the floor [W]
193 real consump(nzcanm) !Consumption for the a.c. in each floor [W]
194 real hsvent(nzcanm) !sensible heat generated by natural ventilation [W]
195 real hlvent(nzcanm) !latent heat generated by natural ventilation [W]
196 real gsrof !heat flux flowing inside the roof [W/m2]
197 real hspv !Sensible heat flux at the roof from the PV panels [W/m2]
198 real gswal(4,nzcanm) !heat flux flowing inside the floors [W/m2]
199 real eppv !Electricity production of PV panels [W]
200 real sfr_indoor,sfpv,tpv_print
201 ! Local:
202 ! -----
203 integer swwal !swich for the physical coefficients calculation
204 integer ilev !index for rooms
205 integer iwal !index for walls
206 integer iflo !index for floors
207 integer ivw !index for vertical walls
208 integer igrd !index for ground
209 integer irof !index for roof
210 real hseqocc(nzcanm) !sensible heat generated by equipments and occupants [W]
211 real hleqocc(nzcanm) !latent heat generated by occupants [W]
212 real hscond(nzcanm) !sensible heat generated by wall conduction [W]
213 real hslev(nzcanm) !sensible heat flux generated inside the room [W]
214 real hllev(nzcanm) !latent heat flux generatd inside the room [W]
215 real surwal(6,nzcanm) !Surface of the walls [m2]
216 real surwal1D(6) !wall surfaces of a generic room [m2]
217 real rsint(6) !short wave radiation reaching the indoor walls[W/m2]
218 real rswalins(6,nzcanm) !internal short wave radiation for the building [W/m2]
219 real twin1D(4) !temperature of windows for a particular room [K]
220 real twal_int(6) !temperature of the first internal layers of a room [K]
221 real rlint(6) !internal wall long wave radiation [w/m2]
222 real rlwalins(6,nzcanm) !internal long wave radiation for the building [W/m2]
223 real hrwalout(4,nzcanm) !external radiative flux to the walls [W/m2]
224 real hrwalins(6,nzcanm) !inside radiative flux to the walls [W/m2]
225 real hrwinout(4,nzcanm) !external radiative flux to the window [W/m2]
226 real hrwinins(4,nzcanm) !inside radiative flux to the window [W/m2]
227 real hrrof !external radiative flux to the roof [W/m2]
228 real hs
229 real hsneed(nzcanm) !sensible heat needed by the room [W]
230 real hlneed(nzcanm) !latent heat needed by the room [W]
231 real hswalins(6,nzcanm) !inside walls sensible heat flux [W/m2]
232 real hswalins1D(6)
233 real hswinins(4,nzcanm) !inside window sensible heat flux [W/m2]
234 real hswinins1D(4)
235 real htot(2) !total heat flux at the wall [W/m2]
236 real twal1D(nwal)
237 real tflo1D(nflo)
238 real tgrd1D(ngrd)
239 real trof1D(nrof)
240 real rswal1D(4)
241 real Qb !Overall heat capacity of the indoor air [J/K]
242 real vollev !volume of the room [m3]
243 real rhoint !density of the internal air [Kg/m3]
244 real cpint !specific heat of the internal air [J/kg.K]
245 real humdry !specific humidiy of dry air [kg water/kg dry air]
246 real radflux !Function to compute the total radiation budget
247 real consumpbuild !Energetic consumption for the entire building [KWh/s]
248 real hsoutbuild !Total sensible heat ejected into the atmosphere[W]
249 !by the air conditioning system and per building
250 real nhourday !number of hours from midnight, local time
251 real hfgrd !Dummy variable to assign hfgr=0 to walls, windows and ground
253 !
254 parameter(hfgrd=0)
255 !--------------------------------------------
256 !Initialization
257 !--------------------------------------------
258 do ilev=1,nzcanm
259 hseqocc(ilev)=0.
260 hleqocc(ilev)=0.
261 hscond(ilev)=0.
262 hslev(ilev)=0.
263 hllev(ilev)=0.
264 enddo
266 !Calculation of the surfaces of the building
267 !--------------------------------------------
270 do ivw=1,6
271 do ilev=1,nzcanm
272 surwal(ivw,ilev)=1. !initialisation
273 end do
274 end do
276 do ilev=1,nlev
277 do ivw=1,2
278 surwal(ivw,ilev)=dzlev*bw
279 end do
280 do ivw=3,4
281 surwal(ivw,ilev)=dzlev*bl
282 end do
283 do ivw=5,6
284 surwal(ivw,ilev)=bw*bl
285 end do
286 end do
289 ! Calculation of the short wave radiations at the internal walls
290 ! ---------------------------------------------------------------
293 do ilev=1,nlev
295 do ivw=1,4
296 rswal1D(ivw)=rswal(ivw,ilev)
297 end do
299 do ivw=1,6
300 surwal1D(ivw)=surwal(ivw,ilev)
301 end do
303 call int_rsrad(albwin,albins,pwin,rswal1D,&
304 surwal1D,bw,bl,dzlev,rsint)
306 do ivw=1,6
307 rswalins(ivw,ilev)=rsint(ivw)
308 end do
310 end do !ilev
314 ! Calculation of the long wave radiation at the internal walls
315 !-------------------------------------------------------------
318 !Intermediate rooms
320 if (nlev.gt.2) then
321 do ilev=2,nlev-1
323 do ivw=1,4
324 twin1D(ivw)=twin(ivw,ilev)
325 twal_int(ivw)=twal(ivw,1,ilev)
326 end do
328 twal_int(5)=tflo(nflo,ilev-1)
329 twal_int(6)=tflo(1,ilev)
331 call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
332 pwin,bw,bl,dzlev,rlint)
335 do ivw=1,6
336 rlwalins(ivw,ilev)=rlint(ivw)
337 end do
339 end do !ilev
340 end if
343 if (nlev.ne.1) then
345 !bottom room
347 do ivw=1,4
348 twin1D(ivw)=twin(ivw,1)
349 twal_int(ivw)=twal(ivw,1,1)
350 end do
352 twal_int(5)=tgrd(ngrd)
353 twal_int(6)=tflo(1,1)
356 call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
357 pwin,bw,bl,dzlev,rlint)
359 do ivw=1,6
360 rlwalins(ivw,1)=rlint(ivw)
361 end do
363 !top room
365 do ivw=1,4
366 twin1D(ivw)=twin(ivw,nlev)
367 twal_int(ivw)=twal(ivw,1,nlev)
368 end do
370 twal_int(5)=tflo(nflo,nlev-1)
371 twal_int(6)=trof(1)
374 call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
375 pwin,bw,bl,dzlev,rlint)
377 do ivw=1,6
378 rlwalins(ivw,nlev)=rlint(ivw)
379 end do
381 else !Top <---> Bottom
383 do ivw=1,4
384 twin1D(ivw)=twin(ivw,1)
385 twal_int(ivw)=twal(ivw,1,1)
386 end do
388 twal_int(5)=tgrd(ngrd)
389 twal_int(6)=trof(1)
391 call int_rlrad(emins,emwin,sigma,twal_int,twin1D, &
392 pwin,bw,bl,dzlev,rlint)
394 do ivw=1,6
395 rlwalins(ivw,1)=rlint(ivw)
396 end do
398 end if
401 ! Calculation of the radiative fluxes
402 ! -----------------------------------
404 !External vertical walls and windows
406 do ilev=1,nlev
407 do ivw=1,4
408 call radfluxs(radflux,albwal,rswal(ivw,ilev), &
409 emwal,rlwal(ivw,ilev),sigma, &
410 twal(ivw,nwal,ilev))
412 hrwalout(ivw,ilev)=radflux
414 hrwinout(ivw,ilev)=emwin*rlwal(ivw,ilev)- &
415 emwin*sigma*(twin(ivw,ilev)**4)
418 end do ! ivw
419 end do ! ilev
421 !Roof and PV panels
423 if (pv_frac_roof.eq.0.) then
424 call radfluxs(radflux,albrof,rs,emrof,rl,sigma,tr_av)
425 hrrof=radflux
426 else
427 call radfluxspv(nzcanm,nlev,albrof,rs,swddif,emrof,rl,tr_av,tout,sigma,radflux,pv_frac_roof,tpv)
428 hrrof=radflux
429 endif
431 !Internal walls for intermediate rooms
433 if(nlev.gt.2) then
435 do ilev=2,nlev-1
436 do ivw=1,4
438 call radfluxs(radflux,albins,rswalins(ivw,ilev), &
439 emins,rlwalins(ivw,ilev),sigma, &
440 twal(ivw,1,ilev))
442 hrwalins(ivw,ilev)=radflux
444 end do !ivw
446 call radfluxs(radflux,albins,rswalins(5,ilev), &
447 emins,rlwalins(5,ilev),sigma,&
448 tflo(nflo,ilev-1))
450 hrwalins(5,ilev)=radflux
452 call radfluxs(radflux,albins,rswalins(6,ilev), &
453 emins,rlwalins(6,ilev),sigma,&
454 tflo(1,ilev))
455 hrwalins(6,ilev)=radflux
457 end do !ilev
459 end if
462 !Internal walls for the bottom and the top room
463 !
464 if (nlev.ne.1) then
466 !bottom floor
468 do ivw=1,4
470 call radfluxs(radflux,albins,rswalins(ivw,1), &
471 emins,rlwalins(ivw,1),sigma, &
472 twal(ivw,1,1))
474 hrwalins(ivw,1)=radflux
476 end do
479 call radfluxs(radflux,albins,rswalins(5,1),&
480 emins,rlwalins(5,1),sigma,& !bottom
481 tgrd(ngrd))
483 hrwalins(5,1)=radflux
486 call radfluxs(radflux,albins,rswalins(6,1),&
487 emins,rlwalins(6,1),sigma,&
488 tflo(1,1))
490 hrwalins(6,1)=radflux
492 !roof floor
494 do ivw=1,4
496 call radfluxs(radflux,albins,rswalins(ivw,nlev), &
497 emins,rlwalins(ivw,nlev),sigma,&
498 twal(ivw,1,nlev))
500 hrwalins(ivw,nlev)=radflux
502 end do !top
505 call radfluxs(radflux,albins,rswalins(5,nlev), &
506 emins,rlwalins(5,nlev),sigma,&
507 tflo(nflo,nlev-1))
509 hrwalins(5,nlev)=radflux
511 call radfluxs(radflux,albins,rswalins(6,nlev), &
512 emins,rlwalins(6,nlev),sigma,&
513 trof(1))
515 hrwalins(6,nlev)=radflux
517 else ! Top <---> Bottom room
519 do ivw=1,4
521 call radfluxs(radflux,albins,rswalins(ivw,1),&
522 emins,rlwalins(ivw,1),sigma, &
523 twal(ivw,1,1))
525 hrwalins(ivw,1)=radflux
527 end do
529 call radfluxs(radflux,albins,rswalins(5,1),&
530 emins,rlwalins(5,1),sigma, &
531 tgrd(ngrd))
533 hrwalins(5,1)=radflux
535 call radfluxs(radflux,albins,rswalins(6,nlev), &
536 emins,rlwalins(6,nlev),sigma,&
537 trof(1))
538 hrwalins(6,1)=radflux
540 end if
543 !Windows
545 do ilev=1,nlev
546 do ivw=1,4
547 hrwinins(ivw,ilev)=emwin*rlwalins(ivw,ilev)- &
548 emwin*sigma*(twin(ivw,ilev)**4)
549 end do
550 end do
553 ! Calculation of the sensible heat fluxes
554 ! ---------------------------------------
556 !Vertical fluxes for walls
558 do ilev=1,nlev
559 do ivw=1,4
561 call hsinsflux (2,2,tlev(ilev),twal(ivw,1,ilev),hs)
563 hswalins(ivw,ilev)=hs
565 end do ! ivw
566 end do ! ilev
569 !Vertical fluxes for windows
571 do ilev=1,nlev
573 do ivw=1,4
575 call hsinsflux (2,1,tlev(ilev),twin(ivw,ilev),hs)
577 hswinins(ivw,ilev)=hs
579 end do ! ivw
581 end do !ilev
583 !Horizontal fluxes
585 if (nlev.gt.2) then
587 do ilev=2,nlev-1
589 call hsinsflux (1,2,tlev(ilev),tflo(nflo,ilev-1),hs)
591 hswalins(5,ilev)=hs
593 call hsinsflux (1,2,tlev(ilev),tflo(1,ilev),hs)
595 hswalins(6,ilev)=hs
597 end do ! ilev
599 end if
601 if (nlev.ne.1) then
603 call hsinsflux (1,2,tlev(1),tgrd(ngrd),hs)
605 hswalins(5,1)=hs !Bottom room
607 call hsinsflux (1,2,tlev(1),tflo(1,1),hs)
609 hswalins(6,1)=hs
611 call hsinsflux (1,2,tlev(nlev),tflo(nflo,nlev-1),hs)
613 hswalins(5,nlev)=hs !Top room
615 call hsinsflux (1,2,tlev(nlev),trof(1),hs)
617 hswalins(6,nlev)=hs
618 sfr_indoor= hswalins(6,nlev)
619 else ! Bottom<--->Top
621 call hsinsflux (1,2,tlev(1),tgrd(ngrd),hs)
623 hswalins(5,1)=hs
625 call hsinsflux (1,2,tlev(nlev),trof(1),hs)
627 hswalins(6,nlev)=hs
629 end if
630 !!
631 !! Calculation of the sensible heat fluxes from the PV panels & electricity producti
634 if(pv_frac_roof.gt.0)then
635 call hsfluxpv(nzcanm,nlev,bl,bw,albrof,rs,swddif,emrof,rl,tr_av,tout,sigma,hspv,eppv,pv_frac_roof,uout,vout,tpv,dt)
636 sfpv=hspv
637 tpv_print=tpv
638 endif
640 !Calculation of the temperature for the different surfaces
641 ! --------------------------------------------------------
643 ! Vertical walls
645 swwal=1
646 do ilev=1,nlev
647 do ivw=1,4
649 htot(1)=hswalins(ivw,ilev)+hrwalins(ivw,ilev)
650 htot(2)=hswalout(ivw,ilev)+hrwalout(ivw,ilev)
651 gswal(ivw,ilev)=htot(2)
653 do iwal=1,nwal
654 twal1D(iwal)=twal(ivw,iwal,ilev)
655 end do
657 call wall(swwal,nwal,dt,dzwal,kwal,cswal,htot,twal1D)
659 do iwal=1,nwal
660 twal(ivw,iwal,ilev)=twal1D(iwal)
661 end do
663 end do ! ivw
664 end do ! ilev
666 ! Windows
668 do ilev=1,nlev
669 do ivw=1,4
671 htot(1)=hswinins(ivw,ilev)+hrwinins(ivw,ilev)
672 htot(2)=hswinout(ivw,ilev)+hrwinout(ivw,ilev)
674 twin(ivw,ilev)=twin(ivw,ilev)+(dt/(cswin*thickwin))* &
675 (htot(1)+htot(2))
677 end do ! ivw
678 end do ! ilev
680 ! Horizontal floors
683 if (nlev.gt.1) then
684 swwal=1
685 do ilev=1,nlev-1
687 htot(1)=hrwalins(6,ilev)+hswalins(6,ilev)
688 htot(2)=hrwalins(5,ilev+1)+hswalins(5,ilev+1)
690 do iflo=1,nflo
691 tflo1D(iflo)=tflo(iflo,ilev)
692 end do
694 call wall(swwal,nflo,dt,dzflo,kflo,csflo,htot,tflo1D)
696 do iflo=1,nflo
697 tflo(iflo,ilev)=tflo1D(iflo)
698 end do
700 end do ! ilev
701 end if
704 ! Ground
706 swwal=1
708 htot(1)=0. !Diriclet b.c. at the internal boundary
709 htot(2)=hswalins(5,1)+hrwalins(5,1)
711 do igrd=1,ngrd
712 tgrd1D(igrd)=tgrd(igrd)
713 end do
715 call wall(swwal,ngrd,dt,dzgrd,kgrd,csgrd,htot,tgrd1D)
717 do igrd=1,ngrd
718 tgrd(igrd)=tgrd1D(igrd)
719 end do
722 ! Roof
724 swwal=1
726 htot(1)=hswalins(6,nlev)+hrwalins(6,nlev)
727 htot(2)=hsrof+hrrof+lsrof
728 gsrof=htot(2)
730 do irof=1,nrof
731 trof1D(irof)=trof(irof)
733 end do
735 if(gr_flag.eq.1)then
736 call wall_gr(hfgr,gr_frac_roof,swwal,nrof,dt,dzrof,krof,csrof,htot,trof1D)
737 else
738 call wall(swwal,nrof,dt,dzrof,krof,csrof,htot,trof1D)
739 endif
740 do irof=1,nrof
741 trof(irof)=trof1D(irof)
743 end do
746 ! Calculation of the heat fluxes and of the temperature of the rooms
747 ! ------------------------------------------------------------------
749 do ilev=1,nlev
751 !Calculation of the heat generated by equipments and occupants
753 call fluxeqocc(nhourday,bw,bl,perflo,hsesf,hsequip,hseqocc(ilev),hleqocc(ilev))
755 !Calculation of the heat generated by natural ventilation
757 vollev=bw*bl*dzlev
758 humdry=shumlev(ilev)/(1.-shumlev(ilev))
759 rhoint=(press(ilev))/(rair*(1.+0.61*humdry)*tlev(ilev))
760 cpint=cp*(1.+0.84*humdry)
763 call fluxvent(cpint,rhoint,vollev,tlev(ilev),tout(ilev), &
764 latent,humout(ilev),rhoout(ilev),shumlev(ilev),&
765 beta,hsvent(ilev),hlvent(ilev))
768 !Calculation of the heat generated by conduction
770 do iwal=1,6
771 hswalins1D(iwal)=hswalins(iwal,ilev)
772 surwal1D(iwal)=surwal(iwal,ilev)
773 end do
775 do iwal=1,4
776 hswinins1D(iwal)=hswinins(iwal,ilev)
777 end do
779 call fluxcond(hswalins1D,hswinins1D,surwal1D,pwin,&
780 hscond(ilev))
782 !Calculation of the heat generated inside the room
784 call fluxroo(hseqocc(ilev),hleqocc(ilev),hsvent(ilev), &
785 hlvent(ilev),hscond(ilev),hslev(ilev),hllev(ilev))
788 !Evolution of the temperature and of the specific humidity
790 Qb=rhoint*cpint*vollev
792 ! temperature regulation
794 call regtemp(sw_cond,nhourday,dt,Qb,hslev(ilev), &
795 tlev(ilev),timeon,timeoff,targtemp,gaptemp,hsneed(ilev))
797 ! humidity regulation
799 call reghum(sw_cond,nhourday,dt,vollev,rhoint,latent, &
800 hllev(ilev),shumlev(ilev),timeon,timeoff,&
801 targhum,gaphum,hlneed(ilev))
802 !
803 !performance of the air conditioning system for the test
804 !
806 call air_cond(hsneed(ilev),hlneed(ilev),dt, &
807 hsout(ilev),hlout(ilev),consump(ilev), cop)
809 tlev(ilev)=tlev(ilev)+(dt/Qb)*(hslev(ilev)-hsneed(ilev))
811 shumlev(ilev)=shumlev(ilev)+(dt/(vollev*rhoint*latent))* &
812 (hllev(ilev)-hlneed(ilev))
814 end do !ilev
816 call consump_total(nzcanm,nlev,consumpbuild,hsoutbuild, &
817 hsout,consump)
819 return
820 end subroutine BEM
821 !====6=8===============================================================72
822 !====6=8===============================================================72
823 subroutine hsfluxpv(nz,n,bl,bw,albr,rs,swddif,emr,rl,tr,tair,sigma,hspv,eppv,pv_frac_roof,uout,vout,tpv,dt)
824 !
825 implicit none
826 !
827 ! Input variables
828 !
829 integer,intent(in) :: nz !Maximum number of vertical levels in the urban grid
830 real,intent(in) :: bl !Building length [m]
831 real,intent(in) :: bw !Building width [m]
832 real,intent(in) :: albr !albedo of the roof (ext.)
833 real,intent(in) :: emr !emissivity of the roof (ext.)
834 real,intent(in) :: rs !external short wave radiation [W/m2]
835 real,intent(in) :: rl !external long wave radiation [W/m2]
836 real,intent(in) :: tr !roof surface temperature [K]
837 real,intent(in) :: pv_frac_roof ! fraction of PV []
838 real,intent(in) :: sigma !Stefan-Boltzmann constant [W/m2.K4]
839 real,intent(in),dimension(1:nz) :: tair !external temperature [K]
840 integer,intent(in) :: n !number of floors in the building
841 real,intent(in), dimension(1:nz) :: uout
842 real,intent(in), dimension(1:nz) :: vout
843 real,intent(in) :: dt
844 real,intent(in) :: swddif
845 ! Output variables
846 !
847 real,intent(inout) :: hspv ! Sensible heat flux from the PV panels to the atmosphere [W/m2]
848 real,intent(inout) :: eppv ! Electricity production from PV panels [W]
849 real,intent(inout) :: tpv !Temperature of the PV panels [K]
850 !
851 ! Local variables
852 !
853 real,parameter :: albpv=0.11 ! albedo of the PV panels
854 real,parameter :: empv_down=0.95 ! emissivity of the PV panels
855 real,parameter :: empv_up=0.79
856 real, parameter :: T_amb=25
857 real, parameter :: tiltangle=0.
858 real, parameter :: a=3.8
859 real, parameter :: b=6.9
861 real, parameter :: r1=2330.
862 real, parameter :: r2=1200.
863 real, parameter :: r3=3000.
864 real, parameter :: c1=677.
865 real, parameter :: c2=1250.
866 real, parameter :: c3=500.
867 real, parameter :: d1=0.0003
868 real, parameter :: d2=0.0005
869 real, parameter :: d3=0.003
870 real, parameter :: F12=1.
871 real :: lwuppv !Long-wave emitted by the PV panels to the sky [W/m2]
872 real :: lwdwr !Long-wave incoming radiation on the roof [W/m2]
873 real :: lwupr !Long-wave coming up from the roof intercepted by the PV panels [W/m2]
874 real :: enerpv !Energy produced by PV panels [W/m2]
875 real :: hc
876 real :: sw_d
877 real :: lw_d
878 real :: lwpv_out
879 real :: tpv_new
880 real :: hdown
881 real :: hup
882 real :: deltat
883 real :: uroof
884 real :: hrad
885 real :: Cm
886 real :: hf
887 Cm=r1*c1*d1+r2*c2*d2+r3*c3*d3
888 hrad=sigma/((1-empv_down)/empv_down+1/F12+(1-emr)/emr)
889 uroof=(uout(n+1)**2+vout(n+1)**2)**0.5
890 deltat=tpv-tair(n+1)
891 hf=2.5*(40./100.*uroof)**(0.5)
892 hc=9.482*abs(deltat)**(1./3.)/(7.238-abs(cos(tiltangle)))
893 hup=sqrt(hc**2.+(hf)**2.)
894 hc=1.810*abs(deltat)**(1./3.)/(1.382+abs(cos(tiltangle)))
895 hdown=sqrt(hc**2.+(hf)**2.)
896 enerpv=effpv*rs*min(1.,1.-0.005*(tpv-(T_amb+273.15)))
897 sw_d=(1-albpv)*(rs)
898 lw_d=empv_up*rl
899 lwpv_out=empv_up*sigma*tpv**4.
900 lwupr=hrad*(tr**4-tpv**4.)
901 hspv=(hup+hdown)*(tpv-tair(n+1))
902 tpv=tpv+(sw_d+lw_d-lwpv_out+lwupr-hspv-enerpv)/Cm*dt
903 eppv=enerpv*pv_frac_roof*bl*bw
904 return
905 end subroutine hsfluxpv
906 !====6=8===============================================================72
907 !====6=8===============================================================72
909 subroutine wall_gr(hfgr,gr_frac_roof,swwall,nz,dt,dz,k,cs,flux,temp)
911 !______________________________________________________________________
913 !The aim of this subroutine is to solve the 1D heat fiffusion equation
914 !for roof, walls and streets:
915 !
916 ! dT/dt=d/dz[K*dT/dz] where:
917 !
918 ! -T is the surface temperature(wall, street, roof)
919 ! -Kz is the heat diffusivity inside the material.
920 !
921 !The resolution is done implicitly with a FV discretisation along the
922 !different layers of the material:
924 ! ____________________________
925 ! n *
926 ! *
927 ! *
928 ! ____________________________
929 ! i+2
930 ! I+1
931 ! ____________________________
932 ! i+1
933 ! I ==> [T(I,n+1)-T(I,n)]/DT=
934 ! ____________________________ [F(i+1)-F(i)]/DZI
935 ! i
936 ! I-1 ==> A*T(n+1)=B where:
937 ! ____________________________
938 ! i-1 * * A is a TRIDIAGONAL matrix.
939 ! * * B=T(n)+S takes into account the sources.
940 ! *
941 ! 1 ____________________________
943 !________________________________________________________________
945 implicit none
947 !Input:
948 !-----
949 real hfgr !Green roof heat flux
950 real gr_frac_roof !Green roof fraction
952 integer nz !Number of layers inside the material
953 real dt !Time step
954 real dz(nz) !Layer sizes [m]
955 real cs(nz) !Specific heat of the material [J/(m3.K)]
956 real k(nz+1) !Thermal conductivity in each layers (face) [W/(m.K)]
957 real flux(2) !Internal and external flux terms.
960 !Input-Output:
961 !-------------
963 integer swwall !swich for the physical coefficients calculation
964 real temp(nz) !Temperature at each layer
966 !Local:
967 !-----
969 real a(-1:1,nz) ! a(-1,*) lower diagonal A(i,i-1)
970 ! a(0,*) principal diagonal A(i,i)
971 ! a(1,*) upper diagonal A(i,i+1).
973 real b(nz) !Coefficients of the second term.
974 real k1(20)
975 real k2(20)
976 real kc(20)
977 save k1,k2,kc
978 integer iz
980 !________________________________________________________________
981 !
982 !Calculation of the coefficients
984 if (swwall.eq.1) then
986 if (nz.gt.20) then
987 write(*,*) 'number of layers in the walls/roofs too big ',nz
988 write(*,*) 'please decrease under of',20
989 stop
990 endif
992 call wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
993 swwall=0
995 end if
997 !Computation of the first value (iz=1) of A and B:
999 a(-1,1)=0.
1000 a(0,1)=1+k2(1)
1001 a(1,1)=-k2(1)
1002 b(1)=temp(1)+flux(1)*kc(1)
1004 do iz=2,nz-1
1005 a(-1,iz)=-k1(iz)
1006 if(iz.eq.5)then
1007 a(-1,iz)=-k1(iz)
1008 a(0,iz)=1+k1(iz)+(1-gr_frac_roof)*k2(iz)
1009 b(iz)=temp(iz)+(gr_frac_roof*hfgr*dt)/dz(iz)
1010 a(1,iz)=-k2(iz)*(1-gr_frac_roof)
1011 else
1012 a(0,iz)=1.+k1(iz)+k2(iz)
1013 b(iz)=temp(iz)
1014 a(1,iz)=-k2(iz)
1015 endif
1018 end do
1020 !Computation of the external value (iz=n) of A and B:
1022 a(-1,nz)=-k1(nz)
1023 a(0,nz)=1.+k1(nz)
1024 a(1,nz)=0.
1025 b(nz)=temp(nz)+kc(nz)*flux(2)
1027 !Resolution of the system A*T(n+1)=B
1029 call tridia(nz,a,b,temp)
1032 return
1033 end subroutine wall_gr
1035 !====6=8===============================================================72
1036 !====6=8===============================================================72
1039 subroutine wall(swwall,nz,dt,dz,k,cs,flux,temp)
1041 !______________________________________________________________________
1043 !The aim of this subroutine is to solve the 1D heat fiffusion equation
1044 !for roof, walls and streets:
1045 !
1046 ! dT/dt=d/dz[K*dT/dz] where:
1047 !
1048 ! -T is the surface temperature(wall, street, roof)
1049 ! -Kz is the heat diffusivity inside the material.
1050 !
1051 !The resolution is done implicitly with a FV discretisation along the
1052 !different layers of the material:
1054 ! ____________________________
1055 ! n *
1056 ! *
1057 ! *
1058 ! ____________________________
1059 ! i+2
1060 ! I+1
1061 ! ____________________________
1062 ! i+1
1063 ! I ==> [T(I,n+1)-T(I,n)]/DT=
1064 ! ____________________________ [F(i+1)-F(i)]/DZI
1065 ! i
1066 ! I-1 ==> A*T(n+1)=B where:
1067 ! ____________________________
1068 ! i-1 * * A is a TRIDIAGONAL matrix.
1069 ! * * B=T(n)+S takes into account the sources.
1070 ! *
1071 ! 1 ____________________________
1073 !________________________________________________________________
1075 implicit none
1077 !Input:
1078 !-----
1079 integer nz !Number of layers inside the material
1080 real dt !Time step
1081 real dz(nz) !Layer sizes [m]
1082 real cs(nz) !Specific heat of the material [J/(m3.K)]
1083 real k(nz+1) !Thermal conductivity in each layers (face) [W/(m.K)]
1084 real flux(2) !Internal and external flux terms.
1087 !Input-Output:
1088 !-------------
1090 integer swwall !swich for the physical coefficients calculation
1091 real temp(nz) !Temperature at each layer
1093 !Local:
1094 !-----
1096 real a(-1:1,nz) ! a(-1,*) lower diagonal A(i,i-1)
1097 ! a(0,*) principal diagonal A(i,i)
1098 ! a(1,*) upper diagonal A(i,i+1).
1100 real b(nz) !Coefficients of the second term.
1101 real k1(20)
1102 real k2(20)
1103 real kc(20)
1104 save k1,k2,kc
1105 integer iz
1107 !________________________________________________________________
1108 !
1109 !Calculation of the coefficients
1111 if (swwall.eq.1) then
1113 if (nz.gt.20) then
1114 write(*,*) 'number of layers in the walls/roofs too big ',nz
1115 write(*,*) 'please decrease under of',20
1116 stop
1117 endif
1119 call wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
1120 swwall=0
1122 end if
1124 !Computation of the first value (iz=1) of A and B:
1126 a(-1,1)=0.
1127 a(0,1)=1+k2(1)
1128 a(1,1)=-k2(1)
1129 b(1)=temp(1)+flux(1)*kc(1)
1131 do iz=2,nz-1
1132 a(-1,iz)=-k1(iz)
1133 a(0,iz)=1+k1(iz)+k2(iz)
1134 b(iz)=temp(iz)
1135 a(1,iz)=-k2(iz)
1136 end do
1138 !Computation of the external value (iz=n) of A and B:
1140 a(-1,nz)=-k1(nz)
1141 a(0,nz)=1+k1(nz)
1142 a(1,nz)=0.
1143 b(nz)=temp(nz)+flux(2)*kc(nz)
1145 !Resolution of the system A*T(n+1)=B
1147 call tridia(nz,a,b,temp)
1150 return
1151 end subroutine wall
1153 !====6=8===============================================================72
1154 !====6=8===============================================================72
1156 subroutine wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
1158 implicit none
1160 !---------------------------------------------------------------------
1161 !Input
1162 !-----
1163 integer nz !Number of layers inside the material
1164 real dt !Time step
1165 real dz(nz) !Layer sizes [m]
1166 real cs(nz) !Specific heat of the material [J/(m3.K)]
1167 real k(nz+1) !Thermal diffusivity in each layers (face) [W/(m.K)]
1170 !Input-Output
1171 !------------
1173 real flux(2) !Internal and external flux terms.
1176 !Output
1177 !------
1178 real k1(20)
1179 real k2(20)
1180 real kc(20)
1182 !Local
1183 !-----
1184 integer iz
1185 real kf(nz)
1187 !------------------------------------------------------------------
1189 do iz=2,nz
1190 kc(iz)=dt/(dz(iz)*cs(iz))
1191 kf(iz)=2*k(iz)/(dz(iz)+dz(iz-1))
1192 end do
1194 kc(1)=dt/(dz(1)*cs(1))
1195 kf(1)=2*k(1)/(dz(1))
1197 do iz=1,nz
1198 k1(iz)=kc(iz)*kf(iz)
1199 end do
1201 do iz=1,nz-1
1202 k2(iz)=kc(iz)*kf(iz+1)*cs(iz)/cs(iz+1)
1203 end do
1205 return
1206 end subroutine wall_coeff
1208 !====6=8===============================================================72
1209 !====6=8===============================================================72
1210 subroutine hsinsflux(swsurf,swwin,tin,tw,hsins)
1212 implicit none
1214 ! --------------------------------------------------------------------
1215 ! This routine computes the internal sensible heat flux.
1216 ! The swsurf, makes rhe difference between a vertical and a
1217 ! horizontal surface.
1218 ! The values of the heat conduction coefficients hc are obtained from the book
1219 ! "Energy Simulation in Building Design". J.A. Clarke.
1220 ! Adam Hilger, Bristol, 362 pp.
1221 ! --------------------------------------------------------------------
1222 !Input
1223 !----
1224 integer swsurf !swich for the type of surface (horizontal/vertical)
1225 integer swwin !swich for the type of surface (window/wall)
1226 real tin !Inside temperature [K]
1227 real tw !Internal wall temperature [K]
1230 !Output
1231 !------
1232 real hsins !internal sensible heat flux [W/m2]
1233 !Local
1234 !-----
1235 real hc !heat conduction coefficient [W/\B0C.m2]
1236 !--------------------------------------------------------------------
1238 if (swsurf.eq.2) then !vertical surface
1239 if (swwin.eq.1) then
1240 hc=5.678*0.99 !window surface (smooth surface)
1241 else
1242 hc=5.678*1.09 !wall surface (rough surface)
1243 endif
1244 hsins=hc*(tin-tw)
1245 endif
1247 if (swsurf.eq.1) then !horizontal surface
1248 if (swwin.eq.1) then
1249 hc=5.678*0.99 !window surface (smooth surface)
1250 else
1251 hc=5.678*1.09 !wall surface (rough surface)
1252 endif
1253 hsins=hc*(tin-tw)
1254 endif
1256 return
1257 end subroutine hsinsflux
1258 !====6=8===============================================================72
1259 !====6=8===============================================================72
1261 subroutine int_rsrad(albwin,albwal,pwin,rswal,&
1262 surwal,bw,bl,zw,rsint)
1264 ! ------------------------------------------------------------------
1265 implicit none
1266 ! ------------------------------------------------------------------
1268 !Input
1269 !-----
1270 real albwin !albedo of the windows
1271 real albwal !albedo of the internal wall
1272 real rswal(4) !incoming short wave radiation [W/m2]
1273 real surwal(6) !surface of the indoor walls [m2]
1274 real bw,bl !width of the walls [m]
1275 real zw !height of the wall [m]
1276 real pwin !window proportion
1278 !Output
1279 !------
1280 real rsint(6) !internal walls short wave radiation [W/m2]
1282 !Local
1283 !-----
1284 real transmit !transmittance of the direct/diffused radiation
1285 real rstr !solar radiation transmitted through the windows
1286 real surtotwal !total indoor surface of the walls in the room
1287 integer iw
1288 real b(6) !second member for the system
1289 real a(6,6) !matrix for the system
1291 !-------------------------------------------------------------------
1294 ! Calculation of the solar radiation transmitted through windows
1296 rstr = 0.
1297 do iw=1,4
1298 transmit=1.-albwin
1299 rstr = rstr+(surwal(iw)*pwin)*(transmit*rswal(iw))
1300 enddo
1302 !We suppose that the radiation is spread isotropically within the
1303 !room when it passes through the windows, so the flux [W/m2] in every
1304 !wall is:
1306 surtotwal=0.
1307 do iw=1,6
1308 surtotwal=surtotwal+surwal(iw)
1309 enddo
1311 rstr=rstr/surtotwal
1313 !Computation of the short wave radiation reaching the internal walls
1315 call algebra_short(rstr,albwal,albwin,bw,bl,zw,pwin,a,b)
1317 call gaussjbem(a,6,b,6)
1319 do iw=1,6
1320 rsint(iw)=b(iw)
1321 enddo
1323 return
1324 end subroutine int_rsrad
1326 !====6=8===============================================================72
1327 !====6=8===============================================================72
1329 subroutine int_rlrad(emwal,emwin,sigma,twal_int,twin,&
1330 pwin,bw,bl,zw,rlint)
1332 ! ------------------------------------------------------------------
1333 implicit none
1334 ! ------------------------------------------------------------------
1336 !Input
1337 !-----
1339 real emwal !emissivity of the internal walls
1340 real emwin !emissivity of the window
1341 real sigma !Stefan-Boltzmann constant [W/m2.K4]
1342 real twal_int(6)!temperature of the first internal layers of a room [K]
1343 real twin(4) !temperature of the windows [K]
1344 real bw !width of the wall
1345 real bl !length of the wall
1346 real zw !height of the wall
1347 real pwin !window proportion
1349 !Output
1350 !------
1352 real rlint(6) !internal walls long wave radiation [W/m2]
1354 !Local
1355 !------
1357 real b(6) !second member vector for the system
1358 real a(6,6) !matrix for the system
1359 integer iw
1360 !----------------------------------------------------------------
1362 !Compute the long wave radiation reachs the internal walls
1364 call algebra_long(emwal,emwin,sigma,twal_int,twin,pwin,&
1365 bw,bl,zw,a,b)
1367 call gaussjbem(a,6,b,6)
1369 do iw=1,6
1370 rlint(iw)=b(iw)
1371 enddo
1373 return
1374 end subroutine int_rlrad
1376 !====6=8===============================================================72
1377 !====6=8===============================================================72
1379 subroutine algebra_short(rstr,albwal,albwin,aw,bw,zw,pwin,a,b)
1381 !--------------------------------------------------------------------
1382 !This routine calculates the algebraic system that will be solved for
1383 !the computation of the total shortwave radiation that reachs every
1384 !indoor wall in a floor.
1385 !Write the matrix system ax=b to solve
1386 !
1387 ! -Rs(1)+a(1,2)Rs(2)+.................+a(1,6)Rs(6)=-Rs=b(1)
1388 !a(2,1)Rs(1)- Rs(2)+.................+a(2,6)Rs(6)=-Rs=b(2)
1389 !a(3,1)Rs(1)+a(3,2)Rs(3)-Rs(3)+...........+a(3,6)Rs(6)=-Rs=b(3)
1390 !a(4,1)Rs(1)+.................-Rs(4)+.....+a(4,6)Rs(6)=-Rs=b(4)
1391 !a(5,1)Rs(1)+.......................-Rs(5)+a(5,6)Rs(6)=-Rs=b(5)
1392 !a(6,1)Rs(1)+....................................-R(6)=-Rs=b(6)
1393 !
1394 !This version suppose the albedo of the indoor walls is the same.
1395 !--------------------------------------------------------------------
1396 implicit none
1397 !--------------------------------------------------------------------
1399 !Input
1400 !-----
1401 real rstr !solar radiation transmitted through the windows
1402 real albwal !albedo of the internal walls
1403 real albwin !albedo of the windows.
1404 real bw !length of the wall
1405 real aw !width of the wall
1406 real zw !height of the wall
1407 real fprl_int !view factor
1408 real fnrm_int !view factor
1409 real pwin !window proportion
1410 !Output
1411 !------
1412 real a(6,6) !Matrix for the system
1413 real b(6) !Second member for the system
1414 !Local
1415 !-----
1416 integer iw,jw
1417 real albm !averaged albedo
1418 !----------------------------------------------------------------
1420 !Initialise the variables
1422 do iw=1,6
1423 b(iw)= 0.
1424 do jw=1,6
1425 a(iw,jw)= 0.
1426 enddo
1427 enddo
1429 !Calculation of the second member b
1431 do iw=1,6
1432 b(iw)=-rstr
1433 end do
1435 !Calculation of the averaged albedo
1437 albm=pwin*albwin+(1-pwin)*albwal
1439 !Calculation of the matrix a
1441 a(1,1)=-1.
1443 call fprl_ints(fprl_int,aw/bw,zw/bw)
1445 a(1,2)=albm*fprl_int
1447 call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1449 a(1,3)=albm*(bw/aw)*fnrm_int
1451 a(1,4)=a(1,3)
1453 call fnrm_ints(fnrm_int,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1455 a(1,5)=albwal*(bw/zw)*fnrm_int
1457 a(1,6)=a(1,5)
1460 a(2,1)=a(1,2)
1461 a(2,2)=-1.
1462 a(2,3)=a(1,3)
1463 a(2,4)=a(1,4)
1464 a(2,5)=a(1,5)
1465 a(2,6)=a(1,6)
1468 call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1470 a(3,1)=albm*(aw/bw)*fnrm_int
1471 a(3,2)=a(3,1)
1472 a(3,3)=-1.
1474 call fprl_ints(fprl_int,zw/aw,bw/aw)
1476 a(3,4)=albm*fprl_int
1478 call fnrm_ints(fnrm_int,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1480 a(3,5)=albwal*(aw/zw)*fnrm_int
1481 a(3,6)=a(3,5)
1484 a(4,1)=a(3,1)
1485 a(4,2)=a(3,2)
1486 a(4,3)=a(3,4)
1487 a(4,4)=-1.
1488 a(4,5)=a(3,5)
1489 a(4,6)=a(3,6)
1491 call fnrm_ints(fnrm_int,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1493 a(5,1)=albm*(zw/bw)*fnrm_int
1495 a(5,2)=a(5,1)
1497 call fnrm_ints(fnrm_int,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1499 a(5,3)=albm*(zw/aw)*fnrm_int
1501 a(5,4)=a(5,3)
1502 a(5,5)=-1.
1504 call fprl_ints(fprl_int,aw/zw,bw/zw)
1506 a(5,6)=albwal*fprl_int
1509 a(6,1)=a(5,1)
1510 a(6,2)=a(5,2)
1511 a(6,3)=a(5,3)
1512 a(6,4)=a(5,4)
1513 a(6,5)=a(5,6)
1514 a(6,6)=-1.
1516 return
1517 end subroutine algebra_short
1519 !====6=8===============================================================72
1520 !====6=8===============================================================72
1522 subroutine algebra_long(emwal,emwin,sigma,twalint,twinint,&
1523 pwin,aw,bw,zw,a,b)
1525 !--------------------------------------------------------------------
1526 !This routine computes the algebraic system that will be solved to
1527 !compute the longwave radiation that reachs the indoor
1528 !walls in a floor.
1529 !Write the matrix system ax=b to solve
1530 !
1531 !a(1,1)Rl(1)+.............................+Rl(6)=b(1)
1532 !a(2,1)Rl(1)+.................+Rl(5)+a(2,6)Rl(6)=b(2)
1533 !a(3,1)Rl(1)+.....+Rl(3)+...........+a(3,6)Rl(6)=b(3)
1534 !a(4,1)Rl(1)+...........+Rl(4)+.....+a(4,6)Rl(6)=b(4)
1535 ! Rl(1)+.......................+a(5,6)Rl(6)=b(5)
1536 !a(6,1)Rl(1)+Rl(2)+.................+a(6,6)Rl(6)=b(6)
1537 !
1538 !--------------------------------------------------------------------
1539 implicit none
1541 !--------------------------------------------------------------------
1543 !Input
1544 !-----
1546 real pwin !window proportion
1547 real emwal !emissivity of the internal walls
1548 real emwin !emissivity of the window
1549 real sigma !Stefan-Boltzmann constant [W/m2.K4]
1550 real twalint(6) !temperature of the first internal layers of a room [K]
1551 real twinint(4) !temperature of the windows [K]
1552 real aw !width of the wall
1553 real bw !length of the wall
1554 real zw !height of the wall
1555 real fprl_int !view factor
1556 real fnrm_int !view factor
1557 real fnrm_intx !view factor
1558 real fnrm_inty !view factor
1560 !Output
1561 !------
1562 real b(6) !second member vector for the system
1563 real a(6,6) !matrix for the system
1564 !Local
1565 !-----
1566 integer iw,jw
1567 real b_wall(6)
1568 real b_wind(6)
1569 real emwal_av !averadge emissivity of the wall
1570 real emwin_av !averadge emissivity of the window
1571 real em_av !averadge emissivity
1572 real twal_int(6) !twalint
1573 real twin(4) !twinint
1574 !------------------------------------------------------------------
1576 !Initialise the variables
1577 !-------------------------
1579 do iw=1,6
1580 b(iw)= 0.
1581 b_wall(iw)=0.
1582 b_wind(iw)=0.
1583 do jw=1,6
1584 a(iw,jw)= 0.
1585 enddo
1586 enddo
1588 do iw=1,6
1589 twal_int(iw)=twalint(iw)
1590 enddo
1592 do iw=1,4
1593 twin(iw)=twinint(iw)
1594 enddo
1596 !Calculation of the averadge emissivities
1597 !-----------------------------------------
1599 emwal_av=(1-pwin)*emwal
1600 emwin_av=pwin*emwin
1601 em_av=emwal_av+emwin_av
1603 !Calculation of the second term for the walls
1604 !-------------------------------------------
1606 call fprl_ints(fprl_int,aw/zw,bw/zw)
1607 call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1608 call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1610 b_wall(1)=(emwal*sigma*(twal_int(5)**4)* &
1611 fprl_int)+ &
1612 (sigma*(emwal_av*(twal_int(3)**4)+ &
1613 emwal_av*(twal_int(4)**4))* &
1614 (zw/aw)*fnrm_intx)+ &
1615 (sigma*(emwal_av*(twal_int(1)**4)+ &
1616 emwal_av*(twal_int(2)**4))* &
1617 (zw/bw)*fnrm_inty)
1619 call fprl_ints(fprl_int,aw/zw,bw/zw)
1620 call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1621 call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1623 b_wall(2)=(emwal*sigma*(twal_int(6)**4)* &
1624 fprl_int)+ &
1625 (sigma*(emwal_av*(twal_int(3)**4)+ &
1626 emwal_av*(twal_int(4)**4))* &
1627 (zw/aw)*fnrm_intx)+ &
1628 (sigma*(emwal_av*(twal_int(1)**4)+ &
1629 emwal_av*(twal_int(2)**4))* &
1630 (zw/bw)*fnrm_inty)
1632 call fprl_ints(fprl_int,zw/aw,bw/aw)
1633 call fnrm_ints(fnrm_intx,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1634 call fnrm_ints(fnrm_inty,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1636 b_wall(3)=(emwal_av*sigma*(twal_int(4)**4)* &
1637 fprl_int)+ &
1638 (sigma*(emwal_av*(twal_int(2)**4)+ &
1639 emwal_av*(twal_int(1)**4))* &
1640 (aw/bw)*fnrm_intx)+ &
1641 (sigma*(emwal*(twal_int(5)**4)+ &
1642 emwal*(twal_int(6)**4))* &
1643 (aw/zw)*fnrm_inty)
1645 call fprl_ints(fprl_int,zw/aw,bw/aw)
1646 call fnrm_ints(fnrm_intx,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1647 call fnrm_ints(fnrm_inty,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1649 b_wall(4)=(emwal_av*sigma*(twal_int(3)**4)* &
1650 fprl_int)+ &
1651 (sigma*(emwal_av*(twal_int(2)**4)+ &
1652 emwal_av*(twal_int(1)**4))* &
1653 (aw/bw)*fnrm_intx)+ &
1654 (sigma*(emwal*(twal_int(5)**4)+ &
1655 emwal*(twal_int(6)**4))* &
1656 (aw/zw)*fnrm_inty)
1658 call fprl_ints(fprl_int,aw/bw,zw/bw)
1659 call fnrm_ints(fnrm_intx,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1660 call fnrm_ints(fnrm_inty,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1662 b_wall(5)=(emwal_av*sigma*(twal_int(2)**4)* &
1663 fprl_int)+ &
1664 (sigma*(emwal_av*(twal_int(3)**4)+ &
1665 emwal_av*(twal_int(4)**4))* &
1666 (bw/aw)*fnrm_intx)+ &
1667 (sigma*(emwal*(twal_int(5)**4)+ &
1668 emwal*(twal_int(6)**4))* &
1669 (bw/zw)*fnrm_inty)
1671 call fprl_ints(fprl_int,aw/bw,zw/bw)
1672 call fnrm_ints(fnrm_intx,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1673 call fnrm_ints(fnrm_inty,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1675 b_wall(6)=(emwal_av*sigma*(twal_int(1)**4)* &
1676 fprl_int)+ &
1677 (sigma*(emwal_av*(twal_int(3)**4)+ &
1678 emwal_av*(twal_int(4)**4))* &
1679 (bw/aw)*fnrm_intx)+ &
1680 (sigma*(emwal*(twal_int(5)**4)+ &
1681 emwal*(twal_int(6)**4))* &
1682 (bw/zw)*fnrm_inty)
1684 !Calculation of the second term for the windows
1685 !---------------------------------------------
1686 call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1687 call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1689 b_wind(1)=(sigma*(emwin_av*(twin(3)**4)+ &
1690 emwin_av*(twin(4)**4))* &
1691 (zw/aw)*fnrm_intx)+ &
1692 (sigma*(emwin_av*(twin(1)**4)+ &
1693 emwin_av*(twin(2)**4))* &
1694 (zw/bw)*fnrm_inty)
1696 call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1697 call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1699 b_wind(2)=(sigma*(emwin_av*(twin(3)**4)+ &
1700 emwin_av*(twin(4)**4))* &
1701 (zw/aw)*fnrm_intx)+ &
1702 (sigma*(emwin_av*(twin(1)**4)+ &
1703 emwin_av*(twin(2)**4))* &
1704 (zw/bw)*fnrm_inty)
1706 call fprl_ints(fprl_int,zw/aw,bw/aw)
1707 call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1709 b_wind(3)=emwin_av*sigma*(twin(4)**4)* &
1710 fprl_int+(sigma*(emwin_av* &
1711 (twin(2)**4)+emwin_av*(twin(1)**4))* &
1712 (aw/bw)*fnrm_int)
1714 call fprl_ints(fprl_int,zw/aw,bw/aw)
1715 call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1717 b_wind(4)=emwin_av*sigma*(twin(3)**4)* &
1718 fprl_int+(sigma*(emwin_av* &
1719 (twin(2)**4)+emwin_av*(twin(1)**4))* &
1720 (aw/bw)*fnrm_int)
1722 call fprl_ints(fprl_int,aw/bw,zw/bw)
1723 call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1725 b_wind(5)=emwin_av*sigma*(twin(2)**4)* &
1726 fprl_int+(sigma*(emwin_av* &
1727 (twin(3)**4)+emwin_av*(twin(4)**4))* &
1728 (bw/aw)*fnrm_int)
1730 call fprl_ints(fprl_int,aw/bw,zw/bw)
1731 call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1733 b_wind(6)=emwin_av*sigma*(twin(1)**4)* &
1734 fprl_int+(sigma*(emwin_av* &
1735 (twin(3)**4)+emwin_av*(twin(4)**4))* &
1736 (bw/aw)*fnrm_int)
1738 !Calculation of the total b term
1739 !-------------------------------
1741 do iw=1,6
1742 b(iw)=b_wall(iw)+b_wind(iw)
1743 end do
1746 !Calculation of the matrix of the system
1747 !----------------------------------------
1749 call fnrm_ints(fnrm_int,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1751 a(1,1)=(em_av-1.)*(zw/bw)*fnrm_int
1753 a(1,2)=a(1,1)
1755 call fnrm_ints(fnrm_int,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1757 a(1,3)=(em_av-1.)*(zw/aw)*fnrm_int
1759 a(1,4)=a(1,3)
1761 call fprl_ints(fprl_int,aw/zw,bw/zw)
1763 a(1,5)=(emwal-1.)*fprl_int
1764 a(1,6)=1.
1766 a(2,1)=a(1,1)
1767 a(2,2)=a(1,2)
1768 a(2,3)=a(1,3)
1769 a(2,4)=a(1,4)
1770 a(2,5)=1.
1771 a(2,6)=a(1,5)
1773 call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1775 a(3,1)=(em_av-1.)*(aw/bw)*fnrm_int
1777 a(3,2)=a(3,1)
1778 a(3,3)=1.
1780 call fprl_ints(fprl_int,zw/aw,bw/aw)
1782 a(3,4)=(em_av-1.)*fprl_int
1784 call fnrm_ints(fnrm_int,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1786 a(3,5)=(emwal-1.)*(aw/zw)*fnrm_int
1788 a(3,6)=a(3,5)
1790 a(4,1)=a(3,1)
1791 a(4,2)=a(3,2)
1792 a(4,3)=a(3,4)
1793 a(4,4)=1.
1794 a(4,5)=a(3,5)
1795 a(4,6)=a(3,6)
1797 a(5,1)=1.
1799 call fprl_ints(fprl_int,aw/bw,zw/bw)
1801 a(5,2)=(em_av-1.)*fprl_int
1803 call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1805 a(5,3)=(em_av-1.)*(bw/aw)*fnrm_int
1807 a(5,4)=a(5,3)
1809 call fnrm_ints(fnrm_int,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1811 a(5,5)=(emwal-1.)*(bw/zw)*fnrm_int
1813 a(5,6)=a(5,5)
1815 a(6,1)=a(5,2)
1816 a(6,2)=1.
1817 a(6,3)=a(5,3)
1818 a(6,4)=a(5,4)
1819 a(6,5)=a(5,5)
1820 a(6,6)=a(6,5)
1822 return
1823 end subroutine algebra_long
1825 !====6=8===============================================================72
1826 !====6=8===============================================================72
1829 subroutine fluxroo(hseqocc,hleqocc,hsvent,hlvent, &
1830 hscond,hslev,hllev)
1832 !-----------------------------------------------------------------------
1833 !This routine calculates the heat flux generated inside the room
1834 !and the heat ejected to the atmosphere.
1835 !----------------------------------------------------------------------
1837 implicit none
1839 !-----------------------------------------------------------------------
1841 !Input
1842 !-----
1843 real hseqocc !sensible heat generated by equipments and occupants [W]
1844 real hleqocc !latent heat generated by occupants [W]
1845 real hsvent !sensible heat generated by natural ventilation [W]
1846 real hlvent !latent heat generated by natural ventilation [W]
1847 real hscond !sensible heat generated by wall conduction
1849 !Output
1850 !------
1851 real hslev !sensible heat flux generated inside the room [W]
1852 real hllev !latent heat flux generatd inside the room
1855 !Calculation of the total sensible heat generated inside the room
1857 hslev=hseqocc+hsvent+hscond
1859 !Calculation of the total latent heat generated inside the room
1861 hllev=hleqocc+hlvent
1863 return
1864 end subroutine fluxroo
1866 !====6=8===============================================================72
1867 !====6=8===============================================================72
1869 subroutine phirat(nhourday,rocc)
1871 !----------------------------------------------------------------------
1872 !This routine calculates the occupation ratio of a floor
1873 !By now we suppose a constant value
1874 !----------------------------------------------------------------------
1876 implicit none
1878 !Input
1879 !-----
1881 real nhourday ! number of hours from midnight (local time)
1883 !Output
1884 !------
1885 real rocc !value between 0 and 1
1887 !!TEST
1888 rocc=1.
1890 return
1891 end subroutine phirat
1893 !====6=8===============================================================72
1894 !====6=8===============================================================72
1896 subroutine phiequ(nhourday,hsesf,hsequip,hsequ)
1898 !----------------------------------------------------------------------
1899 !This routine calculates the sensible heat gain from equipments
1900 !----------------------------------------------------------------------
1901 implicit none
1902 !Input
1903 !-----
1905 real nhourday ! number of hours from midnight, Local time
1906 real, intent(in) :: hsesf
1907 real, intent(in), dimension(24) :: hsequip
1909 !Output
1910 !------
1911 real hsequ !sensible heat gain from equipment [Wm\AF2]
1913 !---------------------------------------------------------------------
1915 hsequ = hsequip(int(nhourday)+1) * hsesf
1917 end subroutine phiequ
1918 !====6=8===============================================================72
1919 !====6=8===============================================================72
1921 subroutine fluxeqocc(nhourday,bw,bl,perflo,hsesf,hsequip,hseqocc,hleqocc)
1923 implicit none
1925 !---------------------------------------------------------------------
1926 !This routine calculates the sensible and the latent heat flux
1927 !generated by equipments and occupants
1928 !---------------------------------------------------------------------
1930 !Input
1931 !-----
1932 real bw !Room width [m]
1933 real bl !Room lengzh [m]
1934 real nhourday !number of hours from the beginning of the day
1935 real, intent(in) :: perflo ! Peak number of occupants per unit floor area
1936 real, intent(in) :: hsesf
1937 real, intent(in), dimension(24) :: hsequip
1939 !Output
1940 !------
1941 real hseqocc !sensible heat generated by equipments and occupants [W]
1942 real hleqocc !latent heat generated by occupants [W]
1943 !Local
1944 !-----
1945 real Af !Air conditioned floor area [m2]
1946 real rocc !Occupation ratio of the floor [0,1]
1947 real hsequ !Heat generated from equipments
1949 real hsocc !Sensible heat generated by a person [W/Person]
1950 !Source Boundary Layer Climates,page 195 (book)
1951 parameter (hsocc=160.)
1953 real hlocc !Latent heat generated by a person [W/Person]
1954 !Source Boundary Layer Climates,page 225 (book)
1955 parameter (hlocc=1.96e6/86400.)
1957 !------------------------------------------------------------------
1958 ! Sensible heat flux
1959 ! ------------------
1961 Af=bw*bl
1963 call phirat(nhourday,rocc)
1965 call phiequ(nhourday,hsesf,hsequip,hsequ)
1967 hseqocc=Af*rocc*perflo*hsocc+Af*hsequ
1969 !
1970 ! Latent heat
1971 ! -----------
1972 !
1974 hleqocc=Af*rocc*perflo*hlocc
1976 return
1977 end subroutine fluxeqocc
1979 !====6=8===============================================================72
1980 !====6=8===============================================================72
1982 subroutine fluxvent(cpint,rhoint,vollev,tlev,tout,latent,&
1983 humout,rhoout,humlev,beta,hsvent,hlvent)
1985 implicit none
1987 !---------------------------------------------------------------------
1988 !This routine calculates the sensible and the latent heat flux
1989 !generated by natural ventilation
1990 !---------------------------------------------------------------------
1992 !Input
1993 !-----
1994 real cpint !specific heat of the indoor air [J/kg.K]
1995 real rhoint !density of the indoor air [Kg/m3]
1996 real vollev !volume of the room [m3]
1997 real tlev !Room temperature [K]
1998 real tout !outside air temperature [K]
1999 real latent !latent heat of evaporation [J/Kg]
2000 real humout !outside absolute humidity [Kgwater/Kgair]
2001 real rhoout !air density [kg/m3]
2002 real humlev !Specific humidity of the indoor air [Kgwater/Kgair]
2003 real, intent(in) :: beta!Thermal efficiency of the heat exchanger
2005 !Output
2006 !------
2007 real hsvent !sensible heat generated by natural ventilation [W]
2008 real hlvent !latent heat generated by natural ventilation [W]
2010 !Local
2011 !-----
2013 !----------------------------------------------------------------------
2015 ! Sensible heat flux
2016 ! ------------------
2018 hsvent=(1.-beta)*cpint*rhoint*(vollev/3600.)* &
2019 (tout-tlev)
2021 ! Latent heat flux
2022 ! ----------------
2024 hlvent=(1.-beta)*latent*rhoint*(vollev/3600.)* &
2025 (humout-humlev)
2028 return
2029 end subroutine fluxvent
2031 !====6=8===============================================================72
2032 !====6=8===============================================================72
2034 subroutine fluxcond(hswalins,hswinins,surwal,pwin,hscond)
2036 implicit none
2038 !---------------------------------------------------------------------
2039 !This routine calculates the sensible heat flux generated by
2040 !wall conduction.
2041 !---------------------------------------------------------------------
2043 !Input
2044 !-----
2045 real hswalins(6) !sensible heat at the internal layers of the wall [W/m2]
2046 real hswinins(4) !internal window sensible heat flux [W/m2]
2047 real surwal(6) !surfaces of the room walls [m2]
2048 real pwin !window proportion
2051 !Output
2052 !------
2054 real hscond !sensible heat generated by wall conduction [W]
2056 !Local
2057 !-----
2059 integer ivw
2061 !----------------------------------------------------------------------
2063 hscond=0.
2065 do ivw=1,4
2066 hscond=hscond+surwal(ivw)*(1-pwin)*hswalins(ivw)+ &
2067 surwal(ivw)*pwin*hswinins(ivw)
2068 end do
2070 do ivw=5,6
2071 hscond=hscond+surwal(ivw)*hswalins(ivw)
2072 end do
2073 !
2074 !Finally we must change the sign in hscond to be proportional
2075 !to the difference (Twall-Tindoor).
2076 !
2077 hscond=(-1)*hscond
2079 return
2080 end subroutine fluxcond
2082 !====6=8===============================================================72
2083 !====6=8===============================================================72
2085 subroutine regtemp(swcond,nhourday,dt,Qb,hsroo, &
2086 tlev,timeon,timeoff,targtemp,gaptemp,hsneed)
2088 implicit none
2090 !---------------------------------------------------------------------
2091 !This routine calculates the sensible heat fluxes,
2092 !after anthropogenic regulation (air conditioning)
2093 !---------------------------------------------------------------------
2095 !Input:
2096 !-----.
2097 integer swcond !swich air conditioning
2098 real nhourday !number of hours from the beginning of the day real
2099 real dt !time step [s]
2100 real Qb !overall heat capacity of the indoor air [J/K]
2101 real hsroo !sensible heat flux generated inside the room [W]
2102 real tlev !room air temperature [K]
2103 real, intent(in) :: timeon ! Initial local time of A/C systems
2104 real, intent(in) :: timeoff ! Ending local time of A/C systems
2105 real, intent(in) :: targtemp! Target temperature of A/C systems
2106 real, intent(in) :: gaptemp ! Comfort range of indoor temperature
2109 !Local:
2110 !-----.
2112 real templev !hipotetical room air temperature [K]
2113 real alpha !variable to control the heating/cooling of
2114 !the air conditining system
2115 !Output:
2116 !-----.
2117 real hsneed !sensible heat extracted to the indoor air [W]
2118 !---------------------------------------------------------------------
2119 !initialize variables
2120 !---------------------
2121 templev = 0.
2122 alpha = 0.
2124 if (swcond.eq.0) then ! there is not air conditioning in the floor
2125 hsneed = 0.
2126 goto 100
2127 else
2128 if ((nhourday.ge.timeon).and.(nhourday.le.timeoff)) then
2129 templev=tlev+(dt/Qb)*hsroo
2130 goto 200
2131 else
2132 hsneed = 0. ! air conditioning is switched off
2133 goto 100
2134 endif
2135 endif
2137 200 continue
2139 if (abs(templev-targtemp).le.gaptemp) then
2140 hsneed = 0.
2141 else
2142 if (templev.gt.(targtemp+gaptemp)) then
2143 hsneed=hsroo-(Qb/dt)*(targtemp+gaptemp-tlev)
2144 alpha=(abs(hsneed-hsroo)/Qb)
2145 if (alpha.gt.temp_rat) then
2146 hsneed=hsroo+temp_rat*Qb
2147 goto 100
2148 else
2149 goto 100
2150 endif
2151 else
2152 hsneed=hsroo-(Qb/dt)*(targtemp-gaptemp-tlev)
2153 alpha=(abs(hsneed-hsroo)/Qb)
2154 if (alpha.gt.temp_rat) then
2155 hsneed=hsroo-temp_rat*Qb
2156 goto 100
2157 else
2158 goto 100
2159 endif
2160 endif
2161 endif
2163 100 continue
2164 return
2165 end subroutine regtemp
2167 !====6=8==============================================================72
2168 !====6=8==============================================================72
2170 subroutine reghum(swcond,nhourday,dt,volroo,rhoint,latent, &
2171 hlroo,shumroo,timeon,timeoff,targhum,gaphum,hlneed)
2173 implicit none
2175 !---------------------------------------------------------------------
2176 !This routine calculates the latent heat fluxes,
2177 !after anthropogenic regulation (air conditioning)
2178 !---------------------------------------------------------------------
2180 !Input:
2181 !-----.
2182 integer swcond !swich air conditioning
2183 real nhourday !number of hours from the beginning of the day real[h]
2184 real dt !time step [s]
2185 real volroo !volume of the room [m3]
2186 real rhoint !density of the internal air [Kg/m3]
2187 real latent !latent heat of evaporation [J/Kg]
2188 real hlroo !latent heat flux generated inside the room [W]
2189 real shumroo !specific humidity of the indoor air [kg/kg]
2190 real, intent(in) :: timeon ! Initial local time of A/C systems
2191 real, intent(in) :: timeoff ! Ending local time of A/C systems
2192 real, intent(in) :: targhum ! Target humidity of the A/C systems
2193 real, intent(in) :: gaphum ! comfort range of the specific humidity
2195 !Local:
2196 !-----.
2198 real humlev !hipotetical specific humidity of the indoor [kg/kg]
2199 real betha !variable to control the drying/moistening of
2200 !the air conditioning system
2201 !Output:
2202 !-----.
2203 real hlneed !latent heat extracted to the indoor air [W]
2204 !------------------------------------------------------------------------
2205 !initialize variables
2206 !---------------------
2207 humlev = 0.
2208 betha = 0.
2210 if (swcond.eq.0) then ! there is not air conditioning in the floor
2211 hlneed = 0.
2212 goto 100
2213 else
2214 if ((nhourday.ge.timeon).and.(nhourday.le.timeoff)) then
2215 humlev=shumroo+(dt/(latent*rhoint*volroo))*hlroo
2216 goto 200
2217 else
2218 hlneed = 0. ! air conditioning is switched off
2219 goto 100
2220 endif
2221 endif
2223 200 continue
2225 if (abs(humlev-targhum).le.gaphum) then
2226 hlneed = 0.
2227 else
2228 if (humlev.gt.(targhum+gaphum)) then
2229 hlneed=hlroo-((latent*rhoint*volroo)/dt)* &
2230 (targhum+gaphum-shumroo)
2231 betha=abs(hlneed-hlroo)/(latent*rhoint*volroo)
2232 if (betha.gt.hum_rat) then
2233 hlneed=hlroo+hum_rat*(latent*rhoint*volroo)
2234 goto 100
2235 else
2236 goto 100
2237 endif
2238 else
2239 hlneed=hlroo-((latent*rhoint*volroo)/dt)* &
2240 (targhum-gaphum-shumroo)
2241 betha=abs(hlneed-hlroo)/(latent*rhoint*volroo)
2242 if (betha.gt.hum_rat) then
2243 hlneed=hlroo-hum_rat*(latent*rhoint*volroo)
2244 goto 100
2245 else
2246 goto 100
2247 endif
2248 endif
2249 endif
2251 100 continue
2252 return
2253 end subroutine reghum
2255 !====6=8==============================================================72
2256 !====6=8==============================================================72
2258 subroutine air_cond(hsneed,hlneed,dt,hsout,hlout,consump,cop)
2260 implicit none
2262 !
2263 !Performance of the air conditioning system
2264 !
2265 !INPUT/OUTPUT VARIABLES
2266 real, intent(in) :: cop
2267 !
2268 !INPUT/OUTPUT VARIABLES
2269 !
2270 real hsneed !sensible heat that is necessary for cooling/heating
2271 !the indoor air temperature [W]
2272 real hlneed !latent heat that is necessary for controling
2273 !the humidity of the indoor air [W]
2274 real dt !time step [s]
2275 !
2276 !OUTPUT VARIABLES
2277 !
2278 real hsout !sensible heat pumped out into the atmosphere [W]
2279 real hlout !latent heat pumped out into the atmosphere [W]
2280 real consump !Electrical consumption of the air conditioning system [W]
2283 !
2284 !Performance of the air conditioning system
2285 !
2286 if (hsneed.gt.0) then ! air conditioning is cooling
2287 ! and the heat is pumped out into the atmosphere
2288 hsout=(1/cop)*(abs(hsneed)+abs(hlneed))+hsneed
2289 hlout=hlneed
2290 consump=(1./cop)*(abs(hsneed)+abs(hlneed))
2291 !! hsout=0.
2292 !! hlout=0.
2294 else if(hsneed.eq.0.) then !air conditioning is not working to regulate the indoor temperature
2295 hlneed=0. !no humidity regulation is considered
2296 hsout=0. !no output into the atmosphere (sensible heat)
2297 hlout=0. !no output into the atmosphere (latent heat)
2298 consump=0. !no electrical consumption
2300 else !! hsneed < 0. !air conditioning is heating
2301 hlneed=0. !no humidity regulation is considered
2302 hlout=0. !no output into the atmosphere (latent heat)
2303 consump=(1./cop)*(abs(hsneed)+abs(hlneed))
2304 !
2305 !!We have two possibilities
2306 !
2307 !! hsout=(1./cop)*(abs(hsneed)+abs(hlneed)) !output into the atmosphere
2308 hsout=0. !no output into the atmosphere
2309 end if
2311 return
2312 end subroutine air_cond
2314 !====6=8==============================================================72
2315 !====6=8==============================================================72
2317 subroutine consump_total(nzcanm,nlev,consumpbuild,hsoutbuild, &
2318 hsout,consump)
2320 implicit none
2322 !-----------------------------------------------------------------------
2323 !Compute the total consumption in kWh/s (1kWh=3.6e+6 J) and sensible heat
2324 !ejected into the atmosphere per building
2325 !------------------------------------------------------------------------
2326 !
2327 !INPUT VARIABLES
2328 !
2329 !
2330 integer nzcanm !Maximum number of vertical levels in the urban grid
2331 real hsout(nzcanm) !sensible heat emitted outside the room [W]
2332 real consump(nzcanm) !Electricity consumption for the a.c. in each floor[W]
2333 !
2334 !OUTPUT VARIABLES
2335 !
2336 real consumpbuild !Energetic consumption for the entire building[kWh/s]
2337 real hsoutbuild !Total sensible heat ejected into the atmosphere
2338 !by the air conditioning systems per building [W]
2339 !
2340 !LOCAL VARIABLES
2341 !
2342 integer ilev
2344 !
2345 !INPUT VARIABLES
2346 !
2347 integer nlev
2349 !
2350 !INITIALIZE VARIABLES
2351 !
2352 consumpbuild=0.
2353 hsoutbuild=0.
2354 !
2355 do ilev=1,nlev
2356 consumpbuild=consumpbuild+consump(ilev)
2357 hsoutbuild=hsoutbuild+hsout(ilev)
2358 enddo !ilev
2360 consumpbuild=consumpbuild/(3.6e+06)
2362 return
2363 end subroutine consump_total
2364 !====6=8==============================================================72
2365 !====6=8==============================================================72
2366 subroutine tridia(n,a,b,x)
2368 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2369 ! + by A. Clappier, EPFL, CH 1015 Lausanne +
2370 ! + phone: ++41-(0)21-693-61-60 +
2371 ! + email:alain.clappier@epfl.ch +
2372 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2374 ! ----------------------------------------------------------------------
2375 ! Resolution of a * x = b where a is a tridiagonal matrix
2376 !
2377 ! ----------------------------------------------------------------------
2379 implicit none
2381 ! Input
2382 integer n
2383 real a(-1:1,n) ! a(-1,*) lower diagonal A(i,i-1)
2384 ! a(0,*) principal diagonal A(i,i)
2385 ! a(1,*) upper diagonal A(i,i+1)
2386 real b(n)
2388 ! Output
2389 real x(n)
2391 ! Local
2392 integer i
2394 ! ----------------------------------------------------------------------
2396 do i=n-1,1,-1
2397 b(i)=b(i)-a(1,i)*b(i+1)/a(0,i+1)
2398 a(0,i)=a(0,i)-a(1,i)*a(-1,i+1)/a(0,i+1)
2399 enddo
2401 do i=2,n
2402 b(i)=b(i)-a(-1,i)*b(i-1)/a(0,i-1)
2403 enddo
2405 do i=1,n
2406 x(i)=b(i)/a(0,i)
2407 enddo
2409 return
2410 end subroutine tridia
2411 !====6=8===============================================================72
2412 !====6=8===============================================================72
2414 subroutine gaussjbem(a,n,b,np)
2416 ! ----------------------------------------------------------------------
2417 ! This routine solve a linear system of n equations of the form
2418 ! A X = B
2419 ! where A is a matrix a(i,j)
2420 ! B a vector and X the solution
2421 ! In output b is replaced by the solution
2422 ! ----------------------------------------------------------------------
2424 implicit none
2426 ! ----------------------------------------------------------------------
2427 ! INPUT:
2428 ! ----------------------------------------------------------------------
2429 integer np
2430 real a(np,np)
2432 ! ----------------------------------------------------------------------
2433 ! OUTPUT:
2434 ! ----------------------------------------------------------------------
2435 real b(np)
2437 ! ----------------------------------------------------------------------
2438 ! LOCAL:
2439 ! ----------------------------------------------------------------------
2440 integer nmax
2441 parameter (nmax=150)
2443 real big,dum
2444 integer i,icol,irow
2445 integer j,k,l,ll,n
2446 integer ipiv(nmax)
2447 real pivinv
2449 ! ----------------------------------------------------------------------
2450 ! END VARIABLES DEFINITIONS
2451 ! ----------------------------------------------------------------------
2453 do j=1,n
2454 ipiv(j)=0.
2455 enddo
2457 do i=1,n
2458 big=0.
2459 do j=1,n
2460 if(ipiv(j).ne.1)then
2461 do k=1,n
2462 if(ipiv(k).eq.0)then
2463 if(abs(a(j,k)).ge.big)then
2464 big=abs(a(j,k))
2465 irow=j
2466 icol=k
2467 endif
2468 elseif(ipiv(k).gt.1)then
2469 FATAL_ERROR('singular matrix in gaussjbem')
2470 endif
2471 enddo
2472 endif
2473 enddo
2475 ipiv(icol)=ipiv(icol)+1
2477 if(irow.ne.icol)then
2478 do l=1,n
2479 dum=a(irow,l)
2480 a(irow,l)=a(icol,l)
2481 a(icol,l)=dum
2482 enddo
2484 dum=b(irow)
2485 b(irow)=b(icol)
2486 b(icol)=dum
2488 endif
2490 if(a(icol,icol).eq.0) FATAL_ERROR('singular matrix in gaussjbem')
2492 pivinv=1./a(icol,icol)
2493 a(icol,icol)=1
2495 do l=1,n
2496 a(icol,l)=a(icol,l)*pivinv
2497 enddo
2499 b(icol)=b(icol)*pivinv
2501 do ll=1,n
2502 if(ll.ne.icol)then
2503 dum=a(ll,icol)
2504 a(ll,icol)=0.
2505 do l=1,n
2506 a(ll,l)=a(ll,l)-a(icol,l)*dum
2507 enddo
2509 b(ll)=b(ll)-b(icol)*dum
2511 endif
2512 enddo
2513 enddo
2515 return
2516 end subroutine gaussjbem
2518 !====6=8===============================================================72
2519 !====6=8===============================================================72
2521 subroutine radfluxs(radflux,alb,rs,em,rl,sigma,twal)
2523 implicit none
2524 !-------------------------------------------------------------------
2525 !This function calculates the radiative fluxe at a surface
2526 !-------------------------------------------------------------------
2529 real alb !albedo of the surface
2530 real rs !shor wave radiation
2531 real em !emissivity of the surface
2532 real rl !lon wave radiation
2533 real sigma !parameter (wall is not black body) [W/m2.K4]
2534 real twal !wall temperature [K]
2535 real radflux
2537 radflux=(1.-alb)*rs+em*rl-em*sigma*twal**4
2539 return
2540 end subroutine radfluxs
2541 !====6=8===============================================================72
2542 !====6=8===============================================================72
2543 subroutine radfluxspv(nz,n,alb,rs,swddif,em,rl,twal,tair,sigma,radflux,pv_frac_roof,tpv)
2544 !
2545 implicit none
2546 !
2547 ! This routine calculates the radiative fluxes at the surfaces
2548 !
2549 ! Integer and real kinds
2550 !
2551 ! integer, parameter :: kind_im = selected_int_kind(6) ! 4 byte integer
2552 ! integer, parameter :: kind_rb = selected_real_kind(12) ! 8 byte real
2553 !
2554 ! Input Variables
2555 !
2556 integer,intent(in) :: nz !Maximum number of vertical levels in the urban grid
2557 real,intent(in) :: alb !albedo of the surface
2558 real,intent(in) :: rs !shortwave radiation [W m-2]
2559 real,intent(in) :: swddif
2560 real,intent(in) :: em !emissivity of the surface
2561 real,intent(in) :: rl !longwave radiation [W m-2]
2562 real,intent(in) :: twal !surface temperature [K]
2563 real,intent(in) :: sigma !Stefan-Boltzmann constant [W/m2.K4]
2564 real,intent(in) :: tpv !Stefan-Boltzmann constant [W/m2.K4]
2565 real,intent(in),dimension(1:nz) :: tair !external temperature [K]
2566 integer,intent(in) :: n !number of floors in the building
2567 real, intent(in) :: pv_frac_roof !
2568 real :: empv
2569 real :: hrad
2570 real :: F12
2571 ! Output variables
2572 !
2573 real,intent(inout) :: radflux !radiative flux at the surface [W m-2]
2574 !
2575 ! Local variables
2576 F12=1.
2577 empv=0.95
2578 hrad=sigma/((1-empv)/empv+1/F12+(1-em)/em)
2579 if ((n+1).gt.nz) then
2580 write(*,*) 'Increase maximum number of vertical levels in the urban grid'
2581 stop
2582 endif
2583 radflux=(1.-alb)*(1.-pv_frac_roof)*rs+em*(1.-pv_frac_roof)*rl+pv_frac_roof*hrad*(tpv**4-twal**4)- &
2584 em*sigma*(1.-pv_frac_roof)*twal**4
2585 return
2586 end subroutine radfluxspv
2587 !====6=8==============================================================72
2588 !====6=8==============================================================72
2589 !
2590 ! we define the view factors fprl and fnrm, which are the angle
2591 ! factors between two equal and parallel planes, fprl, and two
2592 ! equal and orthogonal planes, fnrm, respectively
2593 !
2594 subroutine fprl_ints(fprl_int,vx,vy)
2596 implicit none
2598 real vx,vy
2599 real fprl_int
2601 fprl_int=(2./(3.141592653*vx*vy))* &
2602 (log(sqrt((1.+vx*vx)*(1.+vy*vy)/(1.+vx*vx+vy*vy)))+ &
2603 (vy*sqrt(1.+vx*vx)*atan(vy/sqrt(1.+vx*vx)))+ &
2604 (vx*sqrt(1.+vy*vy)*atan(vx/sqrt(1.+vy*vy)))- &
2605 vy*atan(vy)-vx*atan(vx))
2607 return
2608 end subroutine fprl_ints
2610 !====6=8==============================================================72
2611 !====6=8==============================================================72
2612 !
2613 ! we define the view factors fprl and fnrm, which are the angle
2614 ! factors between two equal and parallel planes, fprl, and two
2615 ! equal and orthogonal planes, fnrm, respectively
2616 !
2618 subroutine fnrm_ints(fnrm_int,wx,wy,wz)
2620 implicit none
2622 real wx,wy,wz
2623 real fnrm_int
2625 fnrm_int=(1./(3.141592653*wy))*(wy*atan(1./wy)+wx*atan(1./wx)- &
2626 (sqrt(wz)*atan(1./sqrt(wz)))+ &
2627 (1./4.)*(log((1.+wx*wx)*(1.+wy*wy)/(1.+wz))+ &
2628 wy*wy*log(wy*wy*(1.+wz)/(wz*(1.+wy*wy)))+ &
2629 wx*wx*log(wx*wx*(1.+wz)/(wz*(1.+wx*wx)))))
2631 return
2632 end subroutine fnrm_ints
2634 !====6=8==============================================================72
2635 !====6=8==============================================================72
2636 END MODULE module_sf_bem