wrfv2_fire/phys/module_sf_bem.F

   1 MODULE module_sf_bem
   2 ! -----------------------------------------------------------------------
   3 !  Variables and constants used in the BEM module
   4 ! -----------------------------------------------------------------------
   5
   6         real emins              !emissivity of the internal walls
   7         parameter (emins=0.9)
   8         real albins             !albedo of the internal walls
   9 !!      parameter (albins=0.5)
  10         parameter (albins=0.3)
  11
  12         real thickwin           !thickness of the window [m]
  13         parameter (thickwin=0.006)
  14         real cswin              !Specific heat of the windows [J/(m3.K)]
  15         parameter(cswin= 2.268e+06)
  16
  17         real temp_rat            !power of the A.C. heating/cooling the indoor air [K/s]
  18         parameter(temp_rat=0.001)
  19
  20         real hum_rat            !power of the A.C. drying/moistening the indoor air [(Kg/kg)/s]
  21         parameter(hum_rat=1.e-06)
  22
  23
  24     CONTAINS
  25
  26 !====6================================================================72
  27 !====6================================================================72
  28
  29         subroutine BEM(nzcanm,nlev,nhourday,dt,bw,bl,dzlev,            &
  30                        nwal,nflo,nrof,ngrd,hswalout,gswal,             &
  31                        hswinout,hsrof,gsrof,                           &
  32                        latent,sigma,albwal,albwin,albrof,              &
  33                        emrof,emwal,emwin,rswal,rlwal,rair,cp,          &
  34                        rhoout,tout,humout,press,                       &
  35                        rs,rl,dzwal,cswal,kwal,pwin,cop,beta,sw_cond,   &
  36                        timeon,timeoff,targtemp,gaptemp,targhum,gaphum, &
  37                        perflo,hsesf,hsequip,dzflo,                     &
  38                        csflo,kflo,dzgrd,csgrd,kgrd,dzrof,csrof,        &
  39                        krof,tlev,shumlev,twal,twin,tflo,tgrd,trof,     &
  40                        hsout,hlout,consump,hsvent,hlvent)
  41
  42
  43 ! ---------------------------------------------------------------------
  44         implicit none
  45
  46 ! ---------------------------------------------------------------------
  47 !                      TOP
  48 !             ---------------------
  49 !             ! ----------------- !--->roof     (-) : level number
  50 !             ! !               ! !             rem: the windows are given
  51 !             ! !---------------! !                  with respect to the
  52 !             ! !---------------! !                  vertical walls-->win(2)
  53 !          (n)! !(1)         (1)!-!(n)
  54 !             ! !---------------! !             2D vision of the building
  55 !   WEST      ! !-------4-------! !     EAST
  56 !           I ! ! 1    ilev    2! ! II               ^
  57 !             ! !-------3--------! !                 !
  58 !             ! !---------------! !--->floor 1       !
  59 !             ! !               ! !                  !
  60 !             ! !               ! !                  !
  61 !             ! ----------------- !          <--------------(n)
  62 !             ------------------------>ground   ------------(1)
  63 !                    BOTTOM
  64 !                               i(6)
  65 !                               i
  66 !                     +---------v-----+
  67 !                    /|              /|    3D vision of a room
  68 !                   / | 4           / |
  69 !                  /  |            /  |
  70 !                 /   |           /   |
  71 !                /    |          /    |
  72 !               +---------------+     |
  73 !               |  1   |        |  2  |
  74 !               |     +---------|-----+
  75 !       dzlev   |    /          |    /
  76 !               |   /    3      |   /
  77 !               |  /bw          |  /
  78 !               | /             | /
  79 !               |/              |/
  80 !               +---------------+
  81 !                     ^ bl
  82 !                     i
  83 !                     i
  84 !                    (5)
  85 !-----------------------------------------------------------------------
  86
  87
  88 ! Input:
  89 ! -----
  90
  91         real dt                         !time step [s]
  92
  93         integer nzcanm                  !Maximum number of vertical levels in the urban grid
  94         integer nlev                    !number of floors in the building
  95         integer nwal                    !number of levels inside the wall
  96         integer nrof                    !number of levels inside the roof
  97         integer nflo                    !number of levels inside the floor
  98         integer ngrd                    !number of levels inside the ground
  99         real dzlev                      !vertical grid resolution [m]
 100         real bl                         !Building length [m]
 101         real bw                         !Building width [m]
 102
 103         real albwal                     !albedo of the walls
 104         real albwin                     !albedo of the windows
 105         real albrof                     !albedo of the roof
 106
 107         real emwal                      !emissivity of the walls
 108
 109         real emrof                      !emissivity of the roof
 110         real emwin                      !emissivity of the windows
 111
 112         real pwin                       !window proportion
 113         real,    intent(in) :: cop      !Coefficient of performance of the A/C systems
 114         real,    intent(in) :: beta     !Thermal efficiency of the heat exchanger
 115         integer, intent(in) :: sw_cond  ! Air Conditioning switch
 116         real,    intent(in) :: timeon   ! Initial local time of A/C systems
 117         real,    intent(in) :: timeoff  ! Ending local time of A/C systems
 118         real,    intent(in) :: targtemp ! Target temperature of A/C systems
 119         real,    intent(in) :: gaptemp  ! Comfort range of indoor temperature
 120         real,    intent(in) :: targhum  ! Target humidity of A/C systems
 121         real,    intent(in) :: gaphum   ! Comfort range of specific humidity
 122         real,    intent(in) :: perflo   ! Peak number of occupants per unit floor area
 123         real,    intent(in) :: hsesf    !
 124         real,    intent(in) :: hsequip(24) !
 125
 126         real cswal(nwal)                !Specific heat of the wall [J/(m3.K)]
 127
 128         real csflo(nflo)                !Specific heat of the floor [J/(m3.K)]
 129         real csrof(nrof)                !Specific heat of the roof [J/(m3.K)]
 130         real csgrd(ngrd)                !Specific heat of the ground [J/(m3.K)]
 131
 132         real kwal(nwal+1)               !Thermal conductivity in each layers of the walls (face) [W/(m.K)]
 133         real kflo(nflo+1)               !Thermal diffusivity in each layers of the floors (face) [W/(m.K)]
 134         real krof(nrof+1)               !Thermal diffusivity in each layers of the roof (face) [W/(m.K)]
 135         real kgrd(ngrd+1)               !Thermal diffusivity in each layers of the ground (face) [W/(m.K)]
 136
 137         real dzwal(nwal)                !Layer sizes of walls [m]
 138         real dzflo(nflo)                !Layer sizes of floors [m]
 139         real dzrof(nrof)                !Layer sizes of roof [m]
 140         real dzgrd(ngrd)                !Layer sizes of ground [m]
 141
 142         real latent                      !latent heat of evaporation [J/Kg]
 143
 144
 145         real rs                         !external short wave radiation [W/m2]
 146         real rl                         !external long wave radiation [W/m2]
 147         real rswal(4,nzcanm)            !short wave radiation reaching the exterior walls [W/m2]
 148         real rlwal(4,nzcanm)            !long wave radiation reaching the walls [W/m2]
 149         real rhoout(nzcanm)             !exterior air density [kg/m3]
 150         real tout(nzcanm)               !external temperature [K]
 151         real humout(nzcanm)             !absolute humidity [Kgwater/Kgair]
 152         real press(nzcanm)              !external air pressure [Pa]
 153
 154         real hswalout(4,nzcanm)         !outside walls sensible heat flux [W/m2]
 155         real hswinout(4,nzcanm)         !outside window sensible heat flux [W/m2]
 156         real hsrof                      !Sensible heat flux at the roof [W/m2]
 157
 158         real rair                       !ideal gas constant  [J.kg-1.K-1]
 159         real sigma                      !parameter (wall is not black body) [W/m2.K4]
 160         real cp                         !specific heat of air [J/kg.K]
 161
 162
 163 !Input-Output
 164 !------------
 165         real tlev(nzcanm)               !temperature of the floors [K]
 166         real shumlev(nzcanm)            !specific humidity of the floor [kg/kg]
 167         real twal(4,nwal,nzcanm)        !walls temperatures [K]
 168         real twin(4,nzcanm)             !windows temperatures [K]
 169         real tflo(nflo,nzcanm-1)        !floor temperatures [K]
 170         real tgrd(ngrd)                 !ground temperature [K]
 171         real trof(nrof)                 !roof temperature [K]
 172         real hsout(nzcanm)              !sensible heat emitted outside the floor [W]
 173         real hlout(nzcanm)              !latent heat emitted outside the floor [W]
 174         real consump(nzcanm)            !Consumption for the a.c. in each floor [W]
 175         real hsvent(nzcanm)             !sensible heat generated by natural ventilation [W]
 176         real hlvent(nzcanm)             !latent heat generated by natural ventilation [W]
 177         real gsrof                      !heat flux flowing inside the roof [W/m²]
 178         real gswal(4,nzcanm)             !heat flux flowing inside the floors [W/m²]
 179
 180 ! Local:
 181 ! -----
 182         integer swwal                   !swich for the physical coefficients calculation
 183         integer ilev                    !index for rooms
 184         integer iwal                    !index for walls
 185         integer iflo                    !index for floors
 186         integer ivw                     !index for vertical walls
 187         integer igrd                    !index for ground
 188         integer irof                    !index for roof
 189         real hseqocc(nzcanm)            !sensible heat generated by equipments and occupants [W]
 190         real hleqocc(nzcanm)            !latent heat generated by occupants [W]
 191         real hscond(nzcanm)             !sensible heat generated by wall conduction [W]
 192         real hslev(nzcanm)              !sensible heat flux generated inside the room [W]
 193         real hllev(nzcanm)              !latent heat flux generatd inside the room [W]
 194         real surwal(6,nzcanm)           !Surface of the walls [m2]
 195         real surwal1D(6)                !wall surfaces of a generic room [m2]
 196         real rsint(6)                   !short wave radiation reaching the indoor walls[W/m2]
 197         real rswalins(6,nzcanm)         !internal short wave radiation for the building [W/m2]
 198         real twin1D(4)                  !temperature of windows for a particular room [K]
 199         real twal_int(6)                !temperature of the first internal layers of a room [K]
 200         real rlint(6)                   !internal wall long wave radiation [w/m2]
 201         real rlwalins(6,nzcanm)         !internal long wave radiation for the building [W/m2]
 202         real hrwalout(4,nzcanm)         !external radiative flux to the walls [W/m2]
 203         real hrwalins(6,nzcanm)         !inside radiative flux to the walls [W/m2]
 204         real hrwinout(4,nzcanm)         !external radiative flux to the window [W/m2]
 205         real hrwinins(4,nzcanm)         !inside radiative flux to the window [W/m2]
 206         real hrrof                      !external radiative flux to the roof [W/m2]
 207         real hs
 208         real hsneed(nzcanm)             !sensible heat needed by the room [W]
 209         real hlneed(nzcanm)             !latent heat needed by the room [W]
 210         real hswalins(6,nzcanm)         !inside walls sensible heat flux [W/m2]
 211         real hswalins1D(6)
 212         real hswinins(4,nzcanm)         !inside window sensible heat flux [W/m2]
 213         real hswinins1D(4)
 214         real htot(2)                    !total heat flux at the wall [W/m2]
 215         real twal1D(nwal)
 216         real tflo1D(nflo)
 217         real tgrd1D(ngrd)
 218         real trof1D(nrof)
 219         real rswal1D(4)
 220         real Qb                         !Overall heat capacity of the indoor air [J/K]
 221         real vollev                     !volume of the room [m3]
 222         real rhoint                     !density of the internal air [Kg/m3]
 223         real cpint                      !specific heat of the internal air [J/kg.K]
 224         real humdry                     !specific humidiy of dry air [kg water/kg dry air]
 225         real radflux                    !Function to compute the total radiation budget
 226         real consumpbuild               !Energetic consumption for the entire building [KWh/s]
 227         real hsoutbuild                 !Total sensible heat ejected into the atmosphere[W]
 228                                         !by the air conditioning system and per building
 229         real nhourday                   !number of hours from midnight, local time
 230 !--------------------------------------------
 231 !Initialization
 232 !--------------------------------------------
 233
 234        do ilev=1,nzcanm
 235           hseqocc(ilev)=0.
 236           hleqocc(ilev)=0.
 237           hscond(ilev)=0.
 238           hslev(ilev)=0.
 239           hllev(ilev)=0.
 240        enddo
 241
 242 !Calculation of the surfaces of the building
 243 !--------------------------------------------
 244
 245
 246         do ivw=1,6
 247         do ilev=1,nzcanm
 248          surwal(ivw,ilev)=1.   !initialisation
 249         end do
 250         end do
 251
 252         do ilev=1,nlev
 253           do ivw=1,2
 254            surwal(ivw,ilev)=dzlev*bw
 255           end do
 256           do ivw=3,4
 257            surwal(ivw,ilev)=dzlev*bl
 258           end do
 259           do ivw=5,6
 260            surwal(ivw,ilev)=bw*bl
 261           end do
 262         end do
 263
 264
 265 ! Calculation of the short wave radiations at the internal walls
 266 ! ---------------------------------------------------------------
 267
 268
 269         do ilev=1,nlev
 270
 271           do ivw=1,4
 272             rswal1D(ivw)=rswal(ivw,ilev)
 273           end do
 274
 275           do ivw=1,6
 276             surwal1D(ivw)=surwal(ivw,ilev)
 277           end do
 278
 279           call int_rsrad(albwin,albins,pwin,rswal1D,&
 280                          surwal1D,bw,bl,dzlev,rsint)
 281
 282           do ivw=1,6
 283             rswalins(ivw,ilev)=rsint(ivw)
 284           end do
 285
 286         end do !ilev
 287
 288
 289
 290 ! Calculation of the long wave radiation at the internal walls
 291 !-------------------------------------------------------------
 292
 293
 294 !Intermediate rooms
 295
 296        if (nlev.gt.2) then
 297         do ilev=2,nlev-1
 298
 299           do ivw=1,4
 300             twin1D(ivw)=twin(ivw,ilev)
 301             twal_int(ivw)=twal(ivw,1,ilev)
 302           end do
 303
 304            twal_int(5)=tflo(nflo,ilev-1)
 305            twal_int(6)=tflo(1,ilev)
 306
 307            call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
 308                           pwin,bw,bl,dzlev,rlint)
 309
 310
 311           do ivw=1,6
 312             rlwalins(ivw,ilev)=rlint(ivw)
 313           end do
 314
 315         end do  !ilev
 316       end if
 317
 318
 319       if (nlev.ne.1) then
 320
 321 !bottom room
 322
 323           do ivw=1,4
 324             twin1D(ivw)=twin(ivw,1)
 325             twal_int(ivw)=twal(ivw,1,1)
 326           end do
 327
 328           twal_int(5)=tgrd(ngrd)
 329           twal_int(6)=tflo(1,1)
 330
 331
 332            call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
 333                           pwin,bw,bl,dzlev,rlint)
 334
 335           do ivw=1,6
 336             rlwalins(ivw,1)=rlint(ivw)
 337           end do
 338
 339 !top room
 340
 341           do ivw=1,4
 342             twin1D(ivw)=twin(ivw,nlev)
 343             twal_int(ivw)=twal(ivw,1,nlev)
 344           end do
 345
 346           twal_int(5)=tflo(nflo,nlev-1)
 347           twal_int(6)=trof(1)
 348
 349
 350            call int_rlrad(emins,emwin,sigma,twal_int,twin1D,&
 351                           pwin,bw,bl,dzlev,rlint)
 352
 353           do ivw=1,6
 354             rlwalins(ivw,nlev)=rlint(ivw)
 355           end do
 356
 357       else   !Top <---> Bottom
 358
 359           do ivw=1,4
 360             twin1D(ivw)=twin(ivw,1)
 361             twal_int(ivw)=twal(ivw,1,1)
 362           end do
 363
 364           twal_int(5)=tgrd(ngrd)
 365           twal_int(6)=trof(1)
 366
 367           call int_rlrad(emins,emwin,sigma,twal_int,twin1D, &
 368                          pwin,bw,bl,dzlev,rlint)
 369
 370           do ivw=1,6
 371             rlwalins(ivw,1)=rlint(ivw)
 372           end do
 373
 374       end if
 375
 376
 377 ! Calculation of the radiative fluxes
 378 ! -----------------------------------
 379
 380 !External vertical walls and windows
 381
 382         do ilev=1,nlev
 383          do ivw=1,4
 384          call radfluxs(radflux,albwal,rswal(ivw,ilev),     &
 385                             emwal,rlwal(ivw,ilev),sigma,   &
 386                             twal(ivw,nwal,ilev))
 387
 388          hrwalout(ivw,ilev)=radflux
 389
 390          hrwinout(ivw,ilev)=emwin*rlwal(ivw,ilev)- &
 391                             emwin*sigma*(twin(ivw,ilev)**4)
 392
 393
 394          end do ! ivw
 395         end do  ! ilev
 396
 397 !Roof
 398
 399         call radfluxs(radflux,albrof,rs,emrof,rl,sigma,trof(nrof))
 400
 401         hrrof=radflux
 402
 403 !Internal walls for intermediate rooms
 404
 405       if(nlev.gt.2) then
 406
 407         do ilev=2,nlev-1
 408          do ivw=1,4
 409
 410          call radfluxs(radflux,albins,rswalins(ivw,ilev),     &
 411                             emins,rlwalins(ivw,ilev),sigma,   &
 412                             twal(ivw,1,ilev))
 413
 414          hrwalins(ivw,ilev)=radflux
 415
 416          end do !ivw
 417
 418          call radfluxs(radflux,albins,rswalins(5,ilev), &
 419                               emins,rlwalins(5,ilev),sigma,&
 420                               tflo(nflo,ilev-1))
 421
 422          hrwalins(5,ilev)=radflux
 423
 424          call radfluxs(radflux,albins,rswalins(6,ilev), &
 425                               emins,rlwalins(6,ilev),sigma,&
 426                               tflo(1,ilev))
 427          hrwalins(6,ilev)=radflux
 428
 429        end do !ilev
 430
 431       end if
 432
 433
 434 !Internal walls for the bottom and the top room
 435 !
 436       if (nlev.ne.1) then
 437
 438 !bottom floor
 439
 440          do ivw=1,4
 441
 442             call radfluxs(radflux,albins,rswalins(ivw,1),  &
 443                             emins,rlwalins(ivw,1),sigma,   &
 444                             twal(ivw,1,1))
 445
 446             hrwalins(ivw,1)=radflux
 447
 448          end do
 449
 450
 451           call radfluxs(radflux,albins,rswalins(5,1),&
 452                            emins,rlwalins(5,1),sigma,&    !bottom
 453                            tgrd(ngrd))
 454
 455           hrwalins(5,1)=radflux
 456
 457
 458           call radfluxs(radflux,albins,rswalins(6,1),&
 459                            emins,rlwalins(6,1),sigma,&
 460                            tflo(1,1))
 461
 462           hrwalins(6,1)=radflux
 463
 464 !roof floor
 465
 466          do ivw=1,4
 467
 468           call radfluxs(radflux,albins,rswalins(ivw,nlev),     &
 469                                 emins,rlwalins(ivw,nlev),sigma,&
 470                                 twal(ivw,1,nlev))
 471
 472           hrwalins(ivw,nlev)=radflux
 473
 474          end do                                          !top
 475
 476
 477          call radfluxs(radflux,albins,rswalins(5,nlev),    &
 478                               emins,rlwalins(5,nlev),sigma,&
 479                               tflo(nflo,nlev-1))
 480
 481          hrwalins(5,nlev)=radflux
 482
 483          call radfluxs(radflux,albins,rswalins(6,nlev), &
 484                               emins,rlwalins(6,nlev),sigma,&
 485                               trof(1))
 486
 487          hrwalins(6,nlev)=radflux
 488
 489       else       ! Top <---> Bottom room
 490
 491          do ivw=1,4
 492
 493             call radfluxs(radflux,albins,rswalins(ivw,1),&
 494                             emins,rlwalins(ivw,1),sigma, &
 495                             twal(ivw,1,1))
 496
 497             hrwalins(ivw,1)=radflux
 498
 499          end do
 500
 501             call radfluxs(radflux,albins,rswalins(5,1),&
 502                            emins,rlwalins(5,1),sigma,  &
 503                            tgrd(ngrd))
 504
 505             hrwalins(5,1)=radflux
 506
 507             call radfluxs(radflux,albins,rswalins(6,nlev),     &
 508                                   emins,rlwalins(6,nlev),sigma,&
 509                                   trof(1))
 510             hrwalins(6,1)=radflux
 511
 512       end if
 513
 514
 515 !Windows
 516
 517          do ilev=1,nlev
 518           do ivw=1,4
 519              hrwinins(ivw,ilev)=emwin*rlwalins(ivw,ilev)-    &
 520                                 emwin*sigma*(twin(ivw,ilev)**4)
 521           end do
 522          end do
 523
 524
 525 ! Calculation of the sensible heat fluxes
 526 ! ---------------------------------------
 527
 528 !Vertical fluxes for walls
 529
 530         do ilev=1,nlev
 531          do ivw=1,4
 532
 533                call hsinsflux (2,2,tlev(ilev),twal(ivw,1,ilev),hs)
 534
 535                hswalins(ivw,ilev)=hs
 536
 537          end do ! ivw
 538         end do ! ilev
 539
 540
 541 !Vertical fluxes for windows
 542
 543         do ilev=1,nlev
 544
 545          do ivw=1,4
 546
 547                call hsinsflux (2,1,tlev(ilev),twin(ivw,ilev),hs)
 548
 549                hswinins(ivw,ilev)=hs
 550
 551          end do ! ivw
 552
 553         end do !ilev
 554
 555 !Horizontal fluxes
 556
 557       if (nlev.gt.2) then
 558
 559         do ilev=2,nlev-1
 560
 561                call hsinsflux (1,2,tlev(ilev),tflo(nflo,ilev-1),hs)
 562
 563                hswalins(5,ilev)=hs
 564
 565                call hsinsflux (1,2,tlev(ilev),tflo(1,ilev),hs)
 566
 567                hswalins(6,ilev)=hs
 568
 569         end do ! ilev
 570
 571       end if
 572
 573       if (nlev.ne.1) then
 574
 575                 call hsinsflux (1,2,tlev(1),tgrd(ngrd),hs)
 576
 577                 hswalins(5,1)=hs                                !Bottom room
 578
 579                 call hsinsflux (1,2,tlev(1),tflo(1,1),hs)
 580
 581                 hswalins(6,1)=hs
 582
 583                 call hsinsflux (1,2,tlev(nlev),tflo(nflo,nlev-1),hs)
 584
 585                 hswalins(5,nlev)=hs                             !Top room
 586
 587                 call hsinsflux (1,2,tlev(nlev),trof(1),hs)
 588
 589                 hswalins(6,nlev)=hs
 590
 591       else  ! Bottom<--->Top
 592
 593                 call hsinsflux (1,2,tlev(1),tgrd(ngrd),hs)
 594
 595                 hswalins(5,1)=hs
 596
 597                 call hsinsflux (1,2,tlev(nlev),trof(1),hs)
 598
 599                 hswalins(6,nlev)=hs
 600
 601       end if
 602
 603
 604 !Calculation of the temperature for the different surfaces
 605 ! --------------------------------------------------------
 606
 607 ! Vertical walls
 608
 609        swwal=1
 610        do ilev=1,nlev
 611         do ivw=1,4
 612
 613            htot(1)=hswalins(ivw,ilev)+hrwalins(ivw,ilev)
 614            htot(2)=hswalout(ivw,ilev)+hrwalout(ivw,ilev)
 615            gswal(ivw,ilev)=htot(2)
 616
 617            do iwal=1,nwal
 618               twal1D(iwal)=twal(ivw,iwal,ilev)
 619            end do
 620
 621            call wall(swwal,nwal,dt,dzwal,kwal,cswal,htot,twal1D)
 622
 623            do iwal=1,nwal
 624               twal(ivw,iwal,ilev)=twal1D(iwal)
 625            end do
 626
 627         end do ! ivw
 628        end do ! ilev
 629
 630 ! Windows
 631
 632        do ilev=1,nlev
 633         do ivw=1,4
 634
 635          htot(1)=hswinins(ivw,ilev)+hrwinins(ivw,ilev)
 636          htot(2)=hswinout(ivw,ilev)+hrwinout(ivw,ilev)
 637
 638          twin(ivw,ilev)=twin(ivw,ilev)+(dt/(cswin*thickwin))* &
 639                         (htot(1)+htot(2))
 640
 641         end do ! ivw
 642        end do ! ilev
 643
 644 ! Horizontal floors
 645
 646
 647       if (nlev.gt.1) then
 648        swwal=1
 649        do ilev=1,nlev-1
 650
 651           htot(1)=hrwalins(6,ilev)+hswalins(6,ilev)
 652           htot(2)=hrwalins(5,ilev+1)+hswalins(5,ilev+1)
 653
 654           do iflo=1,nflo
 655              tflo1D(iflo)=tflo(iflo,ilev)
 656           end do
 657
 658           call wall(swwal,nflo,dt,dzflo,kflo,csflo,htot,tflo1D)
 659
 660          do iflo=1,nflo
 661             tflo(iflo,ilev)=tflo1D(iflo)
 662          end do
 663
 664        end do ! ilev
 665       end if
 666
 667
 668 ! Ground
 669
 670         swwal=1
 671
 672         htot(1)=0.      !Diriclet b.c. at the internal boundary
 673         htot(2)=hswalins(5,1)+hrwalins(5,1)
 674
 675         do igrd=1,ngrd
 676            tgrd1D(igrd)=tgrd(igrd)
 677         end do
 678
 679          call wall(swwal,ngrd,dt,dzgrd,kgrd,csgrd,htot,tgrd1D)
 680
 681         do igrd=1,ngrd
 682            tgrd(igrd)=tgrd1D(igrd)
 683         end do
 684
 685
 686 ! Roof
 687
 688       swwal=1
 689
 690       htot(1)=hswalins(6,nlev)+hrwalins(6,nlev)
 691       htot(2)=hsrof+hrrof
 692       gsrof=htot(2)
 693
 694       do irof=1,nrof
 695          trof1D(irof)=trof(irof)
 696       end do
 697
 698       call wall(swwal,nrof,dt,dzrof,krof,csrof,htot,trof1D)
 699
 700       do irof=1,nrof
 701          trof(irof)=trof1D(irof)
 702       end do
 703
 704 ! Calculation of the heat fluxes and of the temperature of the rooms
 705 ! ------------------------------------------------------------------
 706
 707         do ilev=1,nlev
 708
 709          !Calculation of the heat generated by equipments and occupants
 710
 711          call fluxeqocc(nhourday,bw,bl,perflo,hsesf,hsequip,hseqocc(ilev),hleqocc(ilev))
 712
 713          !Calculation of the heat generated by natural ventilation
 714
 715           vollev=bw*bl*dzlev
 716           humdry=shumlev(ilev)/(1.-shumlev(ilev))
 717           rhoint=(press(ilev))/(rair*(1.+0.61*humdry)*tlev(ilev))
 718           cpint=cp*(1.+0.84*humdry)
 719
 720
 721           call fluxvent(cpint,rhoint,vollev,tlev(ilev),tout(ilev),     &
 722                         latent,humout(ilev),rhoout(ilev),shumlev(ilev),&
 723                         beta,hsvent(ilev),hlvent(ilev))
 724
 725          !Calculation of the heat generated by conduction
 726
 727            do iwal=1,6
 728              hswalins1D(iwal)=hswalins(iwal,ilev)
 729              surwal1D(iwal)=surwal(iwal,ilev)
 730           end do
 731
 732            do iwal=1,4
 733              hswinins1D(iwal)=hswinins(iwal,ilev)
 734            end do
 735
 736           call fluxcond(hswalins1D,hswinins1D,surwal1D,pwin,&
 737                         hscond(ilev))
 738
 739         !Calculation of the heat generated inside the room
 740
 741           call fluxroo(hseqocc(ilev),hleqocc(ilev),hsvent(ilev), &
 742                hlvent(ilev),hscond(ilev),hslev(ilev),hllev(ilev))
 743
 744
 745         !Evolution of the temperature and of the specific humidity
 746
 747           Qb=rhoint*cpint*vollev
 748
 749         ! temperature regulation
 750
 751           call regtemp(sw_cond,nhourday,dt,Qb,hslev(ilev),       &
 752                        tlev(ilev),timeon,timeoff,targtemp,gaptemp,hsneed(ilev))
 753
 754         ! humidity regulation
 755
 756           call reghum(sw_cond,nhourday,dt,vollev,rhoint,latent, &
 757                       hllev(ilev),shumlev(ilev),timeon,timeoff,&
 758                       targhum,gaphum,hlneed(ilev))
 759 !
 760 !performance of the air conditioning system for the test
 761 !
 762
 763           call air_cond(hsneed(ilev),hlneed(ilev),dt, &
 764                         hsout(ilev),hlout(ilev),consump(ilev), cop)
 765
 766           tlev(ilev)=tlev(ilev)+(dt/Qb)*(hslev(ilev)-hsneed(ilev))
 767
 768           shumlev(ilev)=shumlev(ilev)+(dt/(vollev*rhoint*latent))* &
 769                         (hllev(ilev)-hlneed(ilev))
 770
 771         end do !ilev
 772
 773         call consump_total(nzcanm,nlev,consumpbuild,hsoutbuild, &
 774                            hsout,consump)
 775
 776       return
 777       end subroutine BEM
 778
 779 !====6=8===============================================================72
 780 !====6=8===============================================================72
 781
 782         subroutine wall(swwall,nz,dt,dz,k,cs,flux,temp)
 783
 784 !______________________________________________________________________
 785
 786 !The aim of this subroutine is to solve the 1D heat fiffusion equation
 787 !for roof, walls and streets:
 788 !
 789 !       dT/dt=d/dz[K*dT/dz] where:
 790 !
 791 !       -T is the surface temperature(wall, street, roof)
 792 !       -Kz is the heat diffusivity inside the material.
 793 !
 794 !The resolution is done implicitly with a FV discretisation along the
 795 !different layers of the material:
 796
 797 !       ____________________________
 798 !     n             *
 799 !                   *
 800 !                   *
 801 !       ____________________________
 802 !    i+2
 803 !                   I+1
 804 !       ____________________________
 805 !    i+1
 806 !                    I                ==>   [T(I,n+1)-T(I,n)]/DT=
 807 !       ____________________________        [F(i+1)-F(i)]/DZI
 808 !    i
 809 !                   I-1               ==> A*T(n+1)=B where:
 810 !       ____________________________
 811 !   i-1              *                   * A is a TRIDIAGONAL matrix.
 812 !                    *                   * B=T(n)+S takes into account the sources.
 813 !                    *
 814 !     1 ____________________________
 815
 816 !________________________________________________________________
 817
 818         implicit none
 819
 820 !Input:
 821 !-----
 822         integer nz              !Number of layers inside the material
 823         real dt                 !Time step
 824         real dz(nz)             !Layer sizes [m]
 825         real cs(nz)             !Specific heat of the material [J/(m3.K)]
 826         real k(nz+1)            !Thermal conductivity in each layers (face) [W/(m.K)]
 827         real flux(2)            !Internal and external flux terms.
 828
 829 !Input-Output:
 830 !-------------
 831
 832         integer swwall          !swich for the physical coefficients calculation
 833         real temp(nz)           !Temperature at each layer
 834
 835 !Local:
 836 !-----
 837
 838       real a(-1:1,nz)          !  a(-1,*) lower diagonal      A(i,i-1)
 839                                !  a(0,*)  principal diagonal  A(i,i)
 840                                !  a(1,*)  upper diagonal      A(i,i+1).
 841
 842       real b(nz)               !Coefficients of the second term.
 843       real k1(20)
 844       real k2(20)
 845       real kc(20)
 846       save k1,k2,kc
 847       integer iz
 848
 849 !________________________________________________________________
 850 !
 851 !Calculation of the coefficients
 852
 853         if (swwall.eq.1) then
 854
 855            if (nz.gt.20) then
 856               write(*,*) 'number of layers in the walls/roofs too big ',nz
 857               write(*,*) 'please decrease under of',20
 858               stop
 859            endif
 860
 861            call wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
 862            swwall=0
 863
 864         end if
 865
 866 !Computation of the first value (iz=1) of A and B:
 867
 868                  a(-1,1)=0.
 869                  a(0,1)=1+k2(1)
 870                  a(1,1)=-k2(1)
 871
 872                  b(1)=temp(1)+flux(1)*kc(1)
 873
 874 !!
 875 !!We can fixed the internal temperature
 876 !!
 877 !!               a(-1,1)=0.
 878 !!               a(0,1)=1
 879 !!               a(1,1)=0.
 880 !!
 881 !!               b(1)=temp(1)
 882 !!
 883 !Computation of the internal values (iz=2,...,n-1) of A and B:
 884
 885         do iz=2,nz-1
 886                 a(-1,iz)=-k1(iz)
 887                 a(0,iz)=1+k1(iz)+k2(iz)
 888                 a(1,iz)=-k2(iz)
 889                 b(iz)=temp(iz)
 890         end do
 891
 892 !Computation of the external value (iz=n) of A and B:
 893
 894                 a(-1,nz)=-k1(nz)
 895                 a(0,nz)=1+k1(nz)
 896                 a(1,nz)=0.
 897
 898                 b(nz)=temp(nz)+flux(2)*kc(nz)
 899
 900 !Resolution of the system A*T(n+1)=B
 901
 902         call tridia(nz,a,b,temp)
 903
 904         return
 905         end  subroutine wall
 906
 907 !====6=8===============================================================72
 908 !====6=8===============================================================72
 909
 910         subroutine wall_coeff(nz,dt,dz,cs,k,k1,k2,kc)
 911
 912         implicit none
 913
 914 !---------------------------------------------------------------------
 915 !Input
 916 !-----
 917         integer nz              !Number of layers inside the material
 918         real dt                 !Time step
 919         real dz(nz)             !Layer sizes [m]
 920         real cs(nz)             !Specific heat of the material [J/(m3.K)]
 921         real k(nz+1)            !Thermal diffusivity in each layers (face) [W/(m.K)]
 922
 923
 924 !Input-Output
 925 !------------
 926
 927         real flux(2)            !Internal and external flux terms.
 928
 929
 930 !Output
 931 !------
 932         real k1(20)
 933         real k2(20)
 934         real kc(20)
 935
 936 !Local
 937 !-----
 938         integer iz
 939         real kf(nz)
 940
 941 !------------------------------------------------------------------
 942
 943         do iz=2,nz
 944          kc(iz)=dt/(dz(iz)*cs(iz))
 945          kf(iz)=2*k(iz)/(dz(iz)+dz(iz-1))
 946         end do
 947
 948         kc(1)=dt/(dz(1)*cs(1))
 949         kf(1)=2*k(1)/(dz(1))
 950
 951         do iz=1,nz
 952          k1(iz)=kc(iz)*kf(iz)
 953         end do
 954
 955         do iz=1,nz-1
 956          k2(iz)=kc(iz)*kf(iz+1)*cs(iz)/cs(iz+1)
 957         end do
 958
 959         return
 960         end subroutine wall_coeff
 961
 962 !====6=8===============================================================72
 963 !====6=8===============================================================72
 964         subroutine hsinsflux(swsurf,swwin,tin,tw,hsins)
 965
 966         implicit none
 967
 968 ! --------------------------------------------------------------------
 969 ! This routine computes the internal sensible heat flux.
 970 ! The swsurf, makes rhe difference between a vertical and a
 971 ! horizontal surface.
 972 ! The values of the heat conduction coefficients hc are obtained from the book
 973 ! "Energy Simulation in Building Design". J.A. Clarke.
 974 ! Adam Hilger, Bristol, 362 pp.
 975 ! --------------------------------------------------------------------
 976 !Input
 977 !----
 978         integer swsurf  !swich for the type of surface (horizontal/vertical)
 979         integer swwin   !swich for the type of surface (window/wall)
 980         real tin        !Inside temperature [K]
 981         real tw         !Internal wall temperature [K]
 982
 983
 984 !Output
 985 !------
 986         real hsins      !internal sensible heat flux [W/m2]
 987 !Local
 988 !-----
 989         real hc         !heat conduction coefficient [W/°C.m2]
 990 !--------------------------------------------------------------------
 991
 992         if (swsurf.eq.2) then   !vertical surface
 993          if (swwin.eq.1) then
 994             hc=5.678*0.99        !window surface (smooth surface)
 995          else
 996             hc=5.678*1.09        !wall surface (rough surface)
 997          endif
 998          hsins=hc*(tin-tw)
 999         endif
