PageRenderTime 363ms CodeModel.GetById 25ms RepoModel.GetById 1ms 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

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

    

  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.

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