PageRenderTime 56ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_sf_bem.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2345 lines | 1181 code | 600 blank | 564 comment | 45 complexity | eba7fd746f42e2808e01b8d95dc369de MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

  1. MODULE module_sf_bem
  2. ! -----------------------------------------------------------------------
  3. ! Variables and constants used in the BEM module
  4. ! -----------------------------------------------------------------------
  5. real emins !emissivity of the internal walls
  6. parameter (emins=0.9)
  7. real albins !albedo of the internal walls
  8. !! parameter (albins=0.5)
  9. parameter (albins=0.3)
  10. real thickwin !thickness of the window [m]
  11. parameter (thickwin=0.006)
  12. real cswin !Specific heat of the windows [J/(m3.K)]
  13. parameter(cswin= 2.268e+06)
  14. real temp_rat !power of the A.C. heating/cooling the indoor air [K/s]
  15. parameter(temp_rat=0.001)
  16. real hum_rat !power of the A.C. drying/moistening the indoor air [(Kg/kg)/s]
  17. parameter(hum_rat=1.e-06)
  18. CONTAINS
  19. !====6================================================================72
  20. !====6================================================================72
  21. subroutine BEM(nzcanm,nlev,nhourday,dt,bw,bl,dzlev, &
  22. nwal,nflo,nrof,ngrd,hswalout,gswal, &
  23. hswinout,hsrof,gsrof, &
  24. latent,sigma,albwal,albwin,albrof, &
  25. emrof,emwal,emwin,rswal,rlwal,rair,cp, &
  26. rhoout,tout,humout,press, &
  27. rs,rl,dzwal,cswal,kwal,pwin,cop,beta,sw_cond, &
  28. timeon,timeoff,targtemp,gaptemp,targhum,gaphum, &
  29. perflo,hsesf,hsequip,dzflo, &
  30. csflo,kflo,dzgrd,csgrd,kgrd,dzrof,csrof, &
  31. krof,tlev,shumlev,twal,twin,tflo,tgrd,trof, &
  32. hsout,hlout,consump,hsvent,hlvent)
  33. ! ---------------------------------------------------------------------
  34. implicit none
  35. ! ---------------------------------------------------------------------
  36. ! TOP
  37. ! ---------------------
  38. ! ! ----------------- !--->roof (-) : level number
  39. ! ! ! ! ! rem: the windows are given
  40. ! ! !---------------! ! with respect to the
  41. ! ! !---------------! ! vertical walls-->win(2)
  42. ! (n)! !(1) (1)!-!(n)
  43. ! ! !---------------! ! 2D vision of the building
  44. ! WEST ! !-------4-------! ! EAST
  45. ! I ! ! 1 ilev 2! ! II ^
  46. ! ! !-------3--------! ! !
  47. ! ! !---------------! !--->floor 1 !
  48. ! ! ! ! ! !
  49. ! ! ! ! ! !
  50. ! ! ----------------- ! <--------------(n)
  51. ! ------------------------>ground ------------(1)
  52. ! BOTTOM
  53. ! i(6)
  54. ! i
  55. ! +---------v-----+
  56. ! /| /| 3D vision of a room
  57. ! / | 4 / |
  58. ! / | / |
  59. ! / | / |
  60. ! / | / |
  61. ! +---------------+ |
  62. ! | 1 | | 2 |
  63. ! | +---------|-----+
  64. ! dzlev | / | /
  65. ! | / 3 | /
  66. ! | /bw | /
  67. ! | / | /
  68. ! |/ |/
  69. ! +---------------+
  70. ! ^ bl
  71. ! i
  72. ! i
  73. ! (5)
  74. !-----------------------------------------------------------------------
  75. ! Input:
  76. ! -----
  77. real dt !time step [s]
  78. integer nzcanm !Maximum number of vertical levels in the urban grid
  79. integer nlev !number of floors in the building
  80. integer nwal !number of levels inside the wall
  81. integer nrof !number of levels inside the roof
  82. integer nflo !number of levels inside the floor
  83. integer ngrd !number of levels inside the ground
  84. real dzlev !vertical grid resolution [m]
  85. real bl !Building length [m]
  86. real bw !Building width [m]
  87. real albwal !albedo of the walls
  88. real albwin !albedo of the windows
  89. real albrof !albedo of the roof
  90. real emwal !emissivity of the walls
  91. real emrof !emissivity of the roof
  92. real emwin !emissivity of the windows
  93. real pwin !window proportion
  94. real, intent(in) :: cop !Coefficient of performance of the A/C systems
  95. real, intent(in) :: beta !Thermal efficiency of the heat exchanger
  96. integer, intent(in) :: sw_cond ! Air Conditioning switch
  97. real, intent(in) :: timeon ! Initial local time of A/C systems
  98. real, intent(in) :: timeoff ! Ending local time of A/C systems
  99. real, intent(in) :: targtemp ! Target temperature of A/C systems
  100. real, intent(in) :: gaptemp ! Comfort range of indoor temperature
  101. real, intent(in) :: targhum ! Target humidity of A/C systems
  102. real, intent(in) :: gaphum ! Comfort range of specific humidity
  103. real, intent(in) :: perflo ! Peak number of occupants per unit floor area
  104. real, intent(in) :: hsesf !
  105. real, intent(in) :: hsequip(24) !
  106. real cswal(nwal) !Specific heat of the wall [J/(m3.K)]
  107. real csflo(nflo) !Specific heat of the floor [J/(m3.K)]
  108. real csrof(nrof) !Specific heat of the roof [J/(m3.K)]
  109. real csgrd(ngrd) !Specific heat of the ground [J/(m3.K)]
  110. real kwal(nwal+1) !Thermal conductivity in each layers of the walls (face) [W/(m.K)]
  111. real kflo(nflo+1) !Thermal diffusivity in each layers of the floors (face) [W/(m.K)]
  112. real krof(nrof+1) !Thermal diffusivity in each layers of the roof (face) [W/(m.K)]
  113. real kgrd(ngrd+1) !Thermal diffusivity in each layers of the ground (face) [W/(m.K)]
  114. real dzwal(nwal) !Layer sizes of walls [m]
  115. real dzflo(nflo) !Layer sizes of floors [m]
  116. real dzrof(nrof) !Layer sizes of roof [m]
  117. real dzgrd(ngrd) !Layer sizes of ground [m]
  118. real latent !latent heat of evaporation [J/Kg]
  119. real rs !external short wave radiation [W/m2]
  120. real rl !external long wave radiation [W/m2]
  121. real rswal(4,nzcanm) !short wave radiation reaching the exterior walls [W/m2]
  122. real rlwal(4,nzcanm) !long wave radiation reaching the walls [W/m2]
  123. real rhoout(nzcanm) !exterior air density [kg/m3]
  124. real tout(nzcanm) !external temperature [K]
  125. real humout(nzcanm) !absolute humidity [Kgwater/Kgair]
  126. real press(nzcanm) !external air pressure [Pa]
  127. real hswalout(4,nzcanm) !outside walls sensible heat flux [W/m2]
  128. real hswinout(4,nzcanm) !outside window sensible heat flux [W/m2]
  129. real hsrof !Sensible heat flux at the roof [W/m2]
  130. real rair !ideal gas constant [J.kg-1.K-1]
  131. real sigma !parameter (wall is not black body) [W/m2.K4]
  132. real cp !specific heat of air [J/kg.K]
  133. !Input-Output
  134. !------------
  135. real tlev(nzcanm) !temperature of the floors [K]
  136. real shumlev(nzcanm) !specific humidity of the floor [kg/kg]
  137. real twal(4,nwal,nzcanm) !walls temperatures [K]
  138. real twin(4,nzcanm) !windows temperatures [K]
  139. real tflo(nflo,nzcanm-1) !floor temperatures [K]
  140. real tgrd(ngrd) !ground temperature [K]
  141. real trof(nrof) !roof temperature [K]
  142. real hsout(nzcanm) !sensible heat emitted outside the floor [W]
  143. real hlout(nzcanm) !latent heat emitted outside the floor [W]
  144. real consump(nzcanm) !Consumption for the a.c. in each floor [W]
  145. real hsvent(nzcanm) !sensible heat generated by natural ventilation [W]
  146. real hlvent(nzcanm) !latent heat generated by natural ventilation [W]
  147. real gsrof !heat flux flowing inside the roof [W/m˛]
  148. real gswal(4,nzcanm) !heat flux flowing inside the floors [W/m˛]
  149. ! Local:
  150. ! -----
  151. integer swwal !swich for the physical coefficients calculation
  152. integer ilev !index for rooms
  153. integer iwal !index for walls
  154. integer iflo !index for floors
  155. integer ivw !index for vertical walls
  156. integer igrd !index for ground
  157. integer irof !index for roof
  158. real hseqocc(nzcanm) !sensible heat generated by equipments and occupants [W]
  159. real hleqocc(nzcanm) !latent heat generated by occupants [W]
  160. real hscond(nzcanm) !sensible heat generated by wall conduction [W]
  161. real hslev(nzcanm) !sensible heat flux generated inside the room [W]
  162. real hllev(nzcanm) !latent heat flux generatd inside the room [W]
  163. real surwal(6,nzcanm) !Surface of the walls [m2]
  164. real surwal1D(6) !wall surfaces of a generic room [m2]
  165. real rsint(6) !short wave radiation reaching the indoor walls[W/m2]
  166. real rswalins(6,nzcanm) !internal short wave radiation for the building [W/m2]
  167. real twin1D(4) !temperature of windows for a particular room [K]
  168. real twal_int(6) !temperature of the first internal layers of a room [K]
  169. real rlint(6) !internal wall long wave radiation [w/m2]
  170. real rlwalins(6,nzcanm) !internal long wave radiation for the building [W/m2]
  171. real hrwalout(4,nzcanm) !external radiative flux to the walls [W/m2]
  172. real hrwalins(6,nzcanm) !inside radiative flux to the walls [W/m2]
  173. real hrwinout(4,nzcanm) !external radiative flux to the window [W/m2]
  174. real hrwinins(4,nzcanm) !inside radiative flux to the window [W/m2]
  175. real hrrof !external radiative flux to the roof [W/m2]
  176. real hs
  177. real hsneed(nzcanm) !sensible heat needed by the room [W]
  178. real hlneed(nzcanm) !latent heat needed by the room [W]
  179. real hswalins(6,nzcanm) !inside walls sensible heat flux [W/m2]
  180. real hswalins1D(6)
  181. real hswinins(4,nzcanm) !inside window sensible heat flux [W/m2]
  182. real hswinins1D(4)
  183. real htot(2) !total heat flux at the wall [W/m2]
  184. real twal1D(nwal)
  185. real tflo1D(nflo)
  186. real tgrd1D(ngrd)
  187. real trof1D(nrof)
  188. real rswal1D(4)
  189. real Qb !Overall heat capacity of the indoor air [J/K]
  190. real vollev !volume of the room [m3]
  191. real rhoint !density of the internal air [Kg/m3]
  192. real cpint !specific heat of the internal air [J/kg.K]
  193. real humdry !specific humidiy of dry air [kg water/kg dry air]
  194. real radflux !Function to compute the total radiation budget
  195. real consumpbuild !Energetic consumption for the entire building [KWh/s]
  196. real hsoutbuild !Total sensible heat ejected into the atmosphere[W]
  197. !by the air conditioning system and per building
  198. real nhourday !number of hours from midnight, local time
  199. !--------------------------------------------
  200. !Initialization
  201. !--------------------------------------------
  202. do ilev=1,nzcanm
  203. hseqocc(ilev)=0.
  204. hleqocc(ilev)=0.
  205. hscond(ilev)=0.
  206. hslev(ilev)=0.
  207. hllev(ilev)=0.
  208. enddo
  209. !Calculation of the surfaces of the building
  210. !--------------------------------------------
  211. do ivw=1,6
  212. do ilev=1,nzcanm
  213. surwal(ivw,ilev)=1. !initialisation
  214. end do
  215. end do
  216. do ilev=1,nlev
  217. do ivw=1,2
  218. surwal(ivw,ilev)=dzlev*bw
  219. end do
  220. do ivw=3,4
  221. surwal(ivw,ilev)=dzlev*bl
  222. end do
  223. do ivw=5,6
  224. surwal(ivw,ilev)=bw*bl
  225. end do
  226. end do
  227. ! Calculation of the short wave radiations at the internal walls
  228. ! ---------------------------------------------------------------
  229. do ilev=1,nlev
  230. do ivw=1,4
  231. rswal1D(ivw)=rswal(ivw,ilev)
  232. end do
  233. do ivw=1,6
  234. surwal1D(ivw)=surwal(ivw,ilev)
  235. end do
  236. call int_rsrad(albwin,albins,pwin,rswal1D,&
  237. surwal1D,bw,bl,dzlev,rsint)
  238. do ivw=1,6
  239. rswalins(ivw,ilev)=rsint(ivw)
  240. end do
  241. end do !ilev
  242. ! Calculation of the long wave radiation at the internal walls
  243. !-------------------------------------------------------------
  244. !Intermediate rooms
  245. if (nlev.gt.2) then
  246. do ilev=2,nlev-1
  247. do ivw=1,4
  248. twin1D(ivw)=twin(ivw,ilev)
  249. twal_int(ivw)=twal(ivw,1,ilev)
  250. end do
  251. twal_int(5)=tflo(nflo,ilev-1)
  252. twal_int(6)=tflo(1,ilev)
  253. call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
  254. pwin,bw,bl,dzlev,rlint)
  255. do ivw=1,6
  256. rlwalins(ivw,ilev)=rlint(ivw)
  257. end do
  258. end do !ilev
  259. end if
  260. if (nlev.ne.1) then
  261. !bottom room
  262. do ivw=1,4
  263. twin1D(ivw)=twin(ivw,1)
  264. twal_int(ivw)=twal(ivw,1,1)
  265. end do
  266. twal_int(5)=tgrd(ngrd)
  267. twal_int(6)=tflo(1,1)
  268. call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
  269. pwin,bw,bl,dzlev,rlint)
  270. do ivw=1,6
  271. rlwalins(ivw,1)=rlint(ivw)
  272. end do
  273. !top room
  274. do ivw=1,4
  275. twin1D(ivw)=twin(ivw,nlev)
  276. twal_int(ivw)=twal(ivw,1,nlev)
  277. end do
  278. twal_int(5)=tflo(nflo,nlev-1)
  279. twal_int(6)=trof(1)
  280. call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
  281. pwin,bw,bl,dzlev,rlint)
  282. do ivw=1,6
  283. rlwalins(ivw,nlev)=rlint(ivw)
  284. end do
  285. else !Top <---> Bottom
  286. do ivw=1,4
  287. twin1D(ivw)=twin(ivw,1)
  288. twal_int(ivw)=twal(ivw,1,1)
  289. end do
  290. twal_int(5)=tgrd(ngrd)
  291. twal_int(6)=trof(1)
  292. call int_rlrad(emins,emwin,sigma,twal_int,twin1D, &
  293. pwin,bw,bl,dzlev,rlint)
  294. do ivw=1,6
  295. rlwalins(ivw,1)=rlint(ivw)
  296. end do
  297. end if
  298. ! Calculation of the radiative fluxes
  299. ! -----------------------------------
  300. !External vertical walls and windows
  301. do ilev=1,nlev
  302. do ivw=1,4
  303. call radfluxs(radflux,albwal,rswal(ivw,ilev), &
  304. emwal,rlwal(ivw,ilev),sigma, &
  305. twal(ivw,nwal,ilev))
  306. hrwalout(ivw,ilev)=radflux
  307. hrwinout(ivw,ilev)=emwin*rlwal(ivw,ilev)- &
  308. emwin*sigma*(twin(ivw,ilev)**4)
  309. end do ! ivw
  310. end do ! ilev
  311. !Roof
  312. call radfluxs(radflux,albrof,rs,emrof,rl,sigma,trof(nrof))
  313. hrrof=radflux
  314. !Internal walls for intermediate rooms
  315. if(nlev.gt.2) then
  316. do ilev=2,nlev-1
  317. do ivw=1,4
  318. call radfluxs(radflux,albins,rswalins(ivw,ilev), &
  319. emins,rlwalins(ivw,ilev),sigma, &
  320. twal(ivw,1,ilev))
  321. hrwalins(ivw,ilev)=radflux
  322. end do !ivw
  323. call radfluxs(radflux,albins,rswalins(5,ilev), &
  324. emins,rlwalins(5,ilev),sigma,&
  325. tflo(nflo,ilev-1))
  326. hrwalins(5,ilev)=radflux
  327. call radfluxs(radflux,albins,rswalins(6,ilev), &
  328. emins,rlwalins(6,ilev),sigma,&
  329. tflo(1,ilev))
  330. hrwalins(6,ilev)=radflux
  331. end do !ilev
  332. end if
  333. !Internal walls for the bottom and the top room
  334. !
  335. if (nlev.ne.1) then
  336. !bottom floor
  337. do ivw=1,4
  338. call radfluxs(radflux,albins,rswalins(ivw,1), &
  339. emins,rlwalins(ivw,1),sigma, &
  340. twal(ivw,1,1))
  341. hrwalins(ivw,1)=radflux
  342. end do
  343. call radfluxs(radflux,albins,rswalins(5,1),&
  344. emins,rlwalins(5,1),sigma,& !bottom
  345. tgrd(ngrd))
  346. hrwalins(5,1)=radflux
  347. call radfluxs(radflux,albins,rswalins(6,1),&
  348. emins,rlwalins(6,1),sigma,&
  349. tflo(1,1))
  350. hrwalins(6,1)=radflux
  351. !roof floor
  352. do ivw=1,4
  353. call radfluxs(radflux,albins,rswalins(ivw,nlev), &
  354. emins,rlwalins(ivw,nlev),sigma,&
  355. twal(ivw,1,nlev))
  356. hrwalins(ivw,nlev)=radflux
  357. end do !top
  358. call radfluxs(radflux,albins,rswalins(5,nlev), &
  359. emins,rlwalins(5,nlev),sigma,&
  360. tflo(nflo,nlev-1))
  361. hrwalins(5,nlev)=radflux
  362. call radfluxs(radflux,albins,rswalins(6,nlev), &
  363. emins,rlwalins(6,nlev),sigma,&
  364. trof(1))
  365. hrwalins(6,nlev)=radflux
  366. else ! Top <---> Bottom room
  367. do ivw=1,4
  368. call radfluxs(radflux,albins,rswalins(ivw,1),&
  369. emins,rlwalins(ivw,1),sigma, &
  370. twal(ivw,1,1))
  371. hrwalins(ivw,1)=radflux
  372. end do
  373. call radfluxs(radflux,albins,rswalins(5,1),&
  374. emins,rlwalins(5,1),sigma, &
  375. tgrd(ngrd))
  376. hrwalins(5,1)=radflux
  377. call radfluxs(radflux,albins,rswalins(6,nlev), &
  378. emins,rlwalins(6,nlev),sigma,&
  379. trof(1))
  380. hrwalins(6,1)=radflux
  381. end if
  382. !Windows
  383. do ilev=1,nlev
  384. do ivw=1,4
  385. hrwinins(ivw,ilev)=emwin*rlwalins(ivw,ilev)- &
  386. emwin*sigma*(twin(ivw,ilev)**4)
  387. end do
  388. end do
  389. ! Calculation of the sensible heat fluxes
  390. ! ---------------------------------------
  391. !Vertical fluxes for walls
  392. do ilev=1,nlev
  393. do ivw=1,4
  394. call hsinsflux (2,2,tlev(ilev),twal(ivw,1,ilev),hs)
  395. hswalins(ivw,ilev)=hs
  396. end do ! ivw
  397. end do ! ilev
  398. !Vertical fluxes for windows
  399. do ilev=1,nlev
  400. do ivw=1,4
  401. call hsinsflux (2,1,tlev(ilev),twin(ivw,ilev),hs)
  402. hswinins(ivw,ilev)=hs
  403. end do ! ivw
  404. end do !ilev
  405. !Horizontal fluxes
  406. if (nlev.gt.2) then
  407. do ilev=2,nlev-1
  408. call hsinsflux (1,2,tlev(ilev),tflo(nflo,ilev-1),hs)
  409. hswalins(5,ilev)=hs
  410. call hsinsflux (1,2,tlev(ilev),tflo(1,ilev),hs)
  411. hswalins(6,ilev)=hs
  412. end do ! ilev
  413. end if
  414. if (nlev.ne.1) then
  415. call hsinsflux (1,2,tlev(1),tgrd(ngrd),hs)
  416. hswalins(5,1)=hs !Bottom room
  417. call hsinsflux (1,2,tlev(1),tflo(1,1),hs)
  418. hswalins(6,1)=hs
  419. call hsinsflux (1,2,tlev(nlev),tflo(nflo,nlev-1),hs)
  420. hswalins(5,nlev)=hs !Top room
  421. call hsinsflux (1,2,tlev(nlev),trof(1),hs)
  422. hswalins(6,nlev)=hs
  423. else ! Bottom<--->Top
  424. call hsinsflux (1,2,tlev(1),tgrd(ngrd),hs)
  425. hswalins(5,1)=hs
  426. call hsinsflux (1,2,tlev(nlev),trof(1),hs)
  427. hswalins(6,nlev)=hs
  428. end if
  429. !Calculation of the temperature for the different surfaces
  430. ! --------------------------------------------------------
  431. ! Vertical walls
  432. swwal=1
  433. do ilev=1,nlev
  434. do ivw=1,4
  435. htot(1)=hswalins(ivw,ilev)+hrwalins(ivw,ilev)
  436. htot(2)=hswalout(ivw,ilev)+hrwalout(ivw,ilev)
  437. gswal(ivw,ilev)=htot(2)
  438. do iwal=1,nwal
  439. twal1D(iwal)=twal(ivw,iwal,ilev)
  440. end do
  441. call wall(swwal,nwal,dt,dzwal,kwal,cswal,htot,twal1D)
  442. do iwal=1,nwal
  443. twal(ivw,iwal,ilev)=twal1D(iwal)
  444. end do
  445. end do ! ivw
  446. end do ! ilev
  447. ! Windows
  448. do ilev=1,nlev
  449. do ivw=1,4
  450. htot(1)=hswinins(ivw,ilev)+hrwinins(ivw,ilev)
  451. htot(2)=hswinout(ivw,ilev)+hrwinout(ivw,ilev)
  452. twin(ivw,ilev)=twin(ivw,ilev)+(dt/(cswin*thickwin))* &
  453. (htot(1)+htot(2))
  454. end do ! ivw
  455. end do ! ilev
  456. ! Horizontal floors
  457. if (nlev.gt.1) then
  458. swwal=1
  459. do ilev=1,nlev-1
  460. htot(1)=hrwalins(6,ilev)+hswalins(6,ilev)
  461. htot(2)=hrwalins(5,ilev+1)+hswalins(5,ilev+1)
  462. do iflo=1,nflo
  463. tflo1D(iflo)=tflo(iflo,ilev)
  464. end do
  465. call wall(swwal,nflo,dt,dzflo,kflo,csflo,htot,tflo1D)
  466. do iflo=1,nflo
  467. tflo(iflo,ilev)=tflo1D(iflo)
  468. end do
  469. end do ! ilev
  470. end if
  471. ! Ground
  472. swwal=1
  473. htot(1)=0. !Diriclet b.c. at the internal boundary
  474. htot(2)=hswalins(5,1)+hrwalins(5,1)
  475. do igrd=1,ngrd
  476. tgrd1D(igrd)=tgrd(igrd)
  477. end do
  478. call wall(swwal,ngrd,dt,dzgrd,kgrd,csgrd,htot,tgrd1D)
  479. do igrd=1,ngrd
  480. tgrd(igrd)=tgrd1D(igrd)
  481. end do
  482. ! Roof
  483. swwal=1
  484. htot(1)=hswalins(6,nlev)+hrwalins(6,nlev)
  485. htot(2)=hsrof+hrrof
  486. gsrof=htot(2)
  487. do irof=1,nrof
  488. trof1D(irof)=trof(irof)
  489. end do
  490. call wall(swwal,nrof,dt,dzrof,krof,csrof,htot,trof1D)
  491. do irof=1,nrof
  492. trof(irof)=trof1D(irof)
  493. end do
  494. ! Calculation of the heat fluxes and of the temperature of the rooms
  495. ! ------------------------------------------------------------------
  496. do ilev=1,nlev
  497. !Calculation of the heat generated by equipments and occupants
  498. call fluxeqocc(nhourday,bw,bl,perflo,hsesf,hsequip,hseqocc(ilev),hleqocc(ilev))
  499. !Calculation of the heat generated by natural ventilation
  500. vollev=bw*bl*dzlev
  501. humdry=shumlev(ilev)/(1.-shumlev(ilev))
  502. rhoint=(press(ilev))/(rair*(1.+0.61*humdry)*tlev(ilev))
  503. cpint=cp*(1.+0.84*humdry)
  504. call fluxvent(cpint,rhoint,vollev,tlev(ilev),tout(ilev), &
  505. latent,humout(ilev),rhoout(ilev),shumlev(ilev),&
  506. beta,hsvent(ilev),hlvent(ilev))
  507. !Calculation of the heat generated by conduction
  508. do iwal=1,6
  509. hswalins1D(iwal)=hswalins(iwal,ilev)
  510. surwal1D(iwal)=surwal(iwal,ilev)
  511. end do
  512. do iwal=1,4
  513. hswinins1D(iwal)=hswinins(iwal,ilev)
  514. end do
  515. call fluxcond(hswalins1D,hswinins1D,surwal1D,pwin,&
  516. hscond(ilev))
  517. !Calculation of the heat generated inside the room
  518. call fluxroo(hseqocc(ilev),hleqocc(ilev),hsvent(ilev), &
  519. hlvent(ilev),hscond(ilev),hslev(ilev),hllev(ilev))
  520. !Evolution of the temperature and of the specific humidity
  521. Qb=rhoint*cpint*vollev
  522. ! temperature regulation
  523. call regtemp(sw_cond,nhourday,dt,Qb,hslev(ilev), &
  524. tlev(ilev),timeon,timeoff,targtemp,gaptemp,hsneed(ilev))
  525. ! humidity regulation
  526. call reghum(sw_cond,nhourday,dt,vollev,rhoint,latent, &
  527. hllev(ilev),shumlev(ilev),timeon,timeoff,&
  528. targhum,gaphum,hlneed(ilev))
  529. !
  530. !performance of the air conditioning system for the test
  531. !
  532. call air_cond(hsneed(ilev),hlneed(ilev),dt, &
  533. hsout(ilev),hlout(ilev),consump(ilev), cop)
  534. tlev(ilev)=tlev(ilev)+(dt/Qb)*(hslev(ilev)-hsneed(ilev))
  535. shumlev(ilev)=shumlev(ilev)+(dt/(vollev*rhoint*latent))* &
  536. (hllev(ilev)-hlneed(ilev))
  537. end do !ilev
  538. call consump_total(nzcanm,nlev,consumpbuild,hsoutbuild, &
  539. hsout,consump)
  540. return
  541. end subroutine BEM
  542. !====6=8===============================================================72
  543. !====6=8===============================================================72
  544. subroutine wall(swwall,nz,dt,dz,k,cs,flux,temp)
  545. !______________________________________________________________________
  546. !The aim of this subroutine is to solve the 1D heat fiffusion equation
  547. !for roof, walls and streets:
  548. !
  549. ! dT/dt=d/dz[K*dT/dz] where:
  550. !
  551. ! -T is the surface temperature(wall, street, roof)
  552. ! -Kz is the heat diffusivity inside the material.
  553. !
  554. !The resolution is done implicitly with a FV discretisation along the
  555. !different layers of the material:
  556. ! ____________________________
  557. ! n *
  558. ! *
  559. ! *
  560. ! ____________________________
  561. ! i+2
  562. ! I+1
  563. ! ____________________________
  564. ! i+1
  565. ! I ==> [T(I,n+1)-T(I,n)]/DT=
  566. ! ____________________________ [F(i+1)-F(i)]/DZI
  567. ! i
  568. ! I-1 ==> A*T(n+1)=B where:
  569. ! ____________________________
  570. ! i-1 * * A is a TRIDIAGONAL matrix.
  571. ! * * B=T(n)+S takes into account the sources.
  572. ! *
  573. ! 1 ____________________________
  574. !________________________________________________________________
  575. implicit none
  576. !Input:
  577. !-----
  578. integer nz !Number of layers inside the material
  579. real dt !Time step
  580. real dz(nz) !Layer sizes [m]
  581. real cs(nz) !Specific heat of the material [J/(m3.K)]
  582. real k(nz+1) !Thermal conductivity in each layers (face) [W/(m.K)]
  583. real flux(2) !Internal and external flux terms.
  584. !Input-Output:
  585. !-------------
  586. integer swwall !swich for the physical coefficients calculation
  587. real temp(nz) !Temperature at each layer
  588. !Local:
  589. !-----
  590. real a(-1:1,nz) ! a(-1,*) lower diagonal A(i,i-1)
  591. ! a(0,*) principal diagonal A(i,i)
  592. ! a(1,*) upper diagonal A(i,i+1).
  593. real b(nz) !Coefficients of the second term.
  594. real k1(20)
  595. real k2(20)
  596. real kc(20)
  597. save k1,k2,kc
  598. integer iz
  599. !________________________________________________________________
  600. !
  601. !Calculation of the coefficients
  602. if (swwall.eq.1) then
  603. if (nz.gt.20) then
  604. write(*,*) 'number of layers in the walls/roofs too big ',nz
  605. write(*,*) 'please decrease under of',20
  606. stop
  607. endif
  608. call wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
  609. swwall=0
  610. end if
  611. !Computation of the first value (iz=1) of A and B:
  612. a(-1,1)=0.
  613. a(0,1)=1+k2(1)
  614. a(1,1)=-k2(1)
  615. b(1)=temp(1)+flux(1)*kc(1)
  616. !!
  617. !!We can fixed the internal temperature
  618. !!
  619. !! a(-1,1)=0.
  620. !! a(0,1)=1
  621. !! a(1,1)=0.
  622. !!
  623. !! b(1)=temp(1)
  624. !!
  625. !Computation of the internal values (iz=2,...,n-1) of A and B:
  626. do iz=2,nz-1
  627. a(-1,iz)=-k1(iz)
  628. a(0,iz)=1+k1(iz)+k2(iz)
  629. a(1,iz)=-k2(iz)
  630. b(iz)=temp(iz)
  631. end do
  632. !Computation of the external value (iz=n) of A and B:
  633. a(-1,nz)=-k1(nz)
  634. a(0,nz)=1+k1(nz)
  635. a(1,nz)=0.
  636. b(nz)=temp(nz)+flux(2)*kc(nz)
  637. !Resolution of the system A*T(n+1)=B
  638. call tridia(nz,a,b,temp)
  639. return
  640. end subroutine wall
  641. !====6=8===============================================================72
  642. !====6=8===============================================================72
  643. subroutine wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
  644. implicit none
  645. !---------------------------------------------------------------------
  646. !Input
  647. !-----
  648. integer nz !Number of layers inside the material
  649. real dt !Time step
  650. real dz(nz) !Layer sizes [m]
  651. real cs(nz) !Specific heat of the material [J/(m3.K)]
  652. real k(nz+1) !Thermal diffusivity in each layers (face) [W/(m.K)]
  653. !Input-Output
  654. !------------
  655. real flux(2) !Internal and external flux terms.
  656. !Output
  657. !------
  658. real k1(20)
  659. real k2(20)
  660. real kc(20)
  661. !Local
  662. !-----
  663. integer iz
  664. real kf(nz)
  665. !------------------------------------------------------------------
  666. do iz=2,nz
  667. kc(iz)=dt/(dz(iz)*cs(iz))
  668. kf(iz)=2*k(iz)/(dz(iz)+dz(iz-1))
  669. end do
  670. kc(1)=dt/(dz(1)*cs(1))
  671. kf(1)=2*k(1)/(dz(1))
  672. do iz=1,nz
  673. k1(iz)=kc(iz)*kf(iz)
  674. end do
  675. do iz=1,nz-1
  676. k2(iz)=kc(iz)*kf(iz+1)*cs(iz)/cs(iz+1)
  677. end do
  678. return
  679. end subroutine wall_coeff
  680. !====6=8===============================================================72
  681. !====6=8===============================================================72
  682. subroutine hsinsflux(swsurf,swwin,tin,tw,hsins)
  683. implicit none
  684. ! --------------------------------------------------------------------
  685. ! This routine computes the internal sensible heat flux.
  686. ! The swsurf, makes rhe difference between a vertical and a
  687. ! horizontal surface.
  688. ! The values of the heat conduction coefficients hc are obtained from the book
  689. ! "Energy Simulation in Building Design". J.A. Clarke.
  690. ! Adam Hilger, Bristol, 362 pp.
  691. ! --------------------------------------------------------------------
  692. !Input
  693. !----
  694. integer swsurf !swich for the type of surface (horizontal/vertical)
  695. integer swwin !swich for the type of surface (window/wall)
  696. real tin !Inside temperature [K]
  697. real tw !Internal wall temperature [K]
  698. !Output
  699. !------
  700. real hsins !internal sensible heat flux [W/m2]
  701. !Local
  702. !-----
  703. real hc !heat conduction coefficient [W/°C.m2]
  704. !--------------------------------------------------------------------
  705. if (swsurf.eq.2) then !vertical surface
  706. if (swwin.eq.1) then
  707. hc=5.678*0.99 !window surface (smooth surface)
  708. else
  709. hc=5.678*1.09 !wall surface (rough surface)
  710. endif
  711. hsins=hc*(tin-tw)
  712. endif
  713. if (swsurf.eq.1) then !horizontal surface
  714. if (swwin.eq.1) then
  715. hc=5.678*0.99 !window surface (smooth surface)
  716. else
  717. hc=5.678*1.09 !wall surface (rough surface)
  718. endif
  719. hsins=hc*(tin-tw)
  720. endif
  721. return
  722. end subroutine hsinsflux
  723. !====6=8===============================================================72
  724. !====6=8===============================================================72
  725. subroutine int_rsrad(albwin,albwal,pwin,rswal,&
  726. surwal,bw,bl,zw,rsint)
  727. ! ------------------------------------------------------------------
  728. implicit none
  729. ! ------------------------------------------------------------------
  730. !Input
  731. !-----
  732. real albwin !albedo of the windows
  733. real albwal !albedo of the internal wall
  734. real rswal(4) !incoming short wave radiation [W/m2]
  735. real surwal(6) !surface of the indoor walls [m2]
  736. real bw,bl !width of the walls [m]
  737. real zw !height of the wall [m]
  738. real pwin !window proportion
  739. !Output
  740. !------
  741. real rsint(6) !internal walls short wave radiation [W/m2]
  742. !Local
  743. !-----
  744. real transmit !transmittance of the direct/diffused radiation
  745. real rstr !solar radiation transmitted through the windows
  746. real surtotwal !total indoor surface of the walls in the room
  747. integer iw
  748. real b(6) !second member for the system
  749. real a(6,6) !matrix for the system
  750. !-------------------------------------------------------------------
  751. ! Calculation of the solar radiation transmitted through windows
  752. rstr = 0.
  753. do iw=1,4
  754. transmit=1.-albwin
  755. rstr = rstr+(surwal(iw)*pwin)*(transmit*rswal(iw))
  756. enddo
  757. !We suppose that the radiation is spread isotropically within the
  758. !room when it passes through the windows, so the flux [W/m˛] in every
  759. !wall is:
  760. surtotwal=0.
  761. do iw=1,6
  762. surtotwal=surtotwal+surwal(iw)
  763. enddo
  764. rstr=rstr/surtotwal
  765. !Computation of the short wave radiation reaching the internal walls
  766. call algebra_short(rstr,albwal,albwin,bw,bl,zw,pwin,a,b)
  767. call gaussjbem(a,6,b,6)
  768. do iw=1,6
  769. rsint(iw)=b(iw)
  770. enddo
  771. return
  772. end subroutine int_rsrad
  773. !====6=8===============================================================72
  774. !====6=8===============================================================72
  775. subroutine int_rlrad(emwal,emwin,sigma,twal_int,twin,&
  776. pwin,bw,bl,zw,rlint)
  777. ! ------------------------------------------------------------------
  778. implicit none
  779. ! ------------------------------------------------------------------
  780. !Input
  781. !-----
  782. real emwal !emissivity of the internal walls
  783. real emwin !emissivity of the window
  784. real sigma !Stefan-Boltzmann constant [W/m2.K4]
  785. real twal_int(6)!temperature of the first internal layers of a room [K]
  786. real twin(4) !temperature of the windows [K]
  787. real bw !width of the wall
  788. real bl !length of the wall
  789. real zw !height of the wall
  790. real pwin !window proportion
  791. !Output
  792. !------
  793. real rlint(6) !internal walls long wave radiation [W/m2]
  794. !Local
  795. !------
  796. real b(6) !second member vector for the system
  797. real a(6,6) !matrix for the system
  798. integer iw
  799. !----------------------------------------------------------------
  800. !Compute the long wave radiation reachs the internal walls
  801. call algebra_long(emwal,emwin,sigma,twal_int,twin,pwin,&
  802. bw,bl,zw,a,b)
  803. call gaussjbem(a,6,b,6)
  804. do iw=1,6
  805. rlint(iw)=b(iw)
  806. enddo
  807. return
  808. end subroutine int_rlrad
  809. !====6=8===============================================================72
  810. !====6=8===============================================================72
  811. subroutine algebra_short(rstr,albwal,albwin,aw,bw,zw,pwin,a,b)
  812. !--------------------------------------------------------------------
  813. !This routine calculates the algebraic system that will be solved for
  814. !the computation of the total shortwave radiation that reachs every
  815. !indoor wall in a floor.
  816. !Write the matrix system ax=b to solve
  817. !
  818. ! -Rs(1)+a(1,2)Rs(2)+.................+a(1,6)Rs(6)=-Rs=b(1)
  819. !a(2,1)Rs(1)- Rs(2)+.................+a(2,6)Rs(6)=-Rs=b(2)
  820. !a(3,1)Rs(1)+a(3,2)Rs(3)-Rs(3)+...........+a(3,6)Rs(6)=-Rs=b(3)
  821. !a(4,1)Rs(1)+.................-Rs(4)+.....+a(4,6)Rs(6)=-Rs=b(4)
  822. !a(5,1)Rs(1)+.......................-Rs(5)+a(5,6)Rs(6)=-Rs=b(5)
  823. !a(6,1)Rs(1)+....................................-R(6)=-Rs=b(6)
  824. !
  825. !This version suppose the albedo of the indoor walls is the same.
  826. !--------------------------------------------------------------------
  827. implicit none
  828. !--------------------------------------------------------------------
  829. !Input
  830. !-----
  831. real rstr !solar radiation transmitted through the windows
  832. real albwal !albedo of the internal walls
  833. real albwin !albedo of the windows.
  834. real bw !length of the wall
  835. real aw !width of the wall
  836. real zw !height of the wall
  837. real fprl_int !view factor
  838. real fnrm_int !view factor
  839. real pwin !window proportion
  840. !Output
  841. !------
  842. real a(6,6) !Matrix for the system
  843. real b(6) !Second member for the system
  844. !Local
  845. !-----
  846. integer iw,jw
  847. real albm !averaged albedo
  848. !----------------------------------------------------------------
  849. !Initialise the variables
  850. do iw=1,6
  851. b(iw)= 0.
  852. do jw=1,6
  853. a(iw,jw)= 0.
  854. enddo
  855. enddo
  856. !Calculation of the second member b
  857. do iw=1,6
  858. b(iw)=-rstr
  859. end do
  860. !Calculation of the averaged albedo
  861. albm=pwin*albwin+(1-pwin)*albwal
  862. !Calculation of the matrix a
  863. a(1,1)=-1.
  864. call fprl_ints(fprl_int,aw/bw,zw/bw)
  865. a(1,2)=albm*fprl_int
  866. call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
  867. a(1,3)=albm*(bw/aw)*fnrm_int
  868. a(1,4)=a(1,3)
  869. call fnrm_ints(fnrm_int,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
  870. a(1,5)=albwal*(bw/zw)*fnrm_int
  871. a(1,6)=a(1,5)
  872. a(2,1)=a(1,2)
  873. a(2,2)=-1.
  874. a(2,3)=a(1,3)
  875. a(2,4)=a(1,4)
  876. a(2,5)=a(1,5)
  877. a(2,6)=a(1,6)
  878. call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
  879. a(3,1)=albm*(aw/bw)*fnrm_int
  880. a(3,2)=a(3,1)
  881. a(3,3)=-1.
  882. call fprl_ints(fprl_int,zw/aw,bw/aw)
  883. a(3,4)=albm*fprl_int
  884. call fnrm_ints(fnrm_int,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
  885. a(3,5)=albwal*(aw/zw)*fnrm_int
  886. a(3,6)=a(3,5)
  887. a(4,1)=a(3,1)
  888. a(4,2)=a(3,2)
  889. a(4,3)=a(3,4)
  890. a(4,4)=-1.
  891. a(4,5)=a(3,5)
  892. a(4,6)=a(3,6)
  893. call fnrm_ints(fnrm_int,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
  894. a(5,1)=albm*(zw/bw)*fnrm_int
  895. a(5,2)=a(5,1)
  896. call fnrm_ints(fnrm_int,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
  897. a(5,3)=albm*(zw/aw)*fnrm_int
  898. a(5,4)=a(5,3)
  899. a(5,5)=-1.
  900. call fprl_ints(fprl_int,aw/zw,bw/zw)
  901. a(5,6)=albwal*fprl_int
  902. a(6,1)=a(5,1)
  903. a(6,2)=a(5,2)
  904. a(6,3)=a(5,3)
  905. a(6,4)=a(5,4)
  906. a(6,5)=a(5,6)
  907. a(6,6)=-1.
  908. return
  909. end subroutine algebra_short
  910. !====6=8===============================================================72
  911. !====6=8===============================================================72
  912. subroutine algebra_long(emwal,emwin,sigma,twalint,twinint,&
  913. pwin,aw,bw,zw,a,b)
  914. !--------------------------------------------------------------------
  915. !This routine computes the algebraic system that will be solved to
  916. !compute the longwave radiation that reachs the indoor
  917. !walls in a floor.
  918. !Write the matrix system ax=b to solve
  919. !
  920. !a(1,1)Rl(1)+.............................+Rl(6)=b(1)
  921. !a(2,1)Rl(1)+.................+Rl(5)+a(2,6)Rl(6)=b(2)
  922. !a(3,1)Rl(1)+.....+Rl(3)+...........+a(3,6)Rl(6)=b(3)
  923. !a(4,1)Rl(1)+...........+Rl(4)+.....+a(4,6)Rl(6)=b(4)
  924. ! Rl(1)+.......................+a(5,6)Rl(6)=b(5)
  925. !a(6,1)Rl(1)+Rl(2)+.................+a(6,6)Rl(6)=b(6)
  926. !
  927. !--------------------------------------------------------------------
  928. implicit none
  929. !--------------------------------------------------------------------
  930. !Input
  931. !-----
  932. real pwin !window proportion
  933. real emwal !emissivity of the internal walls
  934. real emwin !emissivity of the window
  935. real sigma !Stefan-Boltzmann constant [W/m2.K4]
  936. real twalint(6) !temperature of the first internal layers of a room [K]
  937. real twinint(4) !temperature of the windows [K]
  938. real aw !width of the wall
  939. real bw !length of the wall
  940. real zw !height of the wall
  941. real fprl_int !view factor
  942. real fnrm_int !view factor
  943. real fnrm_intx !view factor
  944. real fnrm_inty !view factor
  945. !Output
  946. !------
  947. real b(6) !second member vector for the system
  948. real a(6,6) !matrix for the system
  949. !Local
  950. !-----
  951. integer iw,jw
  952. real b_wall(6)
  953. real b_wind(6)
  954. real emwal_av !averadge emissivity of the wall
  955. real emwin_av !averadge emissivity of the window
  956. real em_av !averadge emissivity
  957. real twal_int(6) !twalint
  958. real twin(4) !twinint
  959. !------------------------------------------------------------------
  960. !Initialise the variables
  961. !-------------------------
  962. do iw=1,6
  963. b(iw)= 0.
  964. b_wall(iw)=0.
  965. b_wind(iw)=0.
  966. do jw=1,6
  967. a(iw,jw)= 0.
  968. enddo
  969. enddo
  970. do iw=1,6
  971. twal_int(iw)=twalint(iw)
  972. enddo
  973. do iw=1,4
  974. twin(iw)=twinint(iw)
  975. enddo
  976. !Calculation of the averadge emissivities
  977. !-----------------------------------------
  978. emwal_av=(1-pwin)*emwal
  979. emwin_av=pwin*emwin
  980. em_av=emwal_av+emwin_av
  981. !Calculation of the second term for the walls
  982. !-------------------------------------------
  983. call fprl_ints(fprl_int,aw/zw,bw/zw)
  984. call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
  985. call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
  986. b_wall(1)=(emwal*sigma*(twal_int(5)**4)* &
  987. fprl_int)+ &
  988. (sigma*(emwal_av*(twal_int(3)**4)+ &
  989. emwal_av*(twal_int(4)**4))* &
  990. (zw/aw)*fnrm_intx)+ &
  991. (sigma*(emwal_av*(twal_int(1)**4)+ &
  992. emwal_av*(twal_int(2)**4))* &
  993. (zw/bw)*fnrm_inty)
  994. call fprl_ints(fprl_int,aw/zw,bw/zw)
  995. call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
  996. call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
  997. b_wall(2)=(emwal*sigma*(twal_int(6)**4)* &
  998. fprl_int)+ &
  999. (sigma*(emwal_av*(twal_int(3)**4)+ &
  1000. emwal_av*(twal_int(4)**4))* &
  1001. (zw/aw)*fnrm_intx)+ &
  1002. (sigma*(emwal_av*(twal_int(1)**4)+ &
  1003. emwal_av*(twal_int(2)**4))* &
  1004. (zw/bw)*fnrm_inty)
  1005. call fprl_ints(fprl_int,zw/aw,bw/aw)
  1006. call fnrm_ints(fnrm_intx,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
  1007. call fnrm_ints(fnrm_inty,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
  1008. b_wall(3)=(emwal_av*sigma*(twal_int(4)**4)* &
  1009. fprl_int)+ &
  1010. (sigma*(emwal_av*(twal_int(2)**4)+ &
  1011. emwal_av*(twal_int(1)**4))* &
  1012. (aw/bw)*fnrm_intx)+ &
  1013. (sigma*(emwal*(twal_int(5)**4)+ &
  1014. emwal*(twal_int(6)**4))* &
  1015. (aw/zw)*fnrm_inty)
  1016. call fprl_ints(fprl_int,zw/aw,bw/aw)
  1017. call fnrm_ints(fnrm_intx,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
  1018. call fnrm_ints(fnrm_inty,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
  1019. b_wall(4)=(emwal_av*sigma*(twal_int(3)**4)* &
  1020. fprl_int)+ &
  1021. (sigma*(emwal_av*(twal_int(2)**4)+ &
  1022. emwal_av*(twal_int(1)**4))* &
  1023. (aw/bw)*fnrm_intx)+ &
  1024. (sigma*(emwal*(twal_int(5)**4)+ &
  1025. emwal*(twal_int(6)**4))* &
  1026. (aw/zw)*fnrm_inty)
  1027. call fprl_ints(fprl_int,aw/bw,zw/bw)
  1028. call fnrm_ints(fnrm_intx,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
  1029. call fnrm_ints(fnrm_inty,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
  1030. b_wall(5)=(emwal_av*sigma*(twal_int(2)**4)* &
  1031. fprl_int)+ &
  1032. (sigma*(emwal_av*(twal_int(3)**4)+ &
  1033. emwal_av*(twal_int(4)**4))* &
  1034. (bw/aw)*fnrm_intx)+ &
  1035. (sigma*(emwal*(twal_int(5)**4)+ &
  1036. emwal*(twal_int(6)**4))* &
  1037. (bw/zw)*fnrm_inty)
  1038. call fprl_ints(fprl_int,aw/bw,zw/bw)
  1039. call fnrm_ints(fnrm_intx,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
  1040. call fnrm_ints(fnrm_inty,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
  1041. b_wall(6)=(emwal_av*sigma*(twal_int(1)**4)* &
  1042. fprl_int)+ &
  1043. (sigma*(emwal_av*(twal_int(3)**4)+ &
  1044. emwal_av*(twal_int(4)**4))* &
  1045. (bw/aw)*fnrm_intx)+ &
  1046. (sigma*(emwal*(twal_int(5)**4)+ &
  1047. emwal*(twal_int(6)**4))* &
  1048. (bw/zw)*fnrm_inty)
  1049. !Calculation of the second term for the windows
  1050. !---------------------------------------------
  1051. call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
  1052. call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
  1053. b_wind(1)=(sigma*(emwin_av*(twin(3)**4)+ &
  1054. emwin_av*(twin(4)**4))* &
  1055. (zw/aw)*fnrm_intx)+ &
  1056. (sigma*(emwin_av*(twin(1)**4)+ &
  1057. emwin_av*(twin(2)**4))* &
  1058. (zw/bw)*fnrm_inty)
  1059. call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
  1060. call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
  1061. b_wind(2)=(sigma*(emwin_av*(twin(3)**4)+ &
  1062. emwin_av*(twin(4)**4))* &
  1063. (zw/aw)*fnrm_intx)+ &
  1064. (sigma*(emwin_av*(twin(1)**4)+ &
  1065. emwin_av*(twin(2)**4))* &
  1066. (zw/bw)*fnrm_inty)
  1067. call fprl_ints(fprl_int,zw/aw,bw/aw)
  1068. call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
  1069. b_wind(3)=emwin_av*sigma*(twin(4)**4)* &
  1070. fprl_int+(sigma*(emwin_av* &
  1071. (twin(2)**4)+emwin_av*(twin(1)**4))* &
  1072. (aw/bw)*fnrm_int)
  1073. call fprl_ints(fprl_int,zw/aw,bw/aw)
  1074. call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
  1075. b_wind(4)=emwin_av*sigma*(twin(3)**4)* &
  1076. fprl_int+(sigma*(emwin_av* &
  1077. (twin(2)**4)+emwin_av*(twin(1)**4))* &
  1078. (aw/bw)*fnrm_int)
  1079. call fprl_ints(fprl_int,aw/bw,zw/bw)
  1080. call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
  1081. b_wind(5)=emwin_av*sigma*(twin(2)**4)* &
  1082. fprl_int+(sigma*(emwin_av* &
  1083. (twin(3)**4)+emwin_av*(twin(4)**4))* &
  1084. (bw/aw)*fnrm_int)
  1085. call fprl_ints(fprl_int,aw/bw,zw/bw)
  1086. call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
  1087. b_wind(6)=emwin_av*sigma*(twin(1)**4)* &
  1088. fprl_int+(sigma*(emwin_av* &
  1089. (twin(3)**4)+emwin_av*(twin(4)**4))* &
  1090. (bw/aw)*fnrm_int)
  1091. !Calculation of the total b term
  1092. !-------------------------------
  1093. do iw=1,6
  1094. b(iw)=b_wall(iw)+b_wind(iw)
  1095. end do
  1096. !Calculation of the matrix of the system
  1097. !----------------------------------------
  1098. call fnrm_ints(fnrm_int,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
  1099. a(1,1)=(em_av-1.)*(zw/bw)*fnrm_int
  1100. a(1,2)=a(1,1)
  1101. call fnrm_ints(fnrm_int,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
  1102. a(1,3)=(em_av-1.)*(zw/aw)*fnrm_int
  1103. a(1,4)=a(1,3)
  1104. call fprl_ints(fprl_int,aw/zw,bw/zw)
  1105. a(1,5)=(emwal-1.)*fprl_int
  1106. a(1,6)=1.
  1107. a(2,1)=a(1,1)
  1108. a(2,2)=a(1,2)
  1109. a(2,3)=a(1,3)
  1110. a(2,4)=a(1,4)
  1111. a(2,5)=1.
  1112. a(2,6)=a(1,5)
  1113. call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
  1114. a(3,1)=(em_av-1.)*(aw/bw)*fnrm_int
  1115. a(3,2)=a(3,1)
  1116. a(3,3)=1.
  1117. call fprl_ints(fprl_int,zw/aw,bw/aw)
  1118. a(3,4)=(em_av-1.)*fprl_int
  1119. call fnrm_ints(fnrm_int,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
  1120. a(3,5)=(emwal-1.)*(aw/zw)*fnrm_int
  1121. a(3,6)=a(3,5)
  1122. a(4,1)=a(3,1)
  1123. a(4,2)=a(3,2)
  1124. a(4,3)=a(3,4)
  1125. a(4,4)=1.
  1126. a(4,5)=a(3,5)
  1127. a(4,6)=a(3,6)
  1128. a(5,1)=1.
  1129. call fprl_ints(fprl_int,aw/bw,zw/bw)
  1130. a(5,2)=(em_av-1.)*fprl_int
  1131. call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
  1132. a(5,3)=(em_av-1.)*(bw/aw)*fnrm_int
  1133. a(5,4)=a(5,3)
  1134. call fnrm_ints(fnrm_int,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
  1135. a(5,5)=(emwal-1.)*(bw/zw)*fnrm_int
  1136. a(5,6)=a(5,5)
  1137. a(6,1)=a(5,2)
  1138. a(6,2)=1.
  1139. a(6,3)=a(5,3)
  1140. a(6,4)=a(5,4)
  1141. a(6,5)=a(5,5)
  1142. a(6,6)=a(6,5)
  1143. return
  1144. end subroutine algebra_long
  1145. !====6=8===============================================================72
  1146. !====6=8===============================================================72
  1147. subroutine fluxroo(hseqocc,hleqocc,hsvent,hlvent, &
  1148. hscond,hslev,hllev)
  1149. !-----------------------------------------------------------------------
  1150. !This routine calculates the heat flux generated inside the room
  1151. !and the heat ejected to the atmosphere.
  1152. !----------------------------------------------------------------------
  1153. implicit none
  1154. !-----------------------------------------------------------------------
  1155. !Input
  1156. !-----
  1157. real hseqocc !sensible heat generated by equipments and occupants [W]
  1158. real hleqocc !latent heat generated by occupants [W]
  1159. real hsvent !sensible heat generated by natural ventilation [W]
  1160. real hlvent !latent heat generated by natural ventilation [W]
  1161. real hscond !sensible heat generated by wall conduction
  1162. !Output
  1163. !------
  1164. real hslev !sensible heat flux generated inside the room [W]
  1165. real hllev !latent heat flux generatd inside the room
  1166. !Calculation of the total sensible heat generated inside the room
  1167. hslev=hseqocc+hsvent+hscond
  1168. !Calculation of the total latent heat generated inside the room
  1169. hllev=hleqocc+hlvent
  1170. return
  1171. end subroutine fluxroo
  1172. !====6=8===============================================================72
  1173. !====6=8===============================================================72
  1174. subroutine phirat(nhourday,rocc)
  1175. !----------------------------------------------------------------------
  1176. !This routine calculates the occupation ratio of a floor
  1177. !By now we suppose a constant value
  1178. !----------------------------------------------------------------------
  1179. implicit none
  1180. !Input
  1181. !-----
  1182. real nhourday ! number of hours from midnight (local time)
  1183. !Output
  1184. !------
  1185. real rocc !value between 0 and 1
  1186. !!TEST
  1187. rocc=1.
  1188. return
  1189. end subroutine phirat
  1190. !====6=8===============================================================72
  1191. !====6=8===============================================================72
  1192. subroutine phiequ(nhourday,hsesf,hsequip,hsequ)
  1193. !----------------------------------------------------------------------
  1194. !This routine calculates the sensible heat gain from equipments
  1195. !----------------------------------------------------------------------
  1196. implicit none
  1197. !Input
  1198. !-----
  1199. real nhourday ! number of hours from midnight, Local time
  1200. real, intent(in) :: hsesf
  1201. real, intent(in), dimension(24) :: hsequip

Large files files are truncated, but you can click here to view the full file