1000
1001         if (swsurf.eq.1)  then   !horizontal surface
1002          if (swwin.eq.1) then
1003            hc=5.678*0.99        !window surface (smooth surface)
1004          else
1005            hc=5.678*1.09        !wall surface (rough surface)
1006          endif
1007          hsins=hc*(tin-tw)
1008         endif
1009
1010         return
1011         end subroutine hsinsflux
1012 !====6=8===============================================================72
1013 !====6=8===============================================================72
1014
1015         subroutine int_rsrad(albwin,albwal,pwin,rswal,&
1016                              surwal,bw,bl,zw,rsint)
1017
1018 ! ------------------------------------------------------------------
1019         implicit none
1020 ! ------------------------------------------------------------------
1021
1022 !Input
1023 !-----
1024         real albwin             !albedo of the windows
1025         real albwal             !albedo of the internal wall
1026         real rswal(4)           !incoming short wave radiation [W/m2]
1027         real surwal(6)          !surface of the indoor walls [m2]
1028         real bw,bl              !width of the walls [m]
1029         real zw                 !height of the wall [m]
1030         real pwin               !window proportion
1031
1032 !Output
1033 !------
1034         real rsint(6)           !internal walls short wave radiation [W/m2]
1035
1036 !Local
1037 !-----
1038         real transmit   !transmittance of the direct/diffused radiation
1039         real rstr       !solar radiation transmitted through the windows
1040         real surtotwal  !total indoor surface of the walls in the room
1041         integer iw
1042         real b(6)       !second member for the system
1043         real a(6,6)     !matrix for the system
1044
1045 !-------------------------------------------------------------------
1046
1047
1048 ! Calculation of the solar radiation transmitted through windows
1049
1050             rstr = 0.
1051             do iw=1,4
1052                transmit=1.-albwin
1053                rstr = rstr+(surwal(iw)*pwin)*(transmit*rswal(iw))
1054             enddo
1055
1056 !We suppose that the radiation is spread isotropically within the
1057 !room when it passes through the windows, so the flux [W/m²] in every
1058 !wall is:
1059
1060             surtotwal=0.
1061             do iw=1,6
1062                surtotwal=surtotwal+surwal(iw)
1063             enddo
1064
1065             rstr=rstr/surtotwal
1066
1067 !Computation of the short wave radiation reaching the internal walls
1068
1069             call algebra_short(rstr,albwal,albwin,bw,bl,zw,pwin,a,b)
1070
1071             call gaussjbem(a,6,b,6)
1072
1073             do iw=1,6
1074                rsint(iw)=b(iw)
1075             enddo
1076
1077             return
1078             end subroutine int_rsrad
1079
1080 !====6=8===============================================================72
1081 !====6=8===============================================================72
1082
1083         subroutine int_rlrad(emwal,emwin,sigma,twal_int,twin,&
1084                              pwin,bw,bl,zw,rlint)
1085
1086 ! ------------------------------------------------------------------
1087         implicit none
1088 ! ------------------------------------------------------------------
1089
1090 !Input
1091 !-----
1092
1093         real emwal      !emissivity of the internal walls
1094         real emwin      !emissivity of the window
1095         real sigma      !Stefan-Boltzmann constant [W/m2.K4]
1096         real twal_int(6)!temperature of the first internal layers of a room [K]
1097         real twin(4)    !temperature of the windows [K]
1098         real bw         !width of the wall
1099         real bl         !length of the wall
1100         real zw         !height of the wall
1101         real pwin       !window proportion
1102
1103 !Output
1104 !------
1105
1106         real rlint(6)   !internal walls long wave radiation [W/m2]
1107
1108 !Local
1109 !------
1110
1111         real b(6)       !second member vector for the system
1112         real a(6,6)     !matrix for the system
1113         integer iw
1114 !----------------------------------------------------------------
1115
1116 !Compute the long wave radiation reachs the internal walls
1117
1118         call algebra_long(emwal,emwin,sigma,twal_int,twin,pwin,&
1119                           bw,bl,zw,a,b)
1120
1121         call gaussjbem(a,6,b,6)
1122
1123         do iw=1,6
1124            rlint(iw)=b(iw)
1125         enddo
1126
1127         return
1128         end subroutine int_rlrad
1129
1130 !====6=8===============================================================72
1131 !====6=8===============================================================72
1132
1133         subroutine algebra_short(rstr,albwal,albwin,aw,bw,zw,pwin,a,b)
1134
1135 !--------------------------------------------------------------------
1136 !This routine calculates the algebraic system that will be solved for
1137 !the computation of the total shortwave radiation that reachs every
1138 !indoor wall in a floor.
1139 !Write the matrix system ax=b to solve
1140 !
1141 !     -Rs(1)+a(1,2)Rs(2)+.................+a(1,6)Rs(6)=-Rs=b(1)
1142 !a(2,1)Rs(1)-      Rs(2)+.................+a(2,6)Rs(6)=-Rs=b(2)
1143 !a(3,1)Rs(1)+a(3,2)Rs(3)-Rs(3)+...........+a(3,6)Rs(6)=-Rs=b(3)
1144 !a(4,1)Rs(1)+.................-Rs(4)+.....+a(4,6)Rs(6)=-Rs=b(4)
1145 !a(5,1)Rs(1)+.......................-Rs(5)+a(5,6)Rs(6)=-Rs=b(5)
1146 !a(6,1)Rs(1)+....................................-R(6)=-Rs=b(6)
1147 !
1148 !This version suppose the albedo of the indoor walls is the same.
1149 !--------------------------------------------------------------------
1150         implicit none
1151 !--------------------------------------------------------------------
1152
1153 !Input
1154 !-----
1155         real rstr       !solar radiation transmitted through the windows
1156         real albwal     !albedo of the internal walls
1157         real albwin     !albedo of the windows.
1158         real bw         !length of the wall
1159         real aw         !width of the wall
1160         real zw         !height of the wall
1161         real fprl_int   !view factor
1162         real fnrm_int   !view factor
1163         real pwin       !window proportion
1164 !Output
1165 !------
1166         real a(6,6)             !Matrix for the system
1167         real b(6)               !Second member for the system
1168 !Local
1169 !-----
1170         integer iw,jw
1171         real albm               !averaged albedo
1172 !----------------------------------------------------------------
1173
1174 !Initialise the variables
1175
1176         do iw=1,6
1177            b(iw)= 0.
1178           do jw=1,6
1179            a(iw,jw)= 0.
1180           enddo
1181         enddo
1182
1183 !Calculation of the second member b
1184
1185         do iw=1,6
1186          b(iw)=-rstr
1187         end do
1188
1189 !Calculation of the averaged albedo
1190
1191         albm=pwin*albwin+(1-pwin)*albwal
1192
1193 !Calculation of the matrix a
1194
1195             a(1,1)=-1.
1196
1197             call fprl_ints(fprl_int,aw/bw,zw/bw)
1198
1199             a(1,2)=albm*fprl_int
1200
1201             call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1202
1203             a(1,3)=albm*(bw/aw)*fnrm_int
1204
1205             a(1,4)=a(1,3)
1206
1207             call fnrm_ints(fnrm_int,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1208
1209             a(1,5)=albwal*(bw/zw)*fnrm_int
1210
1211             a(1,6)=a(1,5)
1212
1213
1214             a(2,1)=a(1,2)
1215             a(2,2)=-1.
1216             a(2,3)=a(1,3)
1217             a(2,4)=a(1,4)
1218             a(2,5)=a(1,5)
1219             a(2,6)=a(1,6)
1220
1221
1222             call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1223
1224             a(3,1)=albm*(aw/bw)*fnrm_int
1225             a(3,2)=a(3,1)
1226             a(3,3)=-1.
1227
1228             call fprl_ints(fprl_int,zw/aw,bw/aw)
1229
1230             a(3,4)=albm*fprl_int
1231
1232             call fnrm_ints(fnrm_int,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1233
1234             a(3,5)=albwal*(aw/zw)*fnrm_int
1235             a(3,6)=a(3,5)
1236
1237
1238             a(4,1)=a(3,1)
1239             a(4,2)=a(3,2)
1240             a(4,3)=a(3,4)
1241             a(4,4)=-1.
1242             a(4,5)=a(3,5)
1243             a(4,6)=a(3,6)
1244
1245             call fnrm_ints(fnrm_int,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1246
1247             a(5,1)=albm*(zw/bw)*fnrm_int
1248
1249             a(5,2)=a(5,1)
1250
1251             call fnrm_ints(fnrm_int,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1252
1253             a(5,3)=albm*(zw/aw)*fnrm_int
1254
1255             a(5,4)=a(5,3)
1256             a(5,5)=-1.
1257
1258             call fprl_ints(fprl_int,aw/zw,bw/zw)
1259
1260             a(5,6)=albwal*fprl_int
1261
1262
1263             a(6,1)=a(5,1)
1264             a(6,2)=a(5,2)
1265             a(6,3)=a(5,3)
1266             a(6,4)=a(5,4)
1267             a(6,5)=a(5,6)
1268             a(6,6)=-1.
1269
1270         return
1271         end subroutine algebra_short
1272
1273 !====6=8===============================================================72
1274 !====6=8===============================================================72
1275
1276         subroutine algebra_long(emwal,emwin,sigma,twalint,twinint,&
1277                                 pwin,aw,bw,zw,a,b)
1278
1279 !--------------------------------------------------------------------
1280 !This routine computes the algebraic system that will be solved to
1281 !compute the longwave radiation that reachs the indoor
1282 !walls in a floor.
1283 !Write the matrix system ax=b to solve
1284 !
1285 !a(1,1)Rl(1)+.............................+Rl(6)=b(1)
1286 !a(2,1)Rl(1)+.................+Rl(5)+a(2,6)Rl(6)=b(2)
1287 !a(3,1)Rl(1)+.....+Rl(3)+...........+a(3,6)Rl(6)=b(3)
1288 !a(4,1)Rl(1)+...........+Rl(4)+.....+a(4,6)Rl(6)=b(4)
1289 !      Rl(1)+.......................+a(5,6)Rl(6)=b(5)
1290 !a(6,1)Rl(1)+Rl(2)+.................+a(6,6)Rl(6)=b(6)
1291 !
1292 !--------------------------------------------------------------------
1293         implicit none
1294
1295 !--------------------------------------------------------------------
1296
1297 !Input
1298 !-----
1299
1300         real pwin       !window proportion
1301         real emwal      !emissivity of the internal walls
1302         real emwin      !emissivity of the window
1303         real sigma      !Stefan-Boltzmann constant [W/m2.K4]
1304         real twalint(6) !temperature of the first internal layers of a room [K]
1305         real twinint(4) !temperature of the windows [K]
1306         real aw         !width of the wall
1307         real bw         !length of the wall
1308         real zw         !height of the wall
1309         real fprl_int   !view factor
1310         real fnrm_int   !view factor
1311         real fnrm_intx  !view factor
1312         real fnrm_inty  !view factor
1313
1314 !Output
1315 !------
1316         real b(6)       !second member vector for the system
1317         real a(6,6)     !matrix for the system
1318 !Local
1319 !-----
1320         integer iw,jw
1321         real b_wall(6)
1322         real b_wind(6)
1323         real emwal_av           !averadge emissivity of the wall
1324         real emwin_av           !averadge emissivity of the window
1325         real em_av              !averadge emissivity
1326         real twal_int(6)        !twalint
1327         real twin(4)            !twinint
1328 !------------------------------------------------------------------
1329
1330 !Initialise the variables
1331 !-------------------------
1332
1333          do iw=1,6
1334             b(iw)= 0.
1335             b_wall(iw)=0.
1336             b_wind(iw)=0.
1337           do jw=1,6
1338             a(iw,jw)= 0.
1339           enddo
1340          enddo
1341
1342          do iw=1,6
1343             twal_int(iw)=twalint(iw)
1344          enddo
1345
1346          do iw=1,4
1347             twin(iw)=twinint(iw)
1348          enddo
1349
1350 !Calculation of the averadge emissivities
1351 !-----------------------------------------
1352
1353         emwal_av=(1-pwin)*emwal
1354         emwin_av=pwin*emwin
1355         em_av=emwal_av+emwin_av
1356
1357 !Calculation of the second term for the walls
1358 !-------------------------------------------
1359
1360             call fprl_ints(fprl_int,aw/zw,bw/zw)
1361             call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1362             call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1363
1364             b_wall(1)=(emwal*sigma*(twal_int(5)**4)*           &
1365                  fprl_int)+                                    &
1366                  (sigma*(emwal_av*(twal_int(3)**4)+            &
1367                   emwal_av*(twal_int(4)**4))*                  &
1368                  (zw/aw)*fnrm_intx)+                           &
1369                  (sigma*(emwal_av*(twal_int(1)**4)+            &
1370                   emwal_av*(twal_int(2)**4))*                  &
1371                  (zw/bw)*fnrm_inty)
1372
1373             call fprl_ints(fprl_int,aw/zw,bw/zw)
1374             call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1375             call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1376
1377             b_wall(2)=(emwal*sigma*(twal_int(6)**4)*           &
1378                    fprl_int)+                                  &
1379                   (sigma*(emwal_av*(twal_int(3)**4)+           &
1380                   emwal_av*(twal_int(4)**4))*                  &
1381                  (zw/aw)*fnrm_intx)+                           &
1382                  (sigma*(emwal_av*(twal_int(1)**4)+            &
1383                  emwal_av*(twal_int(2)**4))*                   &
1384                  (zw/bw)*fnrm_inty)
1385
1386             call fprl_ints(fprl_int,zw/aw,bw/aw)
1387             call fnrm_ints(fnrm_intx,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1388             call fnrm_ints(fnrm_inty,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1389
1390             b_wall(3)=(emwal_av*sigma*(twal_int(4)**4)*        &
1391                   fprl_int)+                                   &
1392                  (sigma*(emwal_av*(twal_int(2)**4)+            &
1393                   emwal_av*(twal_int(1)**4))*                  &
1394                  (aw/bw)*fnrm_intx)+                           &
1395                  (sigma*(emwal*(twal_int(5)**4)+               &
1396                   emwal*(twal_int(6)**4))*                     &
1397                  (aw/zw)*fnrm_inty)
1398
1399             call fprl_ints(fprl_int,zw/aw,bw/aw)
1400             call fnrm_ints(fnrm_intx,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1401             call fnrm_ints(fnrm_inty,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1402
1403             b_wall(4)=(emwal_av*sigma*(twal_int(3)**4)*        &
1404                   fprl_int)+                                   &
1405                  (sigma*(emwal_av*(twal_int(2)**4)+            &
1406                   emwal_av*(twal_int(1)**4))*                  &
1407                  (aw/bw)*fnrm_intx)+                           &
1408                  (sigma*(emwal*(twal_int(5)**4)+               &
1409                   emwal*(twal_int(6)**4))*                     &
1410                  (aw/zw)*fnrm_inty)
1411
1412             call fprl_ints(fprl_int,aw/bw,zw/bw)
1413             call fnrm_ints(fnrm_intx,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1414             call fnrm_ints(fnrm_inty,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1415
1416             b_wall(5)=(emwal_av*sigma*(twal_int(2)**4)*        &
1417                   fprl_int)+                                   &
1418                  (sigma*(emwal_av*(twal_int(3)**4)+            &
1419                   emwal_av*(twal_int(4)**4))*                  &
1420                  (bw/aw)*fnrm_intx)+                           &
1421                  (sigma*(emwal*(twal_int(5)**4)+               &
1422                   emwal*(twal_int(6)**4))*                     &
1423                  (bw/zw)*fnrm_inty)
1424
1425             call fprl_ints(fprl_int,aw/bw,zw/bw)
1426             call fnrm_ints(fnrm_intx,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1427             call fnrm_ints(fnrm_inty,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1428
1429             b_wall(6)=(emwal_av*sigma*(twal_int(1)**4)*        &
1430                  fprl_int)+                                    &
1431                  (sigma*(emwal_av*(twal_int(3)**4)+            &
1432                   emwal_av*(twal_int(4)**4))*                  &
1433                  (bw/aw)*fnrm_intx)+                           &
1434                  (sigma*(emwal*(twal_int(5)**4)+               &
1435                  emwal*(twal_int(6)**4))*                      &
1436                  (bw/zw)*fnrm_inty)
1437
1438 !Calculation of the second term for the windows
1439 !---------------------------------------------
1440             call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1441             call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1442
1443             b_wind(1)=(sigma*(emwin_av*(twin(3)**4)+          &
1444                   emwin_av*(twin(4)**4))*                     &
1445                  (zw/aw)*fnrm_intx)+                          &
1446                  (sigma*(emwin_av*(twin(1)**4)+               &
1447                   emwin_av*(twin(2)**4))*                     &
1448                  (zw/bw)*fnrm_inty)
1449
1450             call fnrm_ints(fnrm_intx,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1451             call fnrm_ints(fnrm_inty,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1452
1453             b_wind(2)=(sigma*(emwin_av*(twin(3)**4)+          &
1454                   emwin_av*(twin(4)**4))*                     &
1455                  (zw/aw)*fnrm_intx)+                          &
1456                  (sigma*(emwin_av*(twin(1)**4)+               &
1457                   emwin_av*(twin(2)**4))*                     &
1458                  (zw/bw)*fnrm_inty)
1459
1460             call fprl_ints(fprl_int,zw/aw,bw/aw)
1461             call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1462
1463             b_wind(3)=emwin_av*sigma*(twin(4)**4)*            &
1464                  fprl_int+(sigma*(emwin_av*                   &
1465                  (twin(2)**4)+emwin_av*(twin(1)**4))*         &
1466                  (aw/bw)*fnrm_int)
1467
1468             call fprl_ints(fprl_int,zw/aw,bw/aw)
1469             call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1470
1471             b_wind(4)=emwin_av*sigma*(twin(3)**4)*            &
1472                  fprl_int+(sigma*(emwin_av*                   &
1473                   (twin(2)**4)+emwin_av*(twin(1)**4))*        &
1474                  (aw/bw)*fnrm_int)
1475
1476             call fprl_ints(fprl_int,aw/bw,zw/bw)
1477             call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1478
1479             b_wind(5)=emwin_av*sigma*(twin(2)**4)*            &
1480                  fprl_int+(sigma*(emwin_av*                   &
1481                  (twin(3)**4)+emwin_av*(twin(4)**4))*         &
1482                  (bw/aw)*fnrm_int)
1483
1484             call fprl_ints(fprl_int,aw/bw,zw/bw)
1485             call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1486
1487             b_wind(6)=emwin_av*sigma*(twin(1)**4)*            &
1488                  fprl_int+(sigma*(emwin_av*                   &
1489                  (twin(3)**4)+emwin_av*(twin(4)**4))*         &
1490                  (bw/aw)*fnrm_int)
1491
1492 !Calculation of the total b term
1493 !-------------------------------
1494
1495         do iw=1,6
1496          b(iw)=b_wall(iw)+b_wind(iw)
1497         end do
1498
1499
1500 !Calculation of the matrix of the system
1501 !----------------------------------------
1502
1503          call fnrm_ints(fnrm_int,bw/aw,zw/aw,(bw*bw+zw*zw)/(aw*aw))
1504
1505          a(1,1)=(em_av-1.)*(zw/bw)*fnrm_int
1506
1507          a(1,2)=a(1,1)
1508
1509          call fnrm_ints(fnrm_int,aw/bw,zw/bw,(aw*aw+zw*zw)/(bw*bw))
1510
1511          a(1,3)=(em_av-1.)*(zw/aw)*fnrm_int
1512
1513          a(1,4)=a(1,3)
1514
1515          call fprl_ints(fprl_int,aw/zw,bw/zw)
1516
1517          a(1,5)=(emwal-1.)*fprl_int
1518          a(1,6)=1.
1519
1520          a(2,1)=a(1,1)
1521          a(2,2)=a(1,2)
1522          a(2,3)=a(1,3)
1523          a(2,4)=a(1,4)
1524          a(2,5)=1.
1525          a(2,6)=a(1,5)
1526
1527          call fnrm_ints(fnrm_int,bw/zw,aw/zw,(bw*bw+aw*aw)/(zw*zw))
1528
1529          a(3,1)=(em_av-1.)*(aw/bw)*fnrm_int
1530
1531          a(3,2)=a(3,1)
1532          a(3,3)=1.
1533
1534          call fprl_ints(fprl_int,zw/aw,bw/aw)
1535
1536          a(3,4)=(em_av-1.)*fprl_int
1537
1538          call fnrm_ints(fnrm_int,zw/bw,aw/bw,(aw*aw+zw*zw)/(bw*bw))
1539
1540          a(3,5)=(emwal-1.)*(aw/zw)*fnrm_int
1541
1542          a(3,6)=a(3,5)
1543
1544          a(4,1)=a(3,1)
1545          a(4,2)=a(3,2)
1546          a(4,3)=a(3,4)
1547          a(4,4)=1.
1548          a(4,5)=a(3,5)
1549          a(4,6)=a(3,6)
1550
1551          a(5,1)=1.
1552
1553          call fprl_ints(fprl_int,aw/bw,zw/bw)
1554
1555          a(5,2)=(em_av-1.)*fprl_int
1556
1557          call fnrm_ints(fnrm_int,aw/zw,bw/zw,(aw*aw+bw*bw)/(zw*zw))
1558
1559          a(5,3)=(em_av-1.)*(bw/aw)*fnrm_int
1560
1561          a(5,4)=a(5,3)
1562
1563          call fnrm_ints(fnrm_int,zw/aw,bw/aw,(bw*bw+zw*zw)/(aw*aw))
1564
1565          a(5,5)=(emwal-1.)*(bw/zw)*fnrm_int
1566
1567          a(5,6)=a(5,5)
1568
1569          a(6,1)=a(5,2)
1570          a(6,2)=1.
1571          a(6,3)=a(5,3)
1572          a(6,4)=a(5,4)
1573          a(6,5)=a(5,5)
1574          a(6,6)=a(6,5)
1575
1576       return
1577       end subroutine algebra_long
1578
1579 !====6=8===============================================================72
1580 !====6=8===============================================================72
1581
1582
1583         subroutine fluxroo(hseqocc,hleqocc,hsvent,hlvent, &
1584                            hscond,hslev,hllev)
1585
1586 !-----------------------------------------------------------------------
1587 !This routine calculates the heat flux generated inside the room
1588 !and the heat ejected to the atmosphere.
1589 !----------------------------------------------------------------------
1590
1591         implicit none
1592
1593 !-----------------------------------------------------------------------
1594
1595 !Input
1596 !-----
1597         real hseqocc            !sensible heat generated by equipments and occupants [W]
1598         real hleqocc            !latent heat generated by occupants [W]
1599         real hsvent             !sensible heat generated by natural ventilation [W]
1600         real hlvent             !latent heat generated by natural ventilation [W]
1601         real hscond             !sensible heat generated by wall conduction
1602
1603 !Output
1604 !------
1605         real hslev              !sensible heat flux generated inside the room [W]
1606         real hllev              !latent heat flux generatd inside the room
1607
1608
1609 !Calculation of the total sensible heat generated inside the room
1610
1611         hslev=hseqocc+hsvent+hscond
1612
1613 !Calculation of the total latent heat generated inside the room
1614
1615         hllev=hleqocc+hlvent
1616
1617         return
1618         end subroutine fluxroo
1619
1620 !====6=8===============================================================72
1621 !====6=8===============================================================72
1622
1623         subroutine phirat(nhourday,rocc)
1624
1625 !----------------------------------------------------------------------
1626 !This routine calculates the occupation ratio of a floor
1627 !By now we suppose a constant value
1628 !----------------------------------------------------------------------
1629
1630         implicit none
1631
1632 !Input
1633 !-----
1634
1635         real nhourday   ! number of hours from midnight (local time)
1636
1637 !Output
1638 !------
1639         real rocc       !value between 0 and 1
1640
1641 !!TEST
1642         rocc=1.
1643
1644         return
1645         end subroutine phirat
1646
1647 !====6=8===============================================================72
1648 !====6=8===============================================================72
1649
1650         subroutine phiequ(nhourday,hsesf,hsequip,hsequ)
1651
1652 !----------------------------------------------------------------------
1653 !This routine calculates the sensible heat gain from equipments
1654 !----------------------------------------------------------------------
1655         implicit none
1656 !Input
1657 !-----
1658
1659         real nhourday ! number of hours from midnight, Local time
1660         real, intent(in) :: hsesf
1661         real, intent(in), dimension(24) :: hsequip
1662
1663 !Output
1664 !------
1665         real hsequ    !sensible heat gain from equipment [Wm¯2]
1666
1667 !---------------------------------------------------------------------
1668
1669         hsequ = hsequip(int(nhourday)+1) * hsesf
1670
1671         end subroutine phiequ
1672 !====6=8===============================================================72
1673 !====6=8===============================================================72
1674
1675         subroutine fluxeqocc(nhourday,bw,bl,perflo,hsesf,hsequip,hseqocc,hleqocc)
1676
1677         implicit none
1678
1679 !---------------------------------------------------------------------
1680 !This routine calculates the sensible and the latent heat flux
1681 !generated by equipments and occupants
1682 !---------------------------------------------------------------------
1683
1684 !Input
1685 !-----
1686         real bw                 !Room width [m]
1687         real bl                 !Room lengzh [m]
1688         real nhourday           !number of hours from the beginning of the day
1689         real, intent(in) :: perflo ! Peak number of occupants per unit floor area
1690         real, intent(in) :: hsesf
1691         real, intent(in), dimension(24) :: hsequip
1692
1693 !Output
1694 !------
1695         real hseqocc            !sensible heat generated by equipments and occupants [W]
1696         real hleqocc            !latent heat generated by occupants [W]
1697 !Local
1698 !-----
1699         real Af                 !Air conditioned floor area [m2]
1700         real rocc               !Occupation ratio of the floor [0,1]
1701         real hsequ              !Heat generated from equipments
1702
1703         real hsocc              !Sensible heat generated by a person [W/Person]
1704                                 !Source Boundary Layer Climates,page 195 (book)
1705         parameter (hsocc=160.)
1706
1707         real hlocc              !Latent heat generated by a person [W/Person]
1708                                 !Source Boundary Layer Climates,page 225 (book)
1709         parameter (hlocc=1.96e6/86400.)
1710
1711 !------------------------------------------------------------------
1712 !                       Sensible heat flux
1713 !                       ------------------
1714
1715          Af=bw*bl
1716
1717          call phirat(nhourday,rocc)
1718
1719          call phiequ(nhourday,hsesf,hsequip,hsequ)
1720
1721          hseqocc=Af*rocc*perflo*hsocc+Af*hsequ
1722
1723 !
1724 !                       Latent heat
1725 !                       -----------
1726 !
1727
1728          hleqocc=Af*rocc*perflo*hlocc
1729
1730         return
1731         end subroutine fluxeqocc
1732
1733 !====6=8===============================================================72
1734 !====6=8===============================================================72
1735
1736         subroutine fluxvent(cpint,rhoint,vollev,tlev,tout,latent,&
1737                             humout,rhoout,humlev,beta,hsvent,hlvent)
1738
1739         implicit none
1740
1741 !---------------------------------------------------------------------
1742 !This routine calculates the sensible and the latent heat flux
1743 !generated by natural ventilation
1744 !---------------------------------------------------------------------
1745
1746 !Input
1747 !-----
1748         real cpint              !specific heat of the indoor air [J/kg.K]
1749         real rhoint             !density of the indoor air [Kg/m3]
1750         real vollev             !volume of the room [m3]
1751         real tlev               !Room temperature [K]
1752         real tout               !outside air temperature [K]
1753         real latent             !latent heat of evaporation [J/Kg]
1754         real humout             !outside absolute humidity [Kgwater/Kgair]
1755         real rhoout             !air density [kg/m3]
1756         real humlev             !Specific humidity of the indoor air [Kgwater/Kgair]
1757         real, intent(in) :: beta!Thermal efficiency of the heat exchanger
1758
1759 !Output
1760 !------
1761         real hsvent             !sensible heat generated by natural ventilation [W]
1762         real hlvent             !latent heat generated by natural ventilation [W]
1763
1764 !Local
1765 !-----
1766
1767 !----------------------------------------------------------------------
1768
1769 !                       Sensible heat flux
1770 !                       ------------------
1771
1772         hsvent=(1.-beta)*cpint*rhoint*(vollev/3600.)*  &
1773                (tout-tlev)
1774
1775 !                       Latent heat flux
1776 !                       ----------------
1777
1778         hlvent=(1.-beta)*latent*rhoint*(vollev/3600.)* &
1779                (humout-humlev)
1780
1781
1782         return
1783         end subroutine fluxvent
1784
1785 !====6=8===============================================================72
1786 !====6=8===============================================================72
1787
1788         subroutine fluxcond(hswalins,hswinins,surwal,pwin,hscond)
1789
1790         implicit none
1791
1792 !---------------------------------------------------------------------
1793 !This routine calculates the sensible heat flux generated by
1794 !wall conduction.
1795 !---------------------------------------------------------------------
1796
1797 !Input
1798 !-----
1799         real hswalins(6)        !sensible heat at the internal layers of the wall [W/m2]
1800         real hswinins(4)        !internal window sensible heat flux [W/m2]
1801         real surwal(6)          !surfaces of the room walls [m2]
1802         real pwin               !window proportion
1803
1804
1805 !Output
1806 !------
1807
1808         real hscond             !sensible heat generated by wall conduction [W]
1809
1810 !Local
1811 !-----
1812
1813         integer ivw
1814
1815 !----------------------------------------------------------------------
1816
1817           hscond=0.
1818
1819         do ivw=1,4
1820            hscond=hscond+surwal(ivw)*(1-pwin)*hswalins(ivw)+ &
1821                   surwal(ivw)*pwin*hswinins(ivw)
1822         end do
1823
1824         do ivw=5,6
1825            hscond=hscond+surwal(ivw)*hswalins(ivw)
1826         end do
1827 !
1828 !Finally we must change the sign in hscond to be proportional
1829 !to the difference (Twall-Tindoor).
1830 !
1831         hscond=(-1)*hscond
1832
1833         return
1834         end subroutine fluxcond
1835
1836 !====6=8===============================================================72
1837 !====6=8===============================================================72
1838
1839         subroutine regtemp(swcond,nhourday,dt,Qb,hsroo,          &
1840                            tlev,timeon,timeoff,targtemp,gaptemp,hsneed)
1841
1842         implicit none
1843
1844 !---------------------------------------------------------------------
1845 !This routine calculates the sensible heat fluxes,
1846 !after anthropogenic regulation (air conditioning)
1847 !---------------------------------------------------------------------
1848
1849 !Input:
1850 !-----.
1851         integer swcond       !swich air conditioning
1852         real nhourday        !number of hours from the beginning of the day real
1853         real dt              !time step [s]
1854         real Qb              !overall heat capacity of the indoor air [J/K]
1855         real hsroo           !sensible heat flux generated inside the room [W]
1856         real tlev            !room air temperature [K]
1857         real, intent(in) :: timeon  ! Initial local time of A/C systems
1858         real, intent(in) :: timeoff ! Ending local time of A/C systems
1859         real, intent(in) :: targtemp! Target temperature of A/C systems
1860         real, intent(in) :: gaptemp ! Comfort range of indoor temperature
1861
1862
1863 !Local:
1864 !-----.
1865
1866         real templev         !hipotetical room air temperature [K]
1867         real alpha           !variable to control the heating/cooling of
1868                              !the air conditining system
1869 !Output:
1870 !-----.
1871         real hsneed          !sensible heat extracted to the indoor air [W]
1872 !---------------------------------------------------------------------
1873 !initialize variables
1874 !---------------------
1875         templev = 0.
1876         alpha   = 0.
1877
1878         if (swcond.eq.0) then ! there is not air conditioning in the floor
1879             hsneed = 0.
1880             goto 100
1881         else
1882             if ((nhourday.ge.timeon).and.(nhourday.le.timeoff)) then
1883                templev=tlev+(dt/Qb)*hsroo
1884                goto 200
1885             else
1886                hsneed = 0.     ! air conditioning is switched off
1887                goto 100
1888             endif
1889         endif
1890
1891 200     continue
1892
1893         if (abs(templev-targtemp).le.gaptemp) then
1894            hsneed = 0.
1895         else
1896            if (templev.gt.(targtemp+gaptemp)) then
1897               hsneed=hsroo-(Qb/dt)*(targtemp+gaptemp-tlev)
1898               alpha=(abs(hsneed-hsroo)/Qb)
1899               if (alpha.gt.temp_rat) then
1900                   hsneed=hsroo+temp_rat*Qb
1901                   goto 100
1902               else
1903                   goto 100
1904               endif
1905            else
1906               hsneed=hsroo-(Qb/dt)*(targtemp-gaptemp-tlev)
1907               alpha=(abs(hsneed-hsroo)/Qb)
1908               if (alpha.gt.temp_rat) then
1909                   hsneed=hsroo-temp_rat*Qb
1910                   goto 100
1911               else
1912                   goto 100
1913               endif
1914            endif
1915         endif
1916
1917 100     continue
1918         return
1919         end subroutine regtemp
1920
1921 !====6=8==============================================================72
1922 !====6=8==============================================================72
1923
1924          subroutine reghum(swcond,nhourday,dt,volroo,rhoint,latent, &
1925                            hlroo,shumroo,timeon,timeoff,targhum,gaphum,hlneed)
1926
1927          implicit none
1928
1929 !---------------------------------------------------------------------
1930 !This routine calculates the latent heat fluxes,
1931 !after anthropogenic regulation (air conditioning)
1932 !---------------------------------------------------------------------
1933
1934 !Input:
1935 !-----.
1936         integer swcond    !swich air conditioning
1937         real nhourday     !number of hours from the beginning of the day real[h]
1938         real dt           !time step [s]
1939         real volroo       !volume of the room [m3]
1940         real rhoint       !density of the internal air [Kg/m3]
1941         real latent       !latent heat of evaporation [J/Kg]
1942         real hlroo        !latent heat flux generated inside the room [W]
1943         real shumroo      !specific humidity of the indoor air [kg/kg]
1944         real, intent(in) :: timeon  ! Initial local time of A/C systems
1945         real, intent(in) :: timeoff ! Ending local time of A/C systems
1946         real, intent(in) :: targhum ! Target humidity of the A/C systems
1947         real, intent(in) :: gaphum  ! comfort range of the specific humidity
1948
1949 !Local:
1950 !-----.
1951
1952         real humlev       !hipotetical specific humidity of the indoor [kg/kg]
1953         real betha        !variable to control the drying/moistening of
1954                           !the air conditioning system
1955 !Output:
1956 !-----.
1957         real hlneed       !latent heat extracted to the indoor air [W]
1958 !------------------------------------------------------------------------
1959 !initialize variables
1960 !---------------------
1961         humlev = 0.
1962         betha  = 0.
1963
1964         if (swcond.eq.0) then ! there is not air conditioning in the floor
1965             hlneed = 0.
1966             goto 100
1967         else
1968             if ((nhourday.ge.timeon).and.(nhourday.le.timeoff)) then
1969                humlev=shumroo+(dt/(latent*rhoint*volroo))*hlroo
1970                goto 200
1971             else
1972                hlneed = 0.     ! air conditioning is switched off
1973                goto 100
1974             endif
1975         endif
1976
1977 200     continue
1978
1979         if (abs(humlev-targhum).le.gaphum) then
1980            hlneed = 0.
1981         else
1982            if (humlev.gt.(targhum+gaphum)) then
1983               hlneed=hlroo-((latent*rhoint*volroo)/dt)* &
1984                           (targhum+gaphum-shumroo)
1985               betha=abs(hlneed-hlroo)/(latent*rhoint*volroo)
1986               if (betha.gt.hum_rat) then
1987                   hlneed=hlroo+hum_rat*(latent*rhoint*volroo)
1988                   goto 100
1989               else
1990                   goto 100
1991               endif
1992            else
1993               hlneed=hlroo-((latent*rhoint*volroo)/dt)* &
1994                           (targhum-gaphum-shumroo)
1995               betha=abs(hlneed-hlroo)/(latent*rhoint*volroo)
1996               if (betha.gt.hum_rat) then
1997                   hlneed=hlroo-hum_rat*(latent*rhoint*volroo)
1998                   goto 100
1999               else
2000                   goto 100
2001               endif
2002            endif
2003         endif
2004
2005 100     continue
2006         return
2007         end subroutine reghum
2008
2009 !====6=8==============================================================72
2010 !====6=8==============================================================72
2011
2012          subroutine air_cond(hsneed,hlneed,dt,hsout,hlout,consump,cop)
2013
2014          implicit none
2015
2016 !
2017 !Performance of the air conditioning system
2018 !
2019 !INPUT/OUTPUT VARIABLES
2020          real, intent(in) :: cop
2021 !
2022 !INPUT/OUTPUT VARIABLES
2023 !
2024          real hsneed     !sensible heat that is necessary for cooling/heating
2025                          !the indoor air temperature [W]
2026          real hlneed     !latent heat that is necessary for controling
2027                          !the humidity of the indoor air [W]
2028          real dt         !time step [s]
2029 !
2030 !OUTPUT VARIABLES
2031 !
2032          real hsout      !sensible heat pumped out into the atmosphere [W]
2033          real hlout      !latent heat pumped out into the atmosphere [W]
2034          real consump    !Electrical consumption of the air conditioning system [W]
2035
2036
2037 !
2038 !Performance of the air conditioning system
2039 !
2040          if (hsneed.gt.0) then         ! air conditioning is cooling
2041                                        ! and the heat is pumped out into the atmosphere
2042           hsout=(1/cop)*(abs(hsneed)+abs(hlneed))+hsneed
2043           hlout=hlneed
2044           consump=(1./cop)*(abs(hsneed)+abs(hlneed))
2045 !!        hsout=0.
2046 !!        hlout=0.
2047
2048          else if(hsneed.eq.0.) then !air conditioning is not working to regulate the indoor temperature
2049                hlneed=0.       !no humidity regulation is considered
2050                hsout=0.        !no output into the atmosphere (sensible heat)
2051                hlout=0.        !no output into the atmosphere (latent heat)
2052                consump=0.      !no electrical consumption
2053
2054               else  !! hsneed < 0. !air conditioning is heating
2055                hlneed=0.           !no humidity regulation is considered
2056                hlout=0.            !no output into the atmosphere (latent heat)
2057                consump=(1./cop)*(abs(hsneed)+abs(hlneed))
2058 !
2059 !!We have two possibilities
2060 !
2061 !!             hsout=(1./cop)*(abs(hsneed)+abs(hlneed)) !output into the atmosphere
2062                hsout=0.                            !no output into the atmosphere
2063          end if
2064
2065          return
2066          end subroutine air_cond
2067
2068 !====6=8==============================================================72
2069 !====6=8==============================================================72
2070
2071         subroutine consump_total(nzcanm,nlev,consumpbuild,hsoutbuild, &
2072                                  hsout,consump)
2073
2074         implicit none
2075
2076 !-----------------------------------------------------------------------
2077 !Compute the total consumption in kWh/s (1kWh=3.6e+6 J) and sensible heat
2078 !ejected into the atmosphere per building
2079 !------------------------------------------------------------------------
2080 !
2081 !INPUT VARIABLES
2082 !
2083 !
2084         integer nzcanm            !Maximum number of vertical levels in the urban grid
2085         real hsout(nzcanm)        !sensible heat emitted outside the room [W]
2086         real consump(nzcanm)      !Electricity consumption for the a.c. in each floor[W]
2087 !
2088 !OUTPUT VARIABLES
2089 !
2090         real consumpbuild         !Energetic consumption for the entire building[kWh/s]
2091         real hsoutbuild           !Total sensible heat ejected into the atmosphere
2092                                   !by the air conditioning systems per building [W]
2093 !
2094 !LOCAL  VARIABLES
2095 !
2096         integer ilev
2097
2098 !
2099 !INPUT VARIABLES
2100 !
2101         integer nlev
2102
2103 !
2104 !INITIALIZE VARIABLES
2105 !
2106         consumpbuild=0.
2107         hsoutbuild=0.
2108 !
2109         do ilev=1,nlev
2110            consumpbuild=consumpbuild+consump(ilev)
2111            hsoutbuild=hsoutbuild+hsout(ilev)
2112         enddo !ilev
2113
2114         consumpbuild=consumpbuild/(3.6e+06)
2115
2116         return
2117         end subroutine consump_total
2118 !====6=8==============================================================72
2119 !====6=8==============================================================72
2120         subroutine tridia(n,a,b,x)
2121
2122 !     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2123 !     +    by A. Clappier,     EPFL, CH 1015 Lausanne                  +
2124 !     +                        phone: ++41-(0)21-693-61-60             +
2125 !     +                        email:alain.clappier@epfl.ch            +
2126 !     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2127
2128 ! ----------------------------------------------------------------------
2129 !        Resolution of a * x = b    where a is a tridiagonal matrix
2130 !
2131 ! ----------------------------------------------------------------------
2132
2133         implicit none
2134
2135 ! Input
2136         integer n
2137         real a(-1:1,n)           !  a(-1,*) lower diagonal      A(i,i-1)
2138                                !  a(0,*)  principal diagonal  A(i,i)
2139                                !  a(1,*)  upper diagonal      A(i,i+1)
2140         real b(n)
2141
2142 ! Output
2143         real x(n)
2144
2145 ! Local
2146         integer i
2147
2148 ! ----------------------------------------------------------------------
2149
2150         do i=n-1,1,-1
2151            b(i)=b(i)-a(1,i)*b(i+1)/a(0,i+1)
2152            a(0,i)=a(0,i)-a(1,i)*a(-1,i+1)/a(0,i+1)
2153         enddo
2154
2155         do i=2,n
2156            b(i)=b(i)-a(-1,i)*b(i-1)/a(0,i-1)
2157         enddo
2158
2159         do i=1,n
2160            x(i)=b(i)/a(0,i)
2161         enddo
2162
2163         return
2164         end subroutine tridia
2165 !====6=8===============================================================72
2166 !====6=8===============================================================72
2167
2168        subroutine gaussjbem(a,n,b,np)
2169
2170 ! ----------------------------------------------------------------------
2171 ! This routine solve a linear system of n equations of the form
2172 !              A X = B
2173 !  where  A is a matrix a(i,j)
2174 !         B a vector and X the solution
2175 ! In output b is replaced by the solution
2176 ! ----------------------------------------------------------------------
2177
2178        implicit none
2179
2180 ! ----------------------------------------------------------------------
2181 ! INPUT:
2182 ! ----------------------------------------------------------------------
2183        integer np
2184        real a(np,np)
2185
2186 ! ----------------------------------------------------------------------
2187 ! OUTPUT:
2188 ! ----------------------------------------------------------------------
2189        real b(np)
2190
2191 ! ----------------------------------------------------------------------
2192 ! LOCAL:
2193 ! ----------------------------------------------------------------------
2194       integer nmax
2195       parameter (nmax=150)
2196
2197       real big,dum
2198       integer i,icol,irow
2199       integer j,k,l,ll,n
2200       integer ipiv(nmax)
2201       real pivinv
2202
2203 ! ----------------------------------------------------------------------
2204 ! END VARIABLES DEFINITIONS
2205 ! ----------------------------------------------------------------------
2206
2207        do j=1,n
2208           ipiv(j)=0.
2209        enddo
2210
2211       do i=1,n
2212          big=0.
2213          do j=1,n
2214             if(ipiv(j).ne.1)then
2215                do k=1,n
2216                   if(ipiv(k).eq.0)then
2217                      if(abs(a(j,k)).ge.big)then
2218                         big=abs(a(j,k))
2219                         irow=j
2220                         icol=k
2221                      endif
2222                   elseif(ipiv(k).gt.1)then
2223                      pause 'singular matrix in gaussjbem'
2224                   endif
2225                enddo
2226             endif
2227          enddo
2228
2229          ipiv(icol)=ipiv(icol)+1
2230
2231          if(irow.ne.icol)then
2232             do l=1,n
2233                dum=a(irow,l)
2234                a(irow,l)=a(icol,l)
2235                a(icol,l)=dum
2236             enddo
2237
2238             dum=b(irow)
2239             b(irow)=b(icol)
2240             b(icol)=dum
2241
2242          endif
2243
2244          if(a(icol,icol).eq.0)pause 'singular matrix in gaussjbem'
2245
2246          pivinv=1./a(icol,icol)
2247          a(icol,icol)=1
2248
2249          do l=1,n
2250             a(icol,l)=a(icol,l)*pivinv
2251          enddo
2252
2253          b(icol)=b(icol)*pivinv
2254
2255          do ll=1,n
2256             if(ll.ne.icol)then
2257                dum=a(ll,icol)
2258                a(ll,icol)=0.
2259                do l=1,n
2260                   a(ll,l)=a(ll,l)-a(icol,l)*dum
2261                enddo
2262
2263                b(ll)=b(ll)-b(icol)*dum
2264
2265             endif
2266          enddo
2267       enddo
2268
2269       return
2270       end subroutine gaussjbem
2271
2272 !====6=8===============================================================72
2273 !====6=8===============================================================72
2274
2275       subroutine radfluxs(radflux,alb,rs,em,rl,sigma,twal)
2276
2277       implicit none
2278 !-------------------------------------------------------------------
2279 !This function calculates the radiative fluxe at a surface
2280 !-------------------------------------------------------------------
2281
2282
2283         real alb        !albedo of the surface
2284         real rs         !shor wave radiation
2285         real em         !emissivity of the surface
2286         real rl         !lon wave radiation
2287         real sigma      !parameter (wall is not black body) [W/m2.K4]
2288         real twal       !wall temperature [K]
2289         real radflux
2290
2291          radflux=(1.-alb)*rs+em*rl-em*sigma*twal**4
2292
2293       return
2294       end subroutine radfluxs
2295
2296 !====6=8==============================================================72
2297 !====6=8==============================================================72
2298 !
2299 !       we define the view factors fprl and fnrm, which are the angle
2300 !       factors between two equal and parallel planes, fprl, and two
2301 !       equal and orthogonal planes, fnrm, respectively
2302 !
2303         subroutine fprl_ints(fprl_int,vx,vy)
2304
2305         implicit none
2306
2307         real vx,vy
2308         real fprl_int
2309
2310         fprl_int=(2./(3.141592653*vx*vy))*                       &
2311              (log(sqrt((1.+vx*vx)*(1.+vy*vy)/(1.+vx*vx+vy*vy)))+ &
2312               (vy*sqrt(1.+vx*vx)*atan(vy/sqrt(1.+vx*vx)))+       &
2313               (vx*sqrt(1.+vy*vy)*atan(vx/sqrt(1.+vy*vy)))-       &
2314               vy*atan(vy)-vx*atan(vx))
2315
2316         return
2317         end subroutine fprl_ints
2318
2319 !====6=8==============================================================72
2320 !====6=8==============================================================72
2321 !
2322 !       we define the view factors fprl and fnrm, which are the angle
2323 !       factors between two equal and parallel planes, fprl, and two
2324 !       equal and orthogonal planes, fnrm, respectively
2325 !
2326
2327         subroutine fnrm_ints(fnrm_int,wx,wy,wz)
2328
2329         implicit none
2330
2331         real wx,wy,wz
2332         real fnrm_int
2333
2334         fnrm_int=(1./(3.141592653*wy))*(wy*atan(1./wy)+wx*atan(1./wx)- &
2335               (sqrt(wz)*atan(1./sqrt(wz)))+                            &
2336               (1./4.)*(log((1.+wx*wx)*(1.+wy*wy)/(1.+wz))+             &
2337               wy*wy*log(wy*wy*(1.+wz)/(wz*(1.+wy*wy)))+                &
2338               wx*wx*log(wx*wx*(1.+wz)/(wz*(1.+wx*wx)))))
2339
2340         return
2341         end subroutine fnrm_ints
2342
2343 !====6=8==============================================================72
2344 !====6=8==============================================================72
2345 END MODULE module_sf_bem