PageRenderTime 93ms CodeModel.GetById 27ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_sf_bep_bem.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 4101 lines | 2560 code | 730 blank | 811 comment | 59 complexity | 5fab623c5deab0dceb0cde135a58a9af MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. MODULE module_sf_bep_bem
    2. 

  2. 

  3. !USE module_model_constants
    4. 

  4.  USE module_sf_urban
    5. 

  5.  USE module_sf_bem
    6. 

  6. 

  7. ! SGClarke 09/11/2008
    8. 

  8. ! Access urban_param.tbl values through calling urban_param_init in module_physics_init
    9. 

  9. ! for CASE (BEPSCHEME) select sf_urban_physics
    10. 

  10. !
    11. 

  11. ! -----------------------------------------------------------------------
    12. 

  12. !  Dimension for the array used in the BEP module
    13. 

  13. ! -----------------------------------------------------------------------
    14. 

  14. 

  15.       integer nurbm           ! Maximum number of urban classes    
    16. 

  16.       parameter (nurbm=3)
    17. 

  17. 

  18.       integer ndm             ! Maximum number of street directions 
    19. 

  19.       parameter (ndm=2)
    20. 

  20. 

  21.       integer nz_um           ! Maximum number of vertical levels in the urban grid
    22. 

  22.       parameter(nz_um=13)
    23. 

  23. 

  24.       integer ng_u            ! Number of grid levels in the ground
    25. 

  25.       parameter (ng_u=10)
    26. 

  26.       integer nwr_u            ! Number of grid levels in the walls or roofs
    27. 

  27.       parameter (nwr_u=10)
    28. 

  28. 

  29.       integer nf_u             !Number of grid levels in the floors (BEM)
    30. 

  30.       parameter (nf_u=10)
    31. 

  31. 

  32.       integer ngb_u            !Number of grid levels in the ground below building (BEM)
    33. 

  33.       parameter (ngb_u=10)
    34. 

  34. 

  35.       real dz_u                ! Urban grid resolution
    36. 

  36.       parameter (dz_u=5.) 
    37. 

  37. 

  38.       integer nbui_max         !maximum number of types of buildings in an urban class
    39. 

  39.       parameter (nbui_max=4)   !must be less or equal than nz_um 
    40. 

  40. 

  41. !---------------------------------------------------------------------------------
    42. 

  42. !Parameters of the windows. The glasses of windows are considered without films  -
    43. 

  43. !Read the paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour  -
    44. 

  44. !of the total solar energy transmittance of windows".Solar Energy Vol.69,No.4,   -
    45. 

  45. !pp. 321-329, for more details.                                                  -
    46. 

  46. !---------------------------------------------------------------------------------
    47. 

  47.       integer p_num            !number of panes in the windows (1,2 or 3)
    48. 

  48.       parameter (p_num=2)
    49. 

  49.       integer q_num            !category number for the windows (q_num= 4, standard glasses)
    50. 

  50.       parameter(q_num=4)       !Possible values 1,2,...,10
    51. 

  51. 

  52. 

  53. ! The change of ng_u, nwr_u should be done in agreement with the block data
    54. 

  54. !  in the routine "surf_temp" 
    55. 

  55. ! -----------------------------------------------------------------------
    56. 

  56. !  Constant used in the BEP module
    57. 

  57. ! -----------------------------------------------------------------------
    58. 

  58.            
    59. 

  59.       real vk                 ! von Karman constant
    60. 

  60.       real g_u                ! Gravity acceleration
    61. 

  61.       real pi                 !
    62. 

  62.       real r                  ! Perfect gas constant
    63. 

  63.       real cp_u               ! Specific heat at constant pressure
    64. 

  64.       real rcp_u              !
    65. 

  65.       real sigma              !
    66. 

  66.       real p0                 ! Reference pressure at the sea level
    67. 

  67.       real cdrag              ! Drag force constant
    68. 

  68.       real latent             ! Latent heat of vaporization [J/kg] (used in BEM)
    69. 

  69.       parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)        
    70. 

  70.       parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4,latent=2.45e+06)
    71. 

  71. 

  72. ! -----------------------------------------------------------------------     
    73. 

  73. 

  74. 

  75. 

  76. 

  77.    CONTAINS
    78. 

  78.  
    79. 

  79.       subroutine BEP_BEM(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy,      &
    80. 

  80.                       th_phy,rho,p_phy,swdown,glw,                                 &
    81. 

  81.                       gmt,julday,xlong,xlat,                                       &
    82. 

  82.                       declin_urb,cosz_urb2d,omg_urb2d,                             &
    83. 

  83.                       num_urban_layers,                                            &
    84. 

  84.                       trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,                     &
    85. 

  85.                       tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d,             &
    86. 

  86.                       tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d,             &
    87. 

  87.                       cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d,                       &
    88. 

  88.                       sfwin1_urb3d,sfwin2_urb3d,                                   &
    89. 

  89.                       sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,                   &
    90. 

  90.                       a_u,a_v,a_t,a_e,b_u,b_v,                                     &
    91. 

  91.                       b_t,b_e,b_q,dlg,dl_u,sf,vl,                                  &
    92. 

  92.                       rl_up,rs_abs,emiss,grdflx_urb,qv_phy,                        &
    93. 

  93.                       ids,ide, jds,jde, kds,kde,                                   &
    94. 

  94.                       ims,ime, jms,jme, kms,kme,                                   &
    95. 

  95.                       its,ite, jts,jte, kts,kte)                    
    96. 

  96. 

  97.       implicit none
    98. 

  98. 

  99. !------------------------------------------------------------------------
    100. 

  100. !     Input
    101. 

  101. !------------------------------------------------------------------------
    102. 

  102.    INTEGER ::                       ids,ide, jds,jde, kds,kde,  &
    103. 

  103.                                     ims,ime, jms,jme, kms,kme,  &
    104. 

  104.                                     its,ite, jts,jte, kts,kte,  &
    105. 

  105.                                     itimestep
    106. 

  106.  
    107. 

  107. 

  108.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   DZ8W
    109. 

  109.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   P_PHY
    110. 

  110.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   RHO
    111. 

  111.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   TH_PHY
    112. 

  112.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   T_PHY
    113. 

  113.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   U_PHY
    114. 

  114.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   V_PHY
    115. 

  115.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   U
    116. 

  116.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   V
    117. 

  117.    REAL, DIMENSION( ims:ime , jms:jme )        ::   GLW
    118. 

  118.    REAL, DIMENSION( ims:ime , jms:jme )        ::   swdown
    119. 

  119.    REAL, DIMENSION( ims:ime, jms:jme )         ::   UST
    120. 

  120.    INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   UTYPE_URB2D
    121. 

  121.    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   FRC_URB2D
    122. 

  122.    REAL, INTENT(IN  )   ::                                   GMT 
    123. 

  123.    INTEGER, INTENT(IN  ) ::                               JULDAY
    124. 

  124.    REAL, DIMENSION( ims:ime, jms:jme ),                           &
    125. 

  125.          INTENT(IN   )  ::                           XLAT, XLONG
    126. 

  126.    REAL, INTENT(IN) :: DECLIN_URB
    127. 

  127.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
    128. 

  128.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
    129. 

  129.    INTEGER, INTENT(IN  ) :: num_urban_layers
    130. 

  130.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
    131. 

  131.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
    132. 

  132.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
    133. 

  133.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
    134. 

  134. !New variables used for BEM
    135. 

  135.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: qv_phy
    136. 

  136.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tlev_urb3d
    137. 

  137.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: qlev_urb3d
    138. 

  138.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
    139. 

  139.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
    140. 

  140.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tglev_urb3d
    141. 

  141.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tflev_urb3d
    142. 

  142.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
    143. 

  143.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
    144. 

  144.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
    145. 

  145.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
    146. 

  146.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
    147. 

  147.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
    148. 

  148.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
    149. 

  149. !End variables
    150. 

  150.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
    151. 

  151.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
    152. 

  152.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
    153. 

  153.    REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
    154. 

  154. 

  155.    real z(ims:ime,kms:kme,jms:jme)            ! Vertical coordinates
    156. 

  156.    REAL, INTENT(IN )::   DT      ! Time step
    157. 

  157. 

  158. !      
    159. 

  159. !------------------------------------------------------------------------
    160. 

  160. !     Output
    161. 

  161. !------------------------------------------------------------------------ 
    162. 

  162. !
    163. 

  163. !    Implicit and explicit components of the source and sink terms at each levels,
    164. 

  164. !     the fluxes can be computed as follow: FX = A*X + B   example: t_fluxes = a_t * pt + b_t
    165. 

  165.       real a_u(ims:ime,kms:kme,jms:jme)         ! Implicit component for the momemtum in X-direction (center)
    166. 

  166.       real a_v(ims:ime,kms:kme,jms:jme)         ! Implicit component for the momemtum in Y-direction (center)
    167. 

  167.       real a_t(ims:ime,kms:kme,jms:jme)         ! Implicit component for the temperature
    168. 

  168.       real a_e(ims:ime,kms:kme,jms:jme)         ! Implicit component for the TKE
    169. 

  169.       real b_u(ims:ime,kms:kme,jms:jme)         ! Explicit component for the momemtum in X-direction (center)
    170. 

  170.       real b_v(ims:ime,kms:kme,jms:jme)         ! Explicit component for the momemtum in Y-direction (center)
    171. 

  171.       real b_t(ims:ime,kms:kme,jms:jme)         ! Explicit component for the temperature
    172. 

  172.       real b_e(ims:ime,kms:kme,jms:jme)         ! Explicit component for the TKE
    173. 

  173.       real b_q(ims:ime,kms:kme,jms:jme)         ! Explicit component for the Humidity
    174. 

  174.       real dlg(ims:ime,kms:kme,jms:jme)         ! Height above ground (L_ground in formula (24) of the BLM paper). 
    175. 

  175.       real dl_u(ims:ime,kms:kme,jms:jme)        ! Length scale (lb in formula (22) ofthe BLM paper).
    176. 

  176. ! urban surface and volumes        
    177. 

  177.       real sf(ims:ime,kms:kme,jms:jme)           ! surface of the urban grid cells
    178. 

  178.       real vl(ims:ime,kms:kme,jms:jme)             ! volume of the urban grid cells
    179. 

  179. ! urban fluxes
    180. 

  180.       real rl_up(ims:ime,jms:jme) ! upward long wave radiation
    181. 

  181.       real rs_abs(ims:ime,jms:jme) ! absorbed short wave radiation
    182. 

  182.       real emiss(ims:ime,jms:jme)  ! emissivity averaged for urban surfaces
    183. 

  183.       real grdflx_urb(ims:ime,jms:jme)  ! ground heat flux for urban areas
    184. 

  184. !------------------------------------------------------------------------
    185. 

  185. !     Local
    186. 

  186. !------------------------------------------------------------------------
    187. 

  187. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    188. 

  188. !    Building parameters      
    189. 

  189.       real alag_u(nurbm)                      ! Ground thermal diffusivity [m^2 s^-1]
    190. 

  190.       real alaw_u(nurbm)                      ! Wall thermal diffusivity [m^2 s^-1]
    191. 

  191.       real alar_u(nurbm)                      ! Roof thermal diffusivity [m^2 s^-1]
    192. 

  192.       real csg_u(nurbm)                       ! Specific heat of the ground material [J m^3 K^-1]
    193. 

  193.       real csw_u(nurbm)                       ! Specific heat of the wall material [J m^3 K^-1]
    194. 

  194.       real csr_u(nurbm)                       ! Specific heat of the roof material [J m^3 K^-1]
    195. 

  195.       real twini_u(nurbm)                     ! Initial temperature inside the building's wall [K]
    196. 

  196.       real trini_u(nurbm)                     ! Initial temperature inside the building's roof [K]
    197. 

  197.       real tgini_u(nurbm)                     ! Initial road temperature
    198. 

  198. !
    199. 

  199. ! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
    200. 

  200. !
    201. 

  201. !    Radiation paramters
    202. 

  202.       real albg_u(nurbm)                      ! Albedo of the ground
    203. 

  203.       real albw_u(nurbm)                      ! Albedo of the wall
    204. 

  204.       real albr_u(nurbm)                      ! Albedo of the roof
    205. 

  205.       real albwin_u(nurbm)                    ! Albedo of the windows
    206. 

  206.       real emwind_u(nurbm)                    ! Emissivity of windows
    207. 

  207.       real emg_u(nurbm)                       ! Emissivity of ground
    208. 

  208.       real emw_u(nurbm)                       ! Emissivity of wall
    209. 

  209.       real emr_u(nurbm)                       ! Emissivity of roof
    210. 

  210. 

  211. !   fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
    212. 

  212. !   and the short wave radation. 
    213. 

  213.       real fww(nz_um,nz_um,ndm,nurbm)         !  from wall to wall
    214. 

  214.       real fwg(nz_um,ndm,nurbm)               !  from wall to ground
    215. 

  215.       real fgw(nz_um,ndm,nurbm)               !  from ground to wall
    216. 

  216.       real fsw(nz_um,ndm,nurbm)               !  from sky to wall
    217. 

  217.       real fws(nz_um,ndm,nurbm)               !  from sky to wall
    218. 

  218.       real fsg(ndm,nurbm)                     !  from sky to ground
    219. 

  219. 

  220. !    Roughness parameters
    221. 

  221.       real z0g_u(nurbm)         ! The ground's roughness length
    222. 

  222.       real z0r_u(nurbm)         ! The roof's roughness length
    223. 

    224. 

  224. !    Street parameters
    225. 

  225.       integer nd_u(nurbm)       ! Number of street direction for each urban class 
    226. 

  226.       real strd_u(ndm,nurbm)    ! Street length (fix to greater value to the horizontal length of the cells)
    227. 

  227.       real drst_u(ndm,nurbm)    ! Street direction
    228. 

  228.       real ws_u(ndm,nurbm)      ! Street width
    229. 

  229.       real bs_u(ndm,nurbm)      ! Building width
    230. 

  230.       real h_b(nz_um,nurbm)     ! Bulding's heights
    231. 

  231.       real d_b(nz_um,nurbm)     ! Probability that a building has an height h_b
    232. 

  232.       real ss_u(nz_um,nurbm)  ! Probability that a building has an height equal to z
    233. 

  233.       real pb_u(nz_um,nurbm)  ! Probability that a building has an height greater or equal to z
    234. 

  234. 

  235. 

  236. !    Grid parameters
    237. 

  237.       integer nz_u(nurbm)       ! Number of layer in the urban grid
    238. 

  238.       
    239. 

  239.       real z_u(nz_um)         ! Height of the urban grid levels
    240. 

  240. 

  241. !   MT
    242. 

  242.       real cop_u(nurbm)
    243. 

  243.       real pwin_u(nurbm)
    244. 

  244.       real beta_u(nurbm)
    245. 

  245.       integer sw_cond_u(nurbm)
    246. 

  246.       real time_on_u(nurbm)
    247. 

  247.       real time_off_u(nurbm)
    248. 

  248.       real targtemp_u(nurbm)
    249. 

  249.       real gaptemp_u(nurbm)
    250. 

  250.       real targhum_u(nurbm)
    251. 

  251.       real gaphum_u(nurbm)
    252. 

  252.       real perflo_u(nurbm)
    253. 

  253.       real hsesf_u(nurbm)
    254. 

  254.       real hsequip(24)
    255. 

  255. 

  256. ! 1D array used for the input and output of the routine "urban"
    257. 

  257. 

  258.       real z1D(kms:kme)               ! vertical coordinates
    259. 

  259.       real ua1D(kms:kme)                ! wind speed in the x directions
    260. 

  260.       real va1D(kms:kme)                ! wind speed in the y directions
    261. 

  261.       real pt1D(kms:kme)                ! potential temperature
    262. 

  262.       real da1D(kms:kme)                ! air density
    263. 

  263.       real pr1D(kms:kme)                ! air pressure
    264. 

  264.       real pt01D(kms:kme)               ! reference potential temperature
    265. 

  265.       real zr1D                    ! zenith angle
    266. 

  266.       real deltar1D                ! declination of the sun
    267. 

  267.       real ah1D                    ! hour angle (it should come from the radiation routine)
    268. 

  268.       real rs1D                    ! solar radiation
    269. 

  269.       real rld1D                   ! downward flux of the longwave radiation
    270. 

  270. 

  271.       real tw1D(2*ndm,nz_um,nwr_u,nbui_max) ! temperature in each layer of the wall
    272. 

  272.       real tg1D(ndm,ng_u)                   ! temperature in each layer of the ground
    273. 

  273.       real tr1D(ndm,nz_um,nwr_u)   ! temperature in each layer of the roof
    274. 

  274. !
    275. 

  275. !New variable for BEM
    276. 

  276. !
    277. 

  277.       real tlev1D(nz_um,nbui_max)            ! temperature in each floor and in each different type of building
    278. 

  278.       real qlev1D(nz_um,nbui_max)            ! specific humidity in each floor and in each different type of building
    279. 

  279.       real twlev1D(2*ndm,nz_um,nbui_max)     ! temperature in each window in each floor in each different type of building
    280. 

  280.       real tglev1D(ndm,ngb_u,nbui_max)       ! temperature in each layer of the ground below of a type of building
    281. 

  281.       real tflev1D(ndm,nf_u,nz_um-1,nbui_max)! temperature in each layer of the floors in each building
    282. 

  282.       real lflev1D(nz_um,nz_um)           ! latent heat flux due to the air conditioning systems
    283. 

  283.       real sflev1D(nz_um,nz_um)           ! sensible heat flux due to the air conditioning systems
    284. 

  284.       real lfvlev1D(nz_um,nz_um)          ! latent heat flux due to ventilation
    285. 

  285.       real sfvlev1D(nz_um,nz_um)          ! sensible heat flux due to ventilation
    286. 

  286.       real sfwin1D(2*ndm,nz_um,nbui_max)     ! sensible heat flux from windows
    287. 

  287.       real consumlev1D(nz_um,nz_um)       ! consumption due to the air conditioning systems
    288. 

  288.       real qv1D(kms:kme)                  ! specific humidity
    289. 

  289.       real meso_urb                       ! constant to link meso and urban scales [m¯2]
    290. 

  290.       real d_urb(nz_um)    
    291. 

  291.       real sf_ac
    292. 

  292.       integer ibui,nbui
    293. 

  293.       integer nlev(nz_um)
    294. 

  294. !
    295. 

  295. !End new variables
    296. 

  296. !
    297. 

  297. 

  298.       real sfw1D(2*ndm,nz_um,nbui_max)      ! sensible heat flux from walls
    299. 

  299.       real sfg1D(ndm)              ! sensible heat flux from ground (road)
    300. 

  300.       real sfr1D(ndm,nz_um)      ! sensible heat flux from roofs
    301. 

  301.       real sf1D(kms:kme)              ! surface of the urban grid cells
    302. 

  302.       real vl1D(kms:kme)                ! volume of the urban grid cells
    303. 

  303.       real a_u1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the X-direction
    304. 

  304.       real a_v1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the Y-direction
    305. 

  305.       real a_t1D(kms:kme)               ! Implicit component of the heat sources or sinks
    306. 

  306.       real a_e1D(kms:kme)               ! Implicit component of the TKE sources or sinks
    307. 

  307.       real b_u1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the X-direction
    308. 

  308.       real b_v1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the Y-direction
    309. 

  309.       real b_t1D(kms:kme)               ! Explicit component of the heat sources or sinks
    310. 

  310.       real b_ac1D(kms:kme)
    311. 

  311.       real b_e1D(kms:kme)               ! Explicit component of the TKE sources or sinks
    312. 

  312.       real b_q1D(kms:kme)               ! Explicit component of the Humidity sources or sinks
    313. 

  313.       real dlg1D(kms:kme)               ! Height above ground (L_ground in formula (24) of the BLM paper). 
    314. 

  314.       real dl_u1D(kms:kme)              ! Length scale (lb in formula (22) ofthe BLM paper)
    315. 

  315.       real time_bep
    316. 

  316. ! arrays used to collapse indexes
    317. 

  317.       integer ind_zwd(nbui_max,nz_um,nwr_u,ndm)
    318. 

  318.       integer ind_gd(ng_u,ndm)
    319. 

  319.       integer ind_zd(nbui_max,nz_um,ndm)
    320. 

  320.       integer ind_zdf(nz_um,ndm)
    321. 

  321.       integer ind_zrd(nz_um,nwr_u,ndm)
    322. 

  322. !
    323. 

  323.       integer ind_bd(nbui_max,nz_um)
    324. 

  324.       integer ind_wd(nbui_max,nz_um,ndm)
    325. 

  325.       integer ind_gbd(nbui_max,ngb_u,ndm)  
    326. 

  326.       integer ind_fbd(nbui_max,nf_u,nz_um-1,ndm)
    327. 

  327. !
    328. 

  328.       integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
    329. 

  329.       integer it, nint
    330. 

  330.       integer iii
    331. 

  331.       real tempo
    332. 

  332.       logical first
    333. 

  333.       character(len=80) :: text
    334. 

  334.       data first/.true./
    335. 

  335.       save first,time_bep 
    336. 

  336. 

  337.       save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,                      &
    338. 

  338.           albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,fww,fwg,fgw,fsw,fws,fsg, &
    339. 

  339.           z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  &
    340. 

  340.           nz_u,z_u,albwin_u,emwind_u ,                                  &
    341. 

  341.           cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, targtemp_u,  &
    342. 

  342.           gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip
    343. 

  343. 

  344. !------------------------------------------------------------------------
    345. 

  345. !    Calculation of the momentum, heat and turbulent kinetic fluxes
    346. 

  346. !     produced by buildings
    347. 

  347. !
    348. 

  348. ! Reference:
    349. 

  349. ! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
    350. 

  350. ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
    351. 

  351. ! 261-304
    352. 

  352. !
    353. 

  353. ! F. Salamanca and A. Martilli, 2009: 'A new Building Energy Model coupled 
    354. 

  354. ! with an Urban Canopy Parameterization for urban climate simulations_part II. 
    355. 

  355. ! Validation with one dimension off-line simulations'. Theor Appl Climatol
    356. 

  356. ! DOI 10.1007/s00704-009-0143-8 
    357. 

  357. !------------------------------------------------------------------------
    358. 

  358. !prepare the arrays to collapse indexes
    359. 

  359.       if(num_urban_layers.lt.nbui_max*nz_um*ndm*max(nwr_u,ng_u))then
    360. 

  360.         write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um*ndm*max(nwr_u,ng_u)
    361. 

  361.         stop
    362. 

  362.       endif
    363. 

  363. !
    364. 

  364. !New conditions for BEM
    365. 

  365. !
    366. 

  366.       if(num_urban_layers.lt.nbui_max*nz_um)then !limit for indoor temperature and indoor humidity
    367. 

  367.         write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um
    368. 

  368.         stop
    369. 

  369.       endif
    370. 

  370. 

  371.       if(num_urban_layers.lt.nbui_max*nz_um*ndm)then !limit for window temperature
    372. 

  372.         write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um*ndm
    373. 

  373.         stop
    374. 

  374.       endif
    375. 

  375. 

  376.       if(num_urban_layers.lt.nbui_max*ndm*ngb_u)then !limit for ground temperature below a building
    377. 

  377.         write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*ndm*ngb_u
    378. 

  378.         stop
    379. 

  379.       endif
    380. 

  380. 

  381.       if(num_urban_layers.lt.(nz_um-1)*nbui_max*ndm*nf_u)then !limit for floor temperature 
    382. 

  382.         write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*ndm*nf_u*(nz_um-1),num_urban_layers
    383. 

  383.         stop
    384. 

  384.       endif
    385. 

  385. 

  386.       if (ndm.ne.2)then
    387. 

  387.          write(*,*) 'number of directions is not correct',ndm
    388. 

  388.          stop
    389. 

  389.       endif
    390. 

  390. 

  391. !
    392. 

  392. !End of new conditions
    393. 

  393. !
    394. 

  394. !
    395. 

  395. !Initialize collapse indexes
    396. 

  396. !
    397. 

  397.       ind_zwd=0       
    398. 

  398.       ind_gd=0
    399. 

  399.       ind_zd=0
    400. 

  400.       ind_zdf=0
    401. 

  401.       ind_zrd=0
    402. 

  402.       ind_bd=0
    403. 

  403.       ind_wd=0
    404. 

  404.       ind_gbd=0
    405. 

  405.       ind_fbd=0
    406. 

  406. !
    407. 

  407. !End initialization indexes
    408. 

  408. !
    409. 

  409. 

  410.       iii=0
    411. 

  411.       do ibui=1,nbui_max
    412. 

  412.       do iz_u=1,nz_um
    413. 

  413.       do iw=1,nwr_u
    414. 

  414.       do id=1,ndm
    415. 

  415.        iii=iii+1
    416. 

  416.        ind_zwd(ibui,iz_u,iw,id)=iii
    417. 

  417.       enddo
    418. 

  418.       enddo
    419. 

  419.       enddo
    420. 

  420.       enddo
    421. 

  421.       
    422. 

  422.       iii=0
    423. 

  423.       do ig=1,ng_u
    424. 

  424.       do id=1,ndm
    425. 

  425.        iii=iii+1
    426. 

  426.        ind_gd(ig,id)=iii
    427. 

  427.       enddo
    428. 

  428.       enddo
    429. 

  429.       
    430. 

  430.       iii=0
    431. 

  431.       do ibui=1,nbui_max
    432. 

  432.       do iz_u=1,nz_um
    433. 

  433.       do id=1,ndm
    434. 

  434.        iii=iii+1
    435. 

  435.        ind_zd(ibui,iz_u,id)=iii
    436. 

  436.       enddo
    437. 

  437.       enddo
    438. 

  438.       enddo
    439. 

  439.       
    440. 

  440.       iii=0
    441. 

  441.       do iz_u=1,nz_um
    442. 

  442.       do iw=1,nwr_u
    443. 

  443.       do id=1,ndm
    444. 

  444.        iii=iii+1
    445. 

  445.        ind_zrd(iz_u,iw,id)=iii
    446. 

  446.       enddo
    447. 

  447.       enddo
    448. 

  448.       enddo
    449. 

  449.       
    450. 

  450. !
    451. 

  451. !New indexes for BEM
    452. 

  452. !
    453. 

  453.       iii=0
    454. 

  454.       do iz_u=1,nz_um
    455. 

  455.       do id=1,ndm
    456. 

  456.          iii=iii+1
    457. 

  457.          ind_zdf(iz_u,id)=iii
    458. 

  458.       enddo ! id
    459. 

  459.       enddo ! iz_u
    460. 

  460.       
    461. 

  461.       iii=0
    462. 

  462.       do ibui=1,nbui_max  !Type of building
    463. 

  463.       do iz_u=1,nz_um     !vertical levels
    464. 

  464.          iii=iii+1
    465. 

  465.          ind_bd(ibui,iz_u)=iii
    466. 

  466.       enddo !iz_u
    467. 

  467.       enddo !ibui
    468. 

  468.       
    469. 

  469.       
    470. 

  470.       iii=0
    471. 

  471.       do ibui=1,nbui_max !type of building
    472. 

  472.       do iz_u=1,nz_um !vertical levels
    473. 

  473.       do id=1,ndm !direction
    474. 

  474.          iii=iii+1
    475. 

  475.          ind_wd(ibui,iz_u,id)=iii
    476. 

  476.       enddo !id
    477. 

  477.       enddo !iz_u
    478. 

  478.       enddo !ibui
    479. 

  479.       
    480. 

  480.       iii=0
    481. 

  481.       do ibui=1,nbui_max!type of building
    482. 

  482.       do iw=1,ngb_u !layers in the wall (ground below a building)
    483. 

  483.       do id=1,ndm !direction
    484. 

  484.          iii=iii+1
    485. 

  485.          ind_gbd(ibui,iw,id)=iii  
    486. 

  486.       enddo !id
    487. 

  487.       enddo !iw 
    488. 

  488.       enddo !ibui    
    489. 

  489.       
    490. 

  490.       iii=0
    491. 

  491.       do ibui=1,nbui_max !type of building
    492. 

  492.       do iw=1,nf_u !layers in the wall (floor)
    493. 

  493.       do iz_u=1,nz_um-1 !vertical levels
    494. 

  494.       do id=1,ndm  !direction
    495. 

  495.          iii=iii+1
    496. 

  496.          ind_fbd(ibui,iw,iz_u,id)=iii
    497. 

  497.       enddo !id
    498. 

  498.       enddo !iz_u
    499. 

  499.       enddo !iw
    500. 

  500.       enddo !ibui
    501. 

  501.       
    502. 

  502. !
    503. 

  503. !End of new indexes
    504. 

  504. !      
    505. 

  505. 

  506.       do ix=its,ite
    507. 

  507.       do iy=jts,jte
    508. 

  508.        z(ix,kts,iy)=0.
    509. 

  509.        do iz=kts+1,kte+1
    510. 

  510.         z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
    511. 

  511.        enddo
    512. 

  512.       enddo
    513. 

  513.       enddo
    514. 

  514. 

  515.       if (first) then                           ! True only on first call
    516. 

  516.          call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
    517. 

  517.                 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
    518. 

  518.                 emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,&
    519. 

  519.                 cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
    520. 

  520.                 targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)
    521. 

  521. 

  522. !Initialisation of the urban parameters and calculation of the view factor
    523. 

  523.        call icBEP(fww,fwg,fgw,fsw,fws,fsg,                         & 
    524. 

  524.                  z0g_u,z0r_u,                                      & 
    525. 

  525.                  nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,          & 
    526. 

  526.                  nz_u,z_u)                                  
    527. 

  527. 

  528.        first=.false.
    529. 

  529. 

  530.       endif ! first
    531. 

  531.       
    532. 

  532.       do ix=its,ite
    533. 

  533.       do iy=jts,jte
    534. 

  534.         if (FRC_URB2D(ix,iy).gt.0.) then    ! Calling BEP only for existing urban classes.
    535. 

  535. 	
    536. 

  536.          iurb=UTYPE_URB2D(ix,iy)
    537. 

  537.          ibui=0
    538. 

  538.          nlev=0
    539. 

  539.          nbui=0
    540. 

  540.          d_urb=0.
    541. 

  541.          do iz=1,nz_um		   
    542. 

  542.          if(ss_u(iz,iurb).gt.0) then		
    543. 

  543.            ibui=ibui+1		                
    544. 

  544.            nlev(ibui)=iz-1
    545. 

  545.            d_urb(ibui)=ss_u(iz,iurb)
    546. 

  546.            nbui=ibui
    547. 

  547. 	 endif	  
    548. 

  548.          end do  !iz
    549. 

  549.          if (nbui.gt.nbui_max) then
    550. 

  550.             write (*,*) 'nbui_max must be increased to',nbui
    551. 

  551.             stop
    552. 

  552.          endif
    553. 

  553. 

  554.        do iz= kts,kte
    555. 

  555.           ua1D(iz)=u_phy(ix,iz,iy)
    556. 

  556.           va1D(iz)=v_phy(ix,iz,iy)
    557. 

  557. 	  pt1D(iz)=th_phy(ix,iz,iy)
    558. 

  558. 	  da1D(iz)=rho(ix,iz,iy)
    559. 

  559. 	  pr1D(iz)=p_phy(ix,iz,iy)
    560. 

  560. !	  pt01D(iz)=th_phy(ix,iz,iy)
    561. 

  561. 	  pt01D(iz)=300.
    562. 

  562. 	  z1D(iz)=z(ix,iz,iy)
    563. 

  563.           qv1D(iz)=qv_phy(ix,iz,iy)
    564. 

  564.           a_u1D(iz)=0.
    565. 

  565.           a_v1D(iz)=0.
    566. 

  566.           a_t1D(iz)=0.
    567. 

  567.           a_e1D(iz)=0.
    568. 

  568.           b_u1D(iz)=0.
    569. 

  569.           b_v1D(iz)=0.
    570. 

  570.           b_t1D(iz)=0.
    571. 

  571.           b_ac1D(iz)=0.
    572. 

  572.           b_e1D(iz)=0.           
    573. 

  573.          enddo
    574. 

  574. 	 z1D(kte+1)=z(ix,kte+1,iy)
    575. 

  575. 

  576.          do id=1,ndm
    577. 

  577.          do iz_u=1,nz_um
    578. 

  578.          do iw=1,nwr_u
    579. 

  579.          do ibui=1,nbui_max
    580. 

  580.           tw1D(2*id-1,iz_u,iw,ibui)=tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
    581. 

  581.           tw1D(2*id,iz_u,iw,ibui)=tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
    582. 

  582.          enddo
    583. 

  583.          enddo
    584. 

  584.          enddo
    585. 

  585.          enddo
    586. 

  586. 	
    587. 

  587.          do id=1,ndm
    588. 

  588.           do ig=1,ng_u
    589. 

  589.             tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
    590. 

  590.           enddo
    591. 

  591.           do iz_u=1,nz_um
    592. 

  592.           do ir=1,nwr_u
    593. 

  593.             tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)
    594. 

  594.           enddo
    595. 

  595.           enddo
    596. 

  596.          enddo
    597. 

  597. !
    598. 

  598. !Initialize variables for BEM
    599. 

  599. !
    600. 

  600.          tlev1D=0.  !Indoor temperature
    601. 

  601.          qlev1D=0.  !Indoor humidity
    602. 

  602. 

  603.          twlev1D=0. !Window temperature
    604. 

  604.          tglev1D=0. !Ground temperature
    605. 

  605.          tflev1D=0. !Floor temperature
    606. 

  606. 

  607.          sflev1D=0.    !Sensible heat flux from the a.c.
    608. 

  608.          lflev1D=0.    !latent heat flux from the a.c.
    609. 

  609.          consumlev1D=0.!consumption of the a.c.
    610. 

  610.          sfvlev1D=0.   !Sensible heat flux from natural ventilation
    611. 

  611.          lfvlev1D=0.   !Latent heat flux from natural ventilation
    612. 

  612.          sfwin1D=0.    !Sensible heat flux from windows
    613. 

  613.          sfw1D=0.      !Sensible heat flux from walls         
    614. 

  614. 

  615.          do iz_u=1,nz_um    !vertical levels
    616. 

  616.          do ibui=1,nbui_max !Type of building
    617. 

  617.             tlev1D(iz_u,ibui)= tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)  
    618. 

  618.             qlev1D(iz_u,ibui)= qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)  
    619. 

  619.          enddo !ibui
    620. 

  620.          enddo !iz_u
    621. 

  621. 

  622.          do id=1,ndm  !direction
    623. 

  623.             do iz_u=1,nz_um !vertical levels
    624. 

  624.                do ibui=1,nbui_max !type of building
    625. 

  625.                   twlev1D(2*id-1,iz_u,ibui)=tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
    626. 

  626.                   twlev1D(2*id,iz_u,ibui)=tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
    627. 

  627.                   sfwin1D(2*id-1,iz_u,ibui)=sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
    628. 

  628.                   sfwin1D(2*id,iz_u,ibui)=sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
    629. 

  629.                enddo !ibui  
    630. 

  630.             enddo !iz_u
    631. 

  631.          enddo !id
    632. 

  632. 

  633.          do id=1,ndm !direction
    634. 

  634.             do iw=1,ngb_u !layer in the wall
    635. 

  635.                do ibui=1,nbui_max !type of building
    636. 

  636.                   tglev1D(id,iw,ibui)=tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)
    637. 

  637.                enddo !ibui
    638. 

  638.             enddo !iw
    639. 

  639.          enddo !id
    640. 

  640.        
    641. 

  641.          do id=1,ndm !direction
    642. 

  642.             do iw=1,nf_u !layer in the walls
    643. 

  643.                do iz_u=1,nz_um-1 !verticals levels
    644. 

  644.                   do ibui=1,nbui_max !type of building
    645. 

  645.                      tflev1D(id,iw,iz_u,ibui)=tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)
    646. 

  646.                   enddo !ibui
    647. 

  647.                enddo ! iz_u
    648. 

  648.              enddo !iw
    649. 

  649.          enddo !id
    650. 

  650. 

  651. !
    652. 

  652. !End initialization for BEM
    653. 

  653. !         
    654. 

  654. 

  655.          do id=1,ndm
    656. 

  656. 	 do iz=1,nz_um
    657. 

  657.          do ibui=1,nbui_max !type of building
    658. 

  658. 	  sfw1D(2*id-1,iz,ibui)=sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)
    659. 

  659. 	  sfw1D(2*id,iz,ibui)=sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)
    660. 

  660.          enddo
    661. 

  661. 	 enddo
    662. 

  662. 	 enddo
    663. 

  663. 	 
    664. 

  664. 	 do id=1,ndm
    665. 

  665. 	  sfg1D(id)=sfg_urb3d(ix,id,iy)
    666. 

  666. 	 enddo
    667. 

  667. 	 
    668. 

  668. 	 do id=1,ndm
    669. 

  669. 	 do iz=1,nz_um
    670. 

  670. 	  sfr1D(id,iz)=sfr_urb3d(ix,ind_zdf(iz,id),iy)
    671. 

  671. 	 enddo
    672. 

  672. 	 enddo
    673. 

  673.          
    674. 

  674.          rs1D=swdown(ix,iy)
    675. 

  675.          rld1D=glw(ix,iy)
    676. 

  676. 

  677.          zr1D=acos(COSZ_URB2D(ix,iy))
    678. 

  678.          deltar1D=DECLIN_URB
    679. 

  679.          ah1D=OMG_URB2D(ix,iy)       
    680. 

  680. 

  681.          call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D,  &
    682. 

  682.                    zr1D,deltar1D,ah1D,rs1D,rld1D,                   & 
    683. 

  683.                    alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,          & 
    684. 

  684.                    albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, & 
    685. 

  685.                    emwind_u,fww,fwg,fgw,fsw,fws,fsg,                & 
    686. 

  686.                    z0g_u,z0r_u,                                     & 
    687. 

  687.                    nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  & 
    688. 

  688.                    nz_u,z_u,                                        & 
    689. 

  689.                    cop_u,pwin_u,beta_u,sw_cond_u,time_on_u,         &
    690. 

  690.                    time_off_u,targtemp_u,gaptemp_u,targhum_u,       &
    691. 

  691.                    gaphum_u, perflo_u, hsesf_u, hsequip,            &
    692. 

  692.                    tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D,                & 
    693. 

  693.                    a_u1D,a_v1D,a_t1D,a_e1D,                         & 
    694. 

  694.                    b_u1D,b_v1D,b_t1D,b_ac1D,b_e1D,b_q1D,            & 
    695. 

  695.                    dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy),             &
    696. 

  696.                    rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy),    &
    697. 

  697.                    qv1D,tlev1D,qlev1D,sflev1D,lflev1D,consumlev1D,  &
    698. 

  698.                    sfvlev1D,lfvlev1D,twlev1D,tglev1D,tflev1D,sfwin1D,&
    699. 

  699.                    ix,iy)                            
    700. 

  700. 

  701.          do id=1,ndm ! direction 
    702. 

  702. 	    do iz=1,nz_um   !vertical levels
    703. 

  703.                do ibui=1,nbui_max !type of building
    704. 

  704. 	          sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id-1,iz,ibui) 
    705. 

  705. 	          sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id,iz,ibui) 
    706. 

  706. 	       enddo
    707. 

  707. 	    enddo
    708. 

  708.          enddo
    709. 

  709. 	 do id=1,ndm
    710. 

  710. 	  sfg_urb3d(ix,id,iy)=sfg1D(id) 
    711. 

  711. 	 enddo
    712. 

  712.          
    713. 

  713. 	 do id=1,ndm
    714. 

  714. 	 do iz=1,nz_um
    715. 

  715. 	  sfr_urb3d(ix,ind_zdf(iz,id),iy)=sfr1D(id,iz)
    716. 

  716. 	 enddo
    717. 

  717. 	 enddo
    718. 

  718.          
    719. 

  719.          do id=1,ndm
    720. 

  720.          do iz_u=1,nz_um
    721. 

  721.          do iw=1,nwr_u
    722. 

  722.          do ibui=1,nbui_max
    723. 

  723.           tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw,ibui)
    724. 

  724.           tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw,ibui)
    725. 

  725.          enddo
    726. 

  726.          enddo
    727. 

  727.          enddo
    728. 

  728.          enddo
    729. 

  729.          
    730. 

  730.          do id=1,ndm
    731. 

  731.             do ig=1,ng_u
    732. 

  732.                tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
    733. 

  733.             enddo
    734. 

  734.             do iz_u=1,nz_um
    735. 

  735.                do ir=1,nwr_u
    736. 

  736.                   trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
    737. 

  737.                enddo
    738. 

  738.             enddo
    739. 

  739.          enddo
    740. 

  740. !
    741. 

  741. !Outputs of BEM
    742. 

  742. !
    743. 

  743.          do ibui=1,nbui_max !type of building
    744. 

  744.          do iz_u=1,nz_um !vertical levels
    745. 

  745.             tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=tlev1D(iz_u,ibui)  
    746. 

  746.             qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=qlev1D(iz_u,ibui)  
    747. 

  747.          enddo !iz_u
    748. 

  748.          enddo !ibui
    749. 

  749.          do id=1,ndm !direction
    750. 

  750.          do iz_u=1,nz_um !vertical levels
    751. 

  751.          do ibui=1,nbui_max !type of building
    752. 

  752.                tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id-1,iz_u,ibui)
    753. 

  753.                tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id,iz_u,ibui)
    754. 

  754.                sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id-1,iz_u,ibui)
    755. 

  755.                sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id,iz_u,ibui)
    756. 

  756.             enddo !ibui  
    757. 

  757.          enddo !iz_u
    758. 

  758.          enddo !id
    759. 

  759.         
    760. 

  760.          do id=1,ndm !direction
    761. 

  761.             do iw=1,ngb_u !layers in the walls
    762. 

  762.                do ibui=1,nbui_max  !type of building
    763. 

  763.                   tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)=tglev1D(id,iw,ibui)
    764. 

  764.                enddo !ibui
    765. 

  765.             enddo !iw
    766. 

  766.          enddo !id
    767. 

  767. 

  768.          do id=1,ndm   !direction 
    769. 

  769.             do iw=1,nf_u !layer in the walls
    770. 

  770.                do iz_u=1,nz_um-1 !verticals levels
    771. 

  771.                   do ibui=1,nbui_max  !type of building
    772. 

  772.                      tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)=tflev1D(id,iw,iz_u,ibui)
    773. 

  773.                   enddo !ibui
    774. 

  774.                enddo !iz_u
    775. 

  775.              enddo !iw
    776. 

  776.          enddo !id
    777. 

  777.          sf_ac_urb3d(ix,iy)=0.
    778. 

  778.          lf_ac_urb3d(ix,iy)=0.
    779. 

  779.          cm_ac_urb3d(ix,iy)=0.
    780. 

  780.          sfvent_urb3d(ix,iy)=0.
    781. 

  781.          lfvent_urb3d(ix,iy)=0.
    782. 

  782. 

  783.          meso_urb=(1./4.)*FRC_URB2D(ix,iy)/((bs_u(1,iurb)+ws_u(1,iurb))*bs_u(2,iurb))+ &
    784. 

  784.                   (1./4.)*FRC_URB2D(ix,iy)/((bs_u(2,iurb)+ws_u(2,iurb))*bs_u(1,iurb))
    785. 

  785. 

  786.                   
    787. 

  787.          ibui=0
    788. 

  788.          nlev=0
    789. 

  789.          nbui=0
    790. 

  790.          d_urb=0.
    791. 

  791.          do iz=1,nz_um		   
    792. 

  792.          if(ss_u(iz,iurb).gt.0) then		
    793. 

  793.            ibui=ibui+1		                
    794. 

  794.            nlev(ibui)=iz-1
    795. 

  795.            d_urb(ibui)=ss_u(iz,iurb)
    796. 

  796.            nbui=ibui
    797. 

  797. 	 endif	  
    798. 

  798.          end do  !iz
    799. 

  799. 

  800.          do ibui=1,nbui       !type of building   
    801. 

  801.          do iz_u=1,nlev(ibui) !vertical levels                    
    802. 

  802.                sf_ac_urb3d(ix,iy)=sf_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*sflev1D(iz_u,ibui)
    803. 

  803.                lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*lflev1D(iz_u,ibui)
    804. 

  804.                cm_ac_urb3d(ix,iy)=cm_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*consumlev1D(iz_u,ibui)
    805. 

  805.                sfvent_urb3d(ix,iy)=sfvent_urb3d(ix,iy)+meso_urb*d_urb(ibui)*sfvlev1D(iz_u,ibui)
    806. 

  806.                lfvent_urb3d(ix,iy)=lfvent_urb3d(ix,iy)+meso_urb*d_urb(ibui)*lfvlev1D(iz_u,ibui)            
    807. 

  807.          enddo !iz_u
    808. 

  808.          enddo !ibui
    809. 

  809. !
    810. 

  810. !Add the latent heat exchanged throughout the ventilation in the lf_ac_urb3d output variable. 
    811. 

  811. !it is only a print variable
    812. 

  812. !
    813. 

  813. !        lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+lfvent_urb3d(ix,iy)
    814. 

  814. !
    815. 

  815. 

  816.          lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)-lfvent_urb3d(ix,iy)        
    817. 

  817. 

  818. !
    819. 

  819. !End outputs of bem
    820. 

  820. !
    821. 

  821. 

  822.         sf_ac=0.
    823. 

  823.         do iz= kts,kte
    824. 

  824.           sf(ix,iz,iy)=sf1D(iz)
    825. 

  825.           vl(ix,iz,iy)=vl1D(iz)
    826. 

  826.           a_u(ix,iz,iy)=a_u1D(iz)
    827. 

  827.           a_v(ix,iz,iy)=a_v1D(iz)
    828. 

  828.           a_t(ix,iz,iy)=a_t1D(iz)
    829. 

  829.           a_e(ix,iz,iy)=a_e1D(iz)
    830. 

  830.           b_u(ix,iz,iy)=b_u1D(iz)
    831. 

  831.           b_v(ix,iz,iy)=b_v1D(iz)
    832. 

  832.           b_t(ix,iz,iy)=b_t1D(iz)
    833. 

  833.           sf_ac=sf_ac+b_ac1D(iz)*da1D(iz)*cp_u*dz8w(ix,iz,iy)*vl1D(iz)*FRC_URB2D(ix,iy)
    834. 

  834.           b_e(ix,iz,iy)=b_e1D(iz)
    835. 

  835.           b_q(ix,iz,iy)=b_q1D(iz)
    836. 

  836.           dlg(ix,iz,iy)=dlg1D(iz)
    837. 

  837.           dl_u(ix,iz,iy)=dl_u1D(iz)
    838. 

  838.          enddo
    839. 

  839.          sf(ix,kte+1,iy)=sf1D(kte+1)
    840. 

  840. 

  841.          endif ! FRC_URB2D
    842. 

  842.    
    843. 

  843.       enddo  ! iy
    844. 

  844.       enddo  ! ix
    845. 

  845. 

  846.         time_bep=time_bep+dt
    847. 

  847. 

  848.       return
    849. 

  849.       end subroutine BEP_BEM
    850. 

  850.             
    851. 

  851. ! ===6=8===============================================================72
    852. 

  852. 

  853.       subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0,   &  
    854. 

  854.                       zr,deltar,ah,rs,rld,                             & 
    855. 

  855.                       alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,          & 
    856. 

  856.                       albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, & 
    857. 

  857.                       emwind_u,fww,fwg,fgw,fsw,fws,fsg,                & 
    858. 

  858.                       z0g_u,z0r_u,                                     & 
    859. 

  859.                       nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  & 
    860. 

  860.                       nz_u,z_u,                                        & 
    861. 

  861.                       cop_u,pwin_u,beta_u,sw_cond_u,time_on_u,         &
    862. 

  862.                       time_off_u,targtemp_u,gaptemp_u,targhum_u,       &
    863. 

  863.                       gaphum_u, perflo_u, hsesf_u, hsequip,            &
    864. 

  864.                       tw,tg,tr,sfw,sfg,sfr,                            & 
    865. 

  865.                       a_u,a_v,a_t,a_e,                                 &
    866. 

  866.                       b_u,b_v,b_t,b_ac,b_e,b_q,                        & 
    867. 

  867.                       dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb,    &
    868. 

  868.                       qv,tlev,qlev,sflev,lflev,consumlev,              &
    869. 

  869.                       sfvlev,lfvlev,twlev,tglev,tflev,sfwin,ix,iy)                             
    870. 

  870. 

  871. ! ----------------------------------------------------------------------
    872. 

  872. ! This routine computes the effects of buildings on momentum, heat and
    873. 

  873. !  TKE (turbulent kinetic energy) sources or sinks and on the mixing length.
    874. 

  874. ! It provides momentum, heat and TKE sources or sinks at different levels of a
    875. 

  875. !  mesoscale grid defined by the altitude of its cell interfaces "z" and
    876. 

  876. !  its number of levels "nz".
    877. 

  877. ! The meteorological input parameters (wind, temperature, solar radiation)
    878. 

  878. !  are specified on the "mesoscale grid".
    879. 

  879. ! The inputs concerning the building and street charateristics are defined
    880. 

  880. !  on a "urban grid". The "urban grid" is defined with its number of levels
    881. 

  881. !  "nz_u" and its space step "dz_u".
    882. 

  882. ! The input parameters are interpolated on the "urban grid". The sources or sinks
    883. 

  883. !  are calculated on the "urban grid". Finally the sources or sinks are 
    884. 

  884. !  interpolated on the "mesoscale grid".
    885. 

  885.  
    886. 

  886. 

  887. !  Mesoscale grid            Urban grid             Mesoscale grid
    888. 

  888. !  
    889. 

  889. ! z(4)    ---                                               ---
    890. 

  890. !          |                                                 |
    891. 

  891. !          |                                                 |
    892. 

  892. !          |   Interpolation                  Interpolation  |
    893. 

  893. !          |            Sources or sinks calculation         |
    894. 

  894. ! z(3)    ---                                               ---
    895. 

  895. !          | ua               ua_u  ---  uv_a         a_u    |
    896. 

  896. !          | va               va_u   |   uv_b         b_u    |
    897. 

  897. !          | pt               pt_u  ---  uh_b         a_v    |
    898. 

  898. ! z(2)    ---                        |    etc...      etc...---
    899. 

  899. !          |                 z_u(1) ---                      |
    900. 

  900. !          |                         |                       |
    901. 

  901. ! z(1) ------------------------------------------------------------
    902. 

  902. 

  903. !     
    904. 

  904. ! Reference:
    905. 

  905. ! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
    906. 

  906. ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
    907. 

  907. ! 261-304
    908. 

  908.  
    909. 

  909. ! ----------------------------------------------------------------------
    910. 

  910. 

  911.       implicit none
    912. 

  912. 

  913. ! ----------------------------------------------------------------------
    914. 

  914. ! INPUT:
    915. 

  915. ! ----------------------------------------------------------------------
    916. 

  916. 

  917. ! Data relative to the "mesoscale grid"
    918. 

  918. 

  919.       integer kms,kme,kts,kte
    920. 

  920.       real z(kms:kme)               ! Altitude above the ground of the cell interfaces.
    921. 

  921.       real ua(kms:kme)                ! Wind speed in the x direction
    922. 

  922.       real va(kms:kme)                ! Wind speed in the y direction
    923. 

  923.       real pt(kms:kme)                ! Potential temperature
    924. 

  924.       real da(kms:kme)                ! Air density
    925. 

  925.       real pr(kms:kme)                ! Air pressure
    926. 

  926.       real pt0(kms:kme)               ! Reference potential temperature (could be equal to "pt")
    927. 

  927.       real qv(kms:kme)              ! Specific humidity
    928. 

  928.       real dt                    ! Time step
    929. 

  929.       real zr                    ! Zenith angle
    930. 

  930.       real deltar                ! Declination of the sun
    931. 

  931.       real ah                    ! Hour angle
    932. 

  932.       real rs                    ! Solar radiation
    933. 

  933.       real rld                   ! Downward flux of the longwave radiation
    934. 

  934. 

  935. ! Data relative to the "urban grid"
    936. 

  936. 

  937.       integer iurb               ! Current urban class
    938. 

  938. 

  939. !    Building parameters      
    940. 

  940.       real alag_u(nurbm)         ! Ground thermal diffusivity [m^2 s^-1]
    941. 

  941.       real alaw_u(nurbm)         ! Wall thermal diffusivity [m^2 s^-1]
    942. 

  942.       real alar_u(nurbm)         ! Roof thermal diffusivity [m^2 s^-1]
    943. 

  943.       real csg_u(nurbm)          ! Specific heat of the ground material [J m^3 K^-1]
    944. 

  944.       real csw_u(nurbm)          ! Specific heat of the wall material [J m^3 K^-1]
    945. 

  945.       real csr_u(nurbm)          ! Specific heat of the roof material [J m^3 K^-1]
    946. 

  946. 

  947. !    Radiation parameters
    948. 

  948.       real albg_u(nurbm)         ! Albedo of the ground
    949. 

  949.       real albw_u(nurbm)         ! Albedo of the wall
    950. 

  950.       real albr_u(nurbm)         ! Albedo of the roof
    951. 

  951.       real albwin_u(nurbm)       ! Albedo of the windows
    952. 

  952.       real emwind_u(nurbm)       ! Emissivity of windows
    953. 

  953.       real emg_u(nurbm)          ! Emissivity of ground
    954. 

  954.       real emw_u(nurbm)          ! Emissivity of wall
    955. 

  955.       real emr_u(nurbm)          ! Emissivity of roof
    956. 

  956. 

  957. !    fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and 
    958. 

  958. !    short wave radation. 
    959. 

  959. !    The calculation of these factor is explained in the Appendix A of the BLM paper
    960. 

  960.       real fww(nz_um,nz_um,ndm,nurbm)  !  from wall to wall
    961. 

  961.       real fwg(nz_um,ndm,nurbm)        !  from wall to ground
    962. 

  962.       real fgw(nz_um,ndm,nurbm)        !  from ground to wall
    963. 

  963.       real fsw(nz_um,ndm,nurbm)        !  from sky to wall
    964. 

  964.       real fws(nz_um,ndm,nurbm)        !  from wall to sky
    965. 

  965.       real fsg(ndm,nurbm)              !  from sky to ground
    966. 

  966. 

  967. !    Roughness parameters
    968. 

  968.       real z0g_u(nurbm)            ! The ground's roughness length
    969. 

  969.       real z0r_u(nurbm)            ! The roof's roughness length
    970. 

  970.       
    971. 

  971. !    Street parameters
    972. 

  972.       integer nd_u(nurbm)          ! Number of street direction for each urban class 
    973. 

  973.       real strd_u(ndm,nurbm)       ! Street length (set to a greater value then the horizontal length of the cells)
    974. 

  974.       real drst_u(ndm,nurbm)       ! Street direction
    975. 

  975.       real ws_u(ndm,nurbm)         ! Street width
    976. 

  976.       real bs_u(ndm,nurbm)         ! Building width
    977. 

  977.       real h_b(nz_um,nurbm)        ! Bulding's heights
    978. 

  978.       real d_b(nz_um,nurbm)        ! The probability that a building has an height "h_b"
    979. 

  979.       real ss_u(nz_um,nurbm)     ! The probability that a building has an height equal to "z"
    980. 

  980.       real pb_u(nz_um,nurbm)     ! The probability that a building has an height greater or equal to "z"
    981. 

  981.         
    982. 

  982. !    Grid parameters
    983. 

  983.       integer nz_u(nurbm)     ! Number of layer in the urban grid
    984. 

  984.       real z_u(nz_um)       ! Height of the urban grid levels
    985. 

  985. 

  986. !   MT
    987. 

  987.       real cop_u(nurbm)
    988. 

  988.       real pwin_u(nurbm)
    989. 

  989.       real beta_u(nurbm)
    990. 

  990.       integer sw_cond_u(nurbm)
    991. 

  991.       real time_on_u(nurbm)
    992. 

  992.       real time_off_u(nurbm)
    993. 

  993.       real targtemp_u(nurbm)
    994. 

  994.       real gaptemp_u(nurbm)
    995. 

  995.       real targhum_u(nurbm)
    996. 

  996.       real gaphum_u(nurbm)
    997. 

  997.       real perflo_u(nurbm)
    998. 

  998.       real hsesf_u(nurbm)
    999. 

  999.       real hsequip(24)
    1000. 

  1000. 

  1001. ! ----------------------------------------------------------------------
    1002. 

  1002. ! INPUT-OUTPUT
    1003. 

  1003. ! ----------------------------------------------------------------------
    1004. 

  1004. 

  1005. ! Data relative to the "urban grid" which should be stored from the current time step to the next one
    1006. 

  1006. 

  1007.       real tw(2*ndm,nz_um,nwr_u,nbui_max)  ! Temperature in each layer of the wall [K]
    1008. 

  1008.       real tr(ndm,nz_um,nwr_u)  ! Temperature in each layer of the roof [K]
    1009. 

  1009.       real tg(ndm,ng_u)          ! Temperature in each layer of the ground [K]
    1010. 

  1010.       real sfw(2*ndm,nz_um,nbui_max)      ! Sensible heat flux from walls
    1011. 

  1011.       real sfg(ndm)              ! Sensible heat flux from ground (road)
    1012. 

  1012.       real sfr(ndm,nz_um)      ! Sensible heat flux from roofs
    1013. 

  1013.       real gfg(ndm)             ! Heat flux transferred from the surface of the ground (road) towards the interior
    1014. 

  1014.       real gfr(ndm,nz_um)     ! Heat flux transferred from the surface of the roof towards the interior
    1015. 

  1015.       real gfw(2*ndm,nz_um,nbui_max)     ! Heat flux transfered from the surface of the walls towards the interior
    1016. 

  1016. ! ----------------------------------------------------------------------
    1017. 

  1017. ! OUTPUT:
    1018. 

  1018. ! ----------------------------------------------------------------------
    1019. 

  1019. 

  1020. ! Data relative to the "mesoscale grid"
    1021. 

  1021. 

  1022.       real sf(kms:kme)             ! Surface of the "mesoscale grid" cells taking into account the buildings
    1023. 

  1023.       real vl(kms:kme)               ! Volume of the "mesoscale grid" cells taking into account the buildings
    1024. 

  1024.      
    1025. 

  1025. !    Implicit and explicit components of the source and sink terms at each levels,
    1026. 

  1026. !     the fluxes can be computed as follow: FX = A*X + B   example: Heat fluxes = a_t * pt + b_t
    1027. 

  1027.       real a_u(kms:kme)              ! Implicit component of the momentum sources or sinks in the X-direction
    1028. 

  1028.       real a_v(kms:kme)              ! Implicit component of the momentum sources or sinks in the Y-direction
    1029. 

  1029.       real a_t(kms:kme)              ! Implicit component of the heat sources or sinks
    1030. 

  1030.       real a_e(kms:kme)              ! Implicit component of the TKE sources or sinks
    1031. 

  1031.       real b_u(kms:kme)              ! Explicit component of the momentum sources or sinks in the X-direction
    1032. 

  1032.       real b_v(kms:kme)              ! Explicit component of the momentum sources or sinks in the Y-direction
    1033. 

  1033.       real b_t(kms:kme)              ! Explicit component of the heat sources or sinks
    1034. 

  1034.       real b_ac(kms:kme)
    1035. 

  1035.       real b_e(kms:kme)              ! Explicit component of the TKE sources or sinks
    1036. 

  1036.       real b_q(kms:kme)              ! Explicit component of the humidity sources or sinks
    1037. 

  1037.       real dlg(kms:kme)              ! Height above ground (L_ground in formula (24) of the BLM paper). 
    1038. 

  1038.       real dl_u(kms:kme)             ! Length scale (lb in formula (22) ofthe BLM paper).
    1039. 

  1039.       
    1040. 

  1040. ! ----------------------------------------------------------------------
    1041. 

  1041. ! LOCAL:
    1042. 

  1042. ! ----------------------------------------------------------------------
    1043. 

  1043. 

  1044.       real dz(kms:kme)               ! vertical space steps of the "mesoscale grid"
    1045. 

  1045. 

  1046. ! Data interpolated from the "mesoscale grid" to the "urban grid"
    1047. 

  1047. 

  1048.       real ua_u(nz_um)          ! Wind speed in the x direction
    1049. 

  1049.       real va_u(nz_um)          ! Wind speed in the y direction
    1050. 

  1050.       real pt_u(nz_um)          ! Potential temperature
    1051. 

  1051.       real da_u(nz_um)          ! Air density
    1052. 

  1052.       real pt0_u(nz_um)         ! Reference potential temperature
    1053. 

  1053.       real pr_u(nz_um)          ! Air pressure
    1054. 

  1054.       real qv_u(nz_um)          !Specific humidity
    1055. 

  1055. 

  1056. ! Data defining the building and street charateristics
    1057. 

  1057. 

  1058.       real alag(ng_u)           ! Ground thermal diffusivity for the current urban class [m^2 s^-1] 
    1059. 

  1059.       
    1060. 

  1060.       real csg(ng_u)            ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
    1061. 

  1061.       real csr(nwr_u)            ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
    1062. 

  1062.       real csw(nwr_u)            ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
    1063. 

  1063. 

  1064.       real z0(ndm,nz_um)      ! Roughness lengths "profiles"
    1065. 

  1065.       real ws(ndm)              ! Street widths of the current urban class
    1066. 

  1066.       real bs(ndm)              ! Building widths of the current urban class
    1067. 

  1067.       real strd(ndm)            ! Street lengths for the current urban class
    1068. 

  1068.       real drst(ndm)            ! Street directions for the current urban class
    1069. 

  1069.       real ss(nz_um)          ! Probability to have a building with height h
    1070. 

  1070.       real pb(nz_um)          ! Probability to have a building with an height equal
    1071. 

  1071. 

  1072. ! Solar radiation at each level of the "urban grid"
    1073. 

  1073. 

  1074.       real rsg(ndm)             ! Short wave radiation from the ground
    1075. 

  1075.       real rsw(2*ndm,nz_um)     ! Short wave radiation from the walls
    1076. 

  1076.       real rlg(ndm)             ! Long wave radiation from the ground
    1077. 

  1077.       real rlw(2*ndm,nz_um)     ! Long wave radiation from the walls
    1078. 

  1078. 

  1079. ! Potential temperature of the surfaces at each level of the "urban grid"
    1080. 

  1080. 

  1081.       real ptg(ndm)             ! Ground potential temperatures 
    1082. 

  1082.       real ptr(ndm,nz_um)     ! Roof potential temperatures 
    1083. 

  1083.       real ptw(2*ndm,nz_um,nbui_max)     ! Walls potential temperatures 
    1084. 

  1084. 

  1085.  
    1086. 

  1086. ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
    1087. 

  1087. ! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
    1088. 

  1088. ! The fluxes can be computed as follow: Fluxes of X = A*X + B
    1089. 

  1089. ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
    1090. 

  1090. 

  1091.       real uhb_u(ndm,nz_um)   ! U (wind component) Horizontal surfaces, B (explicit) term
    1092. 

  1092.       real uva_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, A (implicit) term
    1093. 

  1093.       real uvb_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, B (explicit) term
    1094. 

  1094.       real vhb_u(ndm,nz_um)   ! V (wind component) Horizontal surfaces, B (explicit) term
    1095. 

  1095.       real vva_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, A (implicit) term
    1096. 

  1096.       real vvb_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, B (explicit) term
    1097. 

  1097.       real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
    1098. 

  1098.       real tva_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, A (implicit) term
    1099. 

  1099.       real tvb_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, B (explicit) term
    1100. 

  1100.       real tvb_ac(2*ndm,nz_um)
    1101. 

  1101.       real ehb_u(ndm,nz_um)   ! Energy (TKE)       Horizontal surfaces, B (explicit) term
    1102. 

  1102.       real evb_u(2*ndm,nz_um)   ! Energy (TKE)         Vertical surfaces, B (explicit) term
    1103. 

  1103.       real qhb_u(ndm,nz_um)     ! Humidity      Horizontal surfaces, B (explicit) term
    1104. 

  1104.       real qvb_u(2*ndm,nz_um)   ! Humidity      Vertical surfaces, B (explicit) term      
    1105. 

  1105. !
    1106. 

  1106.       real rs_abs ! solar radiation absorbed by urban surfaces 
    1107. 

  1107.       real rl_up ! longwave radiation emitted by urban surface to the atmosphere 
    1108. 

  1108.       real emiss ! mean emissivity of the urban surface
    1109. 

  1109.       real grdflx_urb ! ground heat flux
    1110. 

  1110.       real dt_int ! internal time step
    1111. 

  1111.       integer nt_int ! number of internal time step
    1112. 

  1112.       integer iz,id, it_int
    1113. 

  1113.       integer iw,ix,iy
    1114. 

  1114. 

  1115. !---------------------------------------
    1116. 

  1116. !New variables uses in BEM
    1117. 

  1117. !----------------------------------------
    1118. 

  1118.    
    1119. 

  1119.       real tmp_u(nz_um)     !Air Temperature [K]
    1120. 

  1120. 

  1121.       real dzw(nwr_u)       !Layer sizes in the walls
    1122. 

  1122.       real dzr(nwr_u)       !Layer sizes in the roofs
    1123. 

  1123.       real dzf(nf_u)        !Layer sizes in the floors
    1124. 

  1124.       real dzgb(ngb_u)      !Layer sizes in the ground below the buildings
    1125. 

  1125. 

  1126.       real csgb(ngb_u)      !Specific heat of the ground material below the buildings 
    1127. 

  1127.                             !of the current urban class at each ground levels[J m^3 K^-1]
    1128. 

  1128.       real csf(nf_u)        !Specific heat of the floors materials in the buildings 
    1129. 

  1129.                             !of the current urban class at each levels[J m^3 K^-1]
    1130. 

  1130.       real alar(nwr_u+1)    ! Roof thermal diffusivity for the current urban class [W/m K]
    1131. 

  1131.       real alaw(nwr_u+1)    ! Walls thermal diffusivity for the current urban class [W/m K] 
    1132. 

  1132.       real alaf(nf_u+1)     ! Floor thermal diffusivity at each wall layers [W/m K]     
    1133. 

  1133.       real alagb(ngb_u+1)   ! Ground thermal diffusivity below the building at each wall layer [W/m K] 
    1134. 

  1134. 

  1135.       real sfrb(ndm,nbui_max)        ! Sensible heat flux from roofs [W/m²]
    1136. 

  1136.       real gfrb(ndm,nbui_max)        ! Heat flux flowing inside the roofs [W/m²]
    1137. 

  1137.       real sfwb1D(2*ndm,nz_um)    !Sensible heat flux from the walls [W/m²] 
    1138. 

  1138.       real sfwin(2*ndm,nz_um,nbui_max)!Sensible heat flux from windows [W/m²]
    1139. 

  1139.       real sfwinb1D(2*ndm,nz_um)  !Sensible heat flux from windows [W/m²]
    1140. 

  1140.       real gfwb1D(2*ndm,nz_um)    !Heat flux flowing inside the walls [W/m²]
    1141. 

  1141. 

  1142.       real qlev(nz_um,nbui_max)      !specific humidity [kg/kg]
    1143. 

  1143.       real qlevb1D(nz_um)         !specific humidity [kg/kg] 
    1144. 

  1144.       real tlev(nz_um,nbui_max)      !Indoor temperature [K]
    1145. 

  1145.       real tlevb1D(nz_um)         !Indoor temperature [K]
    1146. 

  1146.       real twb1D(2*ndm,nwr_u,nz_um)     !Wall temperature in BEM [K]      
    1147. 

  1147.       real twlev(2*ndm,nz_um,nbui_max)     !Window temperature in BEM [K]
    1148. 

  1148.       real twlevb1D(2*ndm,nz_um)        !Window temperature in BEM [K]
    1149. 

  1149.       real tglev(ndm,ngb_u,nbui_max)        !Ground temperature below a building in BEM [K]
    1150. 

  1150.       real tglevb1D(ngb_u)               !Ground temperature below a building in BEM [K]
    1151. 

  1151.       real tflev(ndm,nf_u,nz_um-1,nbui_max)!Floor temperature in BEM[K]
    1152. 

  1152.       real tflevb1D(nf_u,nz_um-1)       !Floor temperature in BEM[K]
    1153. 

  1153.       real trb(ndm,nwr_u,nbui_max)         !Roof temperature in BEM [K]
    1154. 

  1154.       real trb1D(nwr_u)                 !Roof temperature in BEM [K]
    1155. 

  1155.       
    1156. 

  1156.       real sflev(nz_um,nz_um)     ! sensible heat flux due to the air conditioning systems [W]
    1157. 

  1157.       real lflev(nz_um,nz_um)     ! latent heat flux due to the air conditioning systems  [W]
    1158. 

  1158.       real consumlev(nz_um,nz_um) ! consumption due to the air conditioning systems [W]
    1159. 

  1159.       real sflev1D(nz_um)         ! sensible heat flux due to the air conditioning systems [W]
    1160. 

  1160.       real lflev1D(nz_um)         ! latent heat flux due to the air conditioning systems  [W]
    1161. 

  1161.       real consumlev1D(nz_um)     ! consumption due to the air conditioning systems [W]
    1162. 

  1162.       real sfvlev(nz_um,nz_um)    ! sensible heat flux due to ventilation [W]
    1163. 

  1163.       real lfvlev(nz_um,nz_um)    ! latent heat flux due to ventilation [W]
    1164. 

  1164.       real sfvlev1D(nz_um)        ! sensible heat flux due to ventilation [W]
    1165. 

  1165.       real lfvlev1D(nz_um)        ! Latent heat flux due to ventilation [W]
    1166. 

  1166. 

  1167.       real ptwin(2*ndm,nz_um,nbui_max)  ! window potential temperature
    1168. 

  1168.       real tw_av(2*ndm,nz_um)        ! Averaged temperature of the wall surfaces
    1169. 

  1169.       real twlev_av(2*ndm,nz_um)     ! Averaged temperature of the windows
    1170. 

  1170.       real sfw_av(2*ndm,nz_um)       ! Averaged sensible heat from walls
    1171. 

  1171.       real sfwind_av(2*ndm,nz_um)    ! Averaged sensible heat from windows
    1172. 

  1172. 

  1173.       integer nbui                !Total number of different type of buildings in an urban class
    1174. 

  1174.       integer nlev(nz_um)         !Number of levels in each different type of buildings in an urban class
    1175. 

  1175.       integer ibui,ily  
    1176. 

  1176.       real :: nhourday   ! Number of hours from midnight, local time
    1177. 

  1177. ! ----------------------------------------------------------------------
    1178. 

  1178. ! END VARIABLES DEFINITIONS
    1179. 

  1179. ! ----------------------------------------------------------------------
    1180. 

  1180. 

  1181. ! Fix some usefull parameters for the computation of the sources or sinks
    1182. 

  1182. !
    1183. 

  1183. !initialize inside param
    1184. 

  1184. !
    1185. 

  1185. !      ss=0.
    1186. 

  1186. !      pb=0.
    1187. 

  1187. 

  1188.       do iz=kts,kte
    1189. 

  1189.          dz(iz)=z(iz+1)-z(iz)
    1190. 

  1190.       end do
    1191. 

  1191.       call param(iurb,nz_u(iurb),nd_u(iurb),                          &
    1192. 

  1192.                        csg_u,csg,alag_u,alag,csr_u,csr,               &
    1193. 

  1193.                        alar_u,alar,csw_u,csw,alaw_u,alaw,             &
    1194. 

  1194.                        ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0,                &
    1195. 

  1195.                        strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb,       &
    1196. 

  1196.                        dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)
    1197. 

  1197. 

  1198. ! Interpolation on the "urban grid"
    1199. 

  1199. 

  1200.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,ua,ua_u)
    1201. 

  1201.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,va,va_u)
    1202. 

  1202.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt,pt_u)
    1203. 

  1203.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt0,pt0_u)
    1204. 

  1204.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pr,pr_u)
    1205. 

  1205.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,da,da_u)
    1206. 

  1206.       call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,qv,qv_u)
    1207. 

  1207.                    
    1208. 

  1208. ! Compute the modification of the radiation due to the buildings
    1209. 

  1209. 

  1210.       call averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av, &
    1211. 

  1211.                           sfw_av,sfwind_av,sfw,sfwin)
    1212. 

  1212.       call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws,           &
    1213. 

  1213.                     drst,strd,ss,pb,                              &
    1214. 

  1214.                     tw_av,tg,twlev_av,albg_u(iurb),albw_u(iurb),  &
    1215. 

  1215.                     emw_u(iurb),emg_u(iurb),pwin_u(iurb),albwin_u(iurb),  &
    1216. 

  1216.                     emwind_u(iurb),fww,fwg,fgw,fsw,fsg,           &
    1217. 

  1217.                     zr,deltar,ah,                                 &
    1218. 

  1218.                     rs,rld,rsw,rsg,rlw,rlg)                       
    1219. 

  1219. 

  1220. ! calculation of the urban albedo and the upward long wave radiation
    1221. 

  1221. 

  1222.        call upward_rad(nd_u(iurb),nz_u(iurb),ws,bs,sigma,pb,ss,                &
    1223. 

  1223.                        tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg,                & 
    1224. 

  1224.                        tw_av,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw_av,          & 
    1225. 

  1225.                        tr,emr_u(iurb),albr_u(iurb),emwind_u(iurb),             &
    1226. 

  1226.                        albwin_u(iurb),twlev_av,pwin_u(iurb),sfwind_av,rld,rs,sfr,      & 
    1227. 

  1227.                        rs_abs,rl_up,emiss,grdflx_urb)               
    1228. 

  1228. 

  1229. ! Compute the surface temperatures
    1230. 

  1230.       
    1231. 

  1231.       call surf_temp(nd_u(iurb),pr_u,dt,              & 
    1232. 

  1232.                     rld,rsg,rlg,                      &
    1233. 

  1233.                     tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg)
    1234. 

  1234. 

  1235. 

  1236. ! Call the BEM (Building Energy Model) routine 
    1237. 

  1237.        
    1238. 

  1238.        do iz=1,nz_um !Compute the outdoor temperature 
    1239. 

  1239. 	 tmp_u(iz)=pt_u(iz)*(pr_u(iz)/p0)**(rcp_u) 
    1240. 

  1240.        end do
    1241. 

  1241. 

  1242.        ibui=0
    1243. 

  1243.        nlev=0
    1244. 

  1244.        nbui=0
    1245. 

  1245.      
    1246. 

  1246.        sfrb=0.     !Sensible heat flux from roof
    1247. 

  1247.        gfrb=0.     !Heat flux flowing inside the roof
    1248. 

  1248.        sfwb1D=0.   !Sensible heat flux from walls
    1249. 

  1249.        sfwinb1D=0. !Sensible heat flux from windows
    1250. 

  1250.        gfwb1D=0.   !Heat flux flowing inside the walls[W/m²]
    1251. 

  1251. 

  1252. 

  1253.        twb1D=0.    !Wall temperature
    1254. 

  1254.        twlevb1D=0. !Window temperature
    1255. 

  1255.        tglevb1D=0. !Ground temperature below a building
    1256. 

  1256.        tflevb1D=0. !Floor temperature     
    1257. 

  1257.        trb=0.      !Roof temperature
    1258. 

  1258.        trb1D=0.    !Roof temperature
    1259. 

  1259. 

  1260.        qlevb1D=0. !Indoor humidity
    1261. 

  1261.        tlevb1D=0. !indoor temperature
    1262. 

  1262. 

  1263.        sflev1D=0.    !Sensible heat flux from the a.c.
    1264. 

  1264.        lflev1D=0.    !Latent heat flux from the a.c.
    1265. 

  1265.        consumlev1D=0.!Consumption from the a.c.
    1266. 

  1266.        sfvlev1D=0.   !Sensible heat flux from the natural ventilation
    1267. 

  1267.        lfvlev1D=0.   !Latent heat flux from natural ventilation
    1268. 

  1268. 

  1269.        ptw=0.        !Wall potential temperature
    1270. 

  1270.        ptwin=0.      !Window potential temperature
    1271. 

  1271.        ptr=0.        !Roof potential temperature
    1272. 

  1272.        
    1273. 

  1273.        do iz=1,nz_um		   
    1274. 

  1274.          if(ss(iz).gt.0) then		
    1275. 

  1275.            ibui=ibui+1		                
    1276. 

  1276.            nlev(ibui)=iz-1
    1277. 

  1277.            nbui=ibui
    1278. 

  1278.            do id=1,ndm
    1279. 

  1279.               sfrb(id,ibui)=sfr(id,iz)
    1280. 

  1280.               do ily=1,nwr_u
    1281. 

  1281.                  trb(id,ily,ibui)=tr(id,iz,ily)
    1282. 

  1282.               enddo
    1283. 

  1283.            enddo
    1284. 

  1284. 	 endif	  
    1285. 

  1285.        end do  !iz    
    1286. 

  1286. !--------------------------------------------------------------------------------
    1287. 

  1287. !Loop over BEM  -----------------------------------------------------------------
    1288. 

  1288. !--------------------------------------------------------------------------------
    1289. 

  1289. !--------------------------------------------------------------------------------
    1290. 

  1290.    
    1291. 

  1291.        nhourday=ah/PI*180./15.+12.
    1292. 

  1292.        if (nhourday >= 24) nhourday = nhourday - 24
    1293. 

  1293.        if (nhourday < 0)  nhourday = nhourday + 24
    1294. 

  1294. 

  1295.         do ibui=1,nbui
    1296. 

  1296.          
    1297. 

  1297.           do iz=1,nz_um
    1298. 

  1298.              qlevb1D(iz)=qlev(iz,ibui)
    1299. 

  1299.              tlevb1D(iz)=tlev(iz,ibui) 
    1300. 

  1300.           enddo
    1301. 

  1301.              
    1302. 

  1302.           do id=1,ndm
    1303. 

  1303. 

  1304.              do ily=1,nwr_u
    1305. 

  1305.                 trb1D(ily)=trb(id,ily,ibui)
    1306. 

  1306.              enddo
    1307. 

  1307.              do ily=1,ngb_u
    1308. 

  1308.                 tglevb1D(ily)=tglev(id,ily,ibui) 
    1309. 

  1309.              enddo
    1310. 

  1310.               
    1311. 

  1311.              do ily=1,nf_u
    1312. 

  1312.              do iz=1,nz_um-1
    1313. 

  1313.                tflevb1D(ily,iz)=tflev(id,ily,iz,ibui)
    1314. 

  1314.              enddo
    1315. 

  1315.              enddo
    1316. 

  1316.             
    1317. 

  1317.              do iz=1,nz_um
    1318. 

  1318.                 sfwinb1D(2*id-1,iz)=sfwin(2*id-1,iz,ibui)
    1319. 

  1319.                 sfwinb1D(2*id,iz)=sfwin(2*id,iz,ibui)
    1320. 

  1320.              enddo
    1321. 

  1321.             
    1322. 

  1322.              do iz=1,nz_um
    1323. 

  1323.                 do ily=1,nwr_u
    1324. 

  1324.                    twb1D(2*id-1,ily,iz)=tw(2*id-1,iz,ily,ibui)
    1325. 

  1325.                    twb1D(2*id,ily,iz)=tw(2*id,iz,ily,ibui)
    1326. 

  1326.                 enddo
    1327. 

  1327.                 sfwb1D(2*id-1,iz)=sfw(2*id-1,iz,ibui)
    1328. 

  1328.                 sfwb1D(2*id,iz)=sfw(2*id,iz,ibui)
    1329. 

  1329.                 twlevb1D(2*id-1,iz)=twlev(2*id-1,iz,ibui)
    1330. 

  1330.                 twlevb1D(2*id,iz)=twlev(2*id,iz,ibui)
    1331. 

  1331.              enddo
    1332. 

  1332.           enddo   
    1333. 

  1333. 

  1334.           call BEM(nz_um,nlev(ibui),nhourday,dt,bs_u(1,iurb),                   &
    1335. 

  1335.                    bs_u(2,iurb),dz_u,nwr_u,nf_u,nwr_u,ngb_u,sfwb1D,gfwb1D,      &
    1336. 

  1336.                    sfwinb1D,sfrb(1,ibui),gfrb(1,ibui),                          &
    1337. 

  1337.                    latent,sigma,albw_u(iurb),albwin_u(iurb),albr_u(iurb),       &
    1338. 

  1338.      	           emr_u(iurb),emw_u(iurb),emwind_u(iurb),rsw,rlw,r,cp_u,       &
    1339. 

  1339.      	           da_u,tmp_u,qv_u,pr_u,rs,rld,dzw,csw,alaw,pwin_u(iurb),       &
    1340. 

  1340.                    cop_u(iurb),beta_u(iurb),sw_cond_u(iurb),time_on_u(iurb),    &
    1341. 

  1341.                    time_off_u(iurb),targtemp_u(iurb),gaptemp_u(iurb),           &
    1342. 

  1342.                    targhum_u(iurb),gaphum_u(iurb),perflo_u(iurb),hsesf_u(iurb), &
    1343. 

  1343.                    hsequip,                                                     &
    1344. 

  1344.      	           dzf,csf,alaf,dzgb,csgb,alagb,dzr,csr,                        &
    1345. 

  1345.      	           alar,tlevb1D,qlevb1D,twb1D,twlevb1D,tflevb1D,tglevb1D,       &
    1346. 

  1346.      	           trb1D,sflev1D,lflev1D,consumlev1D,sfvlev1D,lfvlev1D)
    1347. 

  1347.                     
    1348. 

  1348.         
    1349. 

  1349. !
    1350. 

  1350. !Temporal modifications
    1351. 

  1351. !
    1352. 

  1352.          sfrb(2,ibui)=sfrb(1,ibui)
    1353. 

  1353.          gfrb(2,ibui)=gfrb(1,ibui)
    1354. 

  1354. !
    1355. 

  1355. !End temporal modifications  
    1356. 

  1356. !        
    1357. 

  1357.            do iz=1,nz_um
    1358. 

  1358.              qlev(iz,ibui)=qlevb1D(iz)
    1359. 

  1359.              tlev(iz,ibui)=tlevb1D(iz)
    1360. 

  1360.              sflev(iz,ibui)=sflev1D(iz)
    1361. 

  1361.              lflev(iz,ibui)=lflev1D(iz)
    1362. 

  1362.              consumlev(iz,ibui)=consumlev1D(iz)
    1363. 

  1363.              sfvlev(iz,ibui)=sfvlev1D(iz)
    1364. 

  1364.              lfvlev(iz,ibui)=lfvlev1D(iz)
    1365. 

  1365.            enddo
    1366. 

  1366.  
    1367. 

  1367.            do id=1,ndm
    1368. 

  1368.               do ily=1,nwr_u
    1369. 

  1369.                  trb(id,ily,ibui)=trb1D(ily)
    1370. 

  1370.               enddo   
    1371. 

  1371.               do ily=1,ngb_u
    1372. 

  1372.                  tglev(id,ily,ibui)=tglevb1D(ily) 
    1373. 

  1373.               enddo
    1374. 

  1374. 

  1375.               do ily=1,nf_u
    1376. 

  1376.               do iz=1,nz_um-1
    1377. 

  1377.                  tflev(id,ily,iz,ibui)=tflevb1D(ily,iz)
    1378. 

  1378.               enddo
    1379. 

  1379.               enddo
    1380. 

  1380.            
    1381. 

  1381. 

  1382.              do iz=1,nz_um
    1383. 

  1383.                 do ily=1,nwr_u
    1384. 

  1384.                    tw(2*id-1,iz,ily,ibui)=twb1D(2*id-1,ily,iz)
    1385. 

  1385.                    tw(2*id,iz,ily,ibui)=twb1D(2*id,ily,iz)
    1386. 

  1386.                 enddo
    1387. 

  1387.                 gfw(2*id-1,iz,ibui)=gfwb1D(2*id-1,iz)
    1388. 

  1388.                 gfw(2*id,iz,ibui)=gfwb1D(2*id,iz)
    1389. 

  1389.                 twlev(2*id-1,iz,ibui)=twlevb1D(2*id-1,iz)
    1390. 

  1390.                 twlev(2*id,iz,ibui)=twlevb1D(2*id,iz)
    1391. 

  1391.              enddo
    1392. 

  1392.            enddo       
    1393. 

  1393. 

  1394.         enddo !ibui     
    1395. 

  1395. 

  1396. !-----------------------------------------------------------------------------
    1397. 

  1397. !End loop over BEM -----------------------------------------------------------
    1398. 

  1398. !-----------------------------------------------------------------------------
    1399. 

  1399. !-----------------------------------------------------------------------------
    1400. 

  1400. 

  1401.        ibui=0
    1402. 

  1402.        do iz=1,nz_um		   
    1403. 

  1403.          if(ss(iz).gt.0) then		
    1404. 

  1404.            ibui=ibui+1	
    1405. 

  1405.            do id=1,ndm	
    1406. 

  1406.               gfr(id,iz)=gfrb(id,ibui)
    1407. 

  1407.               sfr(id,iz)=sfrb(id,ibui)
    1408. 

  1408.               do ily=1,nwr_u
    1409. 

  1409.                  tr(id,iz,ily)=trb(id,ily,ibui)
    1410. 

  1410.               enddo
    1411. 

  1411.               ptr(id,iz)=tr(id,iz,nwr_u)*(pr_u(iz)/p0)**(-rcp_u)
    1412. 

  1412.            enddo
    1413. 

  1413.          endif
    1414. 

  1414.        enddo !iz
    1415. 

  1415. 

  1416. !Compute the potential temperature for the vertical surfaces of the buildings
    1417. 

  1417.       
    1418. 

  1418.        do id=1,ndm
    1419. 

  1419.           do iz=1,nz_um
    1420. 

  1420.              do ibui=1,nbui
    1421. 

  1421.                 ptw(2*id-1,iz,ibui)=tw(2*id-1,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
    1422. 

  1422.                 ptw(2*id,iz,ibui)=tw(2*id,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
    1423. 

  1423.                 ptwin(2*id-1,iz,ibui)=twlev(2*id-1,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
    1424. 

  1424.                 ptwin(2*id,iz,ibui)=twlev(2*id,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
    1425. 

  1425.              enddo
    1426. 

  1426.           enddo
    1427. 

  1427.        enddo
    1428. 

  1428.              
    1429. 

  1429. ! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
    1430. 

  1430.      
    1431. 

  1431.       call buildings(iurb,nd_u(iurb),nz_u(iurb),z0,ua_u,va_u,                 & 
    1432. 

  1432.                      pt_u,pt0_u,ptg,ptr,da_u,ptw,ptwin,pwin_u(iurb),drst,     &                      
    1433. 

  1433.                      uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,qvb_u,qhb_u,   & 
    1434. 

  1434.                      uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr,               &
    1435. 

  1435.                      sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)  
    1436. 

  1436.      
    1437. 

  1437. ! Calculation of the sensible heat fluxes for the ground, the wall and roof
    1438. 

  1438. ! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
    1439. 

  1439. ! where A and B are the implicit and explicit components of the heat sources or sinks.
    1440. 

  1440. 

  1441.       
    1442. 

  1442. ! Interpolation on the "mesoscale grid"
    1443. 

  1443. 

  1444.       call urban_meso(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z,dz,z_u,pb,ss,bs,ws,sf, & 
    1445. 

  1445.                      vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,                   &
    1446. 

  1446.                      uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u,                            &
    1447. 

  1447.                      a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)                    
    1448. 

  1448.        
    1449. 

  1449. 

  1450. ! Calculation of the length scale taking into account the buildings effects
    1451. 

  1451. 

  1452.       call interp_length(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z_u,z,ss,ws,bs,dlg,dl_u)
    1453. 

  1453. 

  1454.       return
    1455. 

  1455.       end subroutine BEP1D
    1456. 

  1456. 

  1457. ! ===6=8===============================================================72
    1458. 

  1458. ! ===6=8===============================================================72
    1459. 

  1459. 

  1460.        subroutine param(iurb,nz,nd,                                   &
    1461. 

  1461.                        csg_u,csg,alag_u,alag,csr_u,csr,               &
    1462. 

  1462.                        alar_u,alar,csw_u,csw,alaw_u,alaw,             &
    1463. 

  1463.                        ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0,                &  
    1464. 

  1464.                        strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb,       &
    1465. 

  1465.                        dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)        
    1466. 

  1466. 

  1467. ! ----------------------------------------------------------------------
    1468. 

  1468. !    This routine prepare some usefull parameters       
    1469. 

  1469. ! ----------------------------------------------------------------------
    1470. 

  1470. 

  1471.       implicit none
    1472. 

  1472. 

  1473.   
    1474. 

  1474. ! ----------------------------------------------------------------------
    1475. 

  1475. ! INPUT:
    1476. 

  1476. ! ----------------------------------------------------------------------
    1477. 

  1477.       integer iurb                 ! Current urban class
    1478. 

  1478.       integer nz                   ! Number of vertical urban levels in the current class
    1479. 

  1479.       integer nd                   ! Number of street direction for the current urban class
    1480. 

  1480.       real alag_u(nurbm)           ! Ground thermal diffusivity [m^2 s^-1]
    1481. 

  1481.       real alar_u(nurbm)           ! Roof thermal diffusivity [m^2 s^-1]
    1482. 

  1482.       real alaw_u(nurbm)           ! Wall thermal diffusivity [m^2 s^-1]
    1483. 

  1483.       real bs_u(ndm,nurbm)         ! Building width
    1484. 

  1484.       real csg_u(nurbm)            ! Specific heat of the ground material [J m^3 K^-1]
    1485. 

  1485.       real csr_u(nurbm)            ! Specific heat of the roof material [J m^3 K^-1]
    1486. 

  1486.       real csw_u(nurbm)            ! Specific heat of the wall material [J m^3 K^-1]
    1487. 

  1487.       real drst_u(ndm,nurbm)       ! Street direction
    1488. 

  1488.       real strd_u(ndm,nurbm)       ! Street length 
    1489. 

  1489.       real ws_u(ndm,nurbm)         ! Street width
    1490. 

  1490.       real z0g_u(nurbm)            ! The ground's roughness length
    1491. 

  1491.       real z0r_u(nurbm)            ! The roof's roughness length
    1492. 

  1492.       real ss_u(nz_um,nurbm)       ! The probability that a building has an height equal to "z"
    1493. 

  1493.       real pb_u(nz_um,nurbm)       ! The probability that a building has an height greater or equal to "z"
    1494. 

  1494.      
    1495. 

  1495. 

  1496. ! ----------------------------------------------------------------------
    1497. 

  1497. ! OUTPUT:
    1498. 

  1498. ! ----------------------------------------------------------------------
    1499. 

  1499.       real alag(ng_u)           ! Ground thermal diffusivity at each ground levels
    1500. 

  1500.       real csg(ng_u)            ! Specific heat of the ground material at each ground levels
    1501. 

  1501.       real bs(ndm)              ! Building width for the current urban class
    1502. 

  1502.       real drst(ndm)            ! street directions for the current urban class
    1503. 

  1503.       real strd(ndm)            ! Street lengths for the current urban class
    1504. 

  1504.       real ws(ndm)              ! Street widths of the current urban class
    1505. 

  1505.       real z0(ndm,nz_um)      ! Roughness lengths "profiles"
    1506. 

  1506.       real ss(nz_um)          ! Probability to have a building with height h
    1507. 

  1507.       real pb(nz_um)          ! Probability to have a building with an height greater or equal to "z"
    1508. 

  1508. 

  1509. !-----------------------------------------------------------------------------
    1510. 

  1510. !INPUT/OUTPUT
    1511. 

  1511. !-----------------------------------------------------------------------------
    1512. 

  1512. 

  1513.       real dzw(nwr_u)       !Layer sizes in the walls [m]
    1514. 

  1514.       real dzr(nwr_u)       !Layer sizes in the roofs [m]
    1515. 

  1515.       real dzf(nf_u)        !Layer sizes in the floors [m]
    1516. 

  1516.       real dzgb(ngb_u)      !layer sizes in the ground below the buildings [m]
    1517. 

  1517. 

  1518.       real csr(nwr_u)       ! Specific heat of the roof material at each roof levels
    1519. 

  1519.       real csw(nwr_u)       ! Specific heat of the wall material at each wall levels
    1520. 

  1520. 

  1521.       real csf(nf_u)        !Specific heat of the floors materials in the buildings 
    1522. 

  1522.                             !of the current urban class [J m^3 K^-1]
    1523. 

  1523.       real csgb(ngb_u)      !Specific heat of the ground material below the buildings 
    1524. 

  1524.                             !of the current urban class [J m^3 K^-1]
    1525. 

  1525.       real alar(nwr_u+1)    ! Roof thermal diffusivity at each roof levels [W/ m K]
    1526. 

  1526.       real alaw(nwr_u+1)    ! Wall thermal diffusivity at each wall levels [W/ m K]
    1527. 

  1527.       real alaf(nf_u+1)     ! Floor thermal diffusivity at each wall levels [W/m K]
    1528. 

  1528.       real alagb(ngb_u+1)   ! Ground thermal diffusivity below the building at each wall levels [W/m K]
    1529. 

  1529. ! ----------------------------------------------------------------------
    1530. 

  1530. ! LOCAL:
    1531. 

  1531. ! ----------------------------------------------------------------------
    1532. 

  1532.       integer id,ig,ir,iw,iz,iflo
    1533. 

  1533. ! ----------------------------------------------------------------------
    1534. 

  1534. ! END VARIABLES DEFINITIONS
    1535. 

  1535. ! ----------------------------------------------------------------------  
    1536. 

  1536. !Define the layer sizes in the walls
    1537. 

  1537.         
    1538. 

  1538.       ss=0.
    1539. 

  1539.       pb=0.
    1540. 

  1540.       csg=0.
    1541. 

  1541.       alag=0.
    1542. 

  1542.       csr=0.
    1543. 

  1543.       alar=0.
    1544. 

  1544.       csw=0.
    1545. 

  1545.       alaw=0.
    1546. 

  1546.       z0=0.
    1547. 

  1547.       ws=0.
    1548. 

  1548.       bs=0.
    1549. 

  1549.       strd=0.
    1550. 

  1550.       drst=0.
    1551. 

  1551.       csgb=0.
    1552. 

  1552.       alagb=0.
    1553. 

  1553.       csf=0.
    1554. 

  1554.       alaf=0.
    1555. 

  1555.  
    1556. 

  1556.       dzgb=(/0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/)
    1557. 

  1557.       dzr=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)   
    1558. 

  1558.       dzw=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
    1559. 

  1559.       dzf=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02/)   
    1560. 

  1560.       
    1561. 

  1561.       do ig=1,ngb_u
    1562. 

  1562.         csgb(ig) = csg_u(iurb)
    1563. 

  1563.         alagb(ig)= csg_u(iurb)*alag_u(iurb)
    1564. 

  1564.       enddo
    1565. 

  1565.       alagb(ngb_u+1)= csg_u(iurb)*alag_u(iurb)
    1566. 

  1566.       
    1567. 

  1567.       do iflo=1,nf_u
    1568. 

  1568.         csf(iflo) = csw_u(iurb)
    1569. 

  1569.         alaf(iflo)= csw_u(iurb)*alaw_u(iurb) 
    1570. 

  1570.       enddo
    1571. 

  1571.       alaf(nf_u+1)= csw_u(iurb)*alaw_u(iurb) 
    1572. 

  1572.      
    1573. 

  1573.       do ir=1,nwr_u
    1574. 

  1574.         csr(ir) = csr_u(iurb)
    1575. 

  1575.         alar(ir)= csr_u(iurb)*alar_u(iurb)
    1576. 

  1576.       enddo
    1577. 

  1577.       alar(nwr_u+1)= csr_u(iurb)*alar_u(iurb)
    1578. 

  1578. 

  1579.       do iw=1,nwr_u
    1580. 

  1580.         csw(iw) = csw_u(iurb)
    1581. 

  1581.         alaw(iw)= csw_u(iurb)*alaw_u(iurb)
    1582. 

  1582.       enddo
    1583. 

  1583.       alaw(nwr_u+1)=csw_u(iurb)*alaw_u(iurb) 
    1584. 

  1584. 

  1585. !------------------------------------------------------------------------  
    1586. 

  1586.         
    1587. 

  1587.       do iz=1,nz+1
    1588. 

  1588.          ss(iz)=ss_u(iz,iurb)
    1589. 

  1589.          pb(iz)=pb_u(iz,iurb)
    1590. 

  1590.       end do
    1591. 

  1591.                  
    1592. 

  1592.        do ig=1,ng_u
    1593. 

  1593.         csg(ig)=csg_u(iurb)
    1594. 

  1594.         alag(ig)=alag_u(iurb)
    1595. 

  1595.        enddo
    1596. 

  1596.        
    1597. 

  1597.        do id=1,nd
    1598. 

  1598.         z0(id,1)=z0g_u(iurb)
    1599. 

  1599.         do iz=2,nz+1
    1600. 

  1600.          z0(id,iz)=z0r_u(iurb)
    1601. 

  1601.         enddo
    1602. 

  1602.        enddo
    1603. 

  1603.       
    1604. 

  1604.        do id=1,nd
    1605. 

  1605.         ws(id)=ws_u(id,iurb)
    1606. 

  1606.         bs(id)=bs_u(id,iurb)
    1607. 

  1607.         strd(id)=strd_u(id,iurb)
    1608. 

  1608.         drst(id)=drst_u(id,iurb)     
    1609. 

  1609.        enddo
    1610. 

  1610. 

  1611.        
    1612. 

  1612.        return
    1613. 

  1613.        end subroutine param
    1614. 

  1614.        
    1615. 

  1615. ! ===6=8===============================================================72
    1616. 

  1616. ! ===6=8===============================================================72
    1617. 

  1617. 

  1618.       subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
    1619. 

  1619. 

  1620. ! ----------------------------------------------------------------------
    1621. 

  1621. !  This routine interpolate para
    1622. 

  1622. !  meters from the "mesoscale grid" to
    1623. 

  1623. !  the "urban grid".
    1624. 

  1624. !  See p300 Appendix B.1 of the BLM paper.
    1625. 

  1625. ! ----------------------------------------------------------------------
    1626. 

  1626. 

  1627.       implicit none
    1628. 

  1628. 

  1629. ! ----------------------------------------------------------------------
    1630. 

  1630. ! INPUT:
    1631. 

  1631. ! ----------------------------------------------------------------------
    1632. 

  1632. ! Data relative to the "mesoscale grid"
    1633. 

  1633.       integer kts,kte,kms,kme            
    1634. 

  1634.       real z(kms:kme)          ! Altitude of the cell interface
    1635. 

  1635.       real c(kms:kme)            ! Parameter which has to be interpolated
    1636. 

  1636. ! Data relative to the "urban grid"
    1637. 

  1637.       integer nz_u          ! Number of levels
    1638. 

  1638.       real z_u(nz_um)      ! Altitude of the cell interface
    1639. 

  1639. ! ----------------------------------------------------------------------
    1640. 

  1640. ! OUTPUT:
    1641. 

  1641. ! ----------------------------------------------------------------------
    1642. 

  1642.       real c_u(nz_um)        ! Interpolated paramters in the "urban grid"
    1643. 

  1643.        
    1644. 

  1644. ! LOCAL:
    1645. 

  1645. ! ----------------------------------------------------------------------
    1646. 

  1646.       integer iz_u,iz
    1647. 

  1647.       real ctot,dz
    1648. 

  1648. 

  1649. ! ----------------------------------------------------------------------
    1650. 

  1650. ! END VARIABLES DEFINITIONS
    1651. 

  1651. ! ----------------------------------------------------------------------
    1652. 

  1652. 

  1653.        do iz_u=1,nz_u
    1654. 

  1654.         ctot=0.
    1655. 

  1655.         do iz=kts,kte
    1656. 

  1656.          dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
    1657. 

  1657.          ctot=ctot+c(iz)*dz
    1658. 

  1658.         enddo
    1659. 

  1659.         c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
    1660. 

  1660.        enddo
    1661. 

  1661.        
    1662. 

  1662.        return
    1663. 

  1663.        end subroutine interpol
    1664. 

  1664.          
    1665. 

  1665. ! ===6=8===============================================================72       
    1666. 

  1666. ! ===6=8===============================================================72    
    1667. 

  1667. 

  1668.       subroutine  averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av,&
    1669. 

  1669.                                  sfw_av,sfwind_av,sfw,sfwin) 
    1670. 

  1670. 

  1671.       implicit none
    1672. 

  1672. !
    1673. 

  1673. !INPUT VARIABLES
    1674. 

  1674. !
    1675. 

  1675.       real tw(2*ndm,nz_um,nwr_u,nbui_max)        ! Temperature in each layer of the wall [K]
    1676. 

  1676.       real twlev(2*ndm,nz_um,nbui_max)     ! Window temperature in BEM [K]
    1677. 

  1677.       real pb(nz_um)                    ! Probability to have a building with an height equal or greater h
    1678. 

  1678.       real ss(nz_um)                    ! Probability to have a building with height h
    1679. 

  1679.       real sfw(2*ndm,nz_um,nbui_max)             ! Surface fluxes from the walls
    1680. 

  1680.       real sfwin(2*ndm,nz_um,nbui_max)     ! Surface fluxes from the windows
    1681. 

  1681. !
    1682. 

  1682. !OUTPUT VARIABLES
    1683. 

  1683. !
    1684. 

  1684.       real tw_av(2*ndm,nz_um)           ! Averaged temperature of the wall surfaces
    1685. 

  1685.       real twlev_av(2*ndm,nz_um)        ! Averaged temperature of the windows
    1686. 

  1686.       real sfw_av(2*ndm,nz_um)          ! Averaged sensible heat from walls
    1687. 

  1687.       real sfwind_av(2*ndm,nz_um)       ! Averaged sensible heat from windows
    1688. 

  1688. !
    1689. 

  1689. !LOCAL VARIABLES
    1690. 

  1690. !
    1691. 

  1691.       real d_urb(nz_um)    
    1692. 

  1692.       integer nlev(nz_um)            
    1693. 

  1693.       integer id,iz
    1694. 

  1694.       integer nbui,ibui
    1695. 

  1695. !
    1696. 

  1696. !Initialize Variables
    1697. 

  1697. !
    1698. 

  1698.       tw_av=0.
    1699. 

  1699.       twlev_av=0.
    1700. 

  1700.       sfw_av=0.
    1701. 

  1701.       sfwind_av=0.
    1702. 

  1702.       ibui=0
    1703. 

  1703.       nbui=0
    1704. 

  1704.       nlev=0
    1705. 

  1705.       d_urb=0.
    1706. 

  1706. 

  1707.       do iz=1,nz_um		   
    1708. 

  1708.          if(ss(iz).gt.0) then		
    1709. 

  1709.            ibui=ibui+1
    1710. 

  1710.            d_urb(ibui)=ss(iz)
    1711. 

  1711.            nlev(ibui)=iz-1
    1712. 

  1712.            nbui=ibui		               
    1713. 

  1713.          endif
    1714. 

  1714.       enddo
    1715. 

  1715. 

  1716.       do id=1,ndm
    1717. 

  1717.          do iz=1,nz_um-1
    1718. 

  1718.             if (pb(iz+1).gt.0) then
    1719. 

  1719.                 do ibui=1,nbui
    1720. 

  1720.                    if (iz.le.nlev(ibui)) then
    1721. 

  1721.                       tw_av(2*id-1,iz)=tw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
    1722. 

  1722.                                        tw(2*id-1,iz,nwr_u,ibui)**4
    1723. 

  1723.                       tw_av(2*id,iz)=tw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
    1724. 

  1724.                                      tw(2*id,iz,nwr_u,ibui)**4
    1725. 

  1725.                       twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
    1726. 

  1726.                                           twlev(2*id-1,iz,ibui)**4
    1727. 

  1727.                       twlev_av(2*id,iz)=twlev_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
    1728. 

  1728.                                         twlev(2*id,iz,ibui)**4
    1729. 

  1729.                       sfw_av(2*id-1,iz)=sfw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id-1,iz,ibui)
    1730. 

  1730.                       sfw_av(2*id,iz)=sfw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id,iz,ibui)
    1731. 

  1731.                       sfwind_av(2*id-1,iz)=sfwind_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id-1,iz,ibui)
    1732. 

  1732.                       sfwind_av(2*id,iz)=sfwind_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id,iz,ibui)
    1733. 

  1733.                    endif
    1734. 

  1734.                 enddo
    1735. 

  1735.                 tw_av(2*id-1,iz)=tw_av(2*id-1,iz)**(1./4.)
    1736. 

  1736.                 tw_av(2*id,iz)=tw_av(2*id,iz)**(1./4.)
    1737. 

  1737.                 twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)**(1./4.)
    1738. 

  1738.                 twlev_av(2*id,iz)=twlev_av(2*id,iz)**(1./4.)
    1739. 

  1739.             endif
    1740. 

  1740.          enddo !iz         
    1741. 

  1741.       enddo !id
    1742. 

  1742.       return
    1743. 

  1743.       end subroutine averaging_temp
    1744. 

  1744. ! ===6=8===============================================================72       
    1745. 

  1745. ! ===6=8===============================================================72    
    1746. 

  1746. 

  1747.       subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb,    &
    1748. 

  1748.                           tw,tg,twlev,albg,albw,emw,emg,pwin,albwin,   &
    1749. 

  1749.                           emwin,fww,fwg,fgw,fsw,fsg,             &
    1750. 

  1750.                           zr,deltar,ah,                          &    
    1751. 

  1751.                           rs,rl,rsw,rsg,rlw,rlg)                       
    1752. 

  1752.  
    1753. 

  1753. ! ----------------------------------------------------------------------
    1754. 

  1754. ! This routine computes the modification of the short wave and 
    1755. 

  1755. !  long wave radiation due to the buildings.
    1756. 

  1756. ! ----------------------------------------------------------------------
    1757. 

  1757. 

  1758.       implicit none
    1759. 

  1759.  
    1760. 

  1760.  
    1761. 

  1761. ! ----------------------------------------------------------------------
    1762. 

  1762. ! INPUT:
    1763. 

  1763. ! ----------------------------------------------------------------------
    1764. 

  1764.       integer iurb              ! current urban class
    1765. 

  1765.       integer nd                ! Number of street direction for the current urban class
    1766. 

  1766.       integer nz_u              ! Number of layer in the urban grid
    1767. 

  1767.       real z(nz_um)           ! Height of the urban grid levels
    1768. 

  1768.       real ws(ndm)              ! Street widths of the current urban class
    1769. 

  1769.       real drst(ndm)            ! street directions for the current urban class
    1770. 

  1770.       real strd(ndm)            ! Street lengths for the current urban class
    1771. 

  1771.       real ss(nz_um)          ! probability to have a building with height h
    1772. 

  1772.       real pb(nz_um)          ! probability to have a building with an height equal
    1773. 

  1773.       real tw(2*ndm,nz_um)    ! Temperature in each layer of the wall [K]
    1774. 

  1774.       real tg(ndm,ng_u)         ! Temperature in each layer of the ground [K]
    1775. 

  1775.       real albg                 ! Albedo of the ground for the current urban class
    1776. 

  1776.       real albw                 ! Albedo of the wall for the current urban class
    1777. 

  1777.       real emg                  ! Emissivity of ground for the current urban class
    1778. 

  1778.       real emw                  ! Emissivity of wall for the current urban class
    1779. 

  1779.       real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
    1780. 

  1780.       real fsg(ndm,nurbm)             ! View factors from sky to ground
    1781. 

  1781.       real fsw(nz_um,ndm,nurbm)       ! View factors from sky to wall
    1782. 

  1782.       real fws(nz_um,ndm,nurbm)       ! View factors from wall to sky
    1783. 

  1783.       real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
    1784. 

  1784.       real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
    1785. 

  1785.       real ah                   ! Hour angle (it should come from the radiation routine)
    1786. 

  1786.       real zr                   ! zenith angle
    1787. 

  1787.       real deltar               ! Declination of the sun
    1788. 

  1788.       real rs                   ! solar radiation
    1789. 

  1789.       real rl                   ! downward flux of the longwave radiation
    1790. 

  1790. !
    1791. 

  1791. !New variables BEM
    1792. 

  1792. !
    1793. 

  1793.       real twlev(2*ndm,nz_um)         ! Window temperature in BEM [K]
    1794. 

  1794.       real pwin                       ! Coverage area fraction of windows in the walls of the buildings 
    1795. 

  1795.       real albwin                     ! Albedo of the windows for the current urban class
    1796. 

  1796.       real emwin                      ! Emissivity of the windows for the current urban class
    1797. 

  1797.       real alb_av                     ! Averaged albedo (window and wall)
    1798. 

  1798. ! ----------------------------------------------------------------------
    1799. 

  1799. ! OUTPUT:
    1800. 

  1800. ! ----------------------------------------------------------------------
    1801. 

  1801.       real rlg(ndm)             ! Long wave radiation at the ground
    1802. 

  1802.       real rlw(2*ndm,nz_um)     ! Long wave radiation at the walls
    1803. 

  1803.       real rsg(ndm)             ! Short wave radiation at the ground
    1804. 

  1804.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    1805. 

  1805. 

  1806. ! ----------------------------------------------------------------------
    1807. 

  1807. ! LOCAL:
    1808. 

  1808. ! ----------------------------------------------------------------------
    1809. 

  1809. 

  1810.       integer id,iz
    1811. 

  1811. 

  1812. !  Calculation of the shadow effects
    1813. 

  1813. 

  1814.       call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z,  &
    1815. 

  1815.                       rs,rsw,rsg)
    1816. 

  1816. 

  1817. ! Calculation of the reflection effects          
    1818. 

  1818.       do id=1,nd
    1819. 

  1819.          call long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev,      &
    1820. 

  1820.                       fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
    1821. 

  1821. 

  1822.          alb_av=pwin*albwin+(1.-pwin)*albw
    1823. 

  1823.          
    1824. 

  1824.          call short_rad(iurb,nz_u,id,alb_av,albg,fwg,fww,fgw,rsg,rsw,pb)
    1825. 

  1825.   
    1826. 

  1826.       enddo
    1827. 

  1827.       return
    1828. 

  1828.       end subroutine modif_rad
    1829. 

  1829. 

  1830. 

  1831. ! ===6=8===============================================================72  
    1832. 

  1832. ! ===6=8===============================================================72     
    1833. 

  1833. 

  1834.       subroutine surf_temp(nd,pr,dt,rl,rsg,rlg,  &
    1835. 

  1835.                            tg,alag,csg,emg,albg,ptg,sfg,gfg) 
    1836. 

  1836. 

  1837. ! ----------------------------------------------------------------------
    1838. 

  1838. ! Computation of the surface temperatures for walls, ground and roofs 
    1839. 

  1839. ! ----------------------------------------------------------------------
    1840. 

  1840. 

  1841.       implicit none
    1842. 

  1842.                   
    1843. 

  1843. ! ----------------------------------------------------------------------
    1844. 

  1844. ! INPUT:
    1845. 

  1845. ! ----------------------------------------------------------------------
    1846. 

  1846.       
    1847. 

  1847.       integer nd                ! Number of street direction for the current urban class
    1848. 

  1848.       real alag(ng_u)           ! Ground thermal diffusivity for the current urban class [m^2 s^-1] 
    1849. 

  1849. 

  1850.       real albg                 ! Albedo of the ground for the current urban class
    1851. 

  1851. 

  1852.       real csg(ng_u)            ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
    1853. 

  1853. 

  1854.       real dt                   ! Time step
    1855. 

  1855.       real emg                  ! Emissivity of ground for the current urban class
    1856. 

  1856. 

  1857.       real pr(nz_um)            ! Air pressure
    1858. 

  1858.       real rl                   ! Downward flux of the longwave radiation
    1859. 

  1859.       real rlg(ndm)             ! Long wave radiation at the ground
    1860. 

  1860.       real rsg(ndm)             ! Short wave radiation at the ground
    1861. 

  1861.       real sfg(ndm)             ! Sensible heat flux from ground (road)
    1862. 

  1862. 

  1863.       real gfg(ndm)             ! Heat flux transferred from the surface of the ground (road) toward the interior
    1864. 

  1864. 

  1865.       real tg(ndm,ng_u)         ! Temperature in each layer of the ground [K]
    1866. 

  1866. 

  1867.       
    1868. 

  1868. 

  1869. ! ----------------------------------------------------------------------
    1870. 

  1870. ! OUTPUT:
    1871. 

  1871. ! ----------------------------------------------------------------------
    1872. 

  1872.       real ptg(ndm)             ! Ground potential temperatures 
    1873. 

  1873. 

  1874. 

  1875. ! ----------------------------------------------------------------------
    1876. 

  1876. ! LOCAL:
    1877. 

  1877. ! ----------------------------------------------------------------------
    1878. 

  1878.       integer id,ig,ir,iw,iz
    1879. 

  1879. 

  1880.       real rtg(ndm)             ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
    1881. 

  1881. 

  1882.       real tg_tmp(ng_u)
    1883. 

  1883. 

  1884. 

  1885.       real dzg_u(ng_u)          ! Layer sizes in the ground
    1886. 

  1886.       
    1887. 

  1887.       data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
    1888. 

  1888. 

  1889. ! ----------------------------------------------------------------------
    1890. 

  1890. ! END VARIABLES DEFINITIONS
    1891. 

  1891. ! ----------------------------------------------------------------------
    1892. 

  1892. 

  1893.         
    1894. 

  1894.    
    1895. 

  1895.       do id=1,nd
    1896. 

  1896. 

  1897. !      Calculation for the ground surfaces
    1898. 

  1898.        do ig=1,ng_u
    1899. 

  1899.         tg_tmp(ig)=tg(id,ig)
    1900. 

  1900.        end do
    1901. 

  1901. !	        
    1902. 

  1902.        call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg,      &
    1903. 

  1903.                      rsg(id),rlg(id),pr(1),                    &
    1904. 

  1904.                      dt,emg,albg,                              &
    1905. 

  1905.                      rtg(id),sfg(id),gfg(id))    
    1906. 

  1906.        do ig=1,ng_u
    1907. 

  1907.         tg(id,ig)=tg_tmp(ig)
    1908. 

  1908.        end do
    1909. 

  1909. 	
    1910. 

  1910.       end do !id
    1911. 

  1911.       
    1912. 

  1912.       return
    1913. 

  1913.       end subroutine surf_temp
    1914. 

  1914.      
    1915. 

  1915. ! ===6=8===============================================================72     
    1916. 

  1916. ! ===6=8===============================================================72  
    1917. 

  1917. 

  1918.       subroutine buildings(iurb,nd,nz,z0,ua_u,va_u,pt_u,pt0_u,       &
    1919. 

  1919.                         ptg,ptr,da_u,ptw,ptwin,pwin,                 &
    1920. 

  1920.                         drst,uva_u,vva_u,uvb_u,vvb_u,                &
    1921. 

  1921.                         tva_u,tvb_u,evb_u,qvb_u,qhb_u,               &
    1922. 

  1922.                         uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr,   &
    1923. 

  1923.                         sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)                  
    1924. 

  1924. 

  1925. ! ----------------------------------------------------------------------
    1926. 

  1926. ! This routine computes the sources or sinks of the different quantities 
    1927. 

  1927. ! on the urban grid. The actual calculation is done in the subroutines 
    1928. 

  1928. ! called flux_wall, and flux_flat.
    1929. 

  1929. ! ----------------------------------------------------------------------
    1930. 

  1930. 

  1931.       implicit none
    1932. 

  1932. 

  1933.         
    1934. 

  1934. ! ----------------------------------------------------------------------
    1935. 

  1935. ! INPUT:
    1936. 

  1936. ! ----------------------------------------------------------------------
    1937. 

  1937.       integer nd                ! Number of street direction for the current urban class
    1938. 

  1938.       integer nz                ! number of vertical space steps
    1939. 

  1939.       real ua_u(nz_um)          ! Wind speed in the x direction on the urban grid
    1940. 

  1940.       real va_u(nz_um)          ! Wind speed in the y direction on the urban grid
    1941. 

  1941.       real da_u(nz_um)          ! air density on the urban grid
    1942. 

  1942.       real drst(ndm)            ! Street directions for the current urban class
    1943. 

  1943.       real dz
    1944. 

  1944.       real pt_u(nz_um)          ! Potential temperature on the urban grid
    1945. 

  1945.       real pt0_u(nz_um)         ! reference potential temperature on the urban grid
    1946. 

  1946.       real ptg(ndm)             ! Ground potential temperatures 
    1947. 

  1947.       real ptr(ndm,nz_um)       ! Roof potential temperatures 
    1948. 

  1948.       real ptw(2*ndm,nz_um,nbui_max)     ! Walls potential temperatures 
    1949. 

  1949.       real ss(nz_um)            ! probability to have a building with height h
    1950. 

  1950.       real pb(nz_um)
    1951. 

  1951.       real z0(ndm,nz_um)        ! Roughness lengths "profiles"
    1952. 

  1952.       real dt ! time step
    1953. 

  1953.       integer iurb              !Urban class
    1954. 

  1954. 

  1955. !
    1956. 

  1956. !New variables (BEM)
    1957. 

  1957. !
    1958. 

  1958.       real bs_u(ndm,nurbm)    ! Building width [m]
    1959. 

  1959.       real dz_u               ! Urban grid resolution
    1960. 

  1960.       real sflev(nz_um,nz_um)     ! sensible heat flux due to the air conditioning systems  [W]
    1961. 

  1961.       real lflev(nz_um,nz_um)     ! latent heat flux due to the air conditioning systems  [W]
    1962. 

  1962.       real sfvlev(nz_um,nz_um)    ! sensible heat flux due to ventilation [W]
    1963. 

  1963.       real lfvlev(nz_um,nz_um)    ! latent heat flux due to ventilation [W]
    1964. 

  1964.       real qvb_u(2*ndm,nz_um)
    1965. 

  1965.       real qhb_u(ndm,nz_um)
    1966. 

  1966.       real ptwin(2*ndm,nz_um,nbui_max)  ! window potential temperature
    1967. 

  1967.       real pwin
    1968. 

  1968.       real tvb_ac(2*ndm,nz_um)
    1969. 

  1969. ! ----------------------------------------------------------------------
    1970. 

  1970. ! OUTPUT:
    1971. 

  1971. ! ----------------------------------------------------------------------
    1972. 

  1972. ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
    1973. 

  1973. ! vertical surfaces (walls) and horizontal surfaces (roofs and street)
    1974. 

  1974. ! The fluxes can be computed as follow: Fluxes of X = A*X + B
    1975. 

  1975. !  Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
    1976. 

  1976. 

  1977.       real uhb_u(ndm,nz_um)   ! U (wind component) Horizontal surfaces, B (explicit) term
    1978. 

  1978.       real uva_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, A (implicit) term
    1979. 

  1979.       real uvb_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, B (explicit) term
    1980. 

  1980.       real vhb_u(ndm,nz_um)   ! V (wind component) Horizontal surfaces, B (explicit) term
    1981. 

  1981.       real vva_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, A (implicit) term
    1982. 

  1982.       real vvb_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, B (explicit) term
    1983. 

  1983.       real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
    1984. 

  1984.       real tva_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, A (implicit) term
    1985. 

  1985.       real tvb_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, B (explicit) term
    1986. 

  1986.       real ehb_u(ndm,nz_um)   ! Energy (TKE)       Horizontal surfaces, B (explicit) term
    1987. 

  1987.       real evb_u(2*ndm,nz_um)   ! Energy (TKE)         Vertical surfaces, B (explicit) term
    1988. 

  1988.       real sfw(2*ndm,nz_um,nbui_max)   ! sensible heat flux from walls
    1989. 

  1989.       real sfwin(2*ndm,nz_um,nbui_max) ! sensible heat flux form windows
    1990. 

  1990.       real sfr(ndm,nz_um)           ! sensible heat flux from roof
    1991. 

  1991.       real sfg(ndm)                 ! sensible heat flux from street
    1992. 

  1992. 

  1993. ! ----------------------------------------------------------------------
    1994. 

  1994. ! LOCAL:
    1995. 

  1995. ! ----------------------------------------------------------------------
    1996. 

  1996.       real d_urb(nz_um)
    1997. 

  1997.       real uva_tmp
    1998. 

  1998.       real vva_tmp
    1999. 

  1999.       real uvb_tmp
    2000. 

  2000.       real vvb_tmp 
    2001. 

  2001.       real evb_tmp     
    2002. 

  2002.       integer nlev(nz_um)
    2003. 

  2003.       integer id,iz,ibui,nbui
    2004. 

  2004. 

  2005. ! ----------------------------------------------------------------------
    2006. 

  2006. ! END VARIABLES DEFINITIONS
    2007. 

  2007. ! ----------------------------------------------------------------------
    2008. 

  2008.       dz=dz_u
    2009. 

  2009.       ibui=0
    2010. 

  2010.       nbui=0
    2011. 

  2011.       nlev=0
    2012. 

  2012.       d_urb=0.
    2013. 

  2013.       
    2014. 

  2014.       uva_u=0.
    2015. 

  2015.       uvb_u=0.
    2016. 

  2016.       vhb_u=0.
    2017. 

  2017.       vva_u=0.
    2018. 

  2018.       vvb_u=0.
    2019. 

  2019.       thb_u=0.
    2020. 

  2020.       tva_u=0.
    2021. 

  2021.       tvb_u=0.
    2022. 

  2022.       tvb_ac=0.
    2023. 

  2023.       ehb_u=0.
    2024. 

  2024.       evb_u=0.
    2025. 

  2025.       qvb_u=0.
    2026. 

  2026.       qhb_u=0.
    2027. 

  2027. 

  2028.       do iz=1,nz_um		   
    2029. 

  2029.          if(ss(iz).gt.0) then		
    2030. 

  2030.            ibui=ibui+1
    2031. 

  2031.            d_urb(ibui)=ss(iz)
    2032. 

  2032.            nlev(ibui)=iz-1
    2033. 

  2033.            nbui=ibui		               
    2034. 

  2034.          endif
    2035. 

  2035.       enddo
    2036. 

  2036.       if (nbui.gt.nbui_max) then
    2037. 

  2037.          write(*,*) 'nbui_max must be increased to',nbui
    2038. 

  2038.          stop
    2039. 

  2039.       endif
    2040. 

  2040.       do id=1,nd
    2041. 

  2041. 

  2042. !        Calculation at the ground surfaces
    2043. 

  2043. 

  2044.          call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1),  &
    2045. 

  2045.                        ptg(id),uhb_u(id,1),                            & 
    2046. 

  2046.                        vhb_u(id,1),sfg(id),ehb_u(id,1),da_u(1))          
    2047. 

  2047.          thb_u(id,1)=- sfg(id)/(da_u(1)*cp_u)  
    2048. 

  2048.               
    2049. 

  2049. !        Calculation at the roof surfaces    
    2050. 

  2050. 

  2051.          do iz=2,nz
    2052. 

  2052.             if(ss(iz).gt.0)then
    2053. 

  2053.                call flux_flat(dz,z0(id,iz),ua_u(iz),                  &              
    2054. 

  2054.                        va_u(iz),pt_u(iz),pt0_u(iz),                   &   
    2055. 

  2055.                        ptr(id,iz),uhb_u(id,iz),                       &   
    2056. 

  2056.                        vhb_u(id,iz),sfr(id,iz),ehb_u(id,iz),da_u(iz)) 
    2057. 

  2057.                thb_u(id,iz)=- sfr(id,iz)/(da_u(iz)*cp_u) 
    2058. 

  2058.             else
    2059. 

  2059.                uhb_u(id,iz) = 0.0
    2060. 

  2060.                vhb_u(id,iz) = 0.0
    2061. 

  2061.                thb_u(id,iz) = 0.0
    2062. 

  2062.                ehb_u(id,iz) = 0.0
    2063. 

  2063.             endif
    2064. 

  2064.          end do
    2065. 

  2065. 

  2066. !        Calculation at the wall surfaces        
    2067. 

  2067.  
    2068. 

  2068.          do ibui=1,nbui
    2069. 

  2069.          do iz=1,nlev(ibui)  
    2070. 

  2070.                    
    2071. 

  2071.             call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),             &  
    2072. 

  2072.                         ptw(2*id-1,iz,ibui),ptwin(2*id-1,iz,ibui),          &   
    2073. 

  2073.                         uva_tmp,vva_tmp,                                    &   
    2074. 

  2074.                         uvb_tmp,vvb_tmp,                                    &   
    2075. 

  2075.                         sfw(2*id-1,iz,ibui),sfwin(2*id-1,iz,ibui),          &   
    2076. 

  2076.                         evb_tmp,drst(id),dt)      
    2077. 

  2077.    
    2078. 

  2078.             if (pb(iz+1).gt.0.) then
    2079. 

  2079. 

  2080.                     uva_u(2*id-1,iz)=uva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
    2081. 

  2081.                     vva_u(2*id-1,iz)=vva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
    2082. 

  2082.                     uvb_u(2*id-1,iz)=uvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
    2083. 

  2083.                     vvb_u(2*id-1,iz)=vvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
    2084. 

  2084.                     evb_u(2*id-1,iz)=evb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
    2085. 

  2085.                     tvb_u(2*id-1,iz)=tvb_u(2*id-1,iz)-(d_urb(ibui)/pb(iz+1))*                       &
    2086. 

  2086.                                     (sfw(2*id-1,iz,ibui)*(1.-pwin)+sfwin(2*id-1,iz,ibui)*pwin)/     &
    2087. 

  2087.                                     da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
    2088. 

  2088.                                     (dz*bs_u(id,iurb))/da_u(iz)/cp_u
    2089. 

  2089.                     tvb_ac(2*id-1,iz)=tvb_ac(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
    2090. 

  2090.                                     (dz*bs_u(id,iurb))/da_u(iz)/cp_u
    2091. 

  2091.                     qvb_u(2*id-1,iz)=qvb_u(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
    2092. 

  2092.                                     (dz*bs_u(id,iurb))/da_u(iz)/latent
    2093. 

  2093.                                      
    2094. 

  2094.             endif
    2095. 

  2095. 

  2096.             call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),    &   
    2097. 

  2097.                         ptw(2*id,iz,ibui),ptwin(2*id,iz,ibui),     &    
    2098. 

  2098.                         uva_tmp,vva_tmp,                           &    
    2099. 

  2099.                         uvb_tmp,vvb_tmp,                           &    
    2100. 

  2100.                         sfw(2*id,iz,ibui),sfwin(2*id,iz,ibui),     &   
    2101. 

  2101.                         evb_tmp,drst(id),dt) 
    2102. 

  2102. 

  2103.             if (pb(iz+1).gt.0.) then
    2104. 

  2104. 

  2105.                     uva_u(2*id,iz)=uva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
    2106. 

  2106.                     vva_u(2*id,iz)=vva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
    2107. 

  2107.                     uvb_u(2*id,iz)=uvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
    2108. 

  2108.                     vvb_u(2*id,iz)=vvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
    2109. 

  2109.                     evb_u(2*id,iz)=evb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
    2110. 

  2110.                     tvb_u(2*id,iz)=tvb_u(2*id,iz)-(d_urb(ibui)/pb(iz+1))*                    &
    2111. 

  2111.                                     (sfw(2*id,iz,ibui)*(1.-pwin)+sfwin(2*id,iz,ibui)*pwin)/  &
    2112. 

  2112.                                      da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
    2113. 

  2113.                                     (dz*bs_u(id,iurb))/da_u(iz)/cp_u
    2114. 

  2114.                     tvb_ac(2*id,iz)=tvb_ac(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
    2115. 

  2115.                                     (dz*bs_u(id,iurb))/da_u(iz)/cp_u
    2116. 

  2116.                     qvb_u(2*id,iz)=qvb_u(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
    2117. 

  2117.                                     (dz*bs_u(id,iurb))/da_u(iz)/latent
    2118. 

  2118. 

  2119.             endif
    2120. 

  2120. !
    2121. 

  2121.           enddo !iz
    2122. 

  2122.          enddo !ibui
    2123. 

  2123.          
    2124. 

  2124.       end do !id
    2125. 

  2125.                 
    2126. 

  2126.       return
    2127. 

  2127.       end subroutine buildings
    2128. 

  2128.       
    2129. 

  2129. 

  2130. ! ===6=8===============================================================72
    2131. 

  2131. ! ===6=8===============================================================72
    2132. 

  2132. 

  2133.         subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl,    &
    2134. 

  2134.                              uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &       
    2135. 

  2135.                              uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u,       &      
    2136. 

  2136.                              a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)           
    2137. 

  2137. 

  2138. ! ----------------------------------------------------------------------
    2139. 

  2139. !  This routine interpolates the parameters from the "urban grid" to the
    2140. 

  2140. !  "mesoscale grid".
    2141. 

  2141. !  See p300-301 Appendix B.2 of the BLM paper.  
    2142. 

  2142. ! ----------------------------------------------------------------------
    2143. 

  2143. 

  2144.       implicit none
    2145. 

  2145. 

  2146. ! ----------------------------------------------------------------------
    2147. 

  2147. ! INPUT:
    2148. 

  2148. ! ----------------------------------------------------------------------
    2149. 

  2149. ! Data relative to the "mesoscale grid"
    2150. 

  2150.       integer kms,kme,kts,kte               
    2151. 

  2151.       real z(kms:kme)              ! Altitude above the ground of the cell interface
    2152. 

  2152.       real dz(kms:kme)               ! Vertical space steps
    2153. 

  2153. 

  2154. ! Data relative to the "uban grid"
    2155. 

  2155.       integer nz_u              ! Number of layer in the urban grid
    2156. 

  2156.       integer nd                ! Number of street direction for the current urban class
    2157. 

  2157.       real bs(ndm)              ! Building widths of the current urban class
    2158. 

  2158.       real ws(ndm)              ! Street widths of the current urban class
    2159. 

  2159.       real z_u(nz_um)         ! Height of the urban grid levels
    2160. 

  2160.       real pb(nz_um)          ! Probability to have a building with an height equal
    2161. 

  2161.       real ss(nz_um)          ! Probability to have a building with height h
    2162. 

  2162.       real uhb_u(ndm,nz_um)   ! U (x-wind component) Horizontal surfaces, B (explicit) term
    2163. 

  2163.       real uva_u(2*ndm,nz_um)   ! U (x-wind component) Vertical surfaces, A (implicit) term
    2164. 

  2164.       real uvb_u(2*ndm,nz_um)   ! U (x-wind component) Vertical surfaces, B (explicit) term
    2165. 

  2165.       real vhb_u(ndm,nz_um)   ! V (y-wind component) Horizontal surfaces, B (explicit) term
    2166. 

  2166.       real vva_u(2*ndm,nz_um)   ! V (y-wind component) Vertical surfaces, A (implicit) term
    2167. 

  2167.       real vvb_u(2*ndm,nz_um)   ! V (y-wind component) Vertical surfaces, B (explicit) term
    2168. 

  2168.       real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
    2169. 

  2169.       real tva_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, A (implicit) term
    2170. 

  2170.       real tvb_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, B (explicit) term
    2171. 

  2171.       real tvb_ac(2*ndm,nz_um)
    2172. 

  2172.       real ehb_u(ndm,nz_um)   ! Energy (TKE)       Horizontal surfaces, B (explicit) term
    2173. 

  2173.       real evb_u(2*ndm,nz_um)   ! Energy (TKE)         Vertical surfaces, B (explicit) term
    2174. 

  2174. !
    2175. 

  2175. !New variables for BEM
    2176. 

  2176. !
    2177. 

  2177.       real qhb_u(ndm,nz_um)
    2178. 

  2178.       real qvb_u(2*ndm,nz_um)
    2179. 

  2179.      
    2180. 

  2180. 

  2181. ! ----------------------------------------------------------------------
    2182. 

  2182. ! OUTPUT:
    2183. 

  2183. ! ----------------------------------------------------------------------
    2184. 

  2184. ! Data relative to the "mesoscale grid"
    2185. 

  2185.       real sf(kms:kme)             ! Surface of the "mesoscale grid" cells taking into account the buildings
    2186. 

  2186.       real vl(kms:kme)               ! Volume of the "mesoscale grid" cells taking into account the buildings
    2187. 

  2187.       real a_u(kms:kme)              ! Implicit component of the momentum sources or sinks in the X-direction
    2188. 

  2188.       real a_v(kms:kme)              ! Implicit component of the momentum sources or sinks in the Y-direction
    2189. 

  2189.       real a_t(kms:kme)              ! Implicit component of the heat sources or sinks
    2190. 

  2190.       real a_e(kms:kme)              ! Implicit component of the TKE sources or sinks
    2191. 

  2191.       real b_u(kms:kme)              ! Explicit component of the momentum sources or sinks in the X-direction
    2192. 

  2192.       real b_v(kms:kme)              ! Explicit component of the momentum sources or sinks in the Y-direction
    2193. 

  2193.       real b_t(kms:kme)              ! Explicit component of the heat sources or sinks
    2194. 

  2194.       real b_ac(kms:kme)
    2195. 

  2195.       real b_e(kms:kme)              ! Explicit component of the TKE sources or sinks
    2196. 

  2196.       real b_q(kms:kme)               ! Explicit component of the humidity sources or sinks
    2197. 

  2197. ! ----------------------------------------------------------------------
    2198. 

  2198. ! LOCAL:
    2199. 

  2199. ! ----------------------------------------------------------------------
    2200. 

  2200.       real dzz
    2201. 

  2201.       real fact
    2202. 

  2202.       integer id,iz,iz_u
    2203. 

  2203.       real se,sr,st,su,sv,sq
    2204. 

  2204.       real uet(kms:kme)                ! Contribution to TKE due to walls
    2205. 

  2205.       real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb,vqb,vtb_ac
    2206. 

  2206. 

  2207. 

  2208. ! ----------------------------------------------------------------------
    2209. 

  2209. ! END VARIABLES DEFINITIONS
    2210. 

  2210. ! ---------------------------------------------------------------------- 
    2211. 

  2211. 

  2212. ! initialisation
    2213. 

  2213. 

  2214.       do iz=kts,kte
    2215. 

  2215.          a_u(iz)=0.
    2216. 

  2216.          a_v(iz)=0.
    2217. 

  2217.          a_t(iz)=0.
    2218. 

  2218.          a_e(iz)=0.
    2219. 

  2219.          b_u(iz)=0.
    2220. 

  2220.          b_v(iz)=0.
    2221. 

  2221.          b_e(iz)=0.
    2222. 

  2222.          b_t(iz)=0.
    2223. 

  2223.          b_ac(iz)=0.
    2224. 

  2224.          uet(iz)=0.
    2225. 

  2225.          b_q(iz)=0.
    2226. 

  2226.       end do
    2227. 

  2227.             
    2228. 

  2228. ! horizontal surfaces
    2229. 

  2229.       do iz=kts,kte
    2230. 

  2230.          sf(iz)=0.
    2231. 

  2231.          vl(iz)=0.
    2232. 

  2232.       enddo
    2233. 

  2233.       sf(kte+1)=0. 
    2234. 

  2234.       
    2235. 

  2235.       do id=1,nd      
    2236. 

  2236.          do iz=kts+1,kte+1
    2237. 

  2237.             sr=0.
    2238. 

  2238.             do iz_u=2,nz_u
    2239. 

  2239.                if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
    2240. 

  2240.                   sr=pb(iz_u)
    2241. 

  2241.                endif
    2242. 

  2242.             enddo
    2243. 

  2243.             sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
    2244. 

  2244.          enddo
    2245. 

  2245.       enddo
    2246. 

  2246. 

  2247. ! volume      
    2248. 

  2248.       do iz=kts,kte
    2249. 

  2249.          do id=1,nd
    2250. 

  2250.             vtot=0.
    2251. 

  2251.             do iz_u=1,nz_u
    2252. 

  2252.                dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
    2253. 

  2253.                vtot=vtot+pb(iz_u+1)*dzz
    2254. 

  2254.             enddo
    2255. 

  2255.             vtot=vtot/(z(iz+1)-z(iz))
    2256. 

  2256.             vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
    2257. 

  2257.          enddo
    2258. 

  2258.       enddo
    2259. 

  2259.       
    2260. 

  2260. ! horizontal surface impact  
    2261. 

  2261. 

  2262.       do id=1,nd
    2263. 

  2263.       
    2264. 

  2264.          fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
    2265. 

  2265.          b_t(kts)=b_t(kts)+thb_u(id,1)*fact
    2266. 

  2266.          b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
    2267. 

  2267.          b_v(kts)=b_v(kts)+vhb_u(id,1)*fact 
    2268. 

  2268.          b_e(kts)=b_e(kts)+ehb_u(id,1)*fact*(z_u(2)-z_u(1))
    2269. 

  2269.          b_q(kts)=b_q(kts)+qhb_u(id,1)*fact         
    2270. 

  2270. 

  2271.          do iz=kts,kte
    2272. 

  2272.             st=0.
    2273. 

  2273.             su=0.
    2274. 

  2274.             sv=0.
    2275. 

  2275.             se=0.
    2276. 

  2276.             sq=0.
    2277. 

  2277.             do iz_u=2,nz_u
    2278. 

  2278.                if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
    2279. 

  2279.                   st=st+ss(iz_u)*thb_u(id,iz_u)
    2280. 

  2280.                   su=su+ss(iz_u)*uhb_u(id,iz_u)
    2281. 

  2281.                   sv=sv+ss(iz_u)*vhb_u(id,iz_u)          
    2282. 

  2282.                   se=se+ss(iz_u)*ehb_u(id,iz_u)*(z_u(iz_u+1)-z_u(iz_u))
    2283. 

  2283.                   sq=sq+ss(iz_u)*qhb_u(id,iz_u)
    2284. 

  2284.                endif
    2285. 

  2285.             enddo
    2286. 

  2286.       
    2287. 

  2287.             fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
    2288. 

  2288.             b_t(iz)=b_t(iz)+st*fact
    2289. 

  2289.             b_u(iz)=b_u(iz)+su*fact
    2290. 

  2290.             b_v(iz)=b_v(iz)+sv*fact
    2291. 

  2291.             b_e(iz)=b_e(iz)+se*fact
    2292. 

  2292.             b_q(iz)=b_q(iz)+sq*fact
    2293. 

  2293.          enddo
    2294. 

  2294.       enddo              
    2295. 

  2295. 

  2296. ! vertical surface impact
    2297. 

  2297.            
    2298. 

  2298.       do iz=kts,kte 
    2299. 

  2299.          uet(iz)=0.
    2300. 

  2300.          do id=1,nd              
    2301. 

  2301.             vtb=0.
    2302. 

  2302.             vtb_ac=0.
    2303. 

  2303.             vta=0.
    2304. 

  2304.             vua=0.
    2305. 

  2305.             vub=0.
    2306. 

  2306.             vva=0.
    2307. 

  2307.             vvb=0.
    2308. 

  2308.             veb=0.
    2309. 

  2309. 	    vte=0.
    2310. 

  2310.             vqb=0.
    2311. 

  2311.             do iz_u=1,nz_u
    2312. 

  2312.                dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
    2313. 

  2313.                fact=dzz/(ws(id)+bs(id))
    2314. 

  2314.                vtb=vtb+pb(iz_u+1)*                                  &        
    2315. 

  2315.                     (tvb_u(2*id-1,iz_u)+tvb_u(2*id,iz_u))*fact 
    2316. 

  2316.                vtb_ac=vtb_ac+pb(iz_u+1)*                            &        
    2317. 

  2317.                     (tvb_ac(2*id-1,iz_u)+tvb_ac(2*id,iz_u))*fact     
    2318. 

  2318.                vta=vta+pb(iz_u+1)*                                  &        
    2319. 

  2319.                    (tva_u(2*id-1,iz_u)+tva_u(2*id,iz_u))*fact
    2320. 

  2320.                vua=vua+pb(iz_u+1)*                                  &        
    2321. 

  2321.                     (uva_u(2*id-1,iz_u)+uva_u(2*id,iz_u))*fact
    2322. 

  2322.                vva=vva+pb(iz_u+1)*                                  &        
    2323. 

  2323.                     (vva_u(2*id-1,iz_u)+vva_u(2*id,iz_u))*fact
    2324. 

  2324.                vub=vub+pb(iz_u+1)*                                  &        
    2325. 

  2325.                     (uvb_u(2*id-1,iz_u)+uvb_u(2*id,iz_u))*fact
    2326. 

  2326.                vvb=vvb+pb(iz_u+1)*                                  &        
    2327. 

  2327.                     (vvb_u(2*id-1,iz_u)+vvb_u(2*id,iz_u))*fact
    2328. 

  2328.                veb=veb+pb(iz_u+1)*                                  &        
    2329. 

  2329.                     (evb_u(2*id-1,iz_u)+evb_u(2*id,iz_u))*fact
    2330. 

  2330.                vqb=vqb+pb(iz_u+1)*                                  &        
    2331. 

  2331.                     (qvb_u(2*id-1,iz_u)+qvb_u(2*id,iz_u))*fact   
    2332. 

  2332.             enddo
    2333. 

  2333.            
    2334. 

  2334.             fact=1./vl(iz)/dz(iz)/nd
    2335. 

  2335.             b_t(iz)=b_t(iz)+vtb*fact
    2336. 

  2336.             b_ac(iz)=b_ac(iz)+vtb_ac*fact
    2337. 

  2337.             a_t(iz)=a_t(iz)+vta*fact
    2338. 

  2338.             a_u(iz)=a_u(iz)+vua*fact
    2339. 

  2339.             a_v(iz)=a_v(iz)+vva*fact
    2340. 

  2340.             b_u(iz)=b_u(iz)+vub*fact
    2341. 

  2341.             b_v(iz)=b_v(iz)+vvb*fact
    2342. 

  2342.             b_e(iz)=b_e(iz)+veb*fact
    2343. 

  2343.             uet(iz)=uet(iz)+vte*fact
    2344. 

  2344.             b_q(iz)=b_q(iz)+vqb*fact
    2345. 

  2345.          enddo              
    2346. 

  2346.       enddo
    2347. 

  2347.       
    2348. 

  2348. 

  2349.       
    2350. 

  2350.       return
    2351. 

  2351.       end subroutine urban_meso
    2352. 

  2352. 

  2353. ! ===6=8===============================================================72 
    2354. 

  2354. ! ===6=8===============================================================72 
    2355. 

  2355. 

  2356.       subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
    2357. 

  2357.                              dlg,dl_u)
    2358. 

  2358. 

  2359. ! ----------------------------------------------------------------------     
    2360. 

  2360. !    Calculation of the length scales
    2361. 

  2361. !    See p272-274 formula (22) and (24) of the BLM paper    
    2362. 

  2362. ! ----------------------------------------------------------------------     
    2363. 

  2363.      
    2364. 

  2364.       implicit none
    2365. 

  2365. 

  2366. 

  2367. ! ----------------------------------------------------------------------     
    2368. 

  2368. ! INPUT:
    2369. 

  2369. ! ----------------------------------------------------------------------     
    2370. 

  2370.       integer kms,kme,kts,kte                
    2371. 

  2371.       real z(kms:kme)              ! Altitude above the ground of the cell interface
    2372. 

  2372.       integer nd                ! Number of street direction for the current urban class
    2373. 

  2373.       integer nz_u              ! Number of levels in the "urban grid"
    2374. 

  2374.       real z_u(nz_um)         ! Height of the urban grid levels
    2375. 

  2375.       real bs(ndm)              ! Building widths of the current urban class
    2376. 

  2376.       real ss(nz_um)          ! Probability to have a building with height h
    2377. 

  2377.       real ws(ndm)              ! Street widths of the current urban class
    2378. 

  2378. 

  2379. 

  2380. ! ----------------------------------------------------------------------     
    2381. 

  2381. ! OUTPUT:
    2382. 

  2382. ! ----------------------------------------------------------------------     
    2383. 

  2383.       real dlg(kms:kme)              ! Height above ground (L_ground in formula (24) of the BLM paper). 
    2384. 

  2384.       real dl_u(kms:kme)             ! Length scale (lb in formula (22) ofthe BLM paper).
    2385. 

  2385. 

  2386. ! ----------------------------------------------------------------------
    2387. 

  2387. ! LOCAL:
    2388. 

  2388. ! ----------------------------------------------------------------------
    2389. 

  2389.       real dlgtmp
    2390. 

  2390.       integer id,iz,iz_u
    2391. 

  2391.       real sftot
    2392. 

  2392.       real ulu,ssl
    2393. 

  2393. 

  2394. ! ----------------------------------------------------------------------
    2395. 

  2395. ! END VARIABLES DEFINITIONS
    2396. 

  2396. ! ----------------------------------------------------------------------
    2397. 

  2397.    
    2398. 

  2398.         do iz=kts,kte
    2399. 

  2399.          ulu=0.
    2400. 

  2400.          ssl=0.
    2401. 

  2401.          do id=1,nd        
    2402. 

  2402.           do iz_u=2,nz_u
    2403. 

  2403.            if(z_u(iz_u).gt.z(iz))then
    2404. 

  2404.             ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
    2405. 

  2405.             ssl=ssl+ss(iz_u)/nd
    2406. 

  2406.            endif
    2407. 

  2407.           enddo
    2408. 

  2408.          enddo
    2409. 

  2409. 

  2410.         if(ulu.ne.0)then
    2411. 

  2411.           dl_u(iz)=ssl/ulu
    2412. 

  2412.          else
    2413. 

  2413.           dl_u(iz)=0.
    2414. 

  2414.          endif
    2415. 

  2415.         enddo
    2416. 

  2416.        
    2417. 

  2417. 

  2418.         do iz=kts,kte
    2419. 

  2419.          dlg(iz)=0.
    2420. 

  2420.           do id=1,nd
    2421. 

  2421.            sftot=ws(id)  
    2422. 

  2422.            dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
    2423. 

  2423.            do iz_u=1,nz_u
    2424. 

  2424.             if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
    2425. 

  2425.              dlgtmp=dlgtmp+ss(iz_u)*bs(id)/                           &                
    2426. 

  2426.                     ((z(iz)+z(iz+1))/2.-z_u(iz_u))
    2427. 

  2427.              sftot=sftot+ss(iz_u)*bs(id)
    2428. 

  2428.             endif
    2429. 

  2429.            enddo
    2430. 

  2430.            dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
    2431. 

  2431.          enddo
    2432. 

  2432.          dlg(iz)=1./dlg(iz)
    2433. 

  2433.         enddo
    2434. 

  2434.         
    2435. 

  2435.        return
    2436. 

  2436.        end subroutine interp_length
    2437. 

  2437. 

  2438. ! ===6=8===============================================================72
    2439. 

  2439. ! ===6=8===============================================================72   
    2440. 

  2440. 

  2441.       subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z,    &
    2442. 

  2442.                            rs,rsw,rsg)
    2443. 

  2443.         
    2444. 

  2444. ! ----------------------------------------------------------------------
    2445. 

  2445. !         Modification of short wave radiation to take into account
    2446. 

  2446. !         the shadow produced by the buildings
    2447. 

  2447. ! ----------------------------------------------------------------------
    2448. 

  2448. 

  2449.       implicit none
    2450. 

  2450.      
    2451. 

  2451. ! ----------------------------------------------------------------------
    2452. 

  2452. ! INPUT:
    2453. 

  2453. ! ----------------------------------------------------------------------
    2454. 

  2454.       integer nd                ! Number of street direction for the current urban class
    2455. 

  2455.       integer nz_u              ! number of vertical layers defined in the urban grid
    2456. 

  2456.       real ah                   ! Hour angle (it should come from the radiation routine)
    2457. 

  2457.       real deltar               ! Declination of the sun
    2458. 

  2458.       real drst(ndm)            ! street directions for the current urban class
    2459. 

  2459.       real rs                   ! solar radiation
    2460. 

  2460.       real ss(nz_um)          ! probability to have a building with height h
    2461. 

  2461.       real pb(nz_um)          ! Probability that a building has an height greater or equal to h
    2462. 

  2462.       real ws(ndm)              ! Street width of the current urban class
    2463. 

  2463.       real z(nz_um)           ! Height of the urban grid levels
    2464. 

  2464.       real zr                   ! zenith angle
    2465. 

  2465. 

  2466. ! ----------------------------------------------------------------------
    2467. 

  2467. ! OUTPUT:
    2468. 

  2468. ! ----------------------------------------------------------------------
    2469. 

  2469.       real rsg(ndm)             ! Short wave radiation at the ground
    2470. 

  2470.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    2471. 

  2471. 

  2472. ! ----------------------------------------------------------------------
    2473. 

  2473. ! LOCAL:
    2474. 

  2474. ! ----------------------------------------------------------------------
    2475. 

  2475.       integer id,iz,jz
    2476. 

  2476.       real aae,aaw,bbb,phix,rd,rtot,wsd
    2477. 

  2477.       
    2478. 

  2478. ! ----------------------------------------------------------------------
    2479. 

  2479. ! END VARIABLES DEFINITIONS
    2480. 

  2480. ! ----------------------------------------------------------------------
    2481. 

  2481. 

  2482.       if(rs.eq.0.or.sin(zr).eq.1)then
    2483. 

  2483.          do id=1,nd
    2484. 

  2484.             rsg(id)=0.
    2485. 

  2485.             do iz=1,nz_u
    2486. 

  2486.                rsw(2*id-1,iz)=0.
    2487. 

  2487.                rsw(2*id,iz)=0.
    2488. 

  2488.             enddo
    2489. 

  2489.          enddo
    2490. 

  2490.       else            
    2491. 

  2491. !test              
    2492. 

  2492.          if(abs(sin(zr)).gt.1.e-10)then
    2493. 

  2493.           if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
    2494. 

  2494.            bbb=pi/2.
    2495. 

  2495.           elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
    2496. 

  2496.            bbb=-pi/2.
    2497. 

  2497.           else
    2498. 

  2498.            bbb=asin(cos(deltar)*sin(ah)/sin(zr))
    2499. 

  2499.           endif
    2500. 

  2500.          else
    2501. 

  2501.           if(cos(deltar)*sin(ah).ge.0)then
    2502. 

  2502.            bbb=pi/2.
    2503. 

  2503.           elseif(cos(deltar)*sin(ah).lt.0)then
    2504. 

  2504.            bbb=-pi/2.
    2505. 

  2505.           endif
    2506. 

  2506.          endif
    2507. 

  2507. 

  2508.          phix=zr
    2509. 

  2509.            
    2510. 

  2510.          do id=1,nd
    2511. 

  2511.         
    2512. 

  2512.             rsg(id)=0.
    2513. 

  2513.            
    2514. 

  2514.             aae=bbb-drst(id)
    2515. 

  2515.             aaw=bbb-drst(id)+pi
    2516. 

  2516.                     
    2517. 

  2517.             do iz=1,nz_u
    2518. 

  2518.                rsw(2*id-1,iz)=0.
    2519. 

  2519.                rsw(2*id,iz)=0.          
    2520. 

  2520.             if (pb(iz+1).gt.0.) then          
    2521. 

  2521.                do jz=1,nz_u                    
    2522. 

  2522.                 if(abs(sin(aae)).gt.1.e-10)then
    2523. 

  2523.                   call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae,   &    
    2524. 

  2524.                       ws(id),rd)                  
    2525. 

  2525.                  rsw(2*id-1,iz)=rsw(2*id-1,iz)+rs*rd*ss(jz+1)/pb(iz+1)
    2526. 

  2526.                 endif
    2527. 

  2527.               
    2528. 

  2528.                 if(abs(sin(aaw)).gt.1.e-10)then
    2529. 

  2529.                   call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw,   &    
    2530. 

  2530.                       ws(id),rd)
    2531. 

  2531.                  rsw(2*id,iz)=rsw(2*id,iz)+rs*rd*ss(jz+1)/pb(iz+1)                  
    2532. 

  2532.                 endif
    2533. 

  2533.                enddo             
    2534. 

  2534.             endif  
    2535. 

  2535.             enddo
    2536. 

  2536.         if(abs(sin(aae)).gt.1.e-10)then
    2537. 

  2537.             wsd=abs(ws(id)/sin(aae))
    2538. 

  2538.               
    2539. 

  2539.             do jz=1,nz_u           
    2540. 

  2540.                rd=max(0.,wsd-z(jz+1)*tan(phix))
    2541. 

  2541.                rsg(id)=rsg(id)+rs*rd*ss(jz+1)/wsd          
    2542. 

  2542.             enddo
    2543. 

  2543.             rtot=0.
    2544. 

  2544.            
    2545. 

  2545.             do iz=1,nz_u
    2546. 

  2546.                rtot=rtot+(rsw(2*id,iz)+rsw(2*id-1,iz))*            &
    2547. 

  2547.                          (z(iz+1)-z(iz))
    2548. 

  2548.             enddo
    2549. 

  2549.             rtot=rtot+rsg(id)*ws(id)
    2550. 

  2550.         else
    2551. 

  2551.             rsg(id)=rs
    2552. 

  2552.         endif
    2553. 

  2553.             
    2554. 

  2554.          enddo
    2555. 

  2555.       endif
    2556. 

  2556.          
    2557. 

  2557.       return
    2558. 

  2558.       end subroutine shadow_mas
    2559. 

  2559.          
    2560. 

  2560. ! ===6=8===============================================================72     
    2561. 

  2561. ! ===6=8===============================================================72     
    2562. 

  2562. 

  2563.       subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
    2564. 

  2564. 

  2565. ! ----------------------------------------------------------------------
    2566. 

  2566. ! This routine computes the effects of a shadow induced by a building of 
    2567. 

  2567. ! height hu, on a portion of wall between z1 and z2. See equation A10, 
    2568. 

  2568. ! and correction described below formula A11, and figure A1. Basically rd
    2569. 

  2569. ! is the ratio between the horizontal surface illuminated and the portion
    2570. 

  2570. ! of wall. Referring to figure A1, multiplying radiation flux density on 
    2571. 

  2571. ! a horizontal surface (rs) by x1-x2 we have the radiation energy per 
    2572. 

  2572. ! unit time. Dividing this by z2-z1, we obtain the radiation flux 
    2573. 

  2573. ! density reaching the portion of the wall between z2 and z1 
    2574. 

  2574. ! (everything is assumed in 2D)
    2575. 

  2575. ! ----------------------------------------------------------------------
    2576. 

  2576. 

  2577.       implicit none
    2578. 

  2578.       
    2579. 

  2579. ! ----------------------------------------------------------------------
    2580. 

  2580. ! INPUT:
    2581. 

  2581. ! ----------------------------------------------------------------------
    2582. 

  2582.       real aa                   ! Angle between the sun direction and the face of the wall (A12)
    2583. 

  2583.       real hu                   ! Height of the building that generates the shadow
    2584. 

  2584.       real phix                 ! Solar zenith angle
    2585. 

  2585.       real ws                   ! Width of the street
    2586. 

  2586.       real z1                   ! Height of the level z(iz)
    2587. 

  2587.       real z2                   ! Height of the level z(iz+1)
    2588. 

  2588. 

  2589. ! ----------------------------------------------------------------------
    2590. 

  2590. ! OUTPUT:
    2591. 

  2591. ! ----------------------------------------------------------------------
    2592. 

  2592.       real rd                   ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A. 
    2593. 

  2593.                                 ! Multiplying rd by rs (radiation flux 
    2594. 

  2594.                                 ! density on a horizontal surface) gives 
    2595. 

  2595.                                 ! the radiation flux density on the 
    2596. 

  2596.                                 ! portion of wall between z1 and z2. 
    2597. 

  2597. ! ----------------------------------------------------------------------
    2598. 

  2598. ! LOCAL:
    2599. 

  2599. ! ----------------------------------------------------------------------
    2600. 

  2600.       real x1,x2                ! x1,x2 see Fig. A1.
    2601. 

  2601. 

  2602. ! ----------------------------------------------------------------------
    2603. 

  2603. ! END VARIABLES DEFINITIONS
    2604. 

  2604. ! ----------------------------------------------------------------------
    2605. 

  2605. 

  2606.       x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
    2607. 

  2607.       
    2608. 

  2608.       x2=max((hu-z2)*tan(phix),0.)
    2609. 

  2609. 

  2610.       rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
    2611. 

  2611.       
    2612. 

  2612.       return
    2613. 

  2613.       end subroutine shade_wall
    2614. 

  2614. 

  2615. ! ===6=8===============================================================72     
    2616. 

  2616. ! ===6=8===============================================================72     
    2617. 

  2617. 

  2618.       subroutine long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev,&
    2619. 

  2619.                          fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
    2620. 

  2620. 

  2621. ! ----------------------------------------------------------------------
    2622. 

  2622. ! This routine computes the effects of the reflections of long-wave 
    2623. 

  2623. ! radiation in the street canyon by solving the system 
    2624. 

  2624. ! of 2*nz_u+1 eqn. in 2*nz_u+1
    2625. 

  2625. ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
    2626. 

  2626. ! The system is solved by solving A X= B,
    2627. 

  2627. ! with A matrix B vector, and X solution. 
    2628. 

  2628. ! ----------------------------------------------------------------------
    2629. 

  2629. 

  2630.       implicit none
    2631. 

  2631. 

  2632.   
    2633. 

  2633.       
    2634. 

  2634. ! ----------------------------------------------------------------------
    2635. 

  2635. ! INPUT:
    2636. 

  2636. ! ----------------------------------------------------------------------
    2637. 

  2637.       real emg                        ! Emissivity of ground for the current urban class
    2638. 

  2638.       real emw                        ! Emissivity of wall for the current urban class
    2639. 

  2639.       real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
    2640. 

  2640.       real fsg(ndm,nurbm)             ! View factors from sky to ground
    2641. 

  2641.       real fsw(nz_um,ndm,nurbm)       ! View factors from sky to wall
    2642. 

  2642.       real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
    2643. 

  2643.       real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
    2644. 

  2644.       integer id                      ! Current street direction
    2645. 

  2645.       integer iurb                    ! Current urban class
    2646. 

  2646.       integer nz_u                    ! Number of layer in the urban grid
    2647. 

  2647.       real pb(nz_um)                  ! Probability to have a building with an height equal
    2648. 

  2648.       real rl                         ! Downward flux of the longwave radiation
    2649. 

  2649.       real tg(ndm,ng_u)               ! Temperature in each layer of the ground [K]
    2650. 

  2650.       real tw(2*ndm,nz_um)            ! Temperature in each layer of the wall [K]
    2651. 

  2651. !
    2652. 

  2652. !New Variables for BEM
    2653. 

  2653. !
    2654. 

  2654.       real twlev(2*ndm,nz_um)         ! Window temperature in BEM [K]
    2655. 

  2655.       real emwin                      ! Emissivity of windows
    2656. 

  2656.       real pwin                       ! Coverage area fraction of windows in the walls of the buildings (BEM)
    2657. 

  2657.       
    2658. 

  2658. 

  2659. ! ----------------------------------------------------------------------
    2660. 

  2660. ! OUTPUT:
    2661. 

  2661. ! ----------------------------------------------------------------------
    2662. 

  2662.       real rlg(ndm)                   ! Long wave radiation at the ground
    2663. 

  2663.       real rlw(2*ndm,nz_um)           ! Long wave radiation at the walls
    2664. 

  2664. 

  2665. ! ----------------------------------------------------------------------
    2666. 

  2666. ! LOCAL:
    2667. 

  2667. ! ----------------------------------------------------------------------
    2668. 

  2668.       integer i,j
    2669. 

  2669.       real aaa(2*nz_um+1,2*nz_um+1)   ! terms of the matrix
    2670. 

  2670.       real bbb(2*nz_um+1)             ! terms of the vector
    2671. 

  2671. 

  2672. ! ----------------------------------------------------------------------
    2673. 

  2673. ! END VARIABLES DEFINITIONS
    2674. 

  2674. ! ----------------------------------------------------------------------
    2675. 

  2675. 

  2676. 

  2677. ! west wall
    2678. 

  2678.        
    2679. 

  2679.       do i=1,nz_u
    2680. 

  2680.         
    2681. 

  2681.         do j=1,nz_u
    2682. 

  2682.          aaa(i,j)=0.
    2683. 

  2683.         enddo
    2684. 

  2684.         
    2685. 

  2685.         aaa(i,i)=1.        
    2686. 

  2686.        
    2687. 

  2687.         do j=nz_u+1,2*nz_u
    2688. 

  2688.          aaa(i,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
    2689. 

  2689.                   fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
    2690. 

  2690.         enddo
    2691. 

  2691.         
    2692. 

  2692.         aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
    2693. 

  2693.         
    2694. 

  2694.         bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg(id,ng_u)**4
    2695. 

  2695.         do j=1,nz_u
    2696. 

  2696.            bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i,id,iurb)* &
    2697. 

  2697.                  (emw*(1.-pwin)*tw(2*id,j)**4+emwin*pwin*twlev(2*id,j)**4)+ &
    2698. 

  2698.                  fww(j,i,id,iurb)*rl*(1.-pb(j+1))
    2699. 

  2699.         enddo
    2700. 

  2700.         
    2701. 

  2701.        enddo
    2702. 

  2702.        
    2703. 

  2703. ! east wall
    2704. 

  2704. 

  2705.        do i=1+nz_u,2*nz_u
    2706. 

  2706.         
    2707. 

  2707.         do j=1,nz_u
    2708. 

  2708.          aaa(i,j)=-(1.-emw*(1.-pwin)-emwin*pwin)*fww(j,i-nz_u,id,iurb)*pb(j+1)
    2709. 

  2709.         enddo
    2710. 

  2710.         
    2711. 

  2711.         do j=1+nz_u,2*nz_u
    2712. 

  2712.          aaa(i,j)=0.
    2713. 

  2713.         enddo
    2714. 

  2714.         
    2715. 

  2715.         aaa(i,i)=1.
    2716. 

  2716.         
    2717. 

  2717.         aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
    2718. 

  2718.         
    2719. 

  2719.         bbb(i)=fsw(i-nz_u,id,iurb)*rl+  &     
    2720. 

  2720.                emg*fgw(i-nz_u,id,iurb)*sigma*tg(id,ng_u)**4
    2721. 

  2721. 

  2722.         do j=1,nz_u
    2723. 

  2723.          bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i-nz_u,id,iurb)*  &   
    2724. 

  2724.                 (emw*(1.-pwin)*tw(2*id-1,j)**4+emwin*pwin*twlev(2*id-1,j)**4)+&   
    2725. 

  2725.                 fww(j,i-nz_u,id,iurb)*rl*(1.-pb(j+1))
    2726. 

  2726.         enddo
    2727. 

  2727.        
    2728. 

  2728.        enddo
    2729. 

  2729. 

  2730. ! ground
    2731. 

  2731.        do j=1,nz_u
    2732. 

  2732.         aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
    2733. 

  2733.                          fwg(j,id,iurb)*pb(j+1)
    2734. 

  2734.        enddo
    2735. 

  2735.        
    2736. 

  2736.        do j=nz_u+1,2*nz_u
    2737. 

  2737.         aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
    2738. 

  2738.                          fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
    2739. 

  2739.        enddo
    2740. 

  2740.        
    2741. 

  2741.        aaa(2*nz_u+1,2*nz_u+1)=1.
    2742. 

  2742.        
    2743. 

  2743.        bbb(2*nz_u+1)=fsg(id,iurb)*rl
    2744. 

  2744.        
    2745. 

  2745.        do i=1,nz_u
    2746. 

  2746.         bbb(2*nz_u+1)=bbb(2*nz_u+1)+sigma*fwg(i,id,iurb)*pb(i+1)*         &
    2747. 

  2747.                       (emw*(1.-pwin)*(tw(2*id-1,i)**4+tw(2*id,i)**4)+     &
    2748. 

  2748.                       emwin*pwin*(twlev(2*id-1,i)**4+twlev(2*id,i)**4))+  &
    2749. 

  2749.                       2.*fwg(i,id,iurb)*(1.-pb(i+1))*rl                  
    2750. 

  2750.        enddo
    2751. 

  2751.    
    2752. 

  2752. 

  2753.      
    2754. 

  2754.        call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
    2755. 

  2755. 

  2756.        do i=1,nz_u
    2757. 

  2757.         rlw(2*id-1,i)=bbb(i)
    2758. 

  2758.        enddo
    2759. 

  2759.        
    2760. 

  2760.        do i=nz_u+1,2*nz_u
    2761. 

  2761.         rlw(2*id,i-nz_u)=bbb(i)
    2762. 

  2762.        enddo
    2763. 

  2763.        
    2764. 

  2764.        rlg(id)=bbb(2*nz_u+1)
    2765. 

  2765.   
    2766. 

  2766.        return
    2767. 

  2767.        end subroutine long_rad
    2768. 

  2768.              
    2769. 

  2769. ! ===6=8===============================================================72
    2770. 

  2770. ! ===6=8===============================================================72
    2771. 

  2771. 

  2772.        subroutine short_rad(iurb,nz_u,id,albw,                        & 
    2773. 

  2773.                            albg,fwg,fww,fgw,rsg,rsw,pb)
    2774. 

  2774. 

  2775. ! ----------------------------------------------------------------------
    2776. 

  2776. ! This routine computes the effects of the reflections of short-wave 
    2777. 

  2777. ! (solar) radiation in the street canyon by solving the system 
    2778. 

  2778. ! of 2*nz_u+1 eqn. in 2*nz_u+1
    2779. 

  2779. ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
    2780. 

  2780. ! The system is solved by solving A X= B,
    2781. 

  2781. ! with A matrix B vector, and X solution. 
    2782. 

  2782. ! ----------------------------------------------------------------------
    2783. 

  2783. 

  2784.       implicit none
    2785. 

  2785. 

  2786.   
    2787. 

  2787.       
    2788. 

  2788. ! ----------------------------------------------------------------------
    2789. 

  2789. ! INPUT:
    2790. 

  2790. ! ----------------------------------------------------------------------
    2791. 

  2791.       real albg                 ! Albedo of the ground for the current urban class
    2792. 

  2792.       real albw                 ! Albedo of the wall for the current urban class
    2793. 

  2793.       real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
    2794. 

  2794.       real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
    2795. 

  2795.       real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
    2796. 

  2796.       integer id                ! current street direction 
    2797. 

  2797.       integer iurb              ! current urban class
    2798. 

  2798.       integer nz_u              ! Number of layer in the urban grid
    2799. 

  2799.       real pb(nz_um)          ! probability to have a building with an height equal
    2800. 

  2800. 

  2801. ! ----------------------------------------------------------------------
    2802. 

  2802. ! OUTPUT:
    2803. 

  2803. ! ----------------------------------------------------------------------
    2804. 

  2804.       real rsg(ndm)             ! Short wave radiation at the ground
    2805. 

  2805.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    2806. 

  2806. 

  2807. ! ----------------------------------------------------------------------
    2808. 

  2808. ! LOCAL:
    2809. 

  2809. ! ----------------------------------------------------------------------
    2810. 

  2810.       integer i,j
    2811. 

  2811.       real aaa(2*nz_um+1,2*nz_um+1)  ! terms of the matrix
    2812. 

  2812.       real bbb(2*nz_um+1)            ! terms of the vector
    2813. 

  2813. 

  2814. ! ----------------------------------------------------------------------
    2815. 

  2815. ! END VARIABLES DEFINITIONS
    2816. 

  2816. ! ----------------------------------------------------------------------
    2817. 

  2817. 

  2818.       
    2819. 

  2819. ! west wall
    2820. 

  2820.        
    2821. 

  2821.       do i=1,nz_u
    2822. 

  2822.          do j=1,nz_u
    2823. 

  2823.             aaa(i,j)=0.
    2824. 

  2824.          enddo
    2825. 

  2825.          
    2826. 

  2826.          aaa(i,i)=1.        
    2827. 

  2827.          
    2828. 

  2828.          do j=nz_u+1,2*nz_u
    2829. 

  2829.             aaa(i,j)=-albw*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
    2830. 

  2830.          enddo
    2831. 

  2831.          
    2832. 

  2832.          aaa(i,2*nz_u+1)=-albg*fgw(i,id,iurb)
    2833. 

  2833.          bbb(i)=rsw(2*id-1,i)
    2834. 

  2834.          
    2835. 

  2835.       enddo
    2836. 

  2836.        
    2837. 

  2837. ! east wall
    2838. 

  2838. 

  2839.       do i=1+nz_u,2*nz_u
    2840. 

  2840.          do j=1,nz_u
    2841. 

  2841.             aaa(i,j)=-albw*fww(j,i-nz_u,id,iurb)*pb(j+1)
    2842. 

  2842.          enddo
    2843. 

  2843.          
    2844. 

  2844.          do j=1+nz_u,2*nz_u
    2845. 

  2845.             aaa(i,j)=0.
    2846. 

  2846.          enddo
    2847. 

  2847.          
    2848. 

  2848.         aaa(i,i)=1.
    2849. 

  2849.         aaa(i,2*nz_u+1)=-albg*fgw(i-nz_u,id,iurb)
    2850. 

  2850.         bbb(i)=rsw(2*id,i-nz_u)
    2851. 

  2851.         
    2852. 

  2852.       enddo
    2853. 

  2853. 

  2854. ! ground
    2855. 

  2855. 

  2856.       do j=1,nz_u
    2857. 

  2857.          aaa(2*nz_u+1,j)=-albw*fwg(j,id,iurb)*pb(j+1)
    2858. 

  2858.       enddo
    2859. 

  2859.        
    2860. 

  2860.       do j=nz_u+1,2*nz_u
    2861. 

  2861.          aaa(2*nz_u+1,j)=-albw*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
    2862. 

  2862.       enddo
    2863. 

  2863.        
    2864. 

  2864.       aaa(2*nz_u+1,2*nz_u+1)=1.
    2865. 

  2865.       bbb(2*nz_u+1)=rsg(id)
    2866. 

  2866.        
    2867. 

  2867.       call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
    2868. 

  2868. 

  2869.       do i=1,nz_u
    2870. 

  2870.          rsw(2*id-1,i)=bbb(i)
    2871. 

  2871.       enddo
    2872. 

  2872.        
    2873. 

  2873.       do i=nz_u+1,2*nz_u
    2874. 

  2874.          rsw(2*id,i-nz_u)=bbb(i) 
    2875. 

  2875.       enddo
    2876. 

  2876.        
    2877. 

  2877.       rsg(id)=bbb(2*nz_u+1)
    2878. 

  2878.        
    2879. 

  2879.       return
    2880. 

  2880.       end subroutine short_rad
    2881. 

  2881.              
    2882. 

  2882. 

  2883. ! ===6=8===============================================================72     
    2884. 

  2884. ! ===6=8===============================================================72     
    2885. 

  2885.       
    2886. 

  2886.       subroutine gaussj(a,n,b,np)
    2887. 

  2887. 

  2888. ! ----------------------------------------------------------------------
    2889. 

  2889. ! This routine solve a linear system of n equations of the form
    2890. 

  2890. !              A X = B
    2891. 

  2891. !  where  A is a matrix a(i,j)
    2892. 

  2892. !         B a vector and X the solution
    2893. 

  2893. ! In output b is replaced by the solution     
    2894. 

  2894. ! ----------------------------------------------------------------------
    2895. 

  2895. 

  2896.       implicit none
    2897. 

  2897. 

  2898. ! ----------------------------------------------------------------------
    2899. 

  2899. ! INPUT:
    2900. 

  2900. ! ----------------------------------------------------------------------
    2901. 

  2901.       integer np
    2902. 

  2902.       real a(np,np)
    2903. 

  2903. 

  2904. ! ----------------------------------------------------------------------
    2905. 

  2905. ! OUTPUT:
    2906. 

  2906. ! ----------------------------------------------------------------------
    2907. 

  2907.       real b(np)
    2908. 

  2908. 

  2909. ! ----------------------------------------------------------------------
    2910. 

  2910. ! LOCAL:
    2911. 

  2911. ! ----------------------------------------------------------------------
    2912. 

  2912.       integer nmax
    2913. 

  2913.       parameter (nmax=150)
    2914. 

  2914. 

  2915.       real big,dum
    2916. 

  2916.       integer i,icol,irow
    2917. 

  2917.       integer j,k,l,ll,n
    2918. 

  2918.       integer ipiv(nmax)
    2919. 

  2919.       real pivinv
    2920. 

  2920. 

  2921. ! ----------------------------------------------------------------------
    2922. 

  2922. ! END VARIABLES DEFINITIONS
    2923. 

  2923. ! ----------------------------------------------------------------------
    2924. 

  2924.        
    2925. 

  2925.       do j=1,n
    2926. 

  2926.          ipiv(j)=0.
    2927. 

  2927.       enddo
    2928. 

  2928.        
    2929. 

  2929.       do i=1,n
    2930. 

  2930.          big=0.
    2931. 

  2931.          do j=1,n
    2932. 

  2932.             if(ipiv(j).ne.1)then
    2933. 

  2933.                do k=1,n
    2934. 

  2934.                   if(ipiv(k).eq.0)then
    2935. 

  2935.                      if(abs(a(j,k)).ge.big)then
    2936. 

  2936.                         big=abs(a(j,k))
    2937. 

  2937.                         irow=j
    2938. 

  2938.                         icol=k
    2939. 

  2939.                      endif
    2940. 

  2940.                   elseif(ipiv(k).gt.1)then
    2941. 

  2941.                      pause 'singular matrix in gaussj'
    2942. 

  2942.                   endif
    2943. 

  2943.                enddo
    2944. 

  2944.             endif
    2945. 

  2945.          enddo
    2946. 

  2946.          
    2947. 

  2947.          ipiv(icol)=ipiv(icol)+1
    2948. 

  2948.          
    2949. 

  2949.          if(irow.ne.icol)then
    2950. 

  2950.             do l=1,n
    2951. 

  2951.                dum=a(irow,l)
    2952. 

  2952.                a(irow,l)=a(icol,l)
    2953. 

  2953.                a(icol,l)=dum
    2954. 

  2954.             enddo
    2955. 

  2955.             
    2956. 

  2956.             dum=b(irow)
    2957. 

  2957.             b(irow)=b(icol)
    2958. 

  2958.             b(icol)=dum
    2959. 

  2959.           
    2960. 

  2960.          endif
    2961. 

  2961.          
    2962. 

  2962.          if(a(icol,icol).eq.0)pause 'singular matrix in gaussj'
    2963. 

  2963.          
    2964. 

  2964.          pivinv=1./a(icol,icol)
    2965. 

  2965.          a(icol,icol)=1
    2966. 

  2966.          
    2967. 

  2967.          do l=1,n
    2968. 

  2968.             a(icol,l)=a(icol,l)*pivinv
    2969. 

  2969.          enddo
    2970. 

  2970.          
    2971. 

  2971.          b(icol)=b(icol)*pivinv
    2972. 

  2972.          
    2973. 

  2973.          do ll=1,n
    2974. 

  2974.             if(ll.ne.icol)then
    2975. 

  2975.                dum=a(ll,icol)
    2976. 

  2976.                a(ll,icol)=0.
    2977. 

  2977.                do l=1,n
    2978. 

  2978.                   a(ll,l)=a(ll,l)-a(icol,l)*dum
    2979. 

  2979.                enddo
    2980. 

  2980.                
    2981. 

  2981.                b(ll)=b(ll)-b(icol)*dum
    2982. 

  2982.                
    2983. 

  2983.             endif
    2984. 

  2984.          enddo
    2985. 

  2985.       enddo
    2986. 

  2986.       
    2987. 

  2987.       return
    2988. 

  2988.       end subroutine gaussj
    2989. 

  2989.          
    2990. 

  2990. ! ===6=8===============================================================72     
    2991. 

  2991. ! ===6=8===============================================================72     
    2992. 

  2992.        
    2993. 

  2993.       subroutine soil_temp(nz,dz,temp,pt,ala,cs,                       &
    2994. 

  2994.                           rs,rl,press,dt,em,alb,rt,sf,gf)
    2995. 

  2995. 

  2996. ! ----------------------------------------------------------------------
    2997. 

  2997. ! This routine solves the Fourier diffusion equation for heat in 
    2998. 

  2998. ! the material (wall, roof, or ground). Resolution is done implicitely.
    2999. 

  2999. ! Boundary conditions are: 
    3000. 

  3000. ! - fixed temperature at the interior
    3001. 

  3001. ! - energy budget at the surface
    3002. 

  3002. ! ----------------------------------------------------------------------
    3003. 

  3003. 

  3004.       implicit none
    3005. 

  3005. 

  3006.      
    3007. 

  3007.                 
    3008. 

  3008. ! ----------------------------------------------------------------------
    3009. 

  3009. ! INPUT:
    3010. 

  3010. ! ----------------------------------------------------------------------
    3011. 

  3011.       integer nz                ! Number of layers
    3012. 

  3012.       real ala(nz)              ! Thermal diffusivity in each layers [m^2 s^-1] 
    3013. 

  3013.       real alb                  ! Albedo of the surface
    3014. 

  3014.       real cs(nz)               ! Specific heat of the material [J m^3 K^-1]
    3015. 

  3015.       real dt                   ! Time step
    3016. 

  3016.       real em                   ! Emissivity of the surface
    3017. 

  3017.       real press                ! Pressure at ground level
    3018. 

  3018.       real rl                   ! Downward flux of the longwave radiation
    3019. 

  3019.       real rs                   ! Solar radiation
    3020. 

  3020.       real sf                   ! Sensible heat flux at the surface
    3021. 

  3021.       real temp(nz)             ! Temperature in each layer [K]
    3022. 

  3022.       real dz(nz)               ! Layer sizes [m]
    3023. 

  3023. 

  3024. 

  3025. ! ----------------------------------------------------------------------
    3026. 

  3026. ! OUTPUT:
    3027. 

  3027. ! ----------------------------------------------------------------------
    3028. 

  3028.       real gf                   ! Heat flux transferred from the surface toward the interior
    3029. 

  3029.       real pt                   ! Potential temperature at the surface
    3030. 

  3030.       real rt                   ! Total radiation at the surface (solar+incoming long+outgoing long)
    3031. 

  3031. 

  3032. ! ----------------------------------------------------------------------
    3033. 

  3033. ! LOCAL:
    3034. 

  3034. ! ----------------------------------------------------------------------
    3035. 

  3035.       integer iz
    3036. 

  3036.       real a(nz,3)
    3037. 

  3037.       real alpha
    3038. 

  3038.       real c(nz)
    3039. 

  3039.       real cddz(nz+2)
    3040. 

  3040.       real tsig
    3041. 

  3041. 

  3042. ! ----------------------------------------------------------------------
    3043. 

  3043. ! END VARIABLES DEFINITIONS
    3044. 

  3044. ! ----------------------------------------------------------------------
    3045. 

  3045.        
    3046. 

  3046.       tsig=temp(nz)
    3047. 

  3047.       alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
    3048. 

  3048. ! Compute cddz=2*cd/dz  
    3049. 

  3049.         
    3050. 

  3050.       cddz(1)=ala(1)/dz(1)
    3051. 

  3051.       do iz=2,nz
    3052. 

  3052.          cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
    3053. 

  3053.       enddo
    3054. 

  3054.        
    3055. 

  3055.       a(1,1)=0.
    3056. 

  3056.       a(1,2)=1.
    3057. 

  3057.       a(1,3)=0.
    3058. 

  3058.       c(1)=temp(1)
    3059. 

  3059.           
    3060. 

  3060.       do iz=2,nz-1
    3061. 

  3061.          a(iz,1)=-cddz(iz)*dt/dz(iz)
    3062. 

  3062.          a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)          
    3063. 

  3063.          a(iz,3)=-cddz(iz+1)*dt/dz(iz)
    3064. 

  3064.          c(iz)=temp(iz)
    3065. 

  3065.       enddo          
    3066. 

  3066.                      
    3067. 

  3067.       a(nz,1)=-dt*cddz(nz)/dz(nz)
    3068. 

  3068.       a(nz,2)=1.+dt*cddz(nz)/dz(nz)
    3069. 

  3069.       a(nz,3)=0.
    3070. 

  3070.       c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
    3071. 

  3071. 

  3072.       
    3073. 

  3073.       call invert(nz,a,c,temp)
    3074. 

  3074. 

  3075.            
    3076. 

  3076.       pt=temp(nz)*(press/1.e+5)**(-rcp_u)
    3077. 

  3077. 

  3078.       rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)
    3079. 

  3079.                         
    3080. 

  3080.        gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf                                   
    3081. 

  3081.       return
    3082. 

  3082.       end subroutine soil_temp
    3083. 

  3083. 

  3084. ! ===6=8===============================================================72 
    3085. 

  3085. ! ===6=8===============================================================72 
    3086. 

  3086. 

  3087.       subroutine invert(n,a,c,x)
    3088. 

  3088. 

  3089. ! ----------------------------------------------------------------------
    3090. 

  3090. !        Inversion and resolution of a tridiagonal matrix
    3091. 

  3091. !                   A X = C
    3092. 

  3092. ! ----------------------------------------------------------------------
    3093. 

  3093. 

  3094.       implicit none
    3095. 

  3095.                 
    3096. 

  3096. ! ----------------------------------------------------------------------
    3097. 

  3097. ! INPUT:
    3098. 

  3098. ! ----------------------------------------------------------------------
    3099. 

  3099.        integer n
    3100. 

  3100.        real a(n,3)              !  a(*,1) lower diagonal (Ai,i-1)
    3101. 

  3101.                                 !  a(*,2) principal diagonal (Ai,i)
    3102. 

  3102.                                 !  a(*,3) upper diagonal (Ai,i+1)
    3103. 

  3103.        real c(n)
    3104. 

  3104. 

  3105. ! ----------------------------------------------------------------------
    3106. 

  3106. ! OUTPUT:
    3107. 

  3107. ! ----------------------------------------------------------------------
    3108. 

  3108.        real x(n)    
    3109. 

  3109. 

  3110. ! ----------------------------------------------------------------------
    3111. 

  3111. ! LOCAL:
    3112. 

  3112. ! ----------------------------------------------------------------------
    3113. 

  3113.        integer i
    3114. 

  3114. 

  3115. ! ----------------------------------------------------------------------
    3116. 

  3116. ! END VARIABLES DEFINITIONS
    3117. 

  3117. ! ----------------------------------------------------------------------
    3118. 

  3118.                      
    3119. 

  3119.        do i=n-1,1,-1                 
    3120. 

  3120.           c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
    3121. 

  3121.           a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
    3122. 

  3122.        enddo
    3123. 

  3123.        
    3124. 

  3124.        do i=2,n        
    3125. 

  3125.           c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
    3126. 

  3126.        enddo
    3127. 

  3127.         
    3128. 

  3128.        do i=1,n
    3129. 

  3129.           x(i)=c(i)/a(i,2)
    3130. 

  3130.        enddo
    3131. 

  3131. 

  3132.        return
    3133. 

  3133.        end subroutine invert
    3134. 

  3134.   
    3135. 

  3135. 

  3136. ! ===6=8===============================================================72  
    3137. 

  3137. ! ===6=8===============================================================72
    3138. 

  3138.   
    3139. 

  3139.       subroutine flux_wall(ua,va,pt,da,ptw,ptwin,uva,vva,uvb,vvb,  &
    3140. 

  3140.                            sfw,sfwin,evb,drst,dt)         
    3141. 

  3141.        
    3142. 

  3142. ! ----------------------------------------------------------------------
    3143. 

  3143. ! This routine computes the surface sources or sinks of momentum, tke,
    3144. 

  3144. ! and heat from vertical surfaces (walls).   
    3145. 

  3145. ! ----------------------------------------------------------------------
    3146. 

  3146.       implicit none   
    3147. 

  3147.          
    3148. 

  3148. ! INPUT:
    3149. 

  3149. ! -----
    3150. 

  3150.       real drst                 ! street directions for the current urban class
    3151. 

  3151.       real da                   ! air density
    3152. 

  3152.       real pt                   ! potential temperature
    3153. 

  3153.       real ptw                  ! Walls potential temperatures 
    3154. 

  3154.       real ptwin                ! Windows potential temperatures
    3155. 

  3155.       real ua                   ! wind speed
    3156. 

  3156.       real va                   ! wind speed
    3157. 

  3157.       real dt                   !time step
    3158. 

  3158. 

  3159. ! OUTPUT:
    3160. 

  3160. ! ------
    3161. 

  3161. ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
    3162. 

  3162. ! vertical surfaces (walls).
    3163. 

  3163. ! The fluxes can be computed as follow: Fluxes of X = A*X + B
    3164. 

  3164. ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
    3165. 

  3165. 

  3166.       real uva                  ! U (wind component)   Vertical surfaces, A (implicit) term
    3167. 

  3167.       real uvb                  ! U (wind component)   Vertical surfaces, B (explicit) term
    3168. 

  3168.       real vva                  ! V (wind component)   Vertical surfaces, A (implicit) term
    3169. 

  3169.       real vvb                  ! V (wind component)   Vertical surfaces, B (explicit) term
    3170. 

  3170. !     real tva                  ! Temperature          Vertical surfaces, A (implicit) term
    3171. 

  3171. !     real tvb                  ! Temperature          Vertical surfaces, B (explicit) term
    3172. 

  3172.       real evb                  ! Energy (TKE)         Vertical surfaces, B (explicit) term
    3173. 

  3173.       real sfw                  ! Surfaces fluxes from the walls
    3174. 

  3174.       real sfwin                ! Surfaces fluxes from the windows
    3175. 

  3175. 

  3176. ! LOCAL:
    3177. 

  3177. ! -----
    3178. 

  3178.       real hc
    3179. 

  3179.       real hcwin
    3180. 

  3180.       real u_ort
    3181. 

  3181.       real vett
    3182. 

  3182. 

  3183. 

  3184. ! -------------------------
    3185. 

  3185. ! END VARIABLES DEFINITIONS
    3186. 

  3186. ! -------------------------
    3187. 

  3187. 

  3188.       vett=(ua**2+va**2)**.5         
    3189. 

  3189.          
    3190. 

  3190.       u_ort=abs((cos(drst)*ua-sin(drst)*va))
    3191. 

  3191.        
    3192. 

  3192.       uva=-cdrag*u_ort/2.*cos(drst)*cos(drst)
    3193. 

  3193.       vva=-cdrag*u_ort/2.*sin(drst)*sin(drst)
    3194. 

  3194.          
    3195. 

  3195.       uvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*va
    3196. 

  3196.       vvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*ua         
    3197. 

  3197. 

  3198.       if (vett.lt.4.88) then   
    3199. 

  3199.          hc=5.678*(1.09+0.23*(vett/0.3048))  
    3200. 

  3200.       else
    3201. 

  3201.          hc=5.678*0.53*((vett/0.3048)**0.78)
    3202. 

  3202.       endif 
    3203. 

  3203. 

  3204.       if (hc.gt.da*cp_u/dt)then
    3205. 

  3205.          hc=da*cp_u/dt
    3206. 

  3206.       endif
    3207. 

  3207. 

  3208.        if (vett.lt.4.88) then
    3209. 

  3209.           hcwin=5.678*(0.99+0.21*(vett/0.3048))
    3210. 

  3210.        else
    3211. 

  3211.           hcwin=5.678*0.50*((vett/0.3048)**0.78)
    3212. 

  3212.        endif
    3213. 

  3213. 

  3214.        if (hcwin.gt.da*cp_u/dt) then
    3215. 

  3215.            hcwin=da*cp_u/dt
    3216. 

  3216.        endif
    3217. 

  3217.          
    3218. 

  3218. !         tvb=hc*ptw/da/cp_u
    3219. 

  3219. !         tva=-hc/da/cp_u
    3220. 

  3220. !!!!!!!!!!!!!!!!!!!!
    3221. 

  3221. ! explicit 
    3222. 

  3222. 

  3223.       sfw=hc*(pt-ptw)
    3224. 

  3224.       sfwin=hcwin*(pt-ptwin)  
    3225. 

  3225.        
    3226. 

  3226.          
    3227. 

  3227.       evb=cdrag*(abs(u_ort)**3.)/2.
    3228. 

  3228.               
    3229. 

  3229.       return
    3230. 

  3230.       end subroutine flux_wall
    3231. 

  3231.          
    3232. 

  3232. ! ===6=8===============================================================72
    3233. 

  3233. ! ===6=8===============================================================72
    3234. 

  3234. 

  3235.       subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg,                     &
    3236. 

  3236.                           uhb,vhb,sf,ehb,da)
    3237. 

  3237.                                 
    3238. 

  3238. ! ----------------------------------------------------------------------
    3239. 

  3239. !           Calculation of the flux at the ground 
    3240. 

  3240. !           Formulation of Louis (Louis, 1979)       
    3241. 

  3241. ! ----------------------------------------------------------------------
    3242. 

  3242. 

  3243.       implicit none
    3244. 

  3244. 

  3245.       real dz                   ! first vertical level
    3246. 

  3246.       real pt                   ! potential temperature
    3247. 

  3247.       real pt0                  ! reference potential temperature
    3248. 

  3248.       real ptg                  ! ground potential temperature
    3249. 

  3249.       real ua                   ! wind speed
    3250. 

  3250.       real va                   ! wind speed
    3251. 

  3251.       real z0                   ! Roughness length
    3252. 

  3252.       real da                   ! air density
    3253. 

  3253.       
    3254. 

  3254.      
    3255. 

  3255. 

  3256. ! ----------------------------------------------------------------------
    3257. 

  3257. ! OUTPUT:
    3258. 

  3258. ! ----------------------------------------------------------------------
    3259. 

  3259. ! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal 
    3260. 

  3260. !  surfaces (roofs and street)
    3261. 

  3261. ! The fluxes can be computed as follow: Fluxes of X = B
    3262. 

  3262. !  Example: Momentum fluxes on horizontal surfaces =  uhb_u
    3263. 

  3263.       real uhb                  ! U (wind component) Horizontal surfaces, B (explicit) term
    3264. 

  3264.       real vhb                  ! V (wind component) Horizontal surfaces, B (explicit) term
    3265. 

  3265. !     real thb                  ! Temperature        Horizontal surfaces, B (explicit) term
    3266. 

  3266. !     real tva                  ! Temperature          Vertical surfaces, A (implicit) term
    3267. 

  3267. !     real tvb                  ! Temperature          Vertical surfaces, B (explicit) term
    3268. 

  3268.       real ehb                  ! Energy (TKE)       Horizontal surfaces, B (explicit) term
    3269. 

  3269.       real sf
    3270. 

  3270. 

  3271. 

  3272. ! ----------------------------------------------------------------------
    3273. 

  3273. ! LOCAL:
    3274. 

  3274. ! ----------------------------------------------------------------------
    3275. 

  3275.       real aa
    3276. 

  3276.       real al
    3277. 

  3277.       real buu
    3278. 

  3278.       real c
    3279. 

  3279.       real fbuw
    3280. 

  3280.       real fbpt
    3281. 

  3281.       real fh
    3282. 

  3282.       real fm
    3283. 

  3283.       real ric
    3284. 

  3284.       real tstar
    3285. 

  3285.       real ustar
    3286. 

  3286.       real utot
    3287. 

  3287.       real wstar
    3288. 

  3288.       real zz
    3289. 

  3289.       
    3290. 

  3290.       real b,cm,ch,rr,tol
    3291. 

  3291.       parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
    3292. 

  3292. 

  3293. ! ----------------------------------------------------------------------
    3294. 

  3294. ! END VARIABLES DEFINITIONS
    3295. 

  3295. ! ----------------------------------------------------------------------
    3296. 

  3296. 

  3297. 

  3298. ! computation of the ground temperature
    3299. 

  3299.          
    3300. 

  3300.       utot=(ua**2+va**2)**.5
    3301. 

  3301.         
    3302. 

  3302.       
    3303. 

  3303. !!!! Louis formulation
    3304. 

  3304. !
    3305. 

  3305. ! compute the bulk Richardson Number
    3306. 

  3306. 

  3307.       zz=dz/2.
    3308. 

  3308.    
    3309. 

  3309. 

  3310.       utot=max(utot,0.01)
    3311. 

  3311.           
    3312. 

  3312.       ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
    3313. 

  3313.               
    3314. 

  3314.       aa=vk/log(zz/z0)
    3315. 

  3315. 

  3316. ! determine the parameters fm and fh for stable, neutral and unstable conditions
    3317. 

  3317. 

  3318.       if(ric.gt.0)then
    3319. 

  3319.          fm=1/(1+0.5*b*ric)**2
    3320. 

  3320.          fh=fm
    3321. 

  3321.       else
    3322. 

  3322.          c=b*cm*aa*aa*(zz/z0)**.5
    3323. 

  3323.          fm=1-b*ric/(1+c*(-ric)**.5)
    3324. 

  3324.          c=c*ch/cm
    3325. 

  3325.          fh=1-b*ric/(1+c*(-ric)**.5)
    3326. 

  3326.       endif
    3327. 

  3327.       
    3328. 

  3328.       fbuw=-aa*aa*utot*utot*fm
    3329. 

  3329.       fbpt=-aa*aa*utot*(pt-ptg)*fh/rr
    3330. 

  3330.                      
    3331. 

  3331.       ustar=(-fbuw)**.5
    3332. 

  3332.       tstar=-fbpt/ustar
    3333. 

  3333. 

  3334.       al=(vk*g_u*tstar)/(pt*ustar*ustar)                      
    3335. 

  3335.       
    3336. 

  3336.       buu=-g_u/pt0*ustar*tstar
    3337. 

  3337.        
    3338. 

  3338.       uhb=-ustar*ustar*ua/utot
    3339. 

  3339.       vhb=-ustar*ustar*va/utot 
    3340. 

  3340.       sf= ustar*tstar*da*cp_u   
    3341. 

  3341.        
    3342. 

  3342. !     thb= 0.      
    3343. 

  3343.       ehb=buu
    3344. 

  3344. !!!!!!!!!!!!!!!
    3345. 

  3345.          
    3346. 

  3346.       return
    3347. 

  3347.       end subroutine flux_flat
    3348. 

  3348. 

  3349. ! ===6=8===============================================================72
    3350. 

  3350. ! ===6=8===============================================================72
    3351. 

  3351. 

  3352.       subroutine icBEP (fww,fwg,fgw,fsw,fws,fsg,                       &
    3353. 

  3353.                        z0g_u,z0r_u,                                    &
    3354. 

  3354.                        nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,        &
    3355. 

  3355.                        nz_u,z_u)                               
    3356. 

  3356. 

  3357. 

  3358.       implicit none 
    3359. 

  3359.       
    3360. 

  3360.    
    3361. 

  3361. !    Building parameters      
    3362. 

  3362. 

  3363. !    Radiation parameters
    3364. 

  3364. 

  3365. !    Roughness parameters
    3366. 

  3366.       real z0g_u(nurbm)       ! The ground's roughness length      
    3367. 

  3367.       real z0r_u(nurbm)       ! The roof's roughness length     
    3368. 

    3369. 

  3369. !    Street parameters
    3370. 

  3370.       integer nd_u(nurbm)     ! Number of street direction for each urban class
    3371. 

  3371. 

  3372.       real strd_u(ndm,nurbm)  ! Street length (fix to greater value to the horizontal length of the cells)
    3373. 

  3373.       real ws_u(ndm,nurbm)    ! Street width [m]
    3374. 

  3374.       real bs_u(ndm,nurbm)    ! Building width [m]
    3375. 

  3375.       real h_b(nz_um,nurbm)   ! Bulding's heights [m]
    3376. 

  3376.       real d_b(nz_um,nurbm)   ! The probability that a building has an height h_b
    3377. 

  3377. ! -----------------------------------------------------------------------
    3378. 

  3378. !     Output
    3379. 

  3379. !------------------------------------------------------------------------
    3380. 

  3380. 

  3381. 

  3382. 

  3383. !   fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
    3384. 

  3384. !   and the short wave radation. They are the part of radiation from a surface
    3385. 

  3385. !   or from the sky to another surface.
    3386. 

  3386.       real fww(nz_um,nz_um,ndm,nurbm)        !  from wall to wall
    3387. 

  3387.       real fwg(nz_um,ndm,nurbm)              !  from wall to ground
    3388. 

  3388.       real fgw(nz_um,ndm,nurbm)              !  from ground to wall
    3389. 

  3389.       real fsw(nz_um,ndm,nurbm)              !  from sky to wall
    3390. 

  3390.       real fws(nz_um,ndm,nurbm)              !  from wall to sky
    3391. 

  3391.       real fsg(ndm,nurbm)                    !  from sky to ground
    3392. 

  3392. 

  3393.       real ss_u(nz_um,nurbm)     ! The probability that a building has an height equal to z
    3394. 

  3394.       real pb_u(nz_um,nurbm)     ! The probability that a building has an height greater or equal to z
    3395. 

  3395.         
    3396. 

  3396. !    Grid parameters
    3397. 

  3397.       integer nz_u(nurbm)     ! Number of layer in the urban grid
    3398. 

  3398.       real z_u(nz_um)       ! Height of the urban grid levels
    3399. 

  3399. 

  3400. 

  3401. ! -----------------------------------------------------------------------
    3402. 

  3402. !     Local
    3403. 

  3403. !------------------------------------------------------------------------
    3404. 

  3404. 

  3405.       integer iz_u,id,ilu,iurb
    3406. 

  3406. 

  3407.       real dtot
    3408. 

  3408.       real hbmax
    3409. 

  3409. 

  3410. !------------------------------------------------------------------------
    3411. 

  3411.      
    3412. 

  3412. 

  3413. ! -----------------------------------------------------------------------
    3414. 

  3414. !     This routine initialise the urban paramters for the BEP module
    3415. 

  3415. !------------------------------------------------------------------------
    3416. 

  3416. !
    3417. 

  3417. !Initialize some variables
    3418. 

  3418. !
    3419. 

  3419.        nz_u=0
    3420. 

  3420.        z_u=0.    
    3421. 

  3421.        ss_u=0.
    3422. 

  3422.        pb_u=0.
    3423. 

  3423.        fww=0.
    3424. 

  3424.        fwg=0.
    3425. 

  3425.        fgw=0.
    3426. 

  3426.        fsw=0.
    3427. 

  3427.        fws=0.
    3428. 

  3428.        fsg=0.
    3429. 

  3429. 

  3430. ! Computation of the urban levels height
    3431. 

  3431. 

  3432.       z_u(1)=0.
    3433. 

  3433.      
    3434. 

  3434.       do iz_u=1,nz_um-1
    3435. 

  3435.          z_u(iz_u+1)=z_u(iz_u)+dz_u
    3436. 

  3436.       enddo
    3437. 

  3437.       
    3438. 

  3438. ! Normalisation of the building density
    3439. 

  3439. 

  3440.       do iurb=1,nurbm
    3441. 

  3441.          dtot=0.
    3442. 

  3442.          do ilu=1,nz_um
    3443. 

  3443.             dtot=dtot+d_b(ilu,iurb)
    3444. 

  3444.          enddo
    3445. 

  3445.          do ilu=1,nz_um
    3446. 

  3446.             d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
    3447. 

  3447.          enddo
    3448. 

  3448.       enddo      
    3449. 

  3449. 

  3450. ! Compute the view factors, pb and ss 
    3451. 

  3451.       
    3452. 

  3452.       do iurb=1,nurbm         
    3453. 

  3453.          hbmax=0.
    3454. 

  3454.          nz_u(iurb)=0
    3455. 

  3455.          do ilu=1,nz_um
    3456. 

  3456.             if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
    3457. 

  3457.          enddo
    3458. 

  3458.          
    3459. 

  3459.          do iz_u=1,nz_um-1
    3460. 

  3460.             if(z_u(iz_u+1).gt.hbmax)go to 10
    3461. 

  3461.          enddo
    3462. 

  3462.          
    3463. 

  3463.  10      continue
    3464. 

  3464.          nz_u(iurb)=iz_u+1
    3465. 

  3465. 

  3466.          do id=1,nd_u(iurb)
    3467. 

  3467. 

  3468.             call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb),   &    
    3469. 

  3469.                             z_u,ws_u(id,iurb),                      &    
    3470. 

  3470.                             fww,fwg,fgw,fsg,fsw,fws) 
    3471. 

  3471. 

  3472.             do iz_u=1,nz_u(iurb)
    3473. 

  3473.                ss_u(iz_u,iurb)=0.
    3474. 

  3474.                do ilu=1,nz_um
    3475. 

  3475.                   if(z_u(iz_u).le.h_b(ilu,iurb)                      &    
    3476. 

  3476.                     .and.z_u(iz_u+1).gt.h_b(ilu,iurb))then            
    3477. 

  3477.                         ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
    3478. 

  3478.                   endif 
    3479. 

  3479.                enddo
    3480. 

  3480.             enddo
    3481. 

  3481. 

  3482.             pb_u(1,iurb)=1.
    3483. 

  3483.             do iz_u=1,nz_u(iurb)
    3484. 

  3484.                pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
    3485. 

  3485.             enddo
    3486. 

  3486. 

  3487.          enddo
    3488. 

  3488.       end do
    3489. 

  3489.      
    3490. 

  3490.                   
    3491. 

  3491.       return       
    3492. 

  3492.       end subroutine icBEP
    3493. 

  3493. 

  3494. ! ===6=8===============================================================72
    3495. 

  3495. ! ===6=8===============================================================72
    3496. 

  3496. 

  3497.       subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws) 
    3498. 

  3498.      
    3499. 

  3499.       implicit none
    3500. 

  3500. 

  3501.  
    3502. 

  3502. 

  3503. ! -----------------------------------------------------------------------
    3504. 

  3504. !     Input
    3505. 

  3505. !------------------------------------------------------------------------
    3506. 

  3506. 

  3507.       integer iurb            ! Number of the urban class
    3508. 

  3508.       integer nz_u            ! Number of levels in the urban grid
    3509. 

  3509.       integer id              ! Street direction number
    3510. 

  3510.       real ws                 ! Street width
    3511. 

  3511.       real z(nz_um)         ! Height of the urban grid levels
    3512. 

  3512.       real dxy                ! Street lenght
    3513. 

  3513. 

  3514. 

  3515. ! -----------------------------------------------------------------------
    3516. 

  3516. !     Output
    3517. 

  3517. !------------------------------------------------------------------------
    3518. 

  3518. 

  3519. !   fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
    3520. 

  3520. !   and the short wave radation. They are the part of radiation from a surface
    3521. 

  3521. !   or from the sky to another surface.
    3522. 

  3522. 

  3523.       real fww(nz_um,nz_um,ndm,nurbm)            !  from wall to wall
    3524. 

  3524.       real fwg(nz_um,ndm,nurbm)                  !  from wall to ground
    3525. 

  3525.       real fgw(nz_um,ndm,nurbm)                  !  from ground to wall
    3526. 

  3526.       real fsw(nz_um,ndm,nurbm)                  !  from sky to wall
    3527. 

  3527.       real fws(nz_um,ndm,nurbm)                  !  from wall to sky
    3528. 

  3528.       real fsg(ndm,nurbm)                        !  from sky to ground
    3529. 

  3529. 

  3530. 

  3531. ! -----------------------------------------------------------------------
    3532. 

  3532. !     Local
    3533. 

  3533. !------------------------------------------------------------------------
    3534. 

  3534. 

  3535.       integer jz,iz
    3536. 

  3536. 

  3537.       real hut
    3538. 

  3538.       real f1,f2,f12,f23,f123,ftot
    3539. 

  3539.       real fprl,fnrm
    3540. 

  3540.       real a1,a2,a3,a4,a12,a23,a123
    3541. 

  3541. 

  3542. ! -----------------------------------------------------------------------
    3543. 

  3543. !     This routine calculates the view factors
    3544. 

  3544. !------------------------------------------------------------------------
    3545. 

  3545.         
    3546. 

  3546.       hut=z(nz_u+1)
    3547. 

  3547.         
    3548. 

  3548.       do jz=1,nz_u      
    3549. 

  3549.       
    3550. 

  3550. ! radiation from wall to wall
    3551. 

  3551.        
    3552. 

  3552.          do iz=1,nz_u
    3553. 

  3553.      
    3554. 

  3554.             call fprls (fprl,dxy,abs(z(jz+1)-z(iz  )),ws)
    3555. 

  3555.             f123=fprl
    3556. 

  3556.             call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
    3557. 

  3557.             f23=fprl
    3558. 

  3558.             call fprls (fprl,dxy,abs(z(jz  )-z(iz  )),ws)
    3559. 

  3559.             f12=fprl
    3560. 

  3560.             call fprls (fprl,dxy,abs(z(jz  )-z(iz+1)),ws)
    3561. 

  3561.             f2 = fprl
    3562. 

  3562.        
    3563. 

  3563.             a123=dxy*(abs(z(jz+1)-z(iz  )))
    3564. 

  3564.             a12 =dxy*(abs(z(jz  )-z(iz  )))
    3565. 

  3565.             a23 =dxy*(abs(z(jz+1)-z(iz+1)))
    3566. 

  3566.             a1  =dxy*(abs(z(iz+1)-z(iz  )))
    3567. 

  3567.             a2  =dxy*(abs(z(jz  )-z(iz+1)))
    3568. 

  3568.             a3  =dxy*(abs(z(jz+1)-z(jz  )))
    3569. 

  3569.        
    3570. 

  3570.             ftot=0.5*(a123*f123-a23*f23-a12*f12+a2*f2)/a1
    3571. 

  3571.        
    3572. 

  3572.             fww(iz,jz,id,iurb)=ftot*a1/a3
    3573. 

  3573. 

  3574.          enddo 
    3575. 

  3575. 

  3576. ! radiation from ground to wall
    3577. 

  3577.        
    3578. 

  3578.          call fnrms (fnrm,z(jz+1),dxy,ws)
    3579. 

  3579.          f12=fnrm
    3580. 

  3580.          call fnrms (fnrm,z(jz)  ,dxy,ws)
    3581. 

  3581.          f1=fnrm
    3582. 

  3582.        
    3583. 

  3583.          a1 = ws*dxy
    3584. 

  3584.          
    3585. 

  3585.          a12= ws*dxy
    3586. 

  3586.        
    3587. 

  3587.          a4=(z(jz+1)-z(jz))*dxy
    3588. 

  3588.        
    3589. 

  3589.          ftot=(a12*f12-a12*f1)/a1
    3590. 

  3590.                     
    3591. 

  3591.          fgw(jz,id,iurb)=ftot*a1/a4
    3592. 

  3592.      
    3593. 

  3593. !  radiation from sky to wall
    3594. 

  3594.      
    3595. 

  3595.          call fnrms(fnrm,hut-z(jz)  ,dxy,ws)
    3596. 

  3596.          f12 = fnrm
    3597. 

  3597.          call fnrms (fnrm,hut-z(jz+1),dxy,ws)
    3598. 

  3598.          f1 =fnrm
    3599. 

  3599.        
    3600. 

  3600.          a1 = ws*dxy
    3601. 

  3601.        
    3602. 

  3602.          a12= ws*dxy
    3603. 

  3603.               
    3604. 

  3604.          a4 = (z(jz+1)-z(jz))*dxy
    3605. 

  3605.        
    3606. 

  3606.          ftot=(a12*f12-a12*f1)/a1
    3607. 

  3607.         
    3608. 

  3608.          fsw(jz,id,iurb)=ftot*a1/a4       
    3609. 

  3609.       
    3610. 

  3610.       enddo
    3611. 

  3611. 

  3612. ! radiation from wall to sky      
    3613. 

  3613.       do iz=1,nz_u
    3614. 

  3614.        call fnrms(fnrm,ws,dxy,hut-z(iz))
    3615. 

  3615.        f12=fnrm
    3616. 

  3616.        call fnrms(fnrm,ws,dxy,hut-z(iz+1))
    3617. 

  3617.        f1=fnrm
    3618. 

  3618.        a1 = (z(iz+1)-z(iz))*dxy
    3619. 

  3619.        a2 = (hut-z(iz+1))*dxy
    3620. 

  3620.        a12= (hut-z(iz))*dxy
    3621. 

  3621.        a4 = ws*dxy
    3622. 

  3622.        ftot=(a12*f12-a2*f1)/a1
    3623. 

  3623.        fws(iz,id,iurb)=ftot*a1/a4 
    3624. 

  3624.  
    3625. 

  3625.       enddo
    3626. 

  3626. !!!!!!!!!!!!!
    3627. 

  3627. 

  3628. 

  3629.        do iz=1,nz_u
    3630. 

  3630. 

  3631. ! radiation from wall to ground
    3632. 

  3632.       
    3633. 

  3633.          call fnrms (fnrm,ws,dxy,z(iz+1))
    3634. 

  3634.          f12=fnrm
    3635. 

  3635.          call fnrms (fnrm,ws,dxy,z(iz  ))
    3636. 

  3636.          f1 =fnrm
    3637. 

  3637.          
    3638. 

  3638.          a1= (z(iz+1)-z(iz) )*dxy
    3639. 

  3639.        
    3640. 

  3640.          a2 = z(iz)*dxy
    3641. 

  3641.          a12= z(iz+1)*dxy
    3642. 

  3642.          a4 = ws*dxy
    3643. 

  3643. 

  3644.          ftot=(a12*f12-a2*f1)/a1        
    3645. 

  3645.                     
    3646. 

  3646.          fwg(iz,id,iurb)=ftot*a1/a4
    3647. 

  3647.         
    3648. 

  3648.       enddo
    3649. 

  3649. 

  3650. ! radiation from sky to ground
    3651. 

  3651.       
    3652. 

  3652.       call fprls (fprl,dxy,ws,hut)
    3653. 

  3653.       fsg(id,iurb)=fprl
    3654. 

  3654. 

  3655.       return
    3656. 

  3656.       end subroutine view_factors
    3657. 

  3657. 

  3658. ! ===6=8===============================================================72
    3659. 

  3659. ! ===6=8===============================================================72
    3660. 

  3660. 

  3661.       SUBROUTINE fprls (fprl,a,b,c)
    3662. 

  3662. 

  3663.       implicit none
    3664. 

  3664. 

  3665.      
    3666. 

  3666.             
    3667. 

  3667.       real a,b,c
    3668. 

  3668.       real x,y
    3669. 

  3669.       real fprl
    3670. 

  3670. 

  3671. 

  3672.       x=a/c
    3673. 

  3673.       y=b/c
    3674. 

  3674.       
    3675. 

  3675.       if(a.eq.0.or.b.eq.0.)then
    3676. 

  3676.        fprl=0.
    3677. 

  3677.       else
    3678. 

  3678.        fprl=log( ( (1.+x**2)*(1.+y**2)/(1.+x**2+y**2) )**.5)+  &
    3679. 

  3679.            y*((1.+x**2)**.5)*atan(y/((1.+x**2)**.5))+          &  
    3680. 

  3680.            x*((1.+y**2)**.5)*atan(x/((1.+y**2)**.5))-          &   
    3681. 

  3681.            y*atan(y)-x*atan(x)
    3682. 

  3682.        fprl=fprl*2./(pi*x*y)
    3683. 

  3683.       endif
    3684. 

  3684.       
    3685. 

  3685.       return
    3686. 

  3686.       end subroutine fprls
    3687. 

  3687. 

  3688. ! ===6=8===============================================================72     
    3689. 

  3689. ! ===6=8===============================================================72
    3690. 

  3690. 

  3691.       SUBROUTINE fnrms (fnrm,a,b,c)
    3692. 

  3692. 

  3693.       implicit none
    3694. 

  3694. 

  3695. 

  3696. 

  3697.       real a,b,c
    3698. 

  3698.       real x,y,z,a1,a2,a3,a4,a5,a6
    3699. 

  3699.       real fnrm
    3700. 

  3700.       
    3701. 

  3701.       x=a/b
    3702. 

  3702.       y=c/b
    3703. 

  3703.       z=x**2+y**2
    3704. 

  3704.       
    3705. 

  3705.       if(y.eq.0.or.x.eq.0)then
    3706. 

  3706.        fnrm=0.
    3707. 

  3707.       else
    3708. 

  3708.        a1=log( (1.+x*x)*(1.+y*y)/(1.+z) )
    3709. 

  3709.        a2=y*y*log(y*y*(1.+z)/z/(1.+y*y) )
    3710. 

  3710.        a3=x*x*log(x*x*(1.+z)/z/(1.+x*x) )
    3711. 

  3711.        a4=y*atan(1./y)
    3712. 

  3712.        a5=x*atan(1./x)
    3713. 

  3713.        a6=sqrt(z)*atan(1./sqrt(z))
    3714. 

  3714.        fnrm=0.25*(a1+a2+a3)+a4+a5-a6
    3715. 

  3715.        fnrm=fnrm/(pi*y)
    3716. 

  3716.       endif
    3717. 

  3717.       
    3718. 

  3718.       return
    3719. 

  3719.       end subroutine fnrms
    3720. 

  3720.   ! ===6=8===============================================================72  
    3721. 

  3721.      
    3722. 

  3722.       SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
    3723. 

  3723.         twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
    3724. 

  3724.         emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b, &
    3725. 

  3725.         cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
    3726. 

  3726.         targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)
    3727. 

  3727. 

  3728. ! initialization routine, where the variables from the table are read
    3729. 

  3729. 

  3730.       implicit none
    3731. 

  3731.       
    3732. 

  3732.       integer iurb            ! urban class number
    3733. 

  3733. !    Building parameters      
    3734. 

  3734.       real alag_u(nurbm)      ! Ground thermal diffusivity [m^2 s^-1]
    3735. 

  3735.       real alaw_u(nurbm)      ! Wall thermal diffusivity [m^2 s^-1]
    3736. 

  3736.       real alar_u(nurbm)      ! Roof thermal diffusivity [m^2 s^-1]
    3737. 

  3737.       real csg_u(nurbm)       ! Specific heat of the ground material [J m^3 K^-1]
    3738. 

  3738.       real csw_u(nurbm)       ! Specific heat of the wall material [J m^3 K^-1]
    3739. 

  3739.       real csr_u(nurbm)       ! Specific heat of the roof material [J m^3 K^-1]
    3740. 

  3740.       real twini_u(nurbm)     ! Temperature inside the buildings behind the wall [K]
    3741. 

  3741.       real trini_u(nurbm)     ! Temperature inside the buildings behind the roof [K]
    3742. 

  3742.       real tgini_u(nurbm)     ! Initial road temperature
    3743. 

  3743. 

  3744. !    Radiation parameters
    3745. 

  3745.       real albg_u(nurbm)      ! Albedo of the ground
    3746. 

  3746.       real albw_u(nurbm)      ! Albedo of the wall
    3747. 

  3747.       real albr_u(nurbm)      ! Albedo of the roof
    3748. 

  3748.       real albwin_u(nurbm)    ! Albedo of the window
    3749. 

  3749.       real emg_u(nurbm)       ! Emissivity of ground
    3750. 

  3750.       real emw_u(nurbm)       ! Emissivity of wall
    3751. 

  3751.       real emr_u(nurbm)       ! Emissivity of roof
    3752. 

  3752.       real emwind_u(nurbm)    ! Emissivity of windows
    3753. 

  3753. 

  3754. !    Roughness parameters
    3755. 

  3755.       real z0g_u(nurbm)       ! The ground's roughness length      
    3756. 

  3756.       real z0r_u(nurbm)       ! The roof's roughness length
    3757. 

    3758. 

  3758. !    Street parameters
    3759. 

  3759.       integer nd_u(nurbm)     ! Number of street direction for each urban class
    3760. 

  3760. 

  3761.       real strd_u(ndm,nurbm)  ! Street length (fix to greater value to the horizontal length of the cells)
    3762. 

  3762.       real drst_u(ndm,nurbm)  ! Street direction [degree]
    3763. 

  3763.       real ws_u(ndm,nurbm)    ! Street width [m]
    3764. 

  3764.       real bs_u(ndm,nurbm)    ! Building width [m]
    3765. 

  3765.       real h_b(nz_um,nurbm)   ! Bulding's heights [m]
    3766. 

  3766.       real d_b(nz_um,nurbm)   ! The probability that a building has an height h_b
    3767. 

  3767. 

  3768.       integer i,iu
    3769. 

  3769.       integer nurb ! number of urban classes used
    3770. 

  3770. 

  3771.       real, intent(out) :: cop_u(nurbm)
    3772. 

  3772.       real, intent(out) :: pwin_u(nurbm)
    3773. 

  3773.       real, intent(out) :: beta_u(nurbm)
    3774. 

  3774.       integer, intent(out) :: sw_cond_u(nurbm)
    3775. 

  3775.       real, intent(out) :: time_on_u(nurbm)
    3776. 

  3776.       real, intent(out) :: time_off_u(nurbm)
    3777. 

  3777.       real, intent(out) :: targtemp_u(nurbm)
    3778. 

  3778.       real, intent(out) :: gaptemp_u(nurbm)
    3779. 

  3779.       real, intent(out) :: targhum_u(nurbm)
    3780. 

  3780.       real, intent(out) :: gaphum_u(nurbm)
    3781. 

  3781.       real, intent(out) :: perflo_u(nurbm)
    3782. 

  3782.       real, intent(out) :: hsesf_u(nurbm)
    3783. 

  3783.       real, intent(out) :: hsequip(24)
    3784. 

  3784. 

  3785. !
    3786. 

  3786. !We initialize
    3787. 

  3787. !
    3788. 

  3788.         h_b=0.
    3789. 

  3789.         d_b=0.
    3790. 

  3790. 

  3791.        nurb=ICATE
    3792. 

  3792.        do iu=1,nurb                         
    3793. 

  3793.           nd_u(iu)=0
    3794. 

  3794.        enddo
    3795. 

  3795. 

  3796.        csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
    3797. 

  3797.        csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
    3798. 

  3798.        csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
    3799. 

  3799.        do i=1,icate
    3800. 

  3800.          alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
    3801. 

  3801.          alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
    3802. 

  3802.          alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
    3803. 

  3803.        enddo
    3804. 

  3804.        twini_u=TBLEND_TBL
    3805. 

  3805.        trini_u=TRLEND_TBL
    3806. 

  3806.        tgini_u=TGLEND_TBL
    3807. 

  3807.        albw_u=ALBB_TBL
    3808. 

  3808.        albr_u=ALBR_TBL
    3809. 

  3809.        albg_u=ALBG_TBL
    3810. 

  3810.        emw_u=EPSB_TBL
    3811. 

  3811.        emr_u=EPSR_TBL
    3812. 

  3812.        emg_u=EPSG_TBL
    3813. 

  3813.        z0r_u=Z0R_TBL
    3814. 

  3814.        z0g_u=Z0G_TBL
    3815. 

  3815.        nd_u=NUMDIR_TBL
    3816. 

  3816. !MT BEM
    3817. 

  3817.        cop_u = cop_tbl
    3818. 

  3818.        pwin_u = pwin_tbl
    3819. 

  3819.        beta_u = beta_tbl
    3820. 

  3820.        sw_cond_u = sw_cond_tbl
    3821. 

  3821.        time_on_u = time_on_tbl
    3822. 

  3822.        time_off_u = time_off_tbl
    3823. 

  3823.        targtemp_u = targtemp_tbl
    3824. 

  3824.        gaptemp_u = gaptemp_tbl
    3825. 

  3825.        targhum_u = targhum_tbl
    3826. 

  3826.        gaphum_u = gaphum_tbl
    3827. 

  3827.        perflo_u = perflo_tbl
    3828. 

  3828.        hsesf_u = hsesf_tbl
    3829. 

  3829.        hsequip = hsequip_tbl
    3830. 

  3830. 

  3831.        do iu=1,icate
    3832. 

  3832.               if(ndm.lt.nd_u(iu))then
    3833. 

  3833.                 write(*,*)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
    3834. 

  3834.                 write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
    3835. 

  3835.                 stop
    3836. 

  3836.               endif
    3837. 

  3837.          do i=1,nd_u(iu)
    3838. 

  3838.            drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
    3839. 

  3839.            ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
    3840. 

  3840.            bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
    3841. 

  3841.          enddo
    3842. 

  3842.        enddo
    3843. 

  3843.        do iu=1,ICATE
    3844. 

  3844.           if(nz_um.lt.numhgt_tbl(iu)+3)then
    3845. 

  3845.               write(*,*)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
    3846. 

  3846.               write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
    3847. 

  3847.               stop
    3848. 

  3848.           endif
    3849. 

  3849.          do i=1,NUMHGT_TBL(iu)
    3850. 

  3850.            h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
    3851. 

  3851.            d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
    3852. 

  3852.          enddo
    3853. 

  3853.        enddo
    3854. 

  3854. 

  3855.        do i=1,ndm
    3856. 

  3856.         do iu=1,nurbm
    3857. 

  3857.          strd_u(i,iu)=100000.
    3858. 

  3858.         enddo
    3859. 

  3859.        enddo
    3860. 

  3860. 

  3861.        do iu=1,nurb  
    3862. 

  3862.           emwind_u(iu)=0.9                       
    3863. 

  3863.           call albwindow(albwin_u(iu))  
    3864. 

  3864.        enddo
    3865. 

  3865.        
    3866. 

  3866.        return
    3867. 

  3867.        end subroutine init_para
    3868. 

  3868. ! ===6================================================================72
    3869. 

  3869. ! ===6================================================================72
    3870. 

  3870.       subroutine upward_rad(nd_u,nz_u,ws,bs,sigma,pb,ss,                     &
    3871. 

  3871.                        tg,emg_u,albg_u,rlg,rsg,sfg,                          & 
    3872. 

  3872.                        tw,emw_u,albw_u,rlw,rsw,sfw,                          & 
    3873. 

  3873.                        tr,emr_u,albr_u,emwind,albwind,twlev,pwin,            &
    3874. 

  3874.                        sfwind,rld,rs, sfr,                                   & 
    3875. 

  3875.                        rs_abs,rl_up,emiss,grdflx_urb)
    3876. 

  3876. !
    3877. 

  3877. ! IN this surboutine we compute the upward longwave flux, and the albedo
    3878. 

  3878. ! needed for the radiation scheme
    3879. 

  3879. !
    3880. 

  3880.       implicit none
    3881. 

  3881. 

  3882. !
    3883. 

  3883. !INPUT VARIABLES
    3884. 

  3884. !
    3885. 

  3885.       real rsw(2*ndm,nz_um)        ! Short wave radiation at the wall for a given canyon direction [W/m2]
    3886. 

  3886.       real rlw(2*ndm,nz_um)         ! Long wave radiation at the walls for a given canyon direction [W/m2]
    3887. 

  3887.       real rsg(ndm)                   ! Short wave radiation at the canyon for a given canyon direction [W/m2]
    3888. 

  3888.       real rlg(ndm)                   ! Long wave radiation at the ground for a given canyon direction [W/m2]
    3889. 

  3889.       real rs                        ! Short wave radiation at the horizontal surface from the sun [W/m²]  
    3890. 

  3890.       real sfw(2*ndm,nz_um)      ! Sensible heat flux from walls [W/m²]
    3891. 

  3891.       real sfg(ndm)              ! Sensible heat flux from ground (road) [W/m²]
    3892. 

  3892.       real sfr(ndm,nz_um)      ! Sensible heat flux from roofs [W/m²]                      
    3893. 

  3893.       real rld                        ! Long wave radiation from the sky [W/m²]
    3894. 

  3894.       real albg_u                    ! albedo of the ground/street
    3895. 

  3895.       real albw_u                    ! albedo of the walls
    3896. 

  3896.       real albr_u                    ! albedo of the roof 
    3897. 

  3897.       real ws(ndm)                        ! width of the street
    3898. 

  3898.       real bs(ndm)
    3899. 

  3899.                         ! building size
    3900. 

  3900.       real pb(nz_um)                ! Probability to have a building with an height equal or higher   
    3901. 

  3901.       integer nz_u
    3902. 

  3902.       real ss(nz_um)                ! Probability to have a building of a given height
    3903. 

  3903.       real sigma                       
    3904. 

  3904.       real emg_u                       ! emissivity of the street
    3905. 

  3905.       real emw_u                       ! emissivity of the wall
    3906. 

  3906.       real emr_u                       ! emissivity of the roof
    3907. 

  3907.       real tw(2*ndm,nz_um)  ! Temperature in each layer of the wall [K]
    3908. 

  3908.       real tr(ndm,nz_um,nwr_u)  ! Temperature in each layer of the roof [K]
    3909. 

  3909.       real tg(ndm,ng_u)          ! Temperature in each layer of the ground [K]
    3910. 

  3910.       integer id ! street direction
    3911. 

  3911.       integer nd_u ! number of street directions
    3912. 

  3912. !
    3913. 

  3913. !New variables BEM
    3914. 

  3914. !
    3915. 

  3915.       real emwind               !Emissivity of the windows
    3916. 

  3916.       real albwind              !Albedo of the windows
    3917. 

  3917.       real twlev(2*ndm,nz_um)   !Averaged Temperature of the windows 
    3918. 

  3918.       real pwin                 !Coverage area fraction of the windows
    3919. 

  3919.       real gflwin               !Heat stored for the windows
    3920. 

  3920.       real sfwind(2*ndm,nz_um)  !Sensible heat flux from windows [W/m²]
    3921. 

  3921. 

  3922. !OUTPUT/INPUT
    3923. 

  3923.       real rs_abs  ! absrobed solar radiationfor this street direction
    3924. 

  3924.       real rl_up   ! upward longwave radiation for this street direction
    3925. 

  3925.       real emiss ! mean emissivity
    3926. 

  3926.       real grdflx_urb ! ground heat flux 
    3927. 

  3927. !LOCAL
    3928. 

  3928.       integer iz,iw
    3929. 

  3929.       real rl_inc,rl_emit
    3930. 

  3930.       real gfl
    3931. 

  3931.       integer ix,iy,iwrong
    3932. 

  3932. 

  3933.          iwrong=1
    3934. 

  3934.       do iz=1,nz_u+1
    3935. 

  3935.       do id=1,nd_u
    3936. 

  3936.       do iw=1,nwr_u
    3937. 

  3937.         if(tr(id,iz,iw).lt.100.)then
    3938. 

  3938.               write(*,*)'in upward_rad ',iz,id,iw,tr(id,iz,iw) 
    3939. 

  3939.               iwrong=0
    3940. 

  3940.         endif
    3941. 

  3941.       enddo
    3942. 

  3942.       enddo
    3943. 

  3943.       enddo
    3944. 

  3944.            if(iwrong.eq.0)stop
    3945. 

  3945. 

  3946.       rl_up=0.
    3947. 

  3947.  
    3948. 

  3948.       rs_abs=0.
    3949. 

  3949.       rl_inc=0.
    3950. 

  3950.       emiss=0.
    3951. 

  3951.       rl_emit=0.
    3952. 

  3952.       grdflx_urb=0.
    3953. 

  3953.       do id=1,nd_u          
    3954. 

  3954.        rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/nd_u
    3955. 

  3955.        rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u       
    3956. 

  3956.        rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/nd_u
    3957. 

  3957.          gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg(id,ng_u)**4.)+sfg(id)
    3958. 

  3958.          grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u  
    3959. 

  3959.  
    3960. 

  3960.           do iz=2,nz_u
    3961. 

  3961.             rl_emit=rl_emit-(emr_u*sigma*(tr(id,iz,nwr_u)**4.)+(1-emr_u)*rld)*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
    3962. 

  3962.             rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u            
    3963. 

  3963.             rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
    3964. 

  3964.             gfl=(1.-albr_u)*rs+emr_u*rld-emr_u*sigma*(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
    3965. 

  3965.             grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
    3966. 

  3966.          enddo
    3967. 

  3967.            
    3968. 

  3968.          do iz=1,nz_u 
    3969. 

  3969.            
    3970. 

  3970.             rl_emit=rl_emit-(emw_u*(1.-pwin)*sigma*(tw(2*id-1,iz)**4.+tw(2*id,iz)**4.)+ &
    3971. 

  3971.                             (emwind*pwin*sigma*(twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.))+ &
    3972. 

  3972.                 ((1.-emw_u)*(1.-pwin)+pwin*(1.-emwind))*(rlw(2*id-1,iz)+rlw(2*id,iz)))* &
    3973. 

  3973.                             dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
    3974. 

  3974. 

  3975.             rl_inc=rl_inc+((rlw(2*id-1,iz)+rlw(2*id,iz)))*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
    3976. 

  3976. 

  3977.             rs_abs=rs_abs+(((1.-albw_u)*(1.-pwin)+(1.-albwind)*pwin)*(rsw(2*id-1,iz)+rsw(2*id,iz)))*&
    3978. 

  3978.                           dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u 
    3979. 

  3979. 

  3980.             gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) )   &
    3981. 

  3981.              -emw_u*sigma*( tw(2*id-1,iz)**4.+tw(2*id,iz)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))   
    3982. 

  3982. 

  3983.             gflwin=(1.-albwind)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emwind*(rlw(2*id-1,iz)+rlw(2*id,iz))   &
    3984. 

  3984.              -emwind*sigma*( twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.)+(sfwind(2*id-1,iz)+sfwind(2*id,iz)) 
    3985. 

  3985.                
    3986. 

  3986.            
    3987. 

  3987.             grdflx_urb=grdflx_urb-(gfl*(1.-pwin)+pwin*gflwin)*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
    3988. 

  3988. 

  3989.          enddo
    3990. 

  3990.           
    3991. 

  3991.       enddo
    3992. 

  3992.         emiss=(emg_u+emw_u+emr_u)/3.
    3993. 

  3993.         rl_up=(rl_inc+rl_emit)-rld
    3994. 

  3994.        
    3995. 

  3995.          
    3996. 

  3996.       return
    3997. 

  3997. 

  3998.       END SUBROUTINE upward_rad
    3999. 

  3999. 

  4000. !====6=8===============================================================72         
    4001. 

  4001. !====6=8===============================================================72 
    4002. 

  4002. ! ===6================================================================72
    4003. 

  4003. ! ===6================================================================72
    4004. 

  4004. 

  4005.          subroutine albwindow(albwin)
    4006. 

  4006. 		
    4007. 

  4007. !-------------------------------------------------------------------
    4008. 

  4008. 	 implicit none
    4009. 

  4009. 

  4010. 

  4011. ! -------------------------------------------------------------------
    4012. 

  4012. !Based on the 
    4013. 

  4013. !paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour
    4014. 

  4014. !of the total solar energy transmittance of windows"
    4015. 

  4015. !Solar Energy Vol.69,No.4,pp. 321-329.      
    4016. 

  4016. ! -------------------------------------------------------------------
    4017. 

  4017. !Input
    4018. 

  4018. !-----	
    4019. 

  4019.         
    4020. 

  4020. !Output
    4021. 

  4021. !------
    4022. 

  4022.          real albwin	        ! albedo of the window  
    4023. 

  4023. !Local
    4024. 

  4024. !-----
    4025. 

  4025. 	 real a,b,c		!Polynomial coefficients
    4026. 

  4026. 	 real alfa,delta,gama	!Polynomial powers
    4027. 

  4027. 	 real g0	        !transmittance when the angle 
    4028. 

  4028.                                 !of incidence is normal to the surface.
    4029. 

  4029.          real asup,ainf
    4030. 

  4030. 	 real fonc
    4031. 

  4031. 

  4032. !Constants
    4033. 

  4033. !--------------------
    4034. 

  4034.          
    4035. 

  4035.          real epsilon              !accuracy of the integration
    4036. 

  4036.          parameter (epsilon=1.e-07) 
    4037. 

  4037.          real n1,n2                !Index of refraction for glasses and air
    4038. 

  4038.          parameter(n1=1.,n2=1.5)
    4039. 

  4039.          integer intg,k
    4040. 

  4040. !--------------------------------------------------------------------		
    4041. 

  4041.          if (q_num.eq.0) then
    4042. 

  4042.            write(*,*) 'Category parameter of the windows no valid'
    4043. 

  4043.            stop
    4044. 

  4044.          endif
    4045. 

  4045. 

  4046.          g0=4.*n1*n2/((n1+n2)*(n1+n2))
    4047. 

  4047. 	 a=8.
    4048. 

  4048. 	 b=0.25/q_num
    4049. 

  4049.          c=1.-a-b	
    4050. 

  4050. 	 alfa =5.2 + (0.7*q_num)
    4051. 

  4051. 	 delta =2.
    4052. 

  4052. 	 gama =(5.26+0.06*p_num)+(0.73+0.04*p_num)*q_num
    4053. 

  4053. 

  4054.          intg=1
    4055. 

  4055. !----------------------------------------------------------------------
    4056. 

  4056. 

  4057. 

  4058. 100      asup=0.
    4059. 

  4059.          ainf=0.
    4060. 

  4060. 

  4061.          do k=1,intg
    4062. 

  4062.           call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
    4063. 

  4063.           asup=asup+(pi/intg)*fonc
    4064. 

  4064.          enddo
    4065. 

  4065. 

  4066.          intg=intg+1
    4067. 

  4067. 

  4068.          do k=1,intg
    4069. 

  4069.           call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
    4070. 

  4070.           ainf=ainf+(pi/intg)*fonc
    4071. 

  4071.          enddo
    4072. 

  4072. 	 
    4073. 

  4073.          if(abs(asup-ainf).lt.epsilon) then
    4074. 

  4074.            albwin=1-g0+(g0/2.)*asup
    4075. 

  4075.          else
    4076. 

  4076.            goto 100
    4077. 

  4077.          endif
    4078. 

  4078.         
    4079. 

  4079. !---------------------------------------------------------------------- 	
    4080. 

  4080. 	return
    4081. 

  4081. 	end subroutine albwindow
    4082. 

  4082. !====================================================================72
    4083. 

  4083. !====================================================================72
    4084. 

  4084. 

  4085.         subroutine foncs(fonc,x,aa,bb,cc,alf,delt,gam)
    4086. 

  4086. 

  4087.         implicit none
    4088. 

  4088. !
    4089. 

  4089.         real x,aa,bb,cc
    4090. 

  4090.         real alf,delt,gam
    4091. 

  4091.         real fonc
    4092. 

  4092.   
    4093. 

  4093.         fonc=(((aa*(x**alf))/(pi**alf))+   &
    4094. 

  4094.              ((bb*(x**delt))/(pi**delt))+  &
    4095. 

  4095.              ((cc*(x**gam))/(pi**gam)))*sin(x)
    4096. 

  4096.         
    4097. 

  4097.         return
    4098. 

  4098. 	end subroutine foncs
    4099. 

  4099. !====================================================================72
    4100. 

  4100. !====================================================================72  
    4101. 

  4101. END MODULE module_sf_bep_bem
    4102. 

  4102.