PageRenderTime 37ms CodeModel.GetById 31ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_sf_bep.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 3239 lines | 1873 code | 634 blank | 732 comment | 38 complexity | 0c3989f4b2ba6aaa590bea8e69427a24 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. MODULE module_sf_bep
    2. 

  2. 

  3. !USE module_model_constants
    4. 

  4.  USE module_sf_urban
    5. 

  5. 

  6. ! SGClarke 09/11/2008
    7. 

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

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

  9. !
    10. 

  10.  ! -----------------------------------------------------------------------
    11. 

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

  12. ! -----------------------------------------------------------------------
    13. 

  13. 

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

  15.       parameter (nurbm=3)
    16. 

  16. 

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

  18.       parameter (ndm=2)
    19. 

  19. 

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

  21.       parameter(nz_um=13)
    22. 

  22. 

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

  24.       parameter (ng_u=10)
    25. 

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

  26.       parameter (nwr_u=10)
    27. 

  27. 

  28.       real dz_u                 ! Urban grid resolution
    29. 

  29.       parameter (dz_u=5.)
    30. 

  30. 

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

  32. !  in the routine "surf_temp" 
    33. 

  33. ! -----------------------------------------------------------------------
    34. 

  34. !  Constant used in the BEP module
    35. 

  35. ! -----------------------------------------------------------------------
    36. 

  36.            
    37. 

  37.       real vk                 ! von Karman constant
    38. 

  38.       real g_u                  ! Gravity acceleration
    39. 

  39.       real pi                 !
    40. 

  40.       real r                  ! Perfect gas constant
    41. 

  41.       real cp_u                 ! Specific heat at constant pressure
    42. 

  42.       real rcp_u                !
    43. 

  43.       real sigma              !
    44. 

  44.       real p0                 ! Reference pressure at the sea level
    45. 

  45.       real cdrag              ! Drag force constant
    46. 

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

  47.       parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4)
    48. 

  48. 

  49. ! -----------------------------------------------------------------------     
    50. 

  50. 

  51. 

  52. 

  53. 

  54.    CONTAINS
    55. 

  55.  
    56. 

  56.       subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy,      &
    57. 

  57.                       th_phy,rho,p_phy,swdown,glw,                    &
    58. 

  58.                       gmt,julday,xlong,xlat,                          &
    59. 

  59.                       declin_urb,cosz_urb2d,omg_urb2d,                &
    60. 

  60.                       num_urban_layers,                               &
    61. 

  61.                       trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,        &
    62. 

  62.                       sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,      &
    63. 

  63.                       a_u,a_v,a_t,a_e,b_u,b_v,                        &
    64. 

  64.                       b_t,b_e,dlg,dl_u,sf,vl,                         &
    65. 

  65.                       rl_up,rs_abs,emiss,grdflx_urb,                  &
    66. 

  66.                       ids,ide, jds,jde, kds,kde,                      &
    67. 

  67.                       ims,ime, jms,jme, kms,kme,                      &
    68. 

  68.                       its,ite, jts,jte, kts,kte)                    
    69. 

  69. 

  70.       implicit none
    71. 

  71. 

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

  73. !     Input
    74. 

  74. !------------------------------------------------------------------------
    75. 

  75.    INTEGER ::                       ids,ide, jds,jde, kds,kde,  &
    76. 

  76.                                     ims,ime, jms,jme, kms,kme,  &
    77. 

  77.                                     its,ite, jts,jte, kts,kte,  &
    78. 

  78.                                     itimestep
    79. 

  79.  
    80. 

  80. 

  81.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   DZ8W
    82. 

  82.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   P_PHY
    83. 

  83.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   RHO
    84. 

  84.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   TH_PHY
    85. 

  85.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   T_PHY
    86. 

  86.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   U_PHY
    87. 

  87.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   V_PHY
    88. 

  88.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   U
    89. 

  89.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   V
    90. 

  90.    REAL, DIMENSION( ims:ime , jms:jme )        ::   GLW
    91. 

  91.    REAL, DIMENSION( ims:ime , jms:jme )        ::   swdown
    92. 

  92.    REAL, DIMENSION( ims:ime, jms:jme )         ::   UST
    93. 

  93.    INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   UTYPE_URB2D
    94. 

  94.    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   FRC_URB2D
    95. 

  95.    REAL, INTENT(IN  )   ::                                   GMT 
    96. 

  96.    INTEGER, INTENT(IN  ) ::                               JULDAY
    97. 

  97.    REAL, DIMENSION( ims:ime, jms:jme ),                           &
    98. 

  98.          INTENT(IN   )  ::                           XLAT, XLONG
    99. 

  99.    REAL, INTENT(IN) :: DECLIN_URB
    100. 

  100.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
    101. 

  101.    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
    102. 

  102.    INTEGER, INTENT(IN  ) :: num_urban_layers
    103. 

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

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

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

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

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

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

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

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

  111. !      integer nx,ny,nz              ! Number of points in the mesocsale grid
    112. 

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

  113.       REAL, INTENT(IN )::   DT      ! Time step
    114. 

  114. !      real zr(ims:ime,jms:jme)                ! Solar zenith angle
    115. 

  115. !      real deltar(ims:ime,jms:jme)            ! Declination of the sun
    116. 

  116. !      real ah(ims:ime,jms:jme)                ! Hour angle
    117. 

  117. !      real rs(ims:ime,jms:jme)                ! Solar radiation
    118. 

  118. !------------------------------------------------------------------------
    119. 

  119. !     Output
    120. 

  120. !------------------------------------------------------------------------ 
    121. 

  121. !      real tsk(ims:ime,jms:jme)           ! Average of surface temperatures (roads and roofs)
    122. 

  122. !
    123. 

  123. !    Implicit and explicit components of the source and sink terms at each levels,
    124. 

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

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

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

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

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

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

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

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

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

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

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

  135. ! urban surface and volumes        
    136. 

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

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

  138. ! urban fluxes
    139. 

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

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

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

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

  143. !------------------------------------------------------------------------
    144. 

  144. !     Local
    145. 

  145. !------------------------------------------------------------------------
    146. 

  146. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    147. 

  147. !    Building parameters      
    148. 

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

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

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

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

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

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

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

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

  156.       real tgini_u(nurbm)                     ! Initial road temperature
    157. 

  157. !
    158. 

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

  159. !
    160. 

  160. !    Radiation paramters
    161. 

  161.       real albg_u(nurbm)                      ! Albedo of the ground
    162. 

  162.       real albw_u(nurbm)                      ! Albedo of the wall
    163. 

  163.       real albr_u(nurbm)                      ! Albedo of the roof
    164. 

  164.       real emg_u(nurbm)                       ! Emissivity of ground
    165. 

  165.       real emw_u(nurbm)                       ! Emissivity of wall
    166. 

  166.       real emr_u(nurbm)                       ! Emissivity of roof
    167. 

  167. 

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

  169. !   and the short wave radation. 
    170. 

  170.       real fww(nz_um,nz_um,ndm,nurbm)         !  from wall to wall
    171. 

  171.       real fwg(nz_um,ndm,nurbm)               !  from wall to ground
    172. 

  172.       real fgw(nz_um,ndm,nurbm)               !  from ground to wall
    173. 

  173.       real fsw(nz_um,ndm,nurbm)               !  from sky to wall
    174. 

  174.       real fws(nz_um,ndm,nurbm)               !  from sky to wall
    175. 

  175.       real fsg(ndm,nurbm)                     !  from sky to ground
    176. 

  176. 

  177. !    Roughness parameters
    178. 

  178.       real z0g_u(nurbm)         ! The ground's roughness length
    179. 

  179.       real z0r_u(nurbm)         ! The roof's roughness length
    180. 

    181. 

  181. !    Street parameters
    182. 

  182.       integer nd_u(nurbm)       ! Number of street direction for each urban class 
    183. 

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

  184.       real drst_u(ndm,nurbm)    ! Street direction
    185. 

  185.       real ws_u(ndm,nurbm)      ! Street width
    186. 

  186.       real bs_u(ndm,nurbm)      ! Building width
    187. 

  187.       real h_b(nz_um,nurbm)     ! Bulding's heights
    188. 

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

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

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

  191. 

  192. 

  193. !    Grid parameters
    194. 

  194.       integer nz_u(nurbm)       ! Number of layer in the urban grid
    195. 

  195.       
    196. 

  196.       real z_u(nz_um)         ! Height of the urban grid levels
    197. 

  197. 

  198. ! 1D array used for the input and output of the routine "urban"
    199. 

  199. 

  200.       real z1D(kms:kme)               ! vertical coordinates
    201. 

  201.       real ua1D(kms:kme)                ! wind speed in the x directions
    202. 

  202.       real va1D(kms:kme)                ! wind speed in the y directions
    203. 

  203.       real pt1D(kms:kme)                ! potential temperature
    204. 

  204.       real da1D(kms:kme)                ! air density
    205. 

  205.       real pr1D(kms:kme)                ! air pressure
    206. 

  206.       real pt01D(kms:kme)               ! reference potential temperature
    207. 

  207.       real zr1D                    ! zenith angle
    208. 

  208.       real deltar1D                ! declination of the sun
    209. 

  209.       real ah1D                    ! hour angle (it should come from the radiation routine)
    210. 

  210.       real rs1D                    ! solar radiation
    211. 

  211.       real rld1D                   ! downward flux of the longwave radiation
    212. 

  212. 

  213. 

  214.       real tw1D(2*ndm,nz_um,nwr_u)  ! temperature in each layer of the wall
    215. 

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

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

  217.       real sfw1D(2*ndm,nz_um)      ! sensible heat flux from walls
    218. 

  218.       real sfg1D(ndm)              ! sensible heat flux from ground (road)
    219. 

  219.       real sfr1D(ndm,nz_um)      ! sensible heat flux from roofs
    220. 

  220.       real sf1D(kms:kme)              ! surface of the urban grid cells
    221. 

  221.       real vl1D(kms:kme)                ! volume of the urban grid cells
    222. 

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

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

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

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

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

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

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

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

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

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

  232.       real tsk1D                        ! Average of the road surface temperatures
    233. 

  233.       real time_bep
    234. 

  234. ! arrays used to collapse indexes
    235. 

  235.       integer ind_zwd(nz_um,nwr_u,ndm)
    236. 

  236.       integer ind_gd(ng_u,ndm)
    237. 

  237.       integer ind_zd(nz_um,ndm)
    238. 

  238. !
    239. 

  239.       integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
    240. 

  240.       integer it, nint
    241. 

  241.       integer iii
    242. 

  242.       real time_h,tempo,shtot
    243. 

  243.       logical first
    244. 

  244.       character(len=80) :: text
    245. 

  245.       data first/.true./
    246. 

  246.       save first,time_bep 
    247. 

  247. 

  248.       save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,                      &
    249. 

  249.           albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,fww,fwg,fgw,fsw,fws,fsg,   &
    250. 

  250.           z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  &
    251. 

  251.           nz_u,z_u 
    252. 

  252. 

  253. !------------------------------------------------------------------------
    254. 

  254. !    Calculation of the momentum, heat and turbulent kinetic fluxes
    255. 

  255. !     produced by buildings
    256. 

  256. !
    257. 

  257. ! Reference:
    258. 

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

  259. ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
    260. 

  260. ! 261-304
    261. 

  261. !------------------------------------------------------------------------
    262. 

  262. !prepare the arrays to collapse indexes
    263. 

  263. 

  264.       if(num_urban_layers.lt.nz_um*ndm*nwr_u)then
    265. 

  265.         write(*,*)'num_urban_layers too small, please increase to at least ', nz_um*ndm*nwr_u
    266. 

  266.         stop
    267. 

  267.       endif
    268. 

  268.       iii=0
    269. 

  269.       do iz_u=1,nz_um
    270. 

  270.       do iw=1,nwr_u
    271. 

  271.       do id=1,ndm
    272. 

  272.        iii=iii+1
    273. 

  273.        ind_zwd(iz_u,iw,id)=iii
    274. 

  274.       enddo
    275. 

  275.       enddo
    276. 

  276.       enddo
    277. 

  277. 

  278.       iii=0
    279. 

  279.       do ig=1,ng_u
    280. 

  280.       do id=1,ndm
    281. 

  281.        iii=iii+1
    282. 

  282.        ind_gd(ig,id)=iii
    283. 

  283.       enddo
    284. 

  284.       enddo
    285. 

  285. 

  286.       iii=0
    287. 

  287.       do iz_u=1,nz_um
    288. 

  288.       do id=1,ndm
    289. 

  289.        iii=iii+1
    290. 

  290.        ind_zd(iz_u,id)=iii
    291. 

  291.       enddo
    292. 

  292.       enddo
    293. 

  293.       do ix=its,ite
    294. 

  294.       do iy=jts,jte
    295. 

  295.        z(ix,kts,iy)=0.
    296. 

  296.        do iz=kts+1,kte+1
    297. 

  297.         z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
    298. 

  298.        enddo
    299. 

  299.       enddo
    300. 

  300.       enddo
    301. 

  301. 

  302.       if (first) then                           ! True only on first call
    303. 

  303.          call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
    304. 

  304.                 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
    305. 

  305.                 emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
    306. 

  306. 

  307. !      Initialisation of the urban parameters and calculation of the view factors
    308. 

  308.        call icBEP(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,          & 
    309. 

  309.                  albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,           & 
    310. 

  310.                  fww,fwg,fgw,fsw,fws,fsg,                          & 
    311. 

  311.                  z0g_u,z0r_u,                                      & 
    312. 

  312.                  nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,   & 
    313. 

  313.                  nz_u,z_u,                                         & 
    314. 

  314.                  twini_u,trini_u)                                  
    315. 

  315. 

  316.        first=.false.
    317. 

  317. 

  318.       endif ! first
    319. 

  319.       
    320. 

  320.       do ix=its,ite
    321. 

  321.       do iy=jts,jte
    322. 

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

  323. 	
    324. 

  324.          iurb=UTYPE_URB2D(ix,iy)
    325. 

  325. 

  326.        do iz= kts,kte
    327. 

  327.           ua1D(iz)=u_phy(ix,iz,iy)
    328. 

  328.           va1D(iz)=v_phy(ix,iz,iy)
    329. 

  329. 	  pt1D(iz)=th_phy(ix,iz,iy)
    330. 

  330. 	  da1D(iz)=rho(ix,iz,iy)
    331. 

  331. 	  pr1D(iz)=p_phy(ix,iz,iy)
    332. 

  332. !	  pt01D(iz)=th_phy(ix,iz,iy)
    333. 

  333. 	  pt01D(iz)=300.
    334. 

  334. 	  z1D(iz)=z(ix,iz,iy)
    335. 

  335.           a_u1D(iz)=0.
    336. 

  336.           a_v1D(iz)=0.
    337. 

  337.           a_t1D(iz)=0.
    338. 

  338.           a_e1D(iz)=0.
    339. 

  339.           b_u1D(iz)=0.
    340. 

  340.           b_v1D(iz)=0.
    341. 

  341.           b_t1D(iz)=0.
    342. 

  342.           b_e1D(iz)=0.           
    343. 

  343.          enddo
    344. 

  344. 	 z1D(kte+1)=z(ix,kte+1,iy)
    345. 

  345. 

  346.          do id=1,ndm
    347. 

  347.          do iz_u=1,nz_um
    348. 

  348.          do iw=1,nwr_u
    349. 

  349. !          tw1D(2*id-1,iz_u,iw)=tw1_u(ix,iy,ind_zwd(iz_u,iw,id))
    350. 

  350. !          tw1D(2*id,iz_u,iw)=tw2_u(ix,iy,ind_zwd(iz_u,iw,id))
    351. 

  351.             if(ind_zwd(iz_u,iw,id).gt.num_urban_layers)write(*,*)'ind_zwd too big w',ind_zwd(iz_u,iw,id)
    352. 

  352.           tw1D(2*id-1,iz_u,iw)=tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
    353. 

  353.           tw1D(2*id,iz_u,iw)=tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
    354. 

  354.          enddo
    355. 

  355.          enddo
    356. 

  356.          enddo
    357. 

  357. 	
    358. 

  358.          do id=1,ndm
    359. 

  359.           do ig=1,ng_u
    360. 

  360. !           tg1D(id,ig)=tg_u(ix,iy,ind_gd(ig,id))
    361. 

  361.            tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
    362. 

  362.           enddo
    363. 

  363.           do iz_u=1,nz_um
    364. 

  364.           do ir=1,nwr_u
    365. 

  365. !           tr1D(id,iz_u,ir)=tr_u(ix,iy,ind_zwd(iz_u,ir,id))
    366. 

  366.             if(ind_zwd(iz_u,ir,id).gt.num_urban_layers)write(*,*)'ind_zwd too big r',ind_zwd(iz_u,ir,id)
    367. 

  367.            tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)
    368. 

  368.           enddo
    369. 

  369.           enddo
    370. 

  370.          enddo
    371. 

  371. 

  372.          do id=1,ndm
    373. 

  373. 	 do iz=1,nz_um
    374. 

  374. !	  sfw1D(2*id-1,iz)=sfw1(ix,iy,ind_zd(iz,id))
    375. 

  375. !	  sfw1D(2*id,iz)=sfw2(ix,iy,ind_zd(iz,id))
    376. 

  376. 	  sfw1D(2*id-1,iz)=sfw1_urb3d(ix,ind_zd(iz,id),iy)
    377. 

  377. 	  sfw1D(2*id,iz)=sfw2_urb3d(ix,ind_zd(iz,id),iy)
    378. 

  378. 	 enddo
    379. 

  379. 	 enddo
    380. 

  380. 	 
    381. 

  381. 	 do id=1,ndm
    382. 

  382. !	  sfg1D(id)=sfg(ix,iy,id)
    383. 

  383. 	  sfg1D(id)=sfg_urb3d(ix,id,iy)
    384. 

  384. 	 enddo
    385. 

  385. 	 
    386. 

  386. 	 do id=1,ndm
    387. 

  387. 	 do iz=1,nz_um
    388. 

  388. !	  sfr1D(id,iz)=sfr(ix,iy,ind_zd(iz,id))
    389. 

  389. 	  sfr1D(id,iz)=sfr_urb3d(ix,ind_zd(iz,id),iy)
    390. 

  390. 	 enddo
    391. 

  391. 	 enddo
    392. 

  392. 

  393.          
    394. 

  394.          rs1D=swdown(ix,iy)
    395. 

  395.          rld1D=glw(ix,iy)
    396. 

  396.          time_h=(itimestep*dt)/3600.+gmt
    397. 

  397. 

  398.          zr1D=acos(COSZ_URB2D(ix,iy))
    399. 

  399.          deltar1D=DECLIN_URB
    400. 

  400.          ah1D=OMG_URB2D(ix,iy)
    401. 

  401. !	 call angle(xlong(ix,iy),xlat(ix,iy),julday,time_h,zr1D,deltar1D,ah1D)
    402. 

  402. 

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

  404.                    zr1D,deltar1D,ah1D,rs1D,rld1D,                   & 
    405. 

  405.                    alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,          & 
    406. 

  406.                    albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,          & 
    407. 

  407.                    fww,fwg,fgw,fsw,fws,fsg,                         & 
    408. 

  408.                    z0g_u,z0r_u,                                     & 
    409. 

  409.                    nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  & 
    410. 

  410.                    nz_u,z_u,                                        & 
    411. 

  411.                    tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D,                & 
    412. 

  412.                    a_u1D,a_v1D,a_t1D,a_e1D,                         & 
    413. 

  413.                    b_u1D,b_v1D,b_t1D,b_e1D,                         & 
    414. 

  414.                    dlg1D,dl_u1D,tsk1D,sf1D,vl1D,rl_up(ix,iy),       &
    415. 

  415.                    rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy))                            
    416. 

  416. 

  417.          do id=1,ndm
    418. 

  418. 	 do iz=1,nz_um
    419. 

  419. 	  sfw1_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id-1,iz) 
    420. 

  420. 	  sfw2_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id,iz) 
    421. 

  421. 	 enddo
    422. 

  422. 	 enddo
    423. 

  423.  
    424. 

  424. 	 do id=1,ndm
    425. 

  425. 	  sfg_urb3d(ix,id,iy)=sfg1D(id) 
    426. 

  426. 	 enddo
    427. 

  427. 	
    428. 

  428. 	 do id=1,ndm
    429. 

  429. 	 do iz=1,nz_um
    430. 

  430. 	  sfr_urb3d(ix,ind_zd(iz,id),iy)=sfr1D(id,iz)
    431. 

  431. 	 enddo
    432. 

  432. 	 enddo
    433. 

  433. !
    434. 

  434.          do id=1,ndm
    435. 

  435.          do iz_u=1,nz_um
    436. 

  436.          do iw=1,nwr_u
    437. 

  437.           tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw)
    438. 

  438.           tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw)
    439. 

  439.          enddo
    440. 

  440.          enddo
    441. 

  441.          enddo
    442. 

  442.         
    443. 

  443.          do id=1,ndm
    444. 

  444.           do ig=1,ng_u
    445. 

  445.            tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
    446. 

  446.           enddo
    447. 

  447.           do iz_u=1,nz_um
    448. 

  448.           do ir=1,nwr_u
    449. 

  449.            trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
    450. 

  450.           enddo
    451. 

  451.           enddo
    452. 

  452.          enddo
    453. 

  453. 

  454.         do iz= kts,kte
    455. 

  455.           sf(ix,iz,iy)=sf1D(iz)
    456. 

  456.           vl(ix,iz,iy)=vl1D(iz)
    457. 

  457.           a_u(ix,iz,iy)=a_u1D(iz)
    458. 

  458.           a_v(ix,iz,iy)=a_v1D(iz)
    459. 

  459.           a_t(ix,iz,iy)=a_t1D(iz)
    460. 

  460.           a_e(ix,iz,iy)=a_e1D(iz)
    461. 

  461.           b_u(ix,iz,iy)=b_u1D(iz)
    462. 

  462.           b_v(ix,iz,iy)=b_v1D(iz)
    463. 

  463.           b_t(ix,iz,iy)=b_t1D(iz)
    464. 

  464.           b_e(ix,iz,iy)=b_e1D(iz)
    465. 

  465.           dlg(ix,iz,iy)=dlg1D(iz)
    466. 

  466.           dl_u(ix,iz,iy)=dl_u1D(iz)
    467. 

  467.          enddo
    468. 

  468.          sf(ix,kte+1,iy)=sf1D(kte+1)
    469. 

  469. !         tsk(ix,iy)=tsk1D
    470. 

  470. !
    471. 

  471.          endif ! FRC_URB2D
    472. 

  472.    
    473. 

  473.       enddo  ! iy
    474. 

  474.       enddo  ! ix
    475. 

  475. 

  476. 

  477.         time_bep=time_bep+dt
    478. 

  478.          
    479. 

  479.          
    480. 

  480.       return
    481. 

  481.       end subroutine BEP
    482. 

  482.             
    483. 

  483. ! ===6=8===============================================================72
    484. 

  484. 

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

  486.                       zr,deltar,ah,rs,rld,                            & 
    487. 

  487.                       alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,         & 
    488. 

  488.                       albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,         & 
    489. 

  489.                       fww,fwg,fgw,fsw,fws,fsg,                        & 
    490. 

  490.                       z0g_u,z0r_u,                                    & 
    491. 

  491.                       nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, & 
    492. 

  492.                       nz_u,z_u,                                       & 
    493. 

  493.                       tw,tg,tr,sfw,sfg,sfr,                           & 
    494. 

  494.                       a_u,a_v,a_t,a_e,                                &
    495. 

  495.                       b_u,b_v,b_t,b_e,                                & 
    496. 

  496.                       dlg,dl_u,tsk,sf,vl,rl_up,rs_abs,emiss,grdflx_urb)                             
    497. 

  497. 

  498. ! ----------------------------------------------------------------------
    499. 

  499. ! This routine computes the effects of buildings on momentum, heat and
    500. 

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

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

  502. !  mesoscale grid defined by the altitude of its cell interfaces "z" and
    503. 

  503. !  its number of levels "nz".
    504. 

  504. ! The meteorological input parameters (wind, temperature, solar radiation)
    505. 

  505. !  are specified on the "mesoscale grid".
    506. 

  506. ! The inputs concerning the building and street charateristics are defined
    507. 

  507. !  on a "urban grid". The "urban grid" is defined with its number of levels
    508. 

  508. !  "nz_u" and its space step "dz_u".
    509. 

  509. ! The input parameters are interpolated on the "urban grid". The sources or sinks
    510. 

  510. !  are calculated on the "urban grid". Finally the sources or sinks are 
    511. 

  511. !  interpolated on the "mesoscale grid".
    512. 

  512.  
    513. 

  513. 

  514. !  Mesoscale grid            Urban grid             Mesoscale grid
    515. 

  515. !  
    516. 

  516. ! z(4)    ---                                               ---
    517. 

  517. !          |                                                 |
    518. 

  518. !          |                                                 |
    519. 

  519. !          |   Interpolation                  Interpolation  |
    520. 

  520. !          |            Sources or sinks calculation         |
    521. 

  521. ! z(3)    ---                                               ---
    522. 

  522. !          | ua               ua_u  ---  uv_a         a_u    |
    523. 

  523. !          | va               va_u   |   uv_b         b_u    |
    524. 

  524. !          | pt               pt_u  ---  uh_b         a_v    |
    525. 

  525. ! z(2)    ---                        |    etc...      etc...---
    526. 

  526. !          |                 z_u(1) ---                      |
    527. 

  527. !          |                         |                       |
    528. 

  528. ! z(1) ------------------------------------------------------------
    529. 

  529. 

  530. !     
    531. 

  531. ! Reference:
    532. 

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

  533. ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
    534. 

  534. ! 261-304
    535. 

  535.  
    536. 

  536. ! ----------------------------------------------------------------------
    537. 

  537. 

  538.       implicit none
    539. 

  539. 

  540.  
    541. 

  541. 

  542. ! ----------------------------------------------------------------------
    543. 

  543. ! INPUT:
    544. 

  544. ! ----------------------------------------------------------------------
    545. 

  545. 

  546. ! Data relative to the "mesoscale grid"
    547. 

  547. 

  548. !      integer nz                 ! Number of vertical levels
    549. 

  549.       integer kms,kme,kts,kte
    550. 

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

  551.       real ua(kms:kme)                ! Wind speed in the x direction
    552. 

  552.       real va(kms:kme)                ! Wind speed in the y direction
    553. 

  553.       real pt(kms:kme)                ! Potential temperature
    554. 

  554.       real da(kms:kme)                ! Air density
    555. 

  555.       real pr(kms:kme)                ! Air pressure
    556. 

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

  557.       real dt                    ! Time step
    558. 

  558.       real zr                    ! Zenith angle
    559. 

  559.       real deltar                ! Declination of the sun
    560. 

  560.       real ah                    ! Hour angle
    561. 

  561.       real rs                    ! Solar radiation
    562. 

  562.       real rld                   ! Downward flux of the longwave radiation
    563. 

  563. 

  564. ! Data relative to the "urban grid"
    565. 

  565. 

  566.       integer iurb               ! Current urban class
    567. 

  567. 

  568. !    Building parameters      
    569. 

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

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

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

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

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

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

  575. 

  576. !    Radiation parameters
    577. 

  577.       real albg_u(nurbm)         ! Albedo of the ground
    578. 

  578.       real albw_u(nurbm)         ! Albedo of the wall
    579. 

  579.       real albr_u(nurbm)         ! Albedo of the roof
    580. 

  580.       real emg_u(nurbm)          ! Emissivity of ground
    581. 

  581.       real emw_u(nurbm)          ! Emissivity of wall
    582. 

  582.       real emr_u(nurbm)          ! Emissivity of roof
    583. 

  583. 

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

  585. !    short wave radation. 
    586. 

  586. !    The calculation of these factor is explained in the Appendix A of the BLM paper
    587. 

  587.       real fww(nz_um,nz_um,ndm,nurbm)  !  from wall to wall
    588. 

  588.       real fwg(nz_um,ndm,nurbm)        !  from wall to ground
    589. 

  589.       real fgw(nz_um,ndm,nurbm)        !  from ground to wall
    590. 

  590.       real fsw(nz_um,ndm,nurbm)        !  from sky to wall
    591. 

  591.       real fws(nz_um,ndm,nurbm)        !  from wall to sky
    592. 

  592.       real fsg(ndm,nurbm)              !  from sky to ground
    593. 

  593. 

  594. !    Roughness parameters
    595. 

  595.       real z0g_u(nurbm)            ! The ground's roughness length
    596. 

  596.       real z0r_u(nurbm)            ! The roof's roughness length
    597. 

  597.       
    598. 

  598. !    Street parameters
    599. 

  599.       integer nd_u(nurbm)          ! Number of street direction for each urban class 
    600. 

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

  601.       real drst_u(ndm,nurbm)       ! Street direction
    602. 

  602.       real ws_u(ndm,nurbm)         ! Street width
    603. 

  603.       real bs_u(ndm,nurbm)         ! Building width
    604. 

  604.       real h_b(nz_um,nurbm)        ! Bulding's heights
    605. 

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

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

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

  608.         
    609. 

  609. !    Grid parameters
    610. 

  610.       integer nz_u(nurbm)     ! Number of layer in the urban grid
    611. 

  611. !      real dz_u               ! Urban grid resolution
    612. 

  612.       real z_u(nz_um)       ! Height of the urban grid levels
    613. 

  613. 

  614. 

  615. ! ----------------------------------------------------------------------
    616. 

  616. ! INPUT-OUTPUT
    617. 

  617. ! ----------------------------------------------------------------------
    618. 

  618. 

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

  620. 

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

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

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

  624.       real sfw(2*ndm,nz_um)      ! Sensible heat flux from walls
    625. 

  625.       real sfg(ndm)              ! Sensible heat flux from ground (road)
    626. 

  626.       real sfr(ndm,nz_um)      ! Sensible heat flux from roofs
    627. 

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

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

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

  630. ! ----------------------------------------------------------------------
    631. 

  631. ! OUTPUT:
    632. 

  632. ! ----------------------------------------------------------------------
    633. 

  633. 

  634. ! Data relative to the "mesoscale grid"
    635. 

  635. 

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

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

  638.      
    639. 

  639. !    Implicit and explicit components of the source and sink terms at each levels,
    640. 

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

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

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

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

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

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

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

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

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

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

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

  651.       real tsk                  ! Average of the road surface temperatures
    652. 

  652.       
    653. 

  653. ! ----------------------------------------------------------------------
    654. 

  654. ! LOCAL:
    655. 

  655. ! ----------------------------------------------------------------------
    656. 

  656. 

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

  658. 

  659. ! Data interpolated from the "mesoscale grid" to the "urban grid"
    660. 

  660. 

  661.       real ua_u(nz_um)          ! Wind speed in the x direction
    662. 

  662.       real va_u(nz_um)          ! Wind speed in the y direction
    663. 

  663.       real pt_u(nz_um)          ! Potential temperature
    664. 

  664.       real da_u(nz_um)          ! Air density
    665. 

  665.       real pt0_u(nz_um)         ! Reference potential temperature
    666. 

  666.       real pr_u(nz_um)          ! Air pressure
    667. 

  667. 

  668. ! Data defining the building and street charateristics
    669. 

  669. 

  670.       integer nd                ! Number of street direction for the current urban class
    671. 

  671. 

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

  673.       real alar(nwr_u)           ! Roof thermal diffusivity for the current urban class [m^2 s^-1]
    674. 

  674.       real alaw(nwr_u)           ! Walls thermal diffusivity for the current urban class [m^2 s^-1] 
    675. 

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

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

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

  678. 

  679.       real z0(ndm,nz_um)      ! Roughness lengths "profiles"
    680. 

  680.       real ws(ndm)              ! Street widths of the current urban class
    681. 

  681.       real bs(ndm)              ! Building widths of the current urban class
    682. 

  682.       real strd(ndm)            ! Street lengths for the current urban class
    683. 

  683.       real drst(ndm)            ! Street directions for the current urban class
    684. 

  684.       real ss(nz_um)          ! Probability to have a building with height h
    685. 

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

  686. 

  687. ! Solar radiation at each level of the "urban grid"
    688. 

  688. 

  689.       real rsg(ndm)             ! Short wave radiation from the ground
    690. 

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

  691.       real rlg(ndm)             ! Long wave radiation from the ground
    692. 

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

  693. 

  694. ! Potential temperature of the surfaces at each level of the "urban grid"
    695. 

  695. 

  696.       real ptg(ndm)             ! Ground potential temperatures 
    697. 

  697.       real ptr(ndm,nz_um)     ! Roof potential temperatures 
    698. 

  698.       real ptw(2*ndm,nz_um)     ! Walls potential temperatures 
    699. 

  699. 

  700.  
    701. 

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

  702. !  vertical surfaces (walls) ans horizontal surfaces (roofs and street)
    703. 

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

  704. !  Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
    705. 

  705. 

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

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

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

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

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

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

  712.       real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
    713. 

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

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

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

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

  717.       
    718. 

  718. !
    719. 

  719.       real rs_abs ! solar radiation absorbed by urban surfaces 
    720. 

  720.       real rl_up ! longwave radiation emitted by urban surface to the atmosphere 
    721. 

  721.       real emiss ! mean emissivity of the urban surface
    722. 

  722.       real grdflx_urb ! ground heat flux
    723. 

  723.       real shtot,aaa
    724. 

  724.       real dt_int ! internal time step
    725. 

  725.       integer nt_int ! number of internal time step
    726. 

  726.       integer iz,id, it_int
    727. 

  727.       integer iwrong,iw,ix,iy
    728. 

  728. 

  729. ! ----------------------------------------------------------------------
    730. 

  730. ! END VARIABLES DEFINITIONS
    731. 

  731. ! ----------------------------------------------------------------------
    732. 

  732. 

  733. ! Fix some usefull parameters for the computation of the sources or sinks
    734. 

  734. 

  735.       do iz=kts,kte
    736. 

  736.          dz(iz)=z(iz+1)-z(iz)
    737. 

  737.       end do
    738. 

  738.       call param(iurb,nz_u(iurb),nd_u(iurb),         &
    739. 

  739.                        csg_u,csg,alag_u,alag,csr_u,csr,             &
    740. 

  740.                        alar_u,alar,csw_u,csw,alaw_u,alaw,           &
    741. 

  741.                        ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0,              &
    742. 

  742.                        strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb)
    743. 

  743. 

  744. ! Interpolation on the "urban grid"
    745. 

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

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

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

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

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

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

  751. 

  752.                    
    753. 

  753. ! Compute the modification of the radiation due to the buildings
    754. 

  754. 

  755.       call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws, &
    756. 

  756.                     drst,strd,ss,pb,                              &
    757. 

  757.                     tw,tg,albg_u(iurb),albw_u(iurb),              &
    758. 

  758.                     emw_u(iurb),emg_u(iurb),                      &
    759. 

  759.                     fww,fwg,fgw,fsw,fsg,                          &
    760. 

  760.                     zr,deltar,ah,                                 &
    761. 

  761.                     rs,rld,rsw,rsg,rlw,rlg)                       
    762. 

  762.  
    763. 

  763. ! calculation of the urban albedo and the upward long wave radiation
    764. 

  764.        call upward_rad(nd_u(iurb),iurb,nz_u(iurb),ws,bs,sigma,fsw,fsg,pb,ss,   &
    765. 

  765.                        tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg,                & 
    766. 

  766.                        tw,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw,                & 
    767. 

  767.                        tr,emr_u(iurb),albr_u(iurb),rld,rs,sfr,                 & 
    768. 

  768.                        rs_abs,rl_up,emiss,grdflx_urb)               
    769. 

  769.         
    770. 

  770. ! Compute the surface temperatures
    771. 

  771.      
    772. 

  772. 

  773.       call surf_temp(nz_u(iurb),nd_u(iurb),pr_u,dt,ss,                  & 
    774. 

  774.                     rs,rld,rsg,rlg,rsw,rlw,                             &
    775. 

  775.                     tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg,  &
    776. 

  776.                     tr,alar,csr,emr_u(iurb),albr_u(iurb),ptr,sfr,gfr,  &
    777. 

  777.                     tw,alaw,csw,emw_u(iurb),albw_u(iurb),ptw,sfw,gfw)  
    778. 

  778.       
    779. 

  779.       
    780. 

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

  781.        
    782. 

  782.       call buildings(nd_u(iurb),nz_u(iurb),z0,ua_u,va_u,                & 
    783. 

  783.                      pt_u,pt0_u,ptg,ptr,da_u,ptw,drst,                  &                      
    784. 

  784.                      uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,         & 
    785. 

  785.                      uhb_u,vhb_u,thb_u,ehb_u,ss,dt)                        
    786. 

  786.  
    787. 

  787.          
    788. 

  788. ! Calculation of the sensible heat fluxes for the ground, the wall and roof
    789. 

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

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

  791. !  
    792. 

  792. ! 
    793. 

  793.  
    794. 

  794.       do id=1,nd_u(iurb)         
    795. 

  795.          sfg(id)=-da_u(1)*cp_u*thb_u(id,1)
    796. 

  796.          do iz=2,nz_u(iurb)
    797. 

  797.             sfr(id,iz)=-da_u(iz)*cp_u*thb_u(id,iz)
    798. 

  798.          enddo
    799. 

  799.          
    800. 

  800.          do iz=1,nz_u(iurb)
    801. 

  801.             sfw(2*id-1,iz)=-da_u(iz)*cp_u*(tvb_u(2*id-1,iz)+     &
    802. 

  802.                 tva_u(2*id-1,iz)*pt_u(iz))
    803. 

  803.             sfw(2*id,iz)=-da_u(iz)*cp_u*(tvb_u(2*id,iz)+         &
    804. 

  804.                 tva_u(2*id,iz)*pt_u(iz)) 
    805. 

  805.          enddo
    806. 

  806.       enddo
    807. 

  807.       
    808. 

  808. ! calculation of the urban albedo and the upward long wave radiation
    809. 

  809.           
    810. 

  810. !       call upward_rad(nd_u(iurb),iurb,nz_u(iurb),ws,bs,sigma,fsw,fsg,pb,ss,  &
    811. 

  811. !                       tg,emg_u(iurb),albg_u(iurb),rlg,rsg,                & 
    812. 

  812. !                       tw,emw_u(iurb),albw_u(iurb),rlw,rsw,                & 
    813. 

  813. !                       tr,emr_u(iurb),albr_u(iurb),rld,rs,                  & 
    814. 

  814. !                       rs_abs,rl_up,emiss)             
    815. 

  815. 

  816. ! Interpolation on the "mesoscale grid"
    817. 

  817. 

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

  819.                      vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,       &
    820. 

  820.                      uhb_u,vhb_u,thb_u,ehb_u,                            &
    821. 

  821.                      a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)                    
    822. 

  822.        
    823. 

  823. 

  824. ! computation of the mean road temperature tsk (this value could be used 
    825. 

  825. ! to replace the surface temperature in the radiation routines, if needed).
    826. 

  826. 

  827. !      tsk=0.
    828. 

  828. !      do id=1,nd_u(iurb)
    829. 

  829. !         tsk=tsk+tg(id,ng_u)/nd_u(iurb)
    830. 

  830. !      enddo
    831. 

  831. 

  832. ! Calculation of the length scale taking into account the buildings effects
    833. 

  833. 

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

  835. 

  836.       return
    837. 

  837.       end subroutine BEP1D
    838. 

  838. 

  839. ! ===6=8===============================================================72
    840. 

  840. ! ===6=8===============================================================72
    841. 

  841. 

  842.        subroutine param(iurb,nz,nd,                                   &
    843. 

  843.                        csg_u,csg,alag_u,alag,csr_u,csr,               &
    844. 

  844.                        alar_u,alar,csw_u,csw,alaw_u,alaw,             &
    845. 

  845.                        ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0,                &  
    846. 

  846.                        strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb)        
    847. 

  847. 

  848. ! ----------------------------------------------------------------------
    849. 

  849. !    This routine prepare some usefull parameters       
    850. 

  850. ! ----------------------------------------------------------------------
    851. 

  851. 

  852.       implicit none
    853. 

  853. 

  854.   
    855. 

  855. ! ----------------------------------------------------------------------
    856. 

  856. ! INPUT:
    857. 

  857. ! ----------------------------------------------------------------------
    858. 

  858.       integer iurb                 ! Current urban class
    859. 

  859.       integer nz                   ! Number of vertical urban levels in the current class
    860. 

  860.       integer nd                   ! Number of street direction for the current urban class
    861. 

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

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

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

  864.       real bs_u(ndm,nurbm)         ! Building width
    865. 

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

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

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

  868.       real drst_u(ndm,nurbm)       ! Street direction
    869. 

  869.       real strd_u(ndm,nurbm)       ! Street length 
    870. 

  870.       real ws_u(ndm,nurbm)         ! Street width
    871. 

  871.       real z0g_u(nurbm)            ! The ground's roughness length
    872. 

  872.       real z0r_u(nurbm)            ! The roof's roughness length
    873. 

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

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

  875. 

  876. ! ----------------------------------------------------------------------
    877. 

  877. ! OUTPUT:
    878. 

  878. ! ----------------------------------------------------------------------
    879. 

  879.       real alag(ng_u)           ! Ground thermal diffusivity at each ground levels
    880. 

  880.       real alar(nwr_u)           ! Roof thermal diffusivity at each roof levels
    881. 

  881.       real alaw(nwr_u)           ! Wall thermal diffusivity at each wall levels
    882. 

  882.       real csg(ng_u)            ! Specific heat of the ground material at each ground levels
    883. 

  883.       real csr(nwr_u)            ! Specific heat of the roof material at each roof levels
    884. 

  884.       real csw(nwr_u)            ! Specific heat of the wall material at each wall levels
    885. 

  885.       real bs(ndm)              ! Building width for the current urban class
    886. 

  886.       real drst(ndm)            ! street directions for the current urban class
    887. 

  887.       real strd(ndm)            ! Street lengths for the current urban class
    888. 

  888.       real ws(ndm)              ! Street widths of the current urban class
    889. 

  889.       real z0(ndm,nz_um)      ! Roughness lengths "profiles"
    890. 

  890.       real ss(nz_um)          ! Probability to have a building with height h
    891. 

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

  892. 

  893. ! ----------------------------------------------------------------------
    894. 

  894. ! LOCAL:
    895. 

  895. ! ----------------------------------------------------------------------
    896. 

  896.       integer id,ig,ir,iw,iz
    897. 

  897. 

  898. ! ----------------------------------------------------------------------
    899. 

  899. ! END VARIABLES DEFINITIONS
    900. 

  900. ! ----------------------------------------------------------------------     
    901. 

  901. !
    902. 

  902. !Initialize the variables
    903. 

  903. !
    904. 

  904.       ss=0.
    905. 

  905.       pb=0.
    906. 

  906.       csg=0.
    907. 

  907.       alag=0.
    908. 

  908.       csr=0.
    909. 

  909.       alar=0.
    910. 

  910.       csw=0.
    911. 

  911.       alaw=0.
    912. 

  912.       z0=0.
    913. 

  913.       ws=0.
    914. 

  914.       bs=0.
    915. 

  915.       strd=0.
    916. 

  916.       drst=0.
    917. 

  917.       
    918. 

  918.       do iz=1,nz+1
    919. 

  919.          ss(iz)=ss_u(iz,iurb)
    920. 

  920.          pb(iz)=pb_u(iz,iurb)
    921. 

  921.       end do
    922. 

  922.                  
    923. 

  923.        do ig=1,ng_u
    924. 

  924.         csg(ig)=csg_u(iurb)
    925. 

  925.         alag(ig)=alag_u(iurb)
    926. 

  926.        enddo
    927. 

  927.        
    928. 

  928.        do ir=1,nwr_u
    929. 

  929.         csr(ir)=csr_u(iurb)
    930. 

  930.         alar(ir)=alar_u(iurb)
    931. 

  931.        enddo
    932. 

  932.        
    933. 

  933.        do iw=1,nwr_u
    934. 

  934.         csw(iw)=csw_u(iurb)
    935. 

  935.         alaw(iw)=alaw_u(iurb)
    936. 

  936.        enddo
    937. 

  937.        
    938. 

  938.        do id=1,nd
    939. 

  939.         z0(id,1)=z0g_u(iurb)
    940. 

  940.         do iz=2,nz+1
    941. 

  941.          z0(id,iz)=z0r_u(iurb)
    942. 

  942.         enddo
    943. 

  943.        enddo
    944. 

  944.       
    945. 

  945.        do id=1,nd
    946. 

  946.         ws(id)=ws_u(id,iurb)
    947. 

  947.         bs(id)=bs_u(id,iurb)
    948. 

  948.         strd(id)=strd_u(id,iurb)
    949. 

  949.         drst(id)=drst_u(id,iurb)     
    950. 

  950.        enddo
    951. 

  951. 

  952.        
    953. 

  953.        return
    954. 

  954.        end subroutine param
    955. 

  955.        
    956. 

  956. ! ===6=8===============================================================72
    957. 

  957. ! ===6=8===============================================================72
    958. 

  958. 

  959.       subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
    960. 

  960. 

  961. ! ----------------------------------------------------------------------
    962. 

  962. !  This routine interpolate para
    963. 

  963. !  meters from the "mesoscale grid" to
    964. 

  964. !  the "urban grid".
    965. 

  965. !  See p300 Appendix B.1 of the BLM paper.
    966. 

  966. ! ----------------------------------------------------------------------
    967. 

  967. 

  968.       implicit none
    969. 

  969. 

  970. ! ----------------------------------------------------------------------
    971. 

  971. ! INPUT:
    972. 

  972. ! ----------------------------------------------------------------------
    973. 

  973. ! Data relative to the "mesoscale grid"
    974. 

  974.       integer kts,kte,kms,kme            
    975. 

  975.       real z(kms:kme)          ! Altitude of the cell interface
    976. 

  976.       real c(kms:kme)            ! Parameter which has to be interpolated
    977. 

  977. ! Data relative to the "urban grid"
    978. 

  978.       integer nz_u          ! Number of levels
    979. 

  979. !!    real z_u(nz_u+1)      ! Altitude of the cell interface
    980. 

  980.       real z_u(nz_um)      ! Altitude of the cell interface
    981. 

  981. ! ----------------------------------------------------------------------
    982. 

  982. ! OUTPUT:
    983. 

  983. ! ----------------------------------------------------------------------
    984. 

  984. !!    real c_u(nz_u)        ! Interpolated paramters in the "urban grid"
    985. 

  985.       real c_u(nz_um)        ! Interpolated paramters in the "urban grid"
    986. 

  986. 

  987. ! LOCAL:
    988. 

  988. ! ----------------------------------------------------------------------
    989. 

  989.       integer iz_u,iz
    990. 

  990.       real ctot,dz
    991. 

  991. 

  992. ! ----------------------------------------------------------------------
    993. 

  993. ! END VARIABLES DEFINITIONS
    994. 

  994. ! ----------------------------------------------------------------------
    995. 

  995. 

  996.        do iz_u=1,nz_u
    997. 

  997.         ctot=0.
    998. 

  998.         do iz=kts,kte
    999. 

  999.          dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
    1000. 

  1000.          ctot=ctot+c(iz)*dz
    1001. 

  1001.         enddo
    1002. 

  1002.         c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
    1003. 

  1003.        enddo
    1004. 

  1004.        
    1005. 

  1005.        return
    1006. 

  1006.        end subroutine interpol
    1007. 

  1007.          
    1008. 

  1008. ! ===6=8===============================================================72       
    1009. 

  1009. ! ===6=8===============================================================72     
    1010. 

  1010. 

  1011.       subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb,    &
    1012. 

  1012.                           tw,tg,albg,albw,emw,emg,                     &
    1013. 

  1013.                           fww,fwg,fgw,fsw,fsg,                         &
    1014. 

  1014.                           zr,deltar,ah,                                &    
    1015. 

  1015.                           rs,rl,rsw,rsg,rlw,rlg)                       
    1016. 

  1016.  
    1017. 

  1017. ! ----------------------------------------------------------------------
    1018. 

  1018. ! This routine computes the modification of the short wave and 
    1019. 

  1019. !  long wave radiation due to the buildings.
    1020. 

  1020. ! ----------------------------------------------------------------------
    1021. 

  1021. 

  1022.       implicit none
    1023. 

  1023.  
    1024. 

  1024.  
    1025. 

  1025. ! ----------------------------------------------------------------------
    1026. 

  1026. ! INPUT:
    1027. 

  1027. ! ----------------------------------------------------------------------
    1028. 

  1028.       integer iurb              ! current urban class
    1029. 

  1029.       integer nd                ! Number of street direction for the current urban class
    1030. 

  1030.       integer nz_u              ! Number of layer in the urban grid
    1031. 

  1031.       real z(nz_um)           ! Height of the urban grid levels
    1032. 

  1032.       real ws(ndm)              ! Street widths of the current urban class
    1033. 

  1033.       real drst(ndm)            ! street directions for the current urban class
    1034. 

  1034.       real strd(ndm)            ! Street lengths for the current urban class
    1035. 

  1035.       real ss(nz_um)          ! probability to have a building with height h
    1036. 

  1036.       real pb(nz_um)          ! probability to have a building with an height equal
    1037. 

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

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

  1039.       real albg                 ! Albedo of the ground for the current urban class
    1040. 

  1040.       real albw                 ! Albedo of the wall for the current urban class
    1041. 

  1041.       real emg                  ! Emissivity of ground for the current urban class
    1042. 

  1042.       real emw                  ! Emissivity of wall for the current urban class
    1043. 

  1043.       real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
    1044. 

  1044.       real fsg(ndm,nurbm)             ! View factors from sky to ground
    1045. 

  1045.       real fsw(nz_um,ndm,nurbm)       ! View factors from sky to wall
    1046. 

  1046.       real fws(nz_um,ndm,nurbm)       ! View factors from wall to sky
    1047. 

  1047.       real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
    1048. 

  1048.       real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
    1049. 

  1049.       real ah                   ! Hour angle (it should come from the radiation routine)
    1050. 

  1050.       real zr                   ! zenith angle
    1051. 

  1051.       real deltar               ! Declination of the sun
    1052. 

  1052.       real rs                   ! solar radiation
    1053. 

  1053.       real rl                   ! downward flux of the longwave radiation
    1054. 

  1054. ! ----------------------------------------------------------------------
    1055. 

  1055. ! OUTPUT:
    1056. 

  1056. ! ----------------------------------------------------------------------
    1057. 

  1057.       real rlg(ndm)             ! Long wave radiation at the ground
    1058. 

  1058.       real rlw(2*ndm,nz_um)     ! Long wave radiation at the walls
    1059. 

  1059.       real rsg(ndm)             ! Short wave radiation at the ground
    1060. 

  1060.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    1061. 

  1061. 

  1062. ! ----------------------------------------------------------------------
    1063. 

  1063. ! LOCAL:
    1064. 

  1064. ! ----------------------------------------------------------------------
    1065. 

  1065. 

  1066.       integer id,iz
    1067. 

  1067. 

  1068. !  Calculation of the shadow effects
    1069. 

  1069. 

  1070.       call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z,  &
    1071. 

  1071.                      rs,rsw,rsg)
    1072. 

  1072. 

  1073. ! Calculation of the reflection effects          
    1074. 

  1074.       do id=1,nd
    1075. 

  1075.          call long_rad(iurb,nz_u,id,emw,emg,                 &
    1076. 

  1076.                       fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
    1077. 

  1077.          
    1078. 

  1078.          call short_rad(iurb,nz_u,id,albw,albg,fwg,fww,fgw,rsg,rsw,pb)
    1079. 

  1079.   
    1080. 

  1080.       enddo
    1081. 

  1081.       
    1082. 

  1082.       return
    1083. 

  1083.       end subroutine modif_rad
    1084. 

  1084. 

  1085. 

  1086. ! ===6=8===============================================================72  
    1087. 

  1087. ! ===6=8===============================================================72     
    1088. 

  1088. 

  1089.       subroutine surf_temp(nz_u,nd,pr,dt,ss,rs,rl,rsg,rlg,rsw,rlw,     &
    1090. 

  1090.                           tg,alag,csg,emg,albg,ptg,sfg,gfg,             &
    1091. 

  1091.                           tr,alar,csr,emr,albr,ptr,sfr,gfr,             &
    1092. 

  1092.                           tw,alaw,csw,emw,albw,ptw,sfw,gfw)             
    1093. 

  1093.                  
    1094. 

  1094. 

  1095. ! ----------------------------------------------------------------------
    1096. 

  1096. ! Computation of the surface temperatures for walls, ground and roofs 
    1097. 

  1097. ! ----------------------------------------------------------------------
    1098. 

  1098. 

  1099.       implicit none
    1100. 

  1100.       
    1101. 

  1101.       
    1102. 

  1102.       
    1103. 

  1103. ! ----------------------------------------------------------------------
    1104. 

  1104. ! INPUT:
    1105. 

  1105. ! ----------------------------------------------------------------------
    1106. 

  1106.       integer nz_u              ! Number of vertical layers defined in the urban grid
    1107. 

  1107.       integer nd                ! Number of street direction for the current urban class
    1108. 

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

  1109.       real alar(nwr_u)           ! Roof thermal diffusivity for the current urban class [m^2 s^-1]  
    1110. 

  1110.       real alaw(nwr_u)           ! Wall thermal diffusivity for the current urban class [m^2 s^-1]  
    1111. 

  1111.       real albg                 ! Albedo of the ground for the current urban class
    1112. 

  1112.       real albr                 ! Albedo of the roof for the current urban class
    1113. 

  1113.       real albw                 ! Albedo of the wall for the current urban class
    1114. 

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

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

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

  1117.       real dt                   ! Time step
    1118. 

  1118.       real emg                  ! Emissivity of ground for the current urban class
    1119. 

  1119.       real emr                  ! Emissivity of roof for the current urban class
    1120. 

  1120.       real emw                  ! Emissivity of wall for the current urban class
    1121. 

  1121.       real pr(nz_um)            ! Air pressure
    1122. 

  1122.       real rs                   ! Solar radiation
    1123. 

  1123.       real rl                   ! Downward flux of the longwave radiation
    1124. 

  1124.       real rlg(ndm)             ! Long wave radiation at the ground
    1125. 

  1125.       real rlw(2*ndm,nz_um)     ! Long wave radiation at the walls
    1126. 

  1126.       real rsg(ndm)             ! Short wave radiation at the ground
    1127. 

  1127.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    1128. 

  1128.       real sfg(ndm)             ! Sensible heat flux from ground (road)
    1129. 

  1129.       real sfr(ndm,nz_um)     ! Sensible heat flux from roofs
    1130. 

  1130.       real sfw(2*ndm,nz_um)     ! Sensible heat flux from walls
    1131. 

  1131.       real gfg(ndm)             ! Heat flux transferred from the surface of the ground (road) toward the interior
    1132. 

  1132.       real gfr(ndm,nz_um)     ! Heat flux transferred from the surface of the roof toward the interior
    1133. 

  1133.       real gfw(2*ndm,nz_um)     ! Heat flux transfered from the surface of the walls toward the interior
    1134. 

  1134.       real ss(nz_um)          ! Probability to have a building with height h
    1135. 

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

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

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

  1138.       
    1139. 

  1139. 

  1140. ! ----------------------------------------------------------------------
    1141. 

  1141. ! OUTPUT:
    1142. 

  1142. ! ----------------------------------------------------------------------
    1143. 

  1143.       real ptg(ndm)             ! Ground potential temperatures 
    1144. 

  1144.       real ptr(ndm,nz_um)     ! Roof potential temperatures 
    1145. 

  1145.       real ptw(2*ndm,nz_um)     ! Walls potential temperatures 
    1146. 

  1146. 

  1147. ! ----------------------------------------------------------------------
    1148. 

  1148. ! LOCAL:
    1149. 

  1149. ! ----------------------------------------------------------------------
    1150. 

  1150.       integer id,ig,ir,iw,iz
    1151. 

  1151. 

  1152.       real rtg(ndm)             ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
    1153. 

  1153.       real rtr(ndm,nz_um)     ! Total radiation at roof surface (solar+incoming long+outgoing long)
    1154. 

  1154.       real rtw(2*ndm,nz_um)     ! Radiation at walls surface (solar+incoming long+outgoing long)
    1155. 

  1155.       real tg_tmp(ng_u)
    1156. 

  1156.       real tr_tmp(nwr_u)
    1157. 

  1157.       real tw_tmp(nwr_u)
    1158. 

  1158. 

  1159.       real dzg_u(ng_u)          ! Layer sizes in the ground
    1160. 

  1160. 

  1161.       real dzr_u(nwr_u)          ! Layers sizes in the roof
    1162. 

  1162.          
    1163. 

  1163.       real dzw_u(nwr_u)          ! Layer sizes in the wall
    1164. 

  1164.       
    1165. 

  1165. 

  1166.       data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
    1167. 

  1167.       data dzr_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/   
    1168. 

  1168.       data dzw_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
    1169. 

  1169. ! ----------------------------------------------------------------------
    1170. 

  1170. ! END VARIABLES DEFINITIONS
    1171. 

  1171. ! ----------------------------------------------------------------------
    1172. 

  1172. 

  1173.         
    1174. 

  1174.    
    1175. 

  1175.       do id=1,nd
    1176. 

  1176. 

  1177. !      Calculation for the ground surfaces
    1178. 

  1178.        do ig=1,ng_u
    1179. 

  1179.         tg_tmp(ig)=tg(id,ig)
    1180. 

  1180.        end do
    1181. 

  1181. !	        
    1182. 

  1182.        call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg,        &
    1183. 

  1183.                      rsg(id),rlg(id),pr(1),                    &
    1184. 

  1184.                      dt,emg,albg,                              &
    1185. 

  1185.                      rtg(id),sfg(id),gfg(id))    
    1186. 

  1186.        do ig=1,ng_u
    1187. 

  1187.         tg(id,ig)=tg_tmp(ig)
    1188. 

  1188.        end do
    1189. 

  1189. 

  1190. !      Calculation for the roofs surfaces
    1191. 

  1191.          
    1192. 

  1192.        do iz=2,nz_u
    1193. 

  1193.             
    1194. 

  1194.         if(ss(iz).gt.0.)then
    1195. 

  1195.          do ir=1,nwr_u        
    1196. 

  1196.           tr_tmp(ir)=tr(id,iz,ir)
    1197. 

  1197.          end do
    1198. 

  1198.         
    1199. 

  1199.          call soil_temp(nwr_u,dzr_u,tr_tmp,ptr(id,iz),          &
    1200. 

  1200.                        alar,csr,rs,rl,pr(iz),dt,emr,albr,    &
    1201. 

  1201.                        rtr(id,iz),sfr(id,iz),gfr(id,iz))     
    1202. 

  1202.          do ir=1,nwr_u        
    1203. 

  1203.           tr(id,iz,ir)=tr_tmp(ir)
    1204. 

  1204.          end do
    1205. 

  1205.        
    1206. 

  1206.         end if
    1207. 

  1207.             
    1208. 

  1208.        end do !iz
    1209. 

  1209. 

  1210. !      Calculation for the walls surfaces
    1211. 

  1211.          
    1212. 

  1212.        do iz=1,nz_u
    1213. 

  1213.             
    1214. 

  1214.         do iw=1,nwr_u        
    1215. 

  1215.          tw_tmp(iw)=tw(2*id-1,iz,iw)
    1216. 

  1216.         end do
    1217. 

  1217.         call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id-1,iz),alaw,          &
    1218. 

  1218.                       csw,                                           &     
    1219. 

  1219.                       rsw(2*id-1,iz),rlw(2*id-1,iz),                 &     
    1220. 

  1220.                       pr(iz),dt,emw,                                 &    
    1221. 

  1221.                 albw,rtw(2*id-1,iz),sfw(2*id-1,iz),gfw(2*id-1,iz))   
    1222. 

  1222. 

  1223.         do iw=1,nwr_u        
    1224. 

  1224.          tw(2*id-1,iz,iw)=tw_tmp(iw)
    1225. 

  1225.         end do
    1226. 

  1226.             
    1227. 

  1227.         do iw=1,nwr_u        
    1228. 

  1228.          tw_tmp(iw)=tw(2*id,iz,iw)
    1229. 

  1229.         end do
    1230. 

  1230.             
    1231. 

  1231.         call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id,iz),alaw,          &      
    1232. 

  1232.                       csw,                                         &     
    1233. 

  1233.                       rsw(2*id,iz),rlw(2*id,iz),                   &     
    1234. 

  1234.                       pr(iz),dt,emw,                               &     
    1235. 

  1235.                albw,rtw(2*id,iz),sfw(2*id,iz),gfw(2*id,iz))        
    1236. 

  1236.          do iw=1,nwr_u        
    1237. 

  1237.           tw(2*id,iz,iw)=tw_tmp(iw)
    1238. 

  1238.          end do
    1239. 

  1239.                    
    1240. 

  1240.         end do !iz
    1241. 

  1241. 	
    1242. 

  1242.       end do !id
    1243. 

  1243.       
    1244. 

  1244.       return
    1245. 

  1245.       end subroutine surf_temp
    1246. 

  1246.      
    1247. 

  1247. ! ===6=8===============================================================72     
    1248. 

  1248. ! ===6=8===============================================================72  
    1249. 

  1249. 

  1250.       subroutine buildings(nd,nz,z0,ua_u,va_u,pt_u,pt0_u,         &
    1251. 

  1251.                         ptg,ptr,da_u,ptw,                            &
    1252. 

  1252.                         drst,uva_u,vva_u,uvb_u,vvb_u,                &
    1253. 

  1253.                         tva_u,tvb_u,evb_u,                           &
    1254. 

  1254.                         uhb_u,vhb_u,thb_u,ehb_u,ss,dt)                  
    1255. 

  1255. 

  1256. ! ----------------------------------------------------------------------
    1257. 

  1257. ! This routine computes the sources or sinks of the different quantities 
    1258. 

  1258. ! on the urban grid. The actual calculation is done in the subroutines 
    1259. 

  1259. ! called flux_wall, and flux_flat.
    1260. 

  1260. ! ----------------------------------------------------------------------
    1261. 

  1261. 

  1262.       implicit none
    1263. 

  1263. 

  1264.         
    1265. 

  1265. ! ----------------------------------------------------------------------
    1266. 

  1266. ! INPUT:
    1267. 

  1267. ! ----------------------------------------------------------------------
    1268. 

  1268.       integer nd                ! Number of street direction for the current urban class
    1269. 

  1269.       integer nz                ! number of vertical space steps
    1270. 

  1270.       real ua_u(nz_um)          ! Wind speed in the x direction on the urban grid
    1271. 

  1271.       real va_u(nz_um)          ! Wind speed in the y direction on the urban grid
    1272. 

  1272.       real da_u(nz_um)          ! air density on the urban grid
    1273. 

  1273.       real drst(ndm)            ! Street directions for the current urban class
    1274. 

  1274.       real dz
    1275. 

  1275.       real pt_u(nz_um)          ! Potential temperature on the urban grid
    1276. 

  1276.       real pt0_u(nz_um)         ! reference potential temperature on the urban grid
    1277. 

  1277.       real ptg(ndm)             ! Ground potential temperatures 
    1278. 

  1278.       real ptr(ndm,nz_um)     ! Roof potential temperatures 
    1279. 

  1279.       real ptw(2*ndm,nz_um)     ! Walls potential temperatures 
    1280. 

  1280.       real ss(nz_um)          ! probability to have a building with height h
    1281. 

  1281.       real z0(ndm,nz_um)      ! Roughness lengths "profiles"
    1282. 

  1282.       real dt ! time step
    1283. 

  1283. 

  1284. 

  1285. ! ----------------------------------------------------------------------
    1286. 

  1286. ! OUTPUT:
    1287. 

  1287. ! ----------------------------------------------------------------------
    1288. 

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

  1289. ! vertical surfaces (walls) and horizontal surfaces (roofs and street)
    1290. 

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

  1291. !  Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
    1292. 

  1292. 

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

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

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

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

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

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

  1299.       real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
    1300. 

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

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

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

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

  1304.   
    1305. 

  1305. ! ----------------------------------------------------------------------
    1306. 

  1306. ! LOCAL:
    1307. 

  1307. ! ----------------------------------------------------------------------
    1308. 

  1308.       integer id,iz
    1309. 

  1309. 

  1310. ! ----------------------------------------------------------------------
    1311. 

  1311. ! END VARIABLES DEFINITIONS
    1312. 

  1312. ! ----------------------------------------------------------------------
    1313. 

  1313.        dz=dz_u
    1314. 

  1314. 

  1315.       do id=1,nd
    1316. 

  1316. 

  1317. !        Calculation at the ground surfaces
    1318. 

  1318.          call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1),  &
    1319. 

  1319.                        ptg(id),uhb_u(id,1),                            & 
    1320. 

  1320.                        vhb_u(id,1),thb_u(id,1),ehb_u(id,1))            
    1321. 

  1321. 

  1322. !        Calculation at the roof surfaces    
    1323. 

  1323.          do iz=2,nz
    1324. 

  1324.             if(ss(iz).gt.0)then
    1325. 

  1325.                call flux_flat(dz,z0(id,iz),ua_u(iz),                  &              
    1326. 

  1326.                        va_u(iz),pt_u(iz),pt0_u(iz),                   &   
    1327. 

  1327.                        ptr(id,iz),uhb_u(id,iz),                       &   
    1328. 

  1328.                        vhb_u(id,iz),thb_u(id,iz),ehb_u(id,iz))        
    1329. 

  1329.             else
    1330. 

  1330.                uhb_u(id,iz) = 0.0
    1331. 

  1331.                vhb_u(id,iz) = 0.0
    1332. 

  1332.                thb_u(id,iz) = 0.0
    1333. 

  1333.                ehb_u(id,iz) = 0.0
    1334. 

  1334.             endif
    1335. 

  1335.          end do
    1336. 

  1336. 

  1337. !        Calculation at the wall surfaces         
    1338. 

  1338.          do iz=1,nz         
    1339. 

  1339.             call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),     &  
    1340. 

  1340.                         ptw(2*id-1,iz),                             &   
    1341. 

  1341.                         uva_u(2*id-1,iz),vva_u(2*id-1,iz),          &   
    1342. 

  1342.                         uvb_u(2*id-1,iz),vvb_u(2*id-1,iz),          &   
    1343. 

  1343.                         tva_u(2*id-1,iz),tvb_u(2*id-1,iz),          &   
    1344. 

  1344.                         evb_u(2*id-1,iz),drst(id),dt)                  
    1345. 

  1345.                     
    1346. 

  1346.             call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),    &   
    1347. 

  1347.                         ptw(2*id,iz),                              &    
    1348. 

  1348.                         uva_u(2*id,iz),vva_u(2*id,iz),             &    
    1349. 

  1349.                         uvb_u(2*id,iz),vvb_u(2*id,iz),             &    
    1350. 

  1350.                         tva_u(2*id,iz),tvb_u(2*id,iz),             &   
    1351. 

  1351.                         evb_u(2*id,iz),drst(id),dt) 
    1352. 

  1352. !
    1353. 

  1353.        
    1354. 

  1354.          end do
    1355. 

  1355.          
    1356. 

  1356.       end do
    1357. 

  1357.                 
    1358. 

  1358.       return
    1359. 

  1359.       end subroutine buildings
    1360. 

  1360.       
    1361. 

  1361. 

  1362. ! ===6=8===============================================================72
    1363. 

  1363. ! ===6=8===============================================================72
    1364. 

  1364. 

  1365.         subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl,    &
    1366. 

  1366.                              uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &       
    1367. 

  1367.                              uhb_u,vhb_u,thb_u,ehb_u,                   &      
    1368. 

  1368.                              a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)           
    1369. 

  1369. 

  1370. ! ----------------------------------------------------------------------
    1371. 

  1371. !  This routine interpolates the parameters from the "urban grid" to the
    1372. 

  1372. !  "mesoscale grid".
    1373. 

  1373. !  See p300-301 Appendix B.2 of the BLM paper.  
    1374. 

  1374. ! ----------------------------------------------------------------------
    1375. 

  1375. 

  1376.       implicit none
    1377. 

  1377. 

  1378. ! ----------------------------------------------------------------------
    1379. 

  1379. ! INPUT:
    1380. 

  1380. ! ----------------------------------------------------------------------
    1381. 

  1381. ! Data relative to the "mesoscale grid"
    1382. 

  1382.       integer kms,kme,kts,kte               
    1383. 

  1383.       real z(kms:kme)              ! Altitude above the ground of the cell interface
    1384. 

  1384.       real dz(kms:kme)               ! Vertical space steps
    1385. 

  1385. 

  1386. ! Data relative to the "uban grid"
    1387. 

  1387.       integer nz_u              ! Number of layer in the urban grid
    1388. 

  1388.       integer nd                ! Number of street direction for the current urban class
    1389. 

  1389.       real bs(ndm)              ! Building widths of the current urban class
    1390. 

  1390.       real ws(ndm)              ! Street widths of the current urban class
    1391. 

  1391.       real z_u(nz_um)         ! Height of the urban grid levels
    1392. 

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

  1393.       real ss(nz_um)          ! Probability to have a building with height h
    1394. 

  1394.       real uhb_u(ndm,nz_um)   ! U (x-wind component) Horizontal surfaces, B (explicit) term
    1395. 

  1395.       real uva_u(2*ndm,nz_um)   ! U (x-wind component) Vertical surfaces, A (implicit) term
    1396. 

  1396.       real uvb_u(2*ndm,nz_um)   ! U (x-wind component) Vertical surfaces, B (explicit) term
    1397. 

  1397.       real vhb_u(ndm,nz_um)   ! V (y-wind component) Horizontal surfaces, B (explicit) term
    1398. 

  1398.       real vva_u(2*ndm,nz_um)   ! V (y-wind component) Vertical surfaces, A (implicit) term
    1399. 

  1399.       real vvb_u(2*ndm,nz_um)   ! V (y-wind component) Vertical surfaces, B (explicit) term
    1400. 

  1400.       real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
    1401. 

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

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

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

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

  1405. 

  1406. ! ----------------------------------------------------------------------
    1407. 

  1407. ! OUTPUT:
    1408. 

  1408. ! ----------------------------------------------------------------------
    1409. 

  1409. ! Data relative to the "mesoscale grid"
    1410. 

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

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

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

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

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

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

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

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

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

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

  1420.       
    1421. 

  1421. ! ----------------------------------------------------------------------
    1422. 

  1422. ! LOCAL:
    1423. 

  1423. ! ----------------------------------------------------------------------
    1424. 

  1424.       real dzz
    1425. 

  1425.       real fact
    1426. 

  1426.       integer id,iz,iz_u
    1427. 

  1427.       real se,sr,st,su,sv
    1428. 

  1428.       real uet(kms:kme)                ! Contribution to TKE due to walls
    1429. 

  1429.       real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb
    1430. 

  1430. 

  1431. 

  1432. ! ----------------------------------------------------------------------
    1433. 

  1433. ! END VARIABLES DEFINITIONS
    1434. 

  1434. ! ---------------------------------------------------------------------- 
    1435. 

  1435. 

  1436. ! initialisation
    1437. 

  1437. 

  1438.       do iz=kts,kte
    1439. 

  1439.          a_u(iz)=0.
    1440. 

  1440.          a_v(iz)=0.
    1441. 

  1441.          a_t(iz)=0.
    1442. 

  1442.          a_e(iz)=0.
    1443. 

  1443.          b_u(iz)=0.
    1444. 

  1444.          b_v(iz)=0.
    1445. 

  1445.          b_e(iz)=0.
    1446. 

  1446.          b_t(iz)=0.
    1447. 

  1447.          uet(iz)=0.
    1448. 

  1448.       end do
    1449. 

  1449.             
    1450. 

  1450. ! horizontal surfaces
    1451. 

  1451.       do iz=kts,kte
    1452. 

  1452.          sf(iz)=0.
    1453. 

  1453.          vl(iz)=0.
    1454. 

  1454.       enddo
    1455. 

  1455.       sf(kte+1)=0. 
    1456. 

  1456.       
    1457. 

  1457.       do id=1,nd      
    1458. 

  1458.          do iz=kts+1,kte+1
    1459. 

  1459.             sr=0.
    1460. 

  1460.             do iz_u=2,nz_u
    1461. 

  1461.                if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
    1462. 

  1462.                   sr=pb(iz_u)
    1463. 

  1463.                endif
    1464. 

  1464.             enddo
    1465. 

  1465.             sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
    1466. 

  1466.          enddo
    1467. 

  1467.       enddo
    1468. 

  1468. 

  1469. ! volume      
    1470. 

  1470.       do iz=kts,kte
    1471. 

  1471.          do id=1,nd
    1472. 

  1472.             vtot=0.
    1473. 

  1473.             do iz_u=1,nz_u
    1474. 

  1474.                dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
    1475. 

  1475.                vtot=vtot+pb(iz_u+1)*dzz
    1476. 

  1476.             enddo
    1477. 

  1477.             vtot=vtot/(z(iz+1)-z(iz))
    1478. 

  1478.             vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
    1479. 

  1479.          enddo
    1480. 

  1480.       enddo
    1481. 

  1481.       
    1482. 

  1482. ! horizontal surface impact  
    1483. 

  1483. 

  1484.       do id=1,nd
    1485. 

  1485.       
    1486. 

  1486.          fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
    1487. 

  1487.          b_t(kts)=b_t(kts)+thb_u(id,1)*fact
    1488. 

  1488.          b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
    1489. 

  1489.          b_v(kts)=b_v(kts)+vhb_u(id,1)*fact 
    1490. 

  1490.          b_e(kts)=b_e(kts)+ehb_u(id,1)*fact*(z_u(2)-z_u(1))
    1491. 

  1491.          
    1492. 

  1492.          do iz=kts,kte
    1493. 

  1493.             st=0.
    1494. 

  1494.             su=0.
    1495. 

  1495.             sv=0.
    1496. 

  1496.             se=0.
    1497. 

  1497.             do iz_u=2,nz_u
    1498. 

  1498.                if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
    1499. 

  1499.                   st=st+ss(iz_u)*thb_u(id,iz_u)
    1500. 

  1500.                   su=su+ss(iz_u)*uhb_u(id,iz_u)
    1501. 

  1501.                   sv=sv+ss(iz_u)*vhb_u(id,iz_u)          
    1502. 

  1502.                   se=se+ss(iz_u)*ehb_u(id,iz_u)*(z_u(iz_u+1)-z_u(iz_u))
    1503. 

  1503.                endif
    1504. 

  1504.             enddo
    1505. 

  1505.       
    1506. 

  1506.             fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
    1507. 

  1507.             b_t(iz)=b_t(iz)+st*fact
    1508. 

  1508.             b_u(iz)=b_u(iz)+su*fact
    1509. 

  1509.             b_v(iz)=b_v(iz)+sv*fact
    1510. 

  1510.             b_e(iz)=b_e(iz)+se*fact
    1511. 

  1511.          enddo
    1512. 

  1512.       enddo              
    1513. 

  1513. 

  1514. ! vertical surface impact
    1515. 

  1515.            
    1516. 

  1516.       do iz=kts,kte 
    1517. 

  1517.          uet(iz)=0.
    1518. 

  1518.          do id=1,nd              
    1519. 

  1519.             vtb=0.
    1520. 

  1520.             vta=0.
    1521. 

  1521.             vua=0.
    1522. 

  1522.             vub=0.
    1523. 

  1523.             vva=0.
    1524. 

  1524.             vvb=0.
    1525. 

  1525.             veb=0.
    1526. 

  1526. 	    vte=0.
    1527. 

  1527.             do iz_u=1,nz_u
    1528. 

  1528.                dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
    1529. 

  1529.                fact=dzz/(ws(id)+bs(id))
    1530. 

  1530.                vtb=vtb+pb(iz_u+1)*                                  &        
    1531. 

  1531.                     (tvb_u(2*id-1,iz_u)+tvb_u(2*id,iz_u))*fact   
    1532. 

  1532.                vta=vta+pb(iz_u+1)*                                  &        
    1533. 

  1533.                    (tva_u(2*id-1,iz_u)+tva_u(2*id,iz_u))*fact
    1534. 

  1534.                vua=vua+pb(iz_u+1)*                                  &        
    1535. 

  1535.                     (uva_u(2*id-1,iz_u)+uva_u(2*id,iz_u))*fact
    1536. 

  1536.                vva=vva+pb(iz_u+1)*                                  &        
    1537. 

  1537.                     (vva_u(2*id-1,iz_u)+vva_u(2*id,iz_u))*fact
    1538. 

  1538.                vub=vub+pb(iz_u+1)*                                  &        
    1539. 

  1539.                     (uvb_u(2*id-1,iz_u)+uvb_u(2*id,iz_u))*fact
    1540. 

  1540.                vvb=vvb+pb(iz_u+1)*                                  &        
    1541. 

  1541.                     (vvb_u(2*id-1,iz_u)+vvb_u(2*id,iz_u))*fact
    1542. 

  1542.                veb=veb+pb(iz_u+1)*                                  &        
    1543. 

  1543.                     (evb_u(2*id-1,iz_u)+evb_u(2*id,iz_u))*fact
    1544. 

  1544.             enddo
    1545. 

  1545.            
    1546. 

  1546.             fact=1./vl(iz)/dz(iz)/nd
    1547. 

  1547.             b_t(iz)=b_t(iz)+vtb*fact
    1548. 

  1548.             a_t(iz)=a_t(iz)+vta*fact
    1549. 

  1549.             a_u(iz)=a_u(iz)+vua*fact
    1550. 

  1550.             a_v(iz)=a_v(iz)+vva*fact
    1551. 

  1551.             b_u(iz)=b_u(iz)+vub*fact
    1552. 

  1552.             b_v(iz)=b_v(iz)+vvb*fact
    1553. 

  1553.             b_e(iz)=b_e(iz)+veb*fact
    1554. 

  1554.             uet(iz)=uet(iz)+vte*fact
    1555. 

  1555.          enddo              
    1556. 

  1556.       enddo
    1557. 

  1557.       
    1558. 

  1558. 

  1559.       
    1560. 

  1560.       return
    1561. 

  1561.       end subroutine urban_meso
    1562. 

  1562. 

  1563. ! ===6=8===============================================================72 
    1564. 

  1564. 

  1565. ! ===6=8===============================================================72 
    1566. 

  1566. 

  1567.       subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs,              &
    1568. 

  1568.                              dlg,dl_u)
    1569. 

  1569. 

  1570. ! ----------------------------------------------------------------------     
    1571. 

  1571. !    Calculation of the length scales
    1572. 

  1572. !    See p272-274 formula (22) and (24) of the BLM paper    
    1573. 

  1573. ! ----------------------------------------------------------------------     
    1574. 

  1574.      
    1575. 

  1575.       implicit none
    1576. 

  1576. 

  1577. 

  1578. ! ----------------------------------------------------------------------     
    1579. 

  1579. ! INPUT:
    1580. 

  1580. ! ----------------------------------------------------------------------     
    1581. 

  1581.       integer kms,kme,kts,kte                
    1582. 

  1582.       real z(kms:kme)              ! Altitude above the ground of the cell interface
    1583. 

  1583.       integer nd                ! Number of street direction for the current urban class
    1584. 

  1584.       integer nz_u              ! Number of levels in the "urban grid"
    1585. 

  1585.       real z_u(nz_um)         ! Height of the urban grid levels
    1586. 

  1586.       real bs(ndm)              ! Building widths of the current urban class
    1587. 

  1587.       real ss(nz_um)          ! Probability to have a building with height h
    1588. 

  1588.       real ws(ndm)              ! Street widths of the current urban class
    1589. 

  1589. 

  1590. 

  1591. ! ----------------------------------------------------------------------     
    1592. 

  1592. ! OUTPUT:
    1593. 

  1593. ! ----------------------------------------------------------------------     
    1594. 

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

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

  1596. 

  1597. ! ----------------------------------------------------------------------
    1598. 

  1598. ! LOCAL:
    1599. 

  1599. ! ----------------------------------------------------------------------
    1600. 

  1600.       real dlgtmp
    1601. 

  1601.       integer id,iz,iz_u
    1602. 

  1602.       real sftot
    1603. 

  1603.       real ulu,ssl
    1604. 

  1604. 

  1605. ! ----------------------------------------------------------------------
    1606. 

  1606. ! END VARIABLES DEFINITIONS
    1607. 

  1607. ! ----------------------------------------------------------------------
    1608. 

  1608.    
    1609. 

  1609.         do iz=kts,kte
    1610. 

  1610.          ulu=0.
    1611. 

  1611.          ssl=0.
    1612. 

  1612.          do id=1,nd        
    1613. 

  1613.           do iz_u=2,nz_u
    1614. 

  1614.            if(z_u(iz_u).gt.z(iz))then
    1615. 

  1615.             ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
    1616. 

  1616.             ssl=ssl+ss(iz_u)/nd
    1617. 

  1617.            endif
    1618. 

  1618.           enddo
    1619. 

  1619.          enddo
    1620. 

  1620. 

  1621.         if(ulu.ne.0)then
    1622. 

  1622.           dl_u(iz)=ssl/ulu
    1623. 

  1623.          else
    1624. 

  1624.           dl_u(iz)=0.
    1625. 

  1625.          endif
    1626. 

  1626.         enddo
    1627. 

  1627.        
    1628. 

  1628. 

  1629.         do iz=kts,kte
    1630. 

  1630.          dlg(iz)=0.
    1631. 

  1631.           do id=1,nd
    1632. 

  1632.            sftot=ws(id)  
    1633. 

  1633.            dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
    1634. 

  1634.            do iz_u=1,nz_u
    1635. 

  1635.             if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
    1636. 

  1636.              dlgtmp=dlgtmp+ss(iz_u)*bs(id)/                           &                
    1637. 

  1637.                     ((z(iz)+z(iz+1))/2.-z_u(iz_u))
    1638. 

  1638.              sftot=sftot+ss(iz_u)*bs(id)
    1639. 

  1639.             endif
    1640. 

  1640.            enddo
    1641. 

  1641.            dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
    1642. 

  1642.          enddo
    1643. 

  1643.          dlg(iz)=1./dlg(iz)
    1644. 

  1644.         enddo
    1645. 

  1645.         
    1646. 

  1646.        return
    1647. 

  1647.        end subroutine interp_length
    1648. 

  1648. 

  1649. ! ===6=8===============================================================72
    1650. 

  1650. ! ===6=8===============================================================72   
    1651. 

  1651. 

  1652.       subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z,    &
    1653. 

  1653.                            rs,rsw,rsg)
    1654. 

  1654.         
    1655. 

  1655. ! ----------------------------------------------------------------------
    1656. 

  1656. !         Modification of short wave radiation to take into account
    1657. 

  1657. !         the shadow produced by the buildings
    1658. 

  1658. ! ----------------------------------------------------------------------
    1659. 

  1659. 

  1660.       implicit none
    1661. 

  1661.      
    1662. 

  1662. ! ----------------------------------------------------------------------
    1663. 

  1663. ! INPUT:
    1664. 

  1664. ! ----------------------------------------------------------------------
    1665. 

  1665.       integer nd                ! Number of street direction for the current urban class
    1666. 

  1666.       integer nz_u              ! number of vertical layers defined in the urban grid
    1667. 

  1667.       real ah                   ! Hour angle (it should come from the radiation routine)
    1668. 

  1668.       real deltar               ! Declination of the sun
    1669. 

  1669.       real drst(ndm)            ! street directions for the current urban class
    1670. 

  1670.       real rs                   ! solar radiation
    1671. 

  1671.       real ss(nz_um)          ! probability to have a building with height h
    1672. 

  1672.       real pb(nz_um)          ! Probability that a building has an height greater or equal to h
    1673. 

  1673.       real ws(ndm)              ! Street width of the current urban class
    1674. 

  1674.       real z(nz_um)           ! Height of the urban grid levels
    1675. 

  1675.       real zr                   ! zenith angle
    1676. 

  1676. 

  1677. ! ----------------------------------------------------------------------
    1678. 

  1678. ! OUTPUT:
    1679. 

  1679. ! ----------------------------------------------------------------------
    1680. 

  1680.       real rsg(ndm)             ! Short wave radiation at the ground
    1681. 

  1681.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    1682. 

  1682. 

  1683. ! ----------------------------------------------------------------------
    1684. 

  1684. ! LOCAL:
    1685. 

  1685. ! ----------------------------------------------------------------------
    1686. 

  1686.       integer id,iz,jz
    1687. 

  1687.       real aae,aaw,bbb,phix,rd,rtot,wsd
    1688. 

  1688.       
    1689. 

  1689. ! ----------------------------------------------------------------------
    1690. 

  1690. ! END VARIABLES DEFINITIONS
    1691. 

  1691. ! ----------------------------------------------------------------------
    1692. 

  1692. 

  1693.       if(rs.eq.0.or.sin(zr).eq.1)then
    1694. 

  1694.          do id=1,nd
    1695. 

  1695.             rsg(id)=0.
    1696. 

  1696.             do iz=1,nz_u
    1697. 

  1697.                rsw(2*id-1,iz)=0.
    1698. 

  1698.                rsw(2*id,iz)=0.
    1699. 

  1699.             enddo
    1700. 

  1700.          enddo
    1701. 

  1701.       else            
    1702. 

  1702. !test              
    1703. 

  1703.          if(abs(sin(zr)).gt.1.e-10)then
    1704. 

  1704.           if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
    1705. 

  1705.            bbb=pi/2.
    1706. 

  1706.           elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
    1707. 

  1707.            bbb=-pi/2.
    1708. 

  1708.           else
    1709. 

  1709.            bbb=asin(cos(deltar)*sin(ah)/sin(zr))
    1710. 

  1710.           endif
    1711. 

  1711.          else
    1712. 

  1712.           if(cos(deltar)*sin(ah).ge.0)then
    1713. 

  1713.            bbb=pi/2.
    1714. 

  1714.           elseif(cos(deltar)*sin(ah).lt.0)then
    1715. 

  1715.            bbb=-pi/2.
    1716. 

  1716.           endif
    1717. 

  1717.          endif
    1718. 

  1718. 

  1719.          phix=zr
    1720. 

  1720.            
    1721. 

  1721.          do id=1,nd
    1722. 

  1722.         
    1723. 

  1723.             rsg(id)=0.
    1724. 

  1724.            
    1725. 

  1725.             aae=bbb-drst(id)
    1726. 

  1726.             aaw=bbb-drst(id)+pi
    1727. 

  1727.                     
    1728. 

  1728.             do iz=1,nz_u
    1729. 

  1729.                rsw(2*id-1,iz)=0.
    1730. 

  1730.                rsw(2*id,iz)=0.          
    1731. 

  1731.             if(pb(iz+1).gt.0.)then       
    1732. 

  1732.                do jz=1,nz_u                    
    1733. 

  1733.                 if(abs(sin(aae)).gt.1.e-10)then
    1734. 

  1734.                   call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae,   &    
    1735. 

  1735.                       ws(id),rd)                  
    1736. 

  1736.                   rsw(2*id-1,iz)=rsw(2*id-1,iz)+rs*rd*ss(jz+1)/pb(iz+1)
    1737. 

  1737.                   
    1738. 

  1738.                 endif
    1739. 

  1739.               
    1740. 

  1740.                 if(abs(sin(aaw)).gt.1.e-10)then
    1741. 

  1741.                   call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw,   &    
    1742. 

  1742.                       ws(id),rd)
    1743. 

  1743.                   rsw(2*id,iz)=rsw(2*id,iz)+rs*rd*ss(jz+1)/pb(iz+1)   
    1744. 

  1744.                                 
    1745. 

  1745.                 endif
    1746. 

  1746.                enddo             
    1747. 

  1747.             endif  
    1748. 

  1748.             enddo
    1749. 

  1749.         if(abs(sin(aae)).gt.1.e-10)then
    1750. 

  1750.             wsd=abs(ws(id)/sin(aae))
    1751. 

  1751.               
    1752. 

  1752.             do jz=1,nz_u           
    1753. 

  1753.                rd=max(0.,wsd-z(jz+1)*tan(phix))
    1754. 

  1754.                rsg(id)=rsg(id)+rs*rd*ss(jz+1)/wsd          
    1755. 

  1755.             enddo
    1756. 

  1756.             rtot=0.
    1757. 

  1757.            
    1758. 

  1758.             do iz=1,nz_u
    1759. 

  1759.                rtot=rtot+(rsw(2*id,iz)+rsw(2*id-1,iz))*            &
    1760. 

  1760.                          (z(iz+1)-z(iz))
    1761. 

  1761.             enddo
    1762. 

  1762.             rtot=rtot+rsg(id)*ws(id)
    1763. 

  1763.         else
    1764. 

  1764.             rsg(id)=rs
    1765. 

  1765.         endif
    1766. 

  1766.             
    1767. 

  1767.          enddo
    1768. 

  1768.       endif
    1769. 

  1769.          
    1770. 

  1770.       return
    1771. 

  1771.       end subroutine shadow_mas
    1772. 

  1772.          
    1773. 

  1773. ! ===6=8===============================================================72     
    1774. 

  1774. ! ===6=8===============================================================72     
    1775. 

  1775. 

  1776.       subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
    1777. 

  1777. 

  1778. ! ----------------------------------------------------------------------
    1779. 

  1779. ! This routine computes the effects of a shadow induced by a building of 
    1780. 

  1780. ! height hu, on a portion of wall between z1 and z2. See equation A10, 
    1781. 

  1781. ! and correction described below formula A11, and figure A1. Basically rd
    1782. 

  1782. ! is the ratio between the horizontal surface illuminated and the portion
    1783. 

  1783. ! of wall. Referring to figure A1, multiplying radiation flux density on 
    1784. 

  1784. ! a horizontal surface (rs) by x1-x2 we have the radiation energy per 
    1785. 

  1785. ! unit time. Dividing this by z2-z1, we obtain the radiation flux 
    1786. 

  1786. ! density reaching the portion of the wall between z2 and z1 
    1787. 

  1787. ! (everything is assumed in 2D)
    1788. 

  1788. ! ----------------------------------------------------------------------
    1789. 

  1789. 

  1790.       implicit none
    1791. 

  1791.       
    1792. 

  1792. ! ----------------------------------------------------------------------
    1793. 

  1793. ! INPUT:
    1794. 

  1794. ! ----------------------------------------------------------------------
    1795. 

  1795.       real aa                   ! Angle between the sun direction and the face of the wall (A12)
    1796. 

  1796.       real hu                   ! Height of the building that generates the shadow
    1797. 

  1797.       real phix                 ! Solar zenith angle
    1798. 

  1798.       real ws                   ! Width of the street
    1799. 

  1799.       real z1                   ! Height of the level z(iz)
    1800. 

  1800.       real z2                   ! Height of the level z(iz+1)
    1801. 

  1801. 

  1802. ! ----------------------------------------------------------------------
    1803. 

  1803. ! OUTPUT:
    1804. 

  1804. ! ----------------------------------------------------------------------
    1805. 

  1805.       real rd                   ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A. 
    1806. 

  1806.                                 ! Multiplying rd by rs (radiation flux 
    1807. 

  1807.                                 ! density on a horizontal surface) gives 
    1808. 

  1808.                                 ! the radiation flux density on the 
    1809. 

  1809.                                 ! portion of wall between z1 and z2. 
    1810. 

  1810. ! ----------------------------------------------------------------------
    1811. 

  1811. ! LOCAL:
    1812. 

  1812. ! ----------------------------------------------------------------------
    1813. 

  1813.       real x1,x2                ! x1,x2 see Fig. A1.
    1814. 

  1814. 

  1815. ! ----------------------------------------------------------------------
    1816. 

  1816. ! END VARIABLES DEFINITIONS
    1817. 

  1817. ! ----------------------------------------------------------------------
    1818. 

  1818. 

  1819.       x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
    1820. 

  1820.       
    1821. 

  1821.       x2=max((hu-z2)*tan(phix),0.)
    1822. 

  1822. 

  1823.       rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
    1824. 

  1824.       
    1825. 

  1825.       return
    1826. 

  1826.       end subroutine shade_wall
    1827. 

  1827. 

  1828. ! ===6=8===============================================================72     
    1829. 

  1829. ! ===6=8===============================================================72     
    1830. 

  1830. 

  1831.       subroutine long_rad(iurb,nz_u,id,emw,emg,                  &
    1832. 

  1832.                          fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
    1833. 

  1833. 

  1834. ! ----------------------------------------------------------------------
    1835. 

  1835. ! This routine computes the effects of the reflections of long-wave 
    1836. 

  1836. ! radiation in the street canyon by solving the system 
    1837. 

  1837. ! of 2*nz_u+1 eqn. in 2*nz_u+1
    1838. 

  1838. ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
    1839. 

  1839. ! The system is solved by solving A X= B,
    1840. 

  1840. ! with A matrix B vector, and X solution. 
    1841. 

  1841. ! ----------------------------------------------------------------------
    1842. 

  1842. 

  1843.       implicit none
    1844. 

  1844. 

  1845.   
    1846. 

  1846.       
    1847. 

  1847. ! ----------------------------------------------------------------------
    1848. 

  1848. ! INPUT:
    1849. 

  1849. ! ----------------------------------------------------------------------
    1850. 

  1850.       real emg                        ! Emissivity of ground for the current urban class
    1851. 

  1851.       real emw                        ! Emissivity of wall for the current urban class
    1852. 

  1852.       real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
    1853. 

  1853.       real fsg(ndm,nurbm)             ! View factors from sky to ground
    1854. 

  1854.       real fsw(nz_um,ndm,nurbm)       ! View factors from sky to wall
    1855. 

  1855.       real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
    1856. 

  1856.       real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
    1857. 

  1857.       integer id                      ! Current street direction
    1858. 

  1858.       integer iurb                    ! Current urban class
    1859. 

  1859.       integer nz_u                    ! Number of layer in the urban grid
    1860. 

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

  1861.       real rl                         ! Downward flux of the longwave radiation
    1862. 

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

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

  1864.       
    1865. 

  1865. 

  1866. ! ----------------------------------------------------------------------
    1867. 

  1867. ! OUTPUT:
    1868. 

  1868. ! ----------------------------------------------------------------------
    1869. 

  1869.       real rlg(ndm)                   ! Long wave radiation at the ground
    1870. 

  1870.       real rlw(2*ndm,nz_um)           ! Long wave radiation at the walls
    1871. 

  1871. 

  1872. ! ----------------------------------------------------------------------
    1873. 

  1873. ! LOCAL:
    1874. 

  1874. ! ----------------------------------------------------------------------
    1875. 

  1875.       integer i,j
    1876. 

  1876.       real aaa(2*nz_um+1,2*nz_um+1)   ! terms of the matrix
    1877. 

  1877.       real bbb(2*nz_um+1)             ! terms of the vector
    1878. 

  1878. 

  1879. ! ----------------------------------------------------------------------
    1880. 

  1880. ! END VARIABLES DEFINITIONS
    1881. 

  1881. ! ----------------------------------------------------------------------
    1882. 

  1882. 

  1883. 

  1884. ! west wall
    1885. 

  1885.        
    1886. 

  1886.       do i=1,nz_u
    1887. 

  1887.         
    1888. 

  1888.         do j=1,nz_u
    1889. 

  1889.          aaa(i,j)=0.
    1890. 

  1890.         enddo
    1891. 

  1891.         
    1892. 

  1892.         aaa(i,i)=1.        
    1893. 

  1893.        
    1894. 

  1894.         do j=nz_u+1,2*nz_u
    1895. 

  1895.          aaa(i,j)=-(1.-emw)*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
    1896. 

  1896.         enddo
    1897. 

  1897.         
    1898. 

  1898. !!      aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)*pb(i+1)
    1899. 

  1899.         aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
    1900. 

  1900.         
    1901. 

  1901.         bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg(id,ng_u)**4
    1902. 

  1902.         do j=1,nz_u
    1903. 

  1903.          bbb(i)=bbb(i)+pb(j+1)*emw*sigma*fww(j,i,id,iurb)*       &
    1904. 

  1904.                tw(2*id,j,nwr_u)**4+                                 &
    1905. 

  1905.                fww(j,i,id,iurb)*rl*(1.-pb(j+1))
    1906. 

  1906.         enddo
    1907. 

  1907.         
    1908. 

  1908.        enddo
    1909. 

  1909.        
    1910. 

  1910. ! east wall
    1911. 

  1911. 

  1912.        do i=1+nz_u,2*nz_u
    1913. 

  1913.         
    1914. 

  1914.         do j=1,nz_u
    1915. 

  1915.          aaa(i,j)=-(1.-emw)*fww(j,i-nz_u,id,iurb)*pb(j+1)
    1916. 

  1916.         enddo
    1917. 

  1917.         
    1918. 

  1918.         do j=1+nz_u,2*nz_u
    1919. 

  1919.          aaa(i,j)=0.
    1920. 

  1920.         enddo
    1921. 

  1921.         
    1922. 

  1922.         aaa(i,i)=1.
    1923. 

  1923.         
    1924. 

  1924. !!      aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)*pb(i-nz_u+1)
    1925. 

  1925.         aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
    1926. 

  1926.         
    1927. 

  1927.         bbb(i)=fsw(i-nz_u,id,iurb)*rl+                           &     
    1928. 

  1928.               emg*fgw(i-nz_u,id,iurb)*sigma*tg(id,ng_u)**4
    1929. 

  1929. 

  1930.         do j=1,nz_u
    1931. 

  1931.          bbb(i)=bbb(i)+pb(j+1)*emw*sigma*fww(j,i-nz_u,id,iurb)*  &   
    1932. 

  1932.                 tw(2*id-1,j,nwr_u)**4+                              &   
    1933. 

  1933.                 fww(j,i-nz_u,id,iurb)*rl*(1.-pb(j+1))
    1934. 

  1934.         enddo
    1935. 

  1935.        
    1936. 

  1936.        enddo
    1937. 

  1937. 

  1938. ! ground
    1939. 

  1939.        do j=1,nz_u
    1940. 

  1940.         aaa(2*nz_u+1,j)=-(1.-emw)*fwg(j,id,iurb)*pb(j+1)
    1941. 

  1941.        enddo
    1942. 

  1942.        
    1943. 

  1943.        do j=nz_u+1,2*nz_u
    1944. 

  1944.         aaa(2*nz_u+1,j)=-(1.-emw)*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
    1945. 

  1945.        enddo
    1946. 

  1946.        
    1947. 

  1947.        aaa(2*nz_u+1,2*nz_u+1)=1.
    1948. 

  1948.        
    1949. 

  1949.        bbb(2*nz_u+1)=fsg(id,iurb)*rl
    1950. 

  1950.        
    1951. 

  1951.        do i=1,nz_u
    1952. 

  1952.         bbb(2*nz_u+1)=bbb(2*nz_u+1)+emw*sigma*fwg(i,id,iurb)*pb(i+1)*    &
    1953. 

  1953.                       (tw(2*id-1,i,nwr_u)**4+tw(2*id,i,nwr_u)**4)+             &
    1954. 

  1954.                       2.*fwg(i,id,iurb)*(1.-pb(i+1))*rl                  
    1955. 

  1955.        enddo
    1956. 

  1956.    
    1957. 

  1957. 

  1958.      
    1959. 

  1959.        call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
    1960. 

  1960. 

  1961.        do i=1,nz_u
    1962. 

  1962.         rlw(2*id-1,i)=bbb(i)
    1963. 

  1963.        enddo
    1964. 

  1964.        
    1965. 

  1965.        do i=nz_u+1,2*nz_u
    1966. 

  1966.         rlw(2*id,i-nz_u)=bbb(i)
    1967. 

  1967.        enddo
    1968. 

  1968.        
    1969. 

  1969.        rlg(id)=bbb(2*nz_u+1)
    1970. 

  1970.   
    1971. 

  1971.        return
    1972. 

  1972.        end subroutine long_rad
    1973. 

  1973.              
    1974. 

  1974. ! ===6=8===============================================================72
    1975. 

  1975. ! ===6=8===============================================================72
    1976. 

  1976. 

  1977.        subroutine short_rad(iurb,nz_u,id,albw,                        & 
    1978. 

  1978.                            albg,fwg,fww,fgw,rsg,rsw,pb)
    1979. 

  1979. 

  1980. ! ----------------------------------------------------------------------
    1981. 

  1981. ! This routine computes the effects of the reflections of short-wave 
    1982. 

  1982. ! (solar) radiation in the street canyon by solving the system 
    1983. 

  1983. ! of 2*nz_u+1 eqn. in 2*nz_u+1
    1984. 

  1984. ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
    1985. 

  1985. ! The system is solved by solving A X= B,
    1986. 

  1986. ! with A matrix B vector, and X solution. 
    1987. 

  1987. ! ----------------------------------------------------------------------
    1988. 

  1988. 

  1989.       implicit none
    1990. 

  1990. 

  1991.   
    1992. 

  1992.       
    1993. 

  1993. ! ----------------------------------------------------------------------
    1994. 

  1994. ! INPUT:
    1995. 

  1995. ! ----------------------------------------------------------------------
    1996. 

  1996.       real albg                 ! Albedo of the ground for the current urban class
    1997. 

  1997.       real albw                 ! Albedo of the wall for the current urban class
    1998. 

  1998.       real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
    1999. 

  1999.       real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
    2000. 

  2000.       real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
    2001. 

  2001.       integer id                ! current street direction 
    2002. 

  2002.       integer iurb              ! current urban class
    2003. 

  2003.       integer nz_u              ! Number of layer in the urban grid
    2004. 

  2004.       real pb(nz_um)          ! probability to have a building with an height equal
    2005. 

  2005. 

  2006. ! ----------------------------------------------------------------------
    2007. 

  2007. ! OUTPUT:
    2008. 

  2008. ! ----------------------------------------------------------------------
    2009. 

  2009.       real rsg(ndm)             ! Short wave radiation at the ground
    2010. 

  2010.       real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls
    2011. 

  2011. 

  2012. ! ----------------------------------------------------------------------
    2013. 

  2013. ! LOCAL:
    2014. 

  2014. ! ----------------------------------------------------------------------
    2015. 

  2015.       integer i,j
    2016. 

  2016.       real aaa(2*nz_um+1,2*nz_um+1)  ! terms of the matrix
    2017. 

  2017.       real bbb(2*nz_um+1)            ! terms of the vector
    2018. 

  2018. 

  2019. ! ----------------------------------------------------------------------
    2020. 

  2020. ! END VARIABLES DEFINITIONS
    2021. 

  2021. ! ----------------------------------------------------------------------
    2022. 

  2022. 

  2023.       
    2024. 

  2024. ! west wall
    2025. 

  2025.        
    2026. 

  2026.       do i=1,nz_u
    2027. 

  2027.          do j=1,nz_u
    2028. 

  2028.             aaa(i,j)=0.
    2029. 

  2029.          enddo
    2030. 

  2030.          
    2031. 

  2031.          aaa(i,i)=1.        
    2032. 

  2032.          
    2033. 

  2033.          do j=nz_u+1,2*nz_u
    2034. 

  2034.             aaa(i,j)=-albw*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
    2035. 

  2035.          enddo
    2036. 

  2036.          
    2037. 

  2037.          aaa(i,2*nz_u+1)=-albg*fgw(i,id,iurb)
    2038. 

  2038.          bbb(i)=rsw(2*id-1,i)
    2039. 

  2039.          
    2040. 

  2040.       enddo
    2041. 

  2041.        
    2042. 

  2042. ! east wall
    2043. 

  2043. 

  2044.       do i=1+nz_u,2*nz_u
    2045. 

  2045.          do j=1,nz_u
    2046. 

  2046.             aaa(i,j)=-albw*fww(j,i-nz_u,id,iurb)*pb(j+1)
    2047. 

  2047.          enddo
    2048. 

  2048.          
    2049. 

  2049.          do j=1+nz_u,2*nz_u
    2050. 

  2050.             aaa(i,j)=0.
    2051. 

  2051.          enddo
    2052. 

  2052.          
    2053. 

  2053.         aaa(i,i)=1.
    2054. 

  2054.         aaa(i,2*nz_u+1)=-albg*fgw(i-nz_u,id,iurb)
    2055. 

  2055.         bbb(i)=rsw(2*id,i-nz_u)
    2056. 

  2056.         
    2057. 

  2057.       enddo
    2058. 

  2058. 

  2059. ! ground
    2060. 

  2060. 

  2061.       do j=1,nz_u
    2062. 

  2062.          aaa(2*nz_u+1,j)=-albw*fwg(j,id,iurb)*pb(j+1)
    2063. 

  2063.       enddo
    2064. 

  2064.        
    2065. 

  2065.       do j=nz_u+1,2*nz_u
    2066. 

  2066.          aaa(2*nz_u+1,j)=-albw*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
    2067. 

  2067.       enddo
    2068. 

  2068.        
    2069. 

  2069.       aaa(2*nz_u+1,2*nz_u+1)=1.
    2070. 

  2070.       bbb(2*nz_u+1)=rsg(id)
    2071. 

  2071.        
    2072. 

  2072.       call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
    2073. 

  2073. 

  2074.       do i=1,nz_u
    2075. 

  2075.          rsw(2*id-1,i)=bbb(i)
    2076. 

  2076.       enddo
    2077. 

  2077.        
    2078. 

  2078.       do i=nz_u+1,2*nz_u
    2079. 

  2079.          rsw(2*id,i-nz_u)=bbb(i) 
    2080. 

  2080.       enddo
    2081. 

  2081.        
    2082. 

  2082.       rsg(id)=bbb(2*nz_u+1)
    2083. 

  2083.        
    2084. 

  2084.       return
    2085. 

  2085.       end subroutine short_rad
    2086. 

  2086.              
    2087. 

  2087. 

  2088. ! ===6=8===============================================================72     
    2089. 

  2089. ! ===6=8===============================================================72     
    2090. 

  2090.       
    2091. 

  2091.       subroutine gaussj(a,n,b,np)
    2092. 

  2092. 

  2093. ! ----------------------------------------------------------------------
    2094. 

  2094. ! This routine solve a linear system of n equations of the form
    2095. 

  2095. !              A X = B
    2096. 

  2096. !  where  A is a matrix a(i,j)
    2097. 

  2097. !         B a vector and X the solution
    2098. 

  2098. ! In output b is replaced by the solution     
    2099. 

  2099. ! ----------------------------------------------------------------------
    2100. 

  2100. 

  2101.       implicit none
    2102. 

  2102. 

  2103. ! ----------------------------------------------------------------------
    2104. 

  2104. ! INPUT:
    2105. 

  2105. ! ----------------------------------------------------------------------
    2106. 

  2106.       integer np
    2107. 

  2107.       real a(np,np)
    2108. 

  2108. 

  2109. ! ----------------------------------------------------------------------
    2110. 

  2110. ! OUTPUT:
    2111. 

  2111. ! ----------------------------------------------------------------------
    2112. 

  2112.       real b(np)
    2113. 

  2113. 

  2114. ! ----------------------------------------------------------------------
    2115. 

  2115. ! LOCAL:
    2116. 

  2116. ! ----------------------------------------------------------------------
    2117. 

  2117.       integer nmax
    2118. 

  2118.       parameter (nmax=150)
    2119. 

  2119. 

  2120.       real big,dum
    2121. 

  2121.       integer i,icol,irow
    2122. 

  2122.       integer j,k,l,ll,n
    2123. 

  2123.       integer ipiv(nmax)
    2124. 

  2124.       real pivinv
    2125. 

  2125. 

  2126. ! ----------------------------------------------------------------------
    2127. 

  2127. ! END VARIABLES DEFINITIONS
    2128. 

  2128. ! ----------------------------------------------------------------------
    2129. 

  2129.        
    2130. 

  2130.       do j=1,n
    2131. 

  2131.          ipiv(j)=0.
    2132. 

  2132.       enddo
    2133. 

  2133.        
    2134. 

  2134.       do i=1,n
    2135. 

  2135.          big=0.
    2136. 

  2136.          do j=1,n
    2137. 

  2137.             if(ipiv(j).ne.1)then
    2138. 

  2138.                do k=1,n
    2139. 

  2139.                   if(ipiv(k).eq.0)then
    2140. 

  2140.                      if(abs(a(j,k)).ge.big)then
    2141. 

  2141.                         big=abs(a(j,k))
    2142. 

  2142.                         irow=j
    2143. 

  2143.                         icol=k
    2144. 

  2144.                      endif
    2145. 

  2145.                   elseif(ipiv(k).gt.1)then
    2146. 

  2146.                      pause 'singular matrix in gaussj'
    2147. 

  2147.                   endif
    2148. 

  2148.                enddo
    2149. 

  2149.             endif
    2150. 

  2150.          enddo
    2151. 

  2151.          
    2152. 

  2152.          ipiv(icol)=ipiv(icol)+1
    2153. 

  2153.          
    2154. 

  2154.          if(irow.ne.icol)then
    2155. 

  2155.             do l=1,n
    2156. 

  2156.                dum=a(irow,l)
    2157. 

  2157.                a(irow,l)=a(icol,l)
    2158. 

  2158.                a(icol,l)=dum
    2159. 

  2159.             enddo
    2160. 

  2160.             
    2161. 

  2161.             dum=b(irow)
    2162. 

  2162.             b(irow)=b(icol)
    2163. 

  2163.             b(icol)=dum
    2164. 

  2164.           
    2165. 

  2165.          endif
    2166. 

  2166.          
    2167. 

  2167.          if(a(icol,icol).eq.0)pause 'singular matrix in gaussj'
    2168. 

  2168.          
    2169. 

  2169.          pivinv=1./a(icol,icol)
    2170. 

  2170.          a(icol,icol)=1
    2171. 

  2171.          
    2172. 

  2172.          do l=1,n
    2173. 

  2173.             a(icol,l)=a(icol,l)*pivinv
    2174. 

  2174.          enddo
    2175. 

  2175.          
    2176. 

  2176.          b(icol)=b(icol)*pivinv
    2177. 

  2177.          
    2178. 

  2178.          do ll=1,n
    2179. 

  2179.             if(ll.ne.icol)then
    2180. 

  2180.                dum=a(ll,icol)
    2181. 

  2181.                a(ll,icol)=0.
    2182. 

  2182.                do l=1,n
    2183. 

  2183.                   a(ll,l)=a(ll,l)-a(icol,l)*dum
    2184. 

  2184.                enddo
    2185. 

  2185.                
    2186. 

  2186.                b(ll)=b(ll)-b(icol)*dum
    2187. 

  2187.                
    2188. 

  2188.             endif
    2189. 

  2189.          enddo
    2190. 

  2190.       enddo
    2191. 

  2191.       
    2192. 

  2192.       return
    2193. 

  2193.       end subroutine gaussj
    2194. 

  2194.          
    2195. 

  2195. ! ===6=8===============================================================72     
    2196. 

  2196. ! ===6=8===============================================================72     
    2197. 

  2197.        
    2198. 

  2198.       subroutine soil_temp(nz,dz,temp,pt,ala,cs,                       &
    2199. 

  2199.                           rs,rl,press,dt,em,alb,rt,sf,gf)
    2200. 

  2200. 

  2201. ! ----------------------------------------------------------------------
    2202. 

  2202. ! This routine solves the Fourier diffusion equation for heat in 
    2203. 

  2203. ! the material (wall, roof, or ground). Resolution is done implicitely.
    2204. 

  2204. ! Boundary conditions are: 
    2205. 

  2205. ! - fixed temperature at the interior
    2206. 

  2206. ! - energy budget at the surface
    2207. 

  2207. ! ----------------------------------------------------------------------
    2208. 

  2208. 

  2209.       implicit none
    2210. 

  2210. 

  2211.      
    2212. 

  2212.                 
    2213. 

  2213. ! ----------------------------------------------------------------------
    2214. 

  2214. ! INPUT:
    2215. 

  2215. ! ----------------------------------------------------------------------
    2216. 

  2216.       integer nz                ! Number of layers
    2217. 

  2217.       real ala(nz)              ! Thermal diffusivity in each layers [m^2 s^-1] 
    2218. 

  2218.       real alb                  ! Albedo of the surface
    2219. 

  2219.       real cs(nz)               ! Specific heat of the material [J m^3 K^-1]
    2220. 

  2220.       real dt                   ! Time step
    2221. 

  2221.       real em                   ! Emissivity of the surface
    2222. 

  2222.       real press                ! Pressure at ground level
    2223. 

  2223.       real rl                   ! Downward flux of the longwave radiation
    2224. 

  2224.       real rs                   ! Solar radiation
    2225. 

  2225.       real sf                   ! Sensible heat flux at the surface
    2226. 

  2226.       real temp(nz)             ! Temperature in each layer [K]
    2227. 

  2227.       real dz(nz)               ! Layer sizes [m]
    2228. 

  2228. 

  2229. 

  2230. ! ----------------------------------------------------------------------
    2231. 

  2231. ! OUTPUT:
    2232. 

  2232. ! ----------------------------------------------------------------------
    2233. 

  2233.       real gf                   ! Heat flux transferred from the surface toward the interior
    2234. 

  2234.       real pt                   ! Potential temperature at the surface
    2235. 

  2235.       real rt                   ! Total radiation at the surface (solar+incoming long+outgoing long)
    2236. 

  2236. 

  2237. ! ----------------------------------------------------------------------
    2238. 

  2238. ! LOCAL:
    2239. 

  2239. ! ----------------------------------------------------------------------
    2240. 

  2240.       integer iz
    2241. 

  2241.       real a(nz,3)
    2242. 

  2242.       real alpha
    2243. 

  2243.       real c(nz)
    2244. 

  2244.       real cddz(nz+2)
    2245. 

  2245.       real tsig
    2246. 

  2246. 

  2247. ! ----------------------------------------------------------------------
    2248. 

  2248. ! END VARIABLES DEFINITIONS
    2249. 

  2249. ! ----------------------------------------------------------------------
    2250. 

  2250.        
    2251. 

  2251.       tsig=temp(nz)
    2252. 

  2252.       alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
    2253. 

  2253. ! Compute cddz=2*cd/dz  
    2254. 

  2254.         
    2255. 

  2255.       cddz(1)=ala(1)/dz(1)
    2256. 

  2256.       do iz=2,nz
    2257. 

  2257.          cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
    2258. 

  2258.       enddo
    2259. 

  2259. !        cddz(nz+1)=ala(nz+1)/dz(nz)
    2260. 

  2260.        
    2261. 

  2261.       a(1,1)=0.
    2262. 

  2262.       a(1,2)=1.
    2263. 

  2263.       a(1,3)=0.
    2264. 

  2264.       c(1)=temp(1)
    2265. 

  2265.           
    2266. 

  2266.       do iz=2,nz-1
    2267. 

  2267.          a(iz,1)=-cddz(iz)*dt/dz(iz)
    2268. 

  2268.          a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)          
    2269. 

  2269.          a(iz,3)=-cddz(iz+1)*dt/dz(iz)
    2270. 

  2270.          c(iz)=temp(iz)
    2271. 

  2271.       enddo          
    2272. 

  2272.                      
    2273. 

  2273.       a(nz,1)=-dt*cddz(nz)/dz(nz)
    2274. 

  2274.       a(nz,2)=1.+dt*cddz(nz)/dz(nz)
    2275. 

  2275.       a(nz,3)=0.
    2276. 

  2276.       c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
    2277. 

  2277. 

  2278.       
    2279. 

  2279.       call invert(nz,a,c,temp)
    2280. 

  2280. 

  2281.            
    2282. 

  2282.       pt=temp(nz)*(press/1.e+5)**(-rcp_u)
    2283. 

  2283. 

  2284.       rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)
    2285. 

  2285.                         
    2286. 

  2286. !      gf=-cddz(nz)*(temp(nz)-temp(nz-1))*cs(nz)
    2287. 

  2287.        gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf                                   
    2288. 

  2288.       return
    2289. 

  2289.       end subroutine soil_temp
    2290. 

  2290. 

  2291. ! ===6=8===============================================================72 
    2292. 

  2292. ! ===6=8===============================================================72 
    2293. 

  2293. 

  2294.       subroutine invert(n,a,c,x)
    2295. 

  2295. 

  2296. ! ----------------------------------------------------------------------
    2297. 

  2297. !        Inversion and resolution of a tridiagonal matrix
    2298. 

  2298. !                   A X = C
    2299. 

  2299. ! ----------------------------------------------------------------------
    2300. 

  2300. 

  2301.       implicit none
    2302. 

  2302.                 
    2303. 

  2303. ! ----------------------------------------------------------------------
    2304. 

  2304. ! INPUT:
    2305. 

  2305. ! ----------------------------------------------------------------------
    2306. 

  2306.        integer n
    2307. 

  2307.        real a(n,3)              !  a(*,1) lower diagonal (Ai,i-1)
    2308. 

  2308.                                 !  a(*,2) principal diagonal (Ai,i)
    2309. 

  2309.                                 !  a(*,3) upper diagonal (Ai,i+1)
    2310. 

  2310.        real c(n)
    2311. 

  2311. 

  2312. ! ----------------------------------------------------------------------
    2313. 

  2313. ! OUTPUT:
    2314. 

  2314. ! ----------------------------------------------------------------------
    2315. 

  2315.        real x(n)    
    2316. 

  2316. 

  2317. ! ----------------------------------------------------------------------
    2318. 

  2318. ! LOCAL:
    2319. 

  2319. ! ----------------------------------------------------------------------
    2320. 

  2320.        integer i
    2321. 

  2321. 

  2322. ! ----------------------------------------------------------------------
    2323. 

  2323. ! END VARIABLES DEFINITIONS
    2324. 

  2324. ! ----------------------------------------------------------------------
    2325. 

  2325.                      
    2326. 

  2326.        do i=n-1,1,-1                 
    2327. 

  2327.           c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
    2328. 

  2328.           a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
    2329. 

  2329.        enddo
    2330. 

  2330.        
    2331. 

  2331.        do i=2,n        
    2332. 

  2332.           c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
    2333. 

  2333.        enddo
    2334. 

  2334.         
    2335. 

  2335.        do i=1,n
    2336. 

  2336.           x(i)=c(i)/a(i,2)
    2337. 

  2337.        enddo
    2338. 

  2338. 

  2339.        return
    2340. 

  2340.        end subroutine invert
    2341. 

  2341.   
    2342. 

  2342. 

  2343. ! ===6=8===============================================================72  
    2344. 

  2344. ! ===6=8===============================================================72
    2345. 

  2345.   
    2346. 

  2346.       subroutine flux_wall(ua,va,pt,da,ptw,uva,vva,uvb,vvb,            &
    2347. 

  2347.                                   tva,tvb,evb,drst,dt)      
    2348. 

  2348.        
    2349. 

  2349. ! ----------------------------------------------------------------------
    2350. 

  2350. ! This routine computes the surface sources or sinks of momentum, tke,
    2351. 

  2351. ! and heat from vertical surfaces (walls).   
    2352. 

  2352. ! ----------------------------------------------------------------------
    2353. 

  2353. 

  2354.       implicit none
    2355. 

  2355. 

  2356.     
    2357. 

  2357.          
    2358. 

  2358. ! INPUT:
    2359. 

  2359. ! -----
    2360. 

  2360.       real drst                 ! street directions for the current urban class
    2361. 

  2361.       real da                   ! air density
    2362. 

  2362.       real pt                   ! potential temperature
    2363. 

  2363.       real ptw                  ! Walls potential temperatures 
    2364. 

  2364.       real ua                   ! wind speed
    2365. 

  2365.       real va                   ! wind speed
    2366. 

  2366. 

  2367.       real dt                   !time step
    2368. 

  2368. ! OUTPUT:
    2369. 

  2369. ! ------
    2370. 

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

  2371. ! vertical surfaces (walls).
    2372. 

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

  2373. !  Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
    2374. 

  2374.       real uva                  ! U (wind component)   Vertical surfaces, A (implicit) term
    2375. 

  2375.       real uvb                  ! U (wind component)   Vertical surfaces, B (explicit) term
    2376. 

  2376.       real vva                  ! V (wind component)   Vertical surfaces, A (implicit) term
    2377. 

  2377.       real vvb                  ! V (wind component)   Vertical surfaces, B (explicit) term
    2378. 

  2378.       real tva                  ! Temperature          Vertical surfaces, A (implicit) term
    2379. 

  2379.       real tvb                  ! Temperature          Vertical surfaces, B (explicit) term
    2380. 

  2380.       real evb                  ! Energy (TKE)         Vertical surfaces, B (explicit) term
    2381. 

  2381. 

  2382. ! LOCAL:
    2383. 

  2383. ! -----
    2384. 

  2384.       real hc
    2385. 

  2385.       real u_ort
    2386. 

  2386.       real vett
    2387. 

  2387. 

  2388. ! -------------------------
    2389. 

  2389. ! END VARIABLES DEFINITIONS
    2390. 

  2390. ! -------------------------
    2391. 

  2391. 

  2392.       vett=(ua**2+va**2)**.5         
    2393. 

  2393.          
    2394. 

  2394.       u_ort=abs((cos(drst)*ua-sin(drst)*va))
    2395. 

  2395.        
    2396. 

  2396.       uva=-cdrag*u_ort/2.*cos(drst)*cos(drst)
    2397. 

  2397.       vva=-cdrag*u_ort/2.*sin(drst)*sin(drst)
    2398. 

  2398.          
    2399. 

  2399.       uvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*va
    2400. 

  2400.       vvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*ua
    2401. 

  2401.          
    2402. 

  2402.       hc=5.678*(1.09+0.23*(vett/0.3048))  
    2403. 

  2403. 

  2404.       if(hc.gt.da*cp_u/dt)then
    2405. 

  2405.         hc=da*cp_u/dt
    2406. 

  2406.       endif
    2407. 

  2407.          
    2408. 

  2408. !         tvb=hc*ptw/da/cp_u
    2409. 

  2409. !         tva=-hc/da/cp_u
    2410. 

  2410. !!!!!!!!!!!!!!!!!!!!
    2411. 

  2411. ! explicit 
    2412. 

  2412.       tvb=hc*ptw/da/cp_u-hc/da/cp_u*pt !c
    2413. 

  2413.       tva = 0.                  !c
    2414. 

  2414.          
    2415. 

  2415.       evb=cdrag*(abs(u_ort)**3.)/2.
    2416. 

  2416.               
    2417. 

  2417.       return
    2418. 

  2418.       end subroutine flux_wall
    2419. 

  2419.          
    2420. 

  2420. ! ===6=8===============================================================72
    2421. 

  2421. 

  2422. ! ===6=8===============================================================72
    2423. 

  2423. 

  2424.       subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg,                     &
    2425. 

  2425.                           uhb,vhb,thb,ehb)
    2426. 

  2426.                                 
    2427. 

  2427. ! ----------------------------------------------------------------------
    2428. 

  2428. !           Calculation of the flux at the ground 
    2429. 

  2429. !           Formulation of Louis (Louis, 1979)       
    2430. 

  2430. ! ----------------------------------------------------------------------
    2431. 

  2431. 

  2432.       implicit none
    2433. 

  2433. 

  2434.      
    2435. 

  2435. 

  2436. ! ----------------------------------------------------------------------
    2437. 

  2437. ! INPUT:
    2438. 

  2438. ! ----------------------------------------------------------------------
    2439. 

  2439.       real dz                   ! first vertical level
    2440. 

  2440.       real pt                   ! potential temperature
    2441. 

  2441.       real pt0                  ! reference potential temperature
    2442. 

  2442.       real ptg                  ! ground potential temperature
    2443. 

  2443.       real ua                   ! wind speed
    2444. 

  2444.       real va                   ! wind speed
    2445. 

  2445.       real z0                   ! Roughness length
    2446. 

  2446. 

  2447. ! ----------------------------------------------------------------------
    2448. 

  2448. ! OUTPUT:
    2449. 

  2449. ! ----------------------------------------------------------------------
    2450. 

  2450. ! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal 
    2451. 

  2451. !  surfaces (roofs and street)
    2452. 

  2452. ! The fluxes can be computed as follow: Fluxes of X = B
    2453. 

  2453. !  Example: Momentum fluxes on horizontal surfaces =  uhb_u
    2454. 

  2454.       real uhb                  ! U (wind component) Horizontal surfaces, B (explicit) term
    2455. 

  2455.       real vhb                  ! V (wind component) Horizontal surfaces, B (explicit) term
    2456. 

  2456.       real thb                  ! Temperature        Horizontal surfaces, B (explicit) term
    2457. 

  2457.       real tva                  ! Temperature          Vertical surfaces, A (implicit) term
    2458. 

  2458.       real tvb                  ! Temperature          Vertical surfaces, B (explicit) term
    2459. 

  2459.       real ehb                  ! Energy (TKE)       Horizontal surfaces, B (explicit) term
    2460. 

  2460. 

  2461. 

  2462. ! ----------------------------------------------------------------------
    2463. 

  2463. ! LOCAL:
    2464. 

  2464. ! ----------------------------------------------------------------------
    2465. 

  2465.       real aa
    2466. 

  2466.       real al
    2467. 

  2467.       real buu
    2468. 

  2468.       real c
    2469. 

  2469.       real fbuw
    2470. 

  2470.       real fbpt
    2471. 

  2471.       real fh
    2472. 

  2472.       real fm
    2473. 

  2473.       real ric
    2474. 

  2474.       real tstar
    2475. 

  2475.       real ustar
    2476. 

  2476.       real utot
    2477. 

  2477.       real wstar
    2478. 

  2478.       real zz
    2479. 

  2479.       
    2480. 

  2480.       real b,cm,ch,rr,tol
    2481. 

  2481.       parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
    2482. 

  2482. 

  2483. ! ----------------------------------------------------------------------
    2484. 

  2484. ! END VARIABLES DEFINITIONS
    2485. 

  2485. ! ----------------------------------------------------------------------
    2486. 

  2486. 

  2487. 

  2488. ! computation of the ground temperature
    2489. 

  2489.          
    2490. 

  2490.       utot=(ua**2+va**2)**.5
    2491. 

  2491.         
    2492. 

  2492.       
    2493. 

  2493. !!!! Louis formulation
    2494. 

  2494. !
    2495. 

  2495. ! compute the bulk Richardson Number
    2496. 

  2496. 

  2497.       zz=dz/2.
    2498. 

  2498.    
    2499. 

  2499. !        if(tstar.lt.0.)then
    2500. 

  2500. !         wstar=(-ustar*tstar*g*hii/pt)**(1./3.)
    2501. 

  2501. !        else
    2502. 

  2502. !         wstar=0.
    2503. 

  2503. !        endif
    2504. 

  2504. !        
    2505. 

  2505. !      if (utot.le.0.7*wstar) utot=max(0.7*wstar,0.00001)
    2506. 

  2506. 

  2507.       utot=max(utot,0.01)
    2508. 

  2508.           
    2509. 

  2509.       ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
    2510. 

  2510.               
    2511. 

  2511.       aa=vk/log(zz/z0)
    2512. 

  2512. 

  2513. ! determine the parameters fm and fh for stable, neutral and unstable conditions
    2514. 

  2514. 

  2515.       if(ric.gt.0)then
    2516. 

  2516.          fm=1/(1+0.5*b*ric)**2
    2517. 

  2517.          fh=fm
    2518. 

  2518.       else
    2519. 

  2519.          c=b*cm*aa*aa*(zz/z0)**.5
    2520. 

  2520.          fm=1-b*ric/(1+c*(-ric)**.5)
    2521. 

  2521.          c=c*ch/cm
    2522. 

  2522.          fh=1-b*ric/(1+c*(-ric)**.5)
    2523. 

  2523.       endif
    2524. 

  2524.       
    2525. 

  2525.       fbuw=-aa*aa*utot*utot*fm
    2526. 

  2526.       fbpt=-aa*aa*utot*(pt-ptg)*fh/rr
    2527. 

  2527.                      
    2528. 

  2528.       ustar=(-fbuw)**.5
    2529. 

  2529.       tstar=-fbpt/ustar
    2530. 

  2530. 

  2531.       al=(vk*g_u*tstar)/(pt*ustar*ustar)                      
    2532. 

  2532.       
    2533. 

  2533.       buu=-g_u/pt0*ustar*tstar
    2534. 

  2534.        
    2535. 

  2535.       uhb=-ustar*ustar*ua/utot
    2536. 

  2536.       vhb=-ustar*ustar*va/utot 
    2537. 

  2537.       thb=-ustar*tstar       
    2538. 

  2538. !      thb= 0.      
    2539. 

  2539.       ehb=buu
    2540. 

  2540. !!!!!!!!!!!!!!!
    2541. 

  2541.          
    2542. 

  2542.       return
    2543. 

  2543.       end subroutine flux_flat
    2544. 

  2544. 

  2545. ! ===6=8===============================================================72
    2546. 

  2546. ! ===6=8===============================================================72
    2547. 

  2547. 

  2548.       subroutine icBEP (alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,        &
    2549. 

  2549.                        albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,         &
    2550. 

  2550.                        fww,fwg,fgw,fsw,fws,fsg,                        &
    2551. 

  2551.                        z0g_u,z0r_u,                                    &
    2552. 

  2552.                        nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
    2553. 

  2553.                        nz_u,z_u,                                  &
    2554. 

  2554.                        twini_u,trini_u)                               
    2555. 

  2555. 

  2556. 

  2557.       implicit none 
    2558. 

  2558.       
    2559. 

  2559.    
    2560. 

  2560. !    Building parameters      
    2561. 

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

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

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

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

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

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

  2567.       real twini_u(nurbm)     ! Temperature inside the buildings behind the wall [K]
    2568. 

  2568.       real trini_u(nurbm)     ! Temperature inside the buildings behind the roof [K]
    2569. 

  2569. 

  2570. !    Radiation parameters
    2571. 

  2571.       real albg_u(nurbm)      ! Albedo of the ground
    2572. 

  2572.       real albw_u(nurbm)      ! Albedo of the wall
    2573. 

  2573.       real albr_u(nurbm)      ! Albedo of the roof
    2574. 

  2574.       real emg_u(nurbm)       ! Emissivity of ground
    2575. 

  2575.       real emw_u(nurbm)       ! Emissivity of wall
    2576. 

  2576.       real emr_u(nurbm)       ! Emissivity of roof
    2577. 

  2577. 

  2578. !    Roughness parameters
    2579. 

  2579.       real z0g_u(nurbm)       ! The ground's roughness length      
    2580. 

  2580.       real z0r_u(nurbm)       ! The roof's roughness length     
    2581. 

    2582. 

  2582. !    Street parameters
    2583. 

  2583.       integer nd_u(nurbm)     ! Number of street direction for each urban class
    2584. 

  2584. 

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

  2586.       real drst_u(ndm,nurbm)  ! Street direction [degree]
    2587. 

  2587.       real ws_u(ndm,nurbm)    ! Street width [m]
    2588. 

  2588.       real bs_u(ndm,nurbm)    ! Building width [m]
    2589. 

  2589.       real h_b(nz_um,nurbm)   ! Bulding's heights [m]
    2590. 

  2590.       real d_b(nz_um,nurbm)   ! The probability that a building has an height h_b
    2591. 

  2591. ! -----------------------------------------------------------------------
    2592. 

  2592. !     Output
    2593. 

  2593. !------------------------------------------------------------------------
    2594. 

  2594. 

  2595. 

  2596. 

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

  2598. !   and the short wave radation. They are the part of radiation from a surface
    2599. 

  2599. !   or from the sky to another surface.
    2600. 

  2600.       real fww(nz_um,nz_um,ndm,nurbm)        !  from wall to wall
    2601. 

  2601.       real fwg(nz_um,ndm,nurbm)              !  from wall to ground
    2602. 

  2602.       real fgw(nz_um,ndm,nurbm)              !  from ground to wall
    2603. 

  2603.       real fsw(nz_um,ndm,nurbm)              !  from sky to wall
    2604. 

  2604.       real fws(nz_um,ndm,nurbm)              !  from wall to sky
    2605. 

  2605.       real fsg(ndm,nurbm)                    !  from sky to ground
    2606. 

  2606. 

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

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

  2609.         
    2610. 

  2610. !    Grid parameters
    2611. 

  2611.       integer nz_u(nurbm)     ! Number of layer in the urban grid
    2612. 

  2612.       real z_u(nz_um)       ! Height of the urban grid levels
    2613. 

  2613. 

  2614. 

  2615. ! -----------------------------------------------------------------------
    2616. 

  2616. !     Local
    2617. 

  2617. !------------------------------------------------------------------------
    2618. 

  2618. 

  2619.       integer iz_u,id,ilu,iurb
    2620. 

  2620. 

  2621.       real dtot
    2622. 

  2622.       real hbmax
    2623. 

  2623. 

  2624. !------------------------------------------------------------------------
    2625. 

  2625.      
    2626. 

  2626. 

  2627. ! -----------------------------------------------------------------------
    2628. 

  2628. !     This routine initialise the urban paramters for the BEP module
    2629. 

  2629. !------------------------------------------------------------------------
    2630. 

  2630. !
    2631. 

  2631. !Initialize variables
    2632. 

  2632. !
    2633. 

  2633.       nz_u=0
    2634. 

  2634.       z_u=0.
    2635. 

  2635.       ss_u=0.
    2636. 

  2636.       pb_u=0.
    2637. 

  2637.       fww=0.
    2638. 

  2638.       fwg=0.
    2639. 

  2639.       fgw=0.
    2640. 

  2640.       fsw=0.
    2641. 

  2641.       fws=0.
    2642. 

  2642.       fsg=0. 
    2643. 

  2643.         
    2644. 

  2644. ! Computation of the urban levels height
    2645. 

  2645. 

  2646.       z_u(1)=0.
    2647. 

  2647.      
    2648. 

  2648.       do iz_u=1,nz_um-1
    2649. 

  2649.          z_u(iz_u+1)=z_u(iz_u)+dz_u
    2650. 

  2650.       enddo
    2651. 

  2651.       
    2652. 

  2652. ! Normalisation of the building density
    2653. 

  2653. 

  2654.       do iurb=1,nurbm
    2655. 

  2655.          dtot=0.
    2656. 

  2656.          do ilu=1,nz_um
    2657. 

  2657.             dtot=dtot+d_b(ilu,iurb)
    2658. 

  2658.          enddo
    2659. 

  2659.          do ilu=1,nz_um
    2660. 

  2660.             d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
    2661. 

  2661.          enddo
    2662. 

  2662.       enddo      
    2663. 

  2663. 

  2664. ! Compute the view factors, pb and ss 
    2665. 

  2665.       
    2666. 

  2666.       do iurb=1,nurbm         
    2667. 

  2667.          hbmax=0.
    2668. 

  2668.          nz_u(iurb)=0
    2669. 

  2669.          do ilu=1,nz_um
    2670. 

  2670.             if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
    2671. 

  2671.          enddo
    2672. 

  2672.          
    2673. 

  2673.          do iz_u=1,nz_um-1
    2674. 

  2674.             if(z_u(iz_u+1).gt.hbmax)go to 10
    2675. 

  2675.          enddo
    2676. 

  2676.          
    2677. 

  2677.  10      continue
    2678. 

  2678.          nz_u(iurb)=iz_u+1
    2679. 

  2679. 

  2680.          do id=1,nd_u(iurb)
    2681. 

  2681. 

  2682.             call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb),   &    
    2683. 

  2683.                             z_u,ws_u(id,iurb),                      &    
    2684. 

  2684.                             fww,fwg,fgw,fsg,fsw,fws) 
    2685. 

  2685. 

  2686.             do iz_u=1,nz_u(iurb)
    2687. 

  2687.                ss_u(iz_u,iurb)=0.
    2688. 

  2688.                do ilu=1,nz_um
    2689. 

  2689.                   if(z_u(iz_u).le.h_b(ilu,iurb)                      &    
    2690. 

  2690.                     .and.z_u(iz_u+1).gt.h_b(ilu,iurb))then            
    2691. 

  2691.                         ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
    2692. 

  2692.                   endif 
    2693. 

  2693.                enddo
    2694. 

  2694.             enddo
    2695. 

  2695. 

  2696.             pb_u(1,iurb)=1.
    2697. 

  2697.             do iz_u=1,nz_u(iurb)
    2698. 

  2698.                pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
    2699. 

  2699.             enddo
    2700. 

  2700. 

  2701.          enddo
    2702. 

  2702.       end do
    2703. 

  2703.      
    2704. 

  2704.                   
    2705. 

  2705.       return       
    2706. 

  2706.       end subroutine icBEP
    2707. 

  2707. 

  2708. ! ===6=8===============================================================72
    2709. 

  2709. ! ===6=8===============================================================72
    2710. 

  2710. 

  2711.       subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws) 
    2712. 

  2712.      
    2713. 

  2713.       implicit none
    2714. 

  2714. 

  2715.  
    2716. 

  2716. 

  2717. ! -----------------------------------------------------------------------
    2718. 

  2718. !     Input
    2719. 

  2719. !------------------------------------------------------------------------
    2720. 

  2720. 

  2721.       integer iurb            ! Number of the urban class
    2722. 

  2722.       integer nz_u            ! Number of levels in the urban grid
    2723. 

  2723.       integer id              ! Street direction number
    2724. 

  2724.       real ws                 ! Street width
    2725. 

  2725.       real z(nz_um)         ! Height of the urban grid levels
    2726. 

  2726.       real dxy                ! Street lenght
    2727. 

  2727. 

  2728. 

  2729. ! -----------------------------------------------------------------------
    2730. 

  2730. !     Output
    2731. 

  2731. !------------------------------------------------------------------------
    2732. 

  2732. 

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

  2734. !   and the short wave radation. They are the part of radiation from a surface
    2735. 

  2735. !   or from the sky to another surface.
    2736. 

  2736. 

  2737.       real fww(nz_um,nz_um,ndm,nurbm)            !  from wall to wall
    2738. 

  2738.       real fwg(nz_um,ndm,nurbm)                  !  from wall to ground
    2739. 

  2739.       real fgw(nz_um,ndm,nurbm)                  !  from ground to wall
    2740. 

  2740.       real fsw(nz_um,ndm,nurbm)                  !  from sky to wall
    2741. 

  2741.       real fws(nz_um,ndm,nurbm)                  !  from wall to sky
    2742. 

  2742.       real fsg(ndm,nurbm)                        !  from sky to ground
    2743. 

  2743. 

  2744. 

  2745. ! -----------------------------------------------------------------------
    2746. 

  2746. !     Local
    2747. 

  2747. !------------------------------------------------------------------------
    2748. 

  2748. 

  2749.       integer jz,iz
    2750. 

  2750. 

  2751.       real hut
    2752. 

  2752.       real f1,f2,f12,f23,f123,ftot
    2753. 

  2753.       real fprl,fnrm
    2754. 

  2754.       real a1,a2,a3,a4,a12,a23,a123
    2755. 

  2755. 

  2756. ! -----------------------------------------------------------------------
    2757. 

  2757. !     This routine calculates the view factors
    2758. 

  2758. !------------------------------------------------------------------------
    2759. 

  2759.         
    2760. 

  2760.       hut=z(nz_u+1)
    2761. 

  2761.         
    2762. 

  2762.       do jz=1,nz_u      
    2763. 

  2763.       
    2764. 

  2764. ! radiation from wall to wall
    2765. 

  2765.        
    2766. 

  2766.          do iz=1,nz_u
    2767. 

  2767.      
    2768. 

  2768.             call fprls (fprl,dxy,abs(z(jz+1)-z(iz  )),ws)
    2769. 

  2769.             f123=fprl
    2770. 

  2770.             call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
    2771. 

  2771.             f23=fprl
    2772. 

  2772.             call fprls (fprl,dxy,abs(z(jz  )-z(iz  )),ws)
    2773. 

  2773.             f12=fprl
    2774. 

  2774.             call fprls (fprl,dxy,abs(z(jz  )-z(iz+1)),ws)
    2775. 

  2775.             f2 = fprl
    2776. 

  2776.        
    2777. 

  2777.             a123=dxy*(abs(z(jz+1)-z(iz  )))
    2778. 

  2778.             a12 =dxy*(abs(z(jz  )-z(iz  )))
    2779. 

  2779.             a23 =dxy*(abs(z(jz+1)-z(iz+1)))
    2780. 

  2780.             a1  =dxy*(abs(z(iz+1)-z(iz  )))
    2781. 

  2781.             a2  =dxy*(abs(z(jz  )-z(iz+1)))
    2782. 

  2782.             a3  =dxy*(abs(z(jz+1)-z(jz  )))
    2783. 

  2783.        
    2784. 

  2784.             ftot=0.5*(a123*f123-a23*f23-a12*f12+a2*f2)/a1
    2785. 

  2785.        
    2786. 

  2786.             fww(iz,jz,id,iurb)=ftot*a1/a3
    2787. 

  2787. 

  2788.          enddo 
    2789. 

  2789. 

  2790. ! radiation from ground to wall
    2791. 

  2791.        
    2792. 

  2792.          call fnrms (fnrm,z(jz+1),dxy,ws)
    2793. 

  2793.          f12=fnrm
    2794. 

  2794.          call fnrms (fnrm,z(jz)  ,dxy,ws)
    2795. 

  2795.          f1=fnrm
    2796. 

  2796.        
    2797. 

  2797.          a1 = ws*dxy
    2798. 

  2798.          
    2799. 

  2799.          a12= ws*dxy
    2800. 

  2800.        
    2801. 

  2801.          a4=(z(jz+1)-z(jz))*dxy
    2802. 

  2802.        
    2803. 

  2803.          ftot=(a12*f12-a12*f1)/a1
    2804. 

  2804.                     
    2805. 

  2805.          fgw(jz,id,iurb)=ftot*a1/a4
    2806. 

  2806.      
    2807. 

  2807. !  radiation from sky to wall
    2808. 

  2808.      
    2809. 

  2809.          call fnrms(fnrm,hut-z(jz)  ,dxy,ws)
    2810. 

  2810.          f12 = fnrm
    2811. 

  2811.          call fnrms (fnrm,hut-z(jz+1),dxy,ws)
    2812. 

  2812.          f1 =fnrm
    2813. 

  2813.        
    2814. 

  2814.          a1 = ws*dxy
    2815. 

  2815.        
    2816. 

  2816.          a12= ws*dxy
    2817. 

  2817.               
    2818. 

  2818.          a4 = (z(jz+1)-z(jz))*dxy
    2819. 

  2819.        
    2820. 

  2820.          ftot=(a12*f12-a12*f1)/a1
    2821. 

  2821.         
    2822. 

  2822.          fsw(jz,id,iurb)=ftot*a1/a4       
    2823. 

  2823.       
    2824. 

  2824.       enddo
    2825. 

  2825. 

  2826. ! radiation from wall to sky      
    2827. 

  2827.       do iz=1,nz_u
    2828. 

  2828.        call fnrms(fnrm,ws,dxy,hut-z(iz))
    2829. 

  2829.        f12=fnrm
    2830. 

  2830.        call fnrms(fnrm,ws,dxy,hut-z(iz+1))
    2831. 

  2831.        f1=fnrm
    2832. 

  2832.        a1 = (z(iz+1)-z(iz))*dxy
    2833. 

  2833.        a2 = (hut-z(iz+1))*dxy
    2834. 

  2834.        a12= (hut-z(iz))*dxy
    2835. 

  2835.        a4 = ws*dxy
    2836. 

  2836.        ftot=(a12*f12-a2*f1)/a1
    2837. 

  2837.        fws(iz,id,iurb)=ftot*a1/a4 
    2838. 

  2838.  
    2839. 

  2839.       enddo
    2840. 

  2840. !!!!!!!!!!!!!
    2841. 

  2841. 

  2842. 

  2843.        do iz=1,nz_u
    2844. 

  2844. 

  2845. ! radiation from wall to ground
    2846. 

  2846.       
    2847. 

  2847.          call fnrms (fnrm,ws,dxy,z(iz+1))
    2848. 

  2848.          f12=fnrm
    2849. 

  2849.          call fnrms (fnrm,ws,dxy,z(iz  ))
    2850. 

  2850.          f1 =fnrm
    2851. 

  2851.          
    2852. 

  2852.          a1= (z(iz+1)-z(iz) )*dxy
    2853. 

  2853.        
    2854. 

  2854.          a2 = z(iz)*dxy
    2855. 

  2855.          a12= z(iz+1)*dxy
    2856. 

  2856.          a4 = ws*dxy
    2857. 

  2857. 

  2858.          ftot=(a12*f12-a2*f1)/a1        
    2859. 

  2859.                     
    2860. 

  2860.          fwg(iz,id,iurb)=ftot*a1/a4
    2861. 

  2861.         
    2862. 

  2862.       enddo
    2863. 

  2863. 

  2864. ! radiation from sky to ground
    2865. 

  2865.       
    2866. 

  2866.       call fprls (fprl,dxy,ws,hut)
    2867. 

  2867.       fsg(id,iurb)=fprl
    2868. 

  2868. 

  2869.       return
    2870. 

  2870.       end subroutine view_factors
    2871. 

  2871. 

  2872. ! ===6=8===============================================================72
    2873. 

  2873. ! ===6=8===============================================================72
    2874. 

  2874. 

  2875.       SUBROUTINE fprls (fprl,a,b,c)
    2876. 

  2876. 

  2877.       implicit none
    2878. 

  2878. 

  2879.      
    2880. 

  2880.             
    2881. 

  2881.       real a,b,c
    2882. 

  2882.       real x,y
    2883. 

  2883.       real fprl
    2884. 

  2884. 

  2885. 

  2886.       x=a/c
    2887. 

  2887.       y=b/c
    2888. 

  2888.       
    2889. 

  2889.       if(a.eq.0.or.b.eq.0.)then
    2890. 

  2890.        fprl=0.
    2891. 

  2891.       else
    2892. 

  2892.        fprl=log( ( (1.+x**2)*(1.+y**2)/(1.+x**2+y**2) )**.5)+  &
    2893. 

  2893.            y*((1.+x**2)**.5)*atan(y/((1.+x**2)**.5))+          &  
    2894. 

  2894.            x*((1.+y**2)**.5)*atan(x/((1.+y**2)**.5))-          &   
    2895. 

  2895.            y*atan(y)-x*atan(x)
    2896. 

  2896.        fprl=fprl*2./(pi*x*y)
    2897. 

  2897.       endif
    2898. 

  2898.       
    2899. 

  2899.       return
    2900. 

  2900.       end subroutine fprls
    2901. 

  2901. 

  2902. ! ===6=8===============================================================72     
    2903. 

  2903. ! ===6=8===============================================================72
    2904. 

  2904. 

  2905.       SUBROUTINE fnrms (fnrm,a,b,c)
    2906. 

  2906. 

  2907.       implicit none
    2908. 

  2908. 

  2909. 

  2910. 

  2911.       real a,b,c
    2912. 

  2912.       real x,y,z,a1,a2,a3,a4,a5,a6
    2913. 

  2913.       real fnrm
    2914. 

  2914.       
    2915. 

  2915.       x=a/b
    2916. 

  2916.       y=c/b
    2917. 

  2917.       z=x**2+y**2
    2918. 

  2918.       
    2919. 

  2919.       if(y.eq.0.or.x.eq.0)then
    2920. 

  2920.        fnrm=0.
    2921. 

  2921.       else
    2922. 

  2922.        a1=log( (1.+x*x)*(1.+y*y)/(1.+z) )
    2923. 

  2923.        a2=y*y*log(y*y*(1.+z)/z/(1.+y*y) )
    2924. 

  2924.        a3=x*x*log(x*x*(1.+z)/z/(1.+x*x) )
    2925. 

  2925.        a4=y*atan(1./y)
    2926. 

  2926.        a5=x*atan(1./x)
    2927. 

  2927.        a6=sqrt(z)*atan(1./sqrt(z))
    2928. 

  2928.        fnrm=0.25*(a1+a2+a3)+a4+a5-a6
    2929. 

  2929.        fnrm=fnrm/(pi*y)
    2930. 

  2930.       endif
    2931. 

  2931.       
    2932. 

  2932.       return
    2933. 

  2933.       end subroutine fnrms
    2934. 

  2934.   ! ===6=8===============================================================72  
    2935. 

  2935.      
    2936. 

  2936.       SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
    2937. 

  2937.                 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
    2938. 

  2938.                 emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
    2939. 

  2939. 

  2940. ! initialization routine, where the variables from the table are read
    2941. 

  2941. 

  2942.       implicit none
    2943. 

  2943.       
    2944. 

  2944.       integer iurb            ! urban class number
    2945. 

  2945. !    Building parameters      
    2946. 

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

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

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

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

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

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

  2952.       real twini_u(nurbm)     ! Temperature inside the buildings behind the wall [K]
    2953. 

  2953.       real trini_u(nurbm)     ! Temperature inside the buildings behind the roof [K]
    2954. 

  2954.       real tgini_u(nurbm)     ! Initial road temperature
    2955. 

  2955. 

  2956. !    Radiation parameters
    2957. 

  2957.       real albg_u(nurbm)      ! Albedo of the ground
    2958. 

  2958.       real albw_u(nurbm)      ! Albedo of the wall
    2959. 

  2959.       real albr_u(nurbm)      ! Albedo of the roof
    2960. 

  2960.       real emg_u(nurbm)       ! Emissivity of ground
    2961. 

  2961.       real emw_u(nurbm)       ! Emissivity of wall
    2962. 

  2962.       real emr_u(nurbm)       ! Emissivity of roof
    2963. 

  2963. 

  2964. !    Roughness parameters
    2965. 

  2965.       real z0g_u(nurbm)       ! The ground's roughness length      
    2966. 

  2966.       real z0r_u(nurbm)       ! The roof's roughness length
    2967. 

    2968. 

  2968. !    Street parameters
    2969. 

  2969.       integer nd_u(nurbm)     ! Number of street direction for each urban class
    2970. 

  2970. 

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

  2972.       real drst_u(ndm,nurbm)  ! Street direction [degree]
    2973. 

  2973.       real ws_u(ndm,nurbm)    ! Street width [m]
    2974. 

  2974.       real bs_u(ndm,nurbm)    ! Building width [m]
    2975. 

  2975.       real h_b(nz_um,nurbm)   ! Bulding's heights [m]
    2976. 

  2976.       real d_b(nz_um,nurbm)   ! The probability that a building has an height h_b
    2977. 

  2977. 

  2978.       integer i,iu
    2979. 

  2979.       integer nurb ! number of urban classes used
    2980. 

  2980. 

  2981. !
    2982. 

  2982. !Initialize some variables
    2983. 

  2983. !
    2984. 

  2984.        
    2985. 

  2985.        h_b=0.
    2986. 

  2986.        d_b=0.
    2987. 

  2987. 

  2988.        nurb=ICATE
    2989. 

  2989.        do iu=1,nurb                         
    2990. 

  2990.           nd_u(iu)=0
    2991. 

  2991.        enddo
    2992. 

  2992. 

  2993.        csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
    2994. 

  2994.        csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
    2995. 

  2995.        csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
    2996. 

  2996.        do i=1,icate
    2997. 

  2997.          alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
    2998. 

  2998.          alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
    2999. 

  2999.          alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
    3000. 

  3000.        enddo
    3001. 

  3001.        twini_u=TBLEND_TBL
    3002. 

  3002.        trini_u=TRLEND_TBL
    3003. 

  3003.        tgini_u=TGLEND_TBL
    3004. 

  3004.        albw_u=ALBB_TBL
    3005. 

  3005.        albr_u=ALBR_TBL
    3006. 

  3006.        albg_u=ALBG_TBL
    3007. 

  3007.        emw_u=EPSB_TBL
    3008. 

  3008.        emr_u=EPSR_TBL
    3009. 

  3009.        emg_u=EPSG_TBL
    3010. 

  3010.        z0r_u=Z0R_TBL
    3011. 

  3011.        z0g_u=Z0G_TBL
    3012. 

  3012.        nd_u=NUMDIR_TBL
    3013. 

  3013.        do iu=1,icate
    3014. 

  3014.               if(ndm.lt.nd_u(iu))then
    3015. 

  3015.                 write(*,*)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
    3016. 

  3016.                 write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
    3017. 

  3017.                 stop
    3018. 

  3018.               endif
    3019. 

  3019.          do i=1,nd_u(iu)
    3020. 

  3020.            drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
    3021. 

  3021.            ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
    3022. 

  3022.            bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
    3023. 

  3023.          enddo
    3024. 

  3024.        enddo
    3025. 

  3025.        do iu=1,ICATE
    3026. 

  3026.           if(nz_um.lt.numhgt_tbl(iu)+3)then
    3027. 

  3027.               write(*,*)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
    3028. 

  3028.               write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
    3029. 

  3029.               stop
    3030. 

  3030.           endif
    3031. 

  3031.          do i=1,NUMHGT_TBL(iu)
    3032. 

  3032.            h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
    3033. 

  3033.            d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
    3034. 

  3034.          enddo
    3035. 

  3035.        enddo
    3036. 

  3036. 

  3037.        do i=1,ndm
    3038. 

  3038.         do iu=1,nurbm
    3039. 

  3039.          strd_u(i,iu)=100000.
    3040. 

  3040.         enddo
    3041. 

  3041.        enddo
    3042. 

  3042.          
    3043. 

  3043.        return
    3044. 

  3044.       END SUBROUTINE init_para
    3045. 

  3045. !==============================================================
    3046. 

  3046. 

  3047. !==============================================================
    3048. 

  3048.       subroutine angle(along,alat,day,realt,zr,deltar,ah)
    3049. 

  3049. !     ----------------     
    3050. 

  3050. !      
    3051. 

  3051. !         Computation of the solar angles
    3052. 

  3052. !         schayes (1982,atm.  env. , p1407)
    3053. 

  3053. ! Inputs
    3054. 

  3054. !========================
    3055. 

  3055. ! along=longitud
    3056. 

  3056. ! alat=latitude
    3057. 

  3057. ! day=julian day (from the beginning of the year)
    3058. 

  3058. ! realt= time GMT in hours
    3059. 

  3059. ! Outputs
    3060. 

  3060. !============================
    3061. 

  3061. ! zr=solar zenith angle
    3062. 

  3062. ! deltar=declination angle
    3063. 

  3063. ! ah=hour angle
    3064. 

  3064. !===============================
    3065. 

  3065. 

  3066.       implicit none
    3067. 

  3067.       real along,alat, realt, zr, deltar, ah, arg
    3068. 

  3068.       real rad,om,radh,initt, pii, drad, alongt, cphi, sphi
    3069. 

  3069.       real c1, c2, c3, s1, s2, s3, delta, rmsr2, cd, sid 
    3070. 

  3070.       real et, ahor, chor, coznt 
    3071. 

  3071.       integer day 
    3072. 

  3072. 

  3073. 

  3074.       data rad,om,radh,initt/0.0174533,0.0172142,0.26179939,0/
    3075. 

  3075.  
    3076. 

  3076.        zr=0.
    3077. 

  3077.        deltar=0.
    3078. 

  3078.        ah=0.
    3079. 

  3079. 

  3080.        pii = 3.14159265358979312
    3081. 

  3081.        drad = pii/180.
    3082. 

  3082.       
    3083. 

  3083.        alongt=along/15.
    3084. 

  3084.        cphi=cos(alat*drad)
    3085. 

  3085.        sphi=sin(alat*drad)
    3086. 

  3086. !
    3087. 

  3087. !     declination
    3088. 

  3088. !
    3089. 

  3089.        arg=om*day
    3090. 

  3090.        c1=cos(arg)
    3091. 

  3091.        c2=cos(2.*arg)
    3092. 

  3092.        c3=cos(3.*arg)
    3093. 

  3093.        s1=sin(arg)
    3094. 

  3094.        s2=sin(2.*arg)
    3095. 

  3095.        s3=sin(3.*arg)
    3096. 

  3096.        delta=0.33281-22.984*c1-0.3499*c2-0.1398*c3+3.7872*s1+0.03205*s2+0.07187*s3
    3097. 

  3097.        rmsr2=(1./(1.-0.01673*c1))**2
    3098. 

  3098.        deltar=delta*rad
    3099. 

  3099.        cd=cos(deltar)
    3100. 

  3100.        sid=sin(deltar)
    3101. 

  3101. !
    3102. 

  3102. !     time equation in hours
    3103. 

  3103. !
    3104. 

  3104.        et=0.0072*c1-0.0528*c2-0.0012*c3-0.1229*s1-0.1565*s2-0.0041*s3
    3105. 

  3105. !
    3106. 

  3106. !
    3107. 

  3107. !   hour angle  
    3108. 

  3108. !
    3109. 

  3109.       
    3110. 

  3110. !      ifh=0
    3111. 

  3111.       
    3112. 

  3112.  !     ahor=realt-12.+ifh+et+alongt
    3113. 

  3113.       ahor=realt-12.+et+alongt
    3114. 

  3114.       ah=ahor*radh
    3115. 

  3115.       chor=cos(ah)
    3116. 

  3116. !
    3117. 

  3117. !    zenith angle
    3118. 

  3118. !
    3119. 

  3119.       coznt=sphi*sid+cphi*cd*chor
    3120. 

  3120.       
    3121. 

  3121.        zr=acos(coznt)
    3122. 

  3122. 

  3123.       return
    3124. 

  3124. 

  3125.       END SUBROUTINE angle
    3126. 

  3126. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    3127. 

  3127. 

  3128.       subroutine upward_rad(nd_u,iurb,nz_u,ws,bs,sigma,fsw,fsg,pb,ss,           &
    3129. 

  3129.                        tg,emg_u,albg_u,rlg,rsg,sfg,                          & 
    3130. 

  3130.                        tw,emw_u,albw_u,rlw,rsw,sfw,                          & 
    3131. 

  3131.                        tr,emr_u,albr_u,rld,rs, sfr,                            & 
    3132. 

  3132.                        rs_abs,rl_up,emiss,grdflx_urb)
    3133. 

  3133. !
    3134. 

  3134. ! IN this surboutine we compute the upward longwave flux, and the albedo
    3135. 

  3135. ! needed for the radiation scheme
    3136. 

  3136. !
    3137. 

  3137.                        implicit none
    3138. 

  3138. 

  3139. !
    3140. 

  3140. !INPUT VARIABLES
    3141. 

  3141. !
    3142. 

  3142.       real rsw(2*ndm,nz_um)        ! Short wave radiation at the wall for a given canyon direction [W/m2]
    3143. 

  3143.       real rlw(2*ndm,nz_um)         ! Long wave radiation at the walls for a given canyon direction [W/m2]
    3144. 

  3144.       real rsg(ndm)                   ! Short wave radiation at the canyon for a given canyon direction [W/m2]
    3145. 

  3145.       real rlg(ndm)                   ! Long wave radiation at the ground for a given canyon direction [W/m2]
    3146. 

  3146.       real rs                        ! Short wave radiation at the horizontal surface from the sun [W/m²]  
    3147. 

  3147.       real sfw(2*ndm,nz_um)      ! Sensible heat flux from walls [W/m²]
    3148. 

  3148.       real sfg(ndm)              ! Sensible heat flux from ground (road) [W/m²]
    3149. 

  3149.       real sfr(ndm,nz_um)      ! Sensible heat flux from roofs [W/m²]                      
    3150. 

  3150.       real rld                        ! Long wave radiation from the sky [W/m²]
    3151. 

  3151.       real albg_u                    ! albedo of the ground/street
    3152. 

  3152.       real albw_u                    ! albedo of the walls
    3153. 

  3153.       real albr_u                    ! albedo of the roof 
    3154. 

  3154.       real ws(ndm)                        ! width of the street
    3155. 

  3155.       real bs(ndm)
    3156. 

  3156.                         ! building size
    3157. 

  3157.       real pb(nz_um)                ! Probability to have a building with an height equal or higher   
    3158. 

  3158.       integer nz_u
    3159. 

  3159.       real ss(nz_um)                ! Probability to have a building of a given height
    3160. 

  3160.       real sigma                       
    3161. 

  3161.       real emg_u                       ! emissivity of the street
    3162. 

  3162.       real emw_u                       ! emissivity of the wall
    3163. 

  3163.       real emr_u                       ! emissivity of the roof
    3164. 

  3164.       real fsw(nz_um,ndm,nurbm)        ! View factors from sky to wall                
    3165. 

  3165.       real fsg(ndm,nurbm)                      ! groud to sky view factor
    3166. 

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

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

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

  3169.       integer iurb ! urban class
    3170. 

  3170.       integer id ! street direction
    3171. 

  3171.       integer nd_u ! number of street directions
    3172. 

  3172. !OUTPUT/INPUT
    3173. 

  3173.       real rs_abs  ! absrobed solar radiationfor this street direction
    3174. 

  3174.       real rl_up   ! upward longwave radiation for this street direction
    3175. 

  3175.       real emiss ! mean emissivity
    3176. 

  3176.       real grdflx_urb ! ground heat flux 
    3177. 

  3177. !LOCAL
    3178. 

  3178.       integer iz,iw
    3179. 

  3179.       real rl_inc,rl_emit
    3180. 

  3180.       real gfl
    3181. 

  3181.       integer ix,iy,iwrong
    3182. 

  3182. 

  3183.          iwrong=1
    3184. 

  3184.       do iz=1,nz_u+1
    3185. 

  3185.       do id=1,nd_u
    3186. 

  3186.       do iw=1,nwr_u
    3187. 

  3187.         if(tr(id,iz,iw).lt.100.)then
    3188. 

  3188.               write(*,*)'in upward_rad ',iz,id,iw,tr(id,iz,iw) 
    3189. 

  3189.               iwrong=0
    3190. 

  3190.         endif
    3191. 

  3191.       enddo
    3192. 

  3192.       enddo
    3193. 

  3193.       enddo
    3194. 

  3194.            if(iwrong.eq.0)stop
    3195. 

  3195. 

  3196.       rl_up=0.
    3197. 

  3197.  
    3198. 

  3198.       rs_abs=0.
    3199. 

  3199.       rl_inc=0.
    3200. 

  3200.       emiss=0.
    3201. 

  3201.       rl_emit=0.
    3202. 

  3202.       grdflx_urb=0.
    3203. 

  3203.       do id=1,nd_u          
    3204. 

  3204.        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
    3205. 

  3205.        rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u       
    3206. 

  3206.        rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/nd_u
    3207. 

  3207.          gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg(id,ng_u)**4.)+sfg(id)
    3208. 

  3208.          grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u  
    3209. 

  3209.  
    3210. 

  3210.           do iz=2,nz_u
    3211. 

  3211.             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
    3212. 

  3212.             rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u            
    3213. 

  3213.             rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
    3214. 

  3214.             gfl=(1.-albr_u)*rs+emr_u*rld-emr_u*sigma*(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
    3215. 

  3215.             grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
    3216. 

  3216.          enddo
    3217. 

  3217.            
    3218. 

  3218.          do iz=1,nz_u            
    3219. 

  3219.             rl_emit=rl_emit-(emw_u*sigma*( tw(2*id-1,iz,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+          &
    3220. 

  3220.                (1-emw_u)*( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
    3221. 

  3221.             rl_inc=rl_inc+(( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
    3222. 

  3222.             rs_abs=rs_abs+((1.-albw_u)*( rsw(2*id-1,iz)+rsw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u 
    3223. 

  3223.             gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) )   &
    3224. 

  3224.              -emw_u*sigma*( tw(2*id-1,iz,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))            
    3225. 

  3225.             grdflx_urb=grdflx_urb-gfl*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
    3226. 

  3226.          enddo
    3227. 

  3227.           
    3228. 

  3228.       enddo
    3229. 

  3229.         emiss=(emg_u+emw_u+emr_u)/3.
    3230. 

  3230.         rl_up=(rl_inc+rl_emit)-rld
    3231. 

  3231.        
    3232. 

  3232.          
    3233. 

  3233.       return
    3234. 

  3234. 

  3235.       END SUBROUTINE upward_rad
    3236. 

  3236. 

  3237. !====6=8===============================================================72         
    3238. 

  3238. !====6=8===============================================================72 
    3239. 

  3239. END MODULE module_sf_bep
    3240. 

  3240.