| blob | 41b02135117c10c21a780df85d866ef8ff81547f |
1 MODULE module_sf_bep
2 #if defined(mpas)
3 use mpas_atmphys_utilities, only: physics_error_fatal
4 #define FATAL_ERROR(M) call physics_error_fatal( M )
5 #else
6 use module_wrf_error
7 #define FATAL_ERROR(M) call wrf_error_fatal( M)
8 !USE module_model_constants
9 #endif
10 USE module_sf_urban
11 USE module_bep_bem_helper, ONLY: nurbm
13 ! SGClarke 09/11/2008
14 ! Access urban_param.tbl values through calling urban_param_init in module_physics_init
15 ! for CASE (BEPSCHEME) select sf_urban_physics
16 !
17 ! -----------------------------------------------------------------------
18 ! Dimension for the array used in the BEP module
19 ! -----------------------------------------------------------------------
20 integer nurbmax ! Maximum number of urban classes
21 parameter (nurbmax=11)
23 integer ndm ! Maximum number of street directions
24 parameter (ndm=2)
26 integer nz_um ! Maximum number of vertical levels in the urban grid
27 parameter(nz_um=18)
29 integer ng_u ! Number of grid levels in the ground
30 parameter (ng_u=10)
31 integer nwr_u ! Number of grid levels in the walls or roofs
32 parameter (nwr_u=10)
34 real dz_u ! Urban grid resolution
35 parameter (dz_u=5.)
37 ! The change of ng_u, nwr_u should be done in agreement with the block data
38 ! in the routine "surf_temp"
39 ! -----------------------------------------------------------------------
40 ! Constant used in the BEP module
41 ! -----------------------------------------------------------------------
43 real vk ! von Karman constant
44 real g_u ! Gravity acceleration
45 real pi !
46 real r ! Perfect gas constant
47 real cp_u ! Specific heat at constant pressure
48 real rcp_u !
49 real sigma !
50 real p0 ! Reference pressure at the sea level
51 real cdrag ! Drag force constant
52 parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)
53 parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4)
55 ! -----------------------------------------------------------------------
60 CONTAINS
62 subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, &
63 th_phy,rho,p_phy,swdown,glw, &
64 gmt,julday,xlong,xlat, &
65 declin_urb,cosz_urb2d,omg_urb2d, &
66 num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &
67 urban_map_zd, urban_map_zdf, urban_map_bd, urban_map_wd, &
68 urban_map_gbd, urban_map_fbd, &
69 num_urban_hi, &
70 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
71 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
72 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
73 a_u,a_v,a_t,a_e,b_u,b_v, &
74 b_t,b_e,b_q,dlg,dl_u,sf,vl, &
75 rl_up,rs_abs,emiss,grdflx_urb, &
76 ids,ide, jds,jde, kds,kde, &
77 ims,ime, jms,jme, kms,kme, &
78 its,ite, jts,jte, kts,kte)
80 implicit none
82 !------------------------------------------------------------------------
83 ! Input
84 !------------------------------------------------------------------------
85 INTEGER :: ids,ide, jds,jde, kds,kde, &
86 ims,ime, jms,jme, kms,kme, &
87 its,ite, jts,jte, kts,kte, &
88 itimestep
91 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: DZ8W
92 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: P_PHY
93 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: RHO
94 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: TH_PHY
95 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: T_PHY
96 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U_PHY
97 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V_PHY
98 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U
99 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V
100 REAL, DIMENSION( ims:ime , jms:jme ) :: GLW
101 REAL, DIMENSION( ims:ime , jms:jme ) :: swdown
102 REAL, DIMENSION( ims:ime, jms:jme ) :: UST
103 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: UTYPE_URB2D
104 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: FRC_URB2D
105 REAL, INTENT(IN ) :: GMT
106 INTEGER, INTENT(IN ) :: JULDAY
107 REAL, DIMENSION( ims:ime, jms:jme ), &
108 INTENT(IN ) :: XLAT, XLONG
109 REAL, INTENT(IN) :: DECLIN_URB
110 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
111 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
112 INTEGER, INTENT(IN ) :: num_urban_ndm
113 INTEGER, INTENT(IN ) :: urban_map_zrd
114 INTEGER, INTENT(IN ) :: urban_map_zwd
115 INTEGER, INTENT(IN ) :: urban_map_gd
116 INTEGER, INTENT(IN ) :: urban_map_zd
117 INTEGER, INTENT(IN ) :: urban_map_zdf
118 INTEGER, INTENT(IN ) :: urban_map_bd
119 INTEGER, INTENT(IN ) :: urban_map_wd
120 INTEGER, INTENT(IN ) :: urban_map_gbd
121 INTEGER, INTENT(IN ) :: urban_map_fbd
122 INTEGER, INTENT(IN ) :: num_urban_hi
123 REAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
124 REAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
125 REAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
126 REAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
127 REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
128 REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
129 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
130 REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
131 REAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
132 REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lp_urb2d
133 REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lb_urb2d
134 REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: hgt_urb2d
136 ! integer nx,ny,nz ! Number of points in the mesocsale grid
137 real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates
138 REAL, INTENT(IN ):: DT ! Time step
139 ! real zr(ims:ime,jms:jme) ! Solar zenith angle
140 ! real deltar(ims:ime,jms:jme) ! Declination of the sun
141 ! real ah(ims:ime,jms:jme) ! Hour angle
142 ! real rs(ims:ime,jms:jme) ! Solar radiation
143 !------------------------------------------------------------------------
144 ! Output
145 !------------------------------------------------------------------------
146 ! real tsk(ims:ime,jms:jme) ! Average of surface temperatures (roads and roofs)
147 !
148 ! Implicit and explicit components of the source and sink terms at each levels,
149 ! the fluxes can be computed as follow: FX = A*X + B example: t_fluxes = a_t * pt + b_t
150 real a_u(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in X-direction (center)
151 real a_v(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in Y-direction (center)
152 real a_t(ims:ime,kms:kme,jms:jme) ! Implicit component for the temperature
153 real a_e(ims:ime,kms:kme,jms:jme) ! Implicit component for the TKE
154 real b_u(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in X-direction (center)
155 real b_v(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in Y-direction (center)
156 real b_t(ims:ime,kms:kme,jms:jme) ! Explicit component for the temperature
157 real b_e(ims:ime,kms:kme,jms:jme) ! Explicit component for the TKE
158 real b_q(ims:ime,kms:kme,jms:jme) ! Explicit component for the humidity
159 real dlg(ims:ime,kms:kme,jms:jme) ! Height above ground (L_ground in formula (24) of the BLM paper).
160 real dl_u(ims:ime,kms:kme,jms:jme) ! Length scale (lb in formula (22) ofthe BLM paper).
161 ! urban surface and volumes
162 real sf(ims:ime,kms:kme,jms:jme) ! surface of the urban grid cells
163 real vl(ims:ime,kms:kme,jms:jme) ! volume of the urban grid cells
164 ! urban fluxes
165 real rl_up(its:ite,jts:jte) ! upward long wave radiation
166 real rs_abs(its:ite,jts:jte) ! absorbed short wave radiation
167 real emiss(its:ite,jts:jte) ! emissivity averaged for urban surfaces
168 real grdflx_urb(its:ite,jts:jte) ! ground heat flux for urban areas
169 !------------------------------------------------------------------------
170 ! Local
171 !------------------------------------------------------------------------
172 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174 real hi_urb(its:ite,1:nz_um,jts:jte) ! Height histograms of buildings
175 real hi_urb1D(nz_um) ! Height histograms of buildings
176 real hb_u(nz_um) ! Bulding's heights
177 real ss_urb(nz_um) ! Probability that a building has an height equal to z
178 real pb_urb(nz_um) ! Probability that a building has an height greater or equal to z
179 integer nz_urb(nurbmax) ! Number of layer in the urban grid
180 integer nzurban(nurbmax)
182 ! Building parameters
183 real alag_u(nurbmax) ! Ground thermal diffusivity [m^2 s^-1]
184 real alaw_u(nurbmax) ! Wall thermal diffusivity [m^2 s^-1]
185 real alar_u(nurbmax) ! Roof thermal diffusivity [m^2 s^-1]
186 real csg_u(nurbmax) ! Specific heat of the ground material [J m^3 K^-1]
187 real csw_u(nurbmax) ! Specific heat of the wall material [J m^3 K^-1]
188 real csr_u(nurbmax) ! Specific heat of the roof material [J m^3 K^-1]
189 real twini_u(nurbmax) ! Initial temperature inside the building's wall [K]
190 real trini_u(nurbmax) ! Initial temperature inside the building's roof [K]
191 real tgini_u(nurbmax) ! Initial road temperature
192 !
193 ! Building materials
194 !
195 real csg(ng_u) ! Specific heat of the ground material [J m^3 K^-1]
196 real csr(nwr_u) ! Specific heat of the roof material [J m^3 K^-1]
197 real csw(nwr_u) ! Specific heat of the wall material [J m^3 K^-1]
198 real alag(ng_u) ! Ground thermal diffusivity [m^2 s^-1]
199 real alaw(nwr_u) ! Wall thermal diffusivity [m^2 s^-1]
200 real alar(nwr_u) ! Roof thermal diffusivity [m^2 s^-1]
201 !
202 ! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
203 !
204 ! Radiation parameters
205 real albg_u(nurbmax) ! Albedo of the ground
206 real albw_u(nurbmax) ! Albedo of the wall
207 real albr_u(nurbmax) ! Albedo of the roof
208 real emg_u(nurbmax) ! Emissivity of ground
209 real emw_u(nurbmax) ! Emissivity of wall
210 real emr_u(nurbmax) ! Emissivity of roof
212 ! fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
213 ! and the short wave radation.
214 real fww_u(nz_um,nz_um,ndm,nurbmax) ! from wall to wall
215 real fwg_u(nz_um,ndm,nurbmax) ! from wall to ground
216 real fgw_u(nz_um,ndm,nurbmax) ! from ground to wall
217 real fsw_u(nz_um,ndm,nurbmax) ! from sky to wall
218 real fws_u(nz_um,ndm,nurbmax) ! from sky to wall
219 real fsg_u(ndm,nurbmax) ! from sky to ground
221 ! Roughness parameters
222 real z0g_u(nurbmax) ! The ground's roughness length
223 real z0r_u(nurbmax) ! The roof's roughness length
225 ! Roughness parameters
226 real z0(ndm,nz_um) ! Roughness lengths "profiles"
228 ! Street parameters
229 integer nd_u(nurbmax) ! Number of street direction for each urban class
230 real strd_u(ndm,nurbmax) ! Street length (fix to greater value to the horizontal length of the cells)
231 real drst_u(ndm,nurbmax) ! Street direction
232 real ws_u(ndm,nurbmax) ! Street width
233 real bs_u(ndm,nurbmax) ! Building width
234 real h_b(nz_um,nurbmax) ! Bulding's heights
235 real d_b(nz_um,nurbmax) ! Probability that a building has an height h_b
236 real ss_u(nz_um,nurbmax) ! Probability that a building has an height equal to z
237 real pb_u(nz_um,nurbmax) ! Probability that a building has an height greater or equal to z
238 !
239 ! Street parameters
240 !
241 real bs(ndm) ! Building width
242 real ws(ndm) ! Street width
243 real drst(ndm) ! street directions
244 real strd(ndm) ! Street lengths
245 real ss(nz_um) ! Probability to have a building with height h
246 real pb(nz_um) ! Probability to have a building with an height equal
248 ! Grid parameters
250 integer nz_u(nurbmax) ! Number of layer in the urban grid
252 real z_u(nz_um) ! Height of the urban grid levels
254 ! 1D array used for the input and output of the routine "urban"
256 real z1D(kms:kme) ! vertical coordinates
257 real ua1D(kms:kme) ! wind speed in the x directions
258 real va1D(kms:kme) ! wind speed in the y directions
259 real pt1D(kms:kme) ! potential temperature
260 real da1D(kms:kme) ! air density
261 real pr1D(kms:kme) ! air pressure
262 real pt01D(kms:kme) ! reference potential temperature
263 real zr1D ! zenith angle
264 real deltar1D ! declination of the sun
265 real ah1D ! hour angle (it should come from the radiation routine)
266 real rs1D ! solar radiation
267 real rld1D ! downward flux of the longwave radiation
270 real tw1D(2*ndm,nz_um,nwr_u) ! temperature in each layer of the wall
271 real tg1D(ndm,ng_u) ! temperature in each layer of the ground
272 real tr1D(ndm,nz_um,nwr_u) ! temperature in each layer of the roof
273 real sfw1D(2*ndm,nz_um) ! sensible heat flux from walls
274 real sfg1D(ndm) ! sensible heat flux from ground (road)
275 real sfr1D(ndm,nz_um) ! sensible heat flux from roofs
276 real sf1D(kms:kme) ! surface of the urban grid cells
277 real vl1D(kms:kme) ! volume of the urban grid cells
278 real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
279 real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
280 real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
281 real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
282 real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
283 real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
284 real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
285 real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
286 real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
287 real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
288 real time_bep
289 ! arrays used to collapse indexes
290 integer ind_zwd(nz_um,nwr_u,ndm)
291 integer ind_gd(ng_u,ndm)
292 integer ind_zd(nz_um,ndm)
293 !
294 integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
295 integer it, nint
296 integer iii
297 real time_h,tempo
298 logical first
299 character(len=80) :: text
300 data first/.true./
301 save first,time_bep
303 save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
304 albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
305 z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
306 nz_u,z_u
308 !------------------------------------------------------------------------
309 ! Calculation of the momentum, heat and turbulent kinetic fluxes
310 ! produced by builgings
311 !
312 ! Reference:
313 ! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
314 ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
315 ! 261-304
316 !------------------------------------------------------------------------
317 !prepare the arrays to collapse indexes
319 if(urban_map_zrd.lt.nz_um*ndm*nwr_u)then
320 write(*,*)'urban_map_zrd too small, please increase to at least ', nz_um*ndm*nwr_u
321 stop
322 endif
323 iii=0
324 do iz_u=1,nz_um
325 do iw=1,nwr_u
326 do id=1,ndm
327 iii=iii+1
328 ind_zwd(iz_u,iw,id)=iii
329 enddo
330 enddo
331 enddo
333 iii=0
334 do ig=1,ng_u
335 do id=1,ndm
336 iii=iii+1
337 ind_gd(ig,id)=iii
338 enddo
339 enddo
341 iii=0
342 do iz_u=1,nz_um
343 do id=1,ndm
344 iii=iii+1
345 ind_zd(iz_u,id)=iii
346 enddo
347 enddo
349 if (num_urban_hi.ge.nz_um)then
350 write(*,*)'nz_um too small, please increase to at least ', num_urban_hi+1
351 stop
352 endif
354 do ix=its,ite
355 do iy=jts,jte
356 do iz_u=1,nz_um
357 hi_urb(ix,iz_u,iy)=0.
358 enddo
359 enddo
360 enddo
362 do ix=its,ite
363 do iy=jts,jte
364 z(ix,kts,iy)=0.
365 do iz=kts+1,kte+1
366 z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
367 enddo !iz
368 do iz_u=1,num_urban_hi
369 hi_urb(ix,iz_u,iy)= hi_urb2d(ix,iz_u,iy)
370 enddo !iz_u
371 enddo
372 enddo
375 if (first) then ! True only on first call
377 call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
378 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
379 emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
381 ! Initialisation of the urban parameters and calculation of the view factors
383 call icBEP(nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u)
385 first=.false.
387 endif ! first
389 do ix=its,ite
390 do iy=jts,jte
391 if (FRC_URB2D(ix,iy).gt.0.) then ! Calling BEP only for existing urban classes.
392 iurb=UTYPE_URB2D(ix,iy)
394 hi_urb1D=0.
395 do iz_u=1,nz_um
396 hi_urb1D(iz_u)=hi_urb(ix,iz_u,iy)
397 enddo
399 call icBEPHI_XY(hb_u,hi_urb1D,ss_urb,pb_urb, &
400 nz_urb(iurb),z_u)
402 call param(iurb,nz_u(iurb),nz_urb(iurb),nzurban(iurb), &
403 nd_u(iurb),csg_u,csg,alag_u,alag,csr_u,csr, &
404 alar_u,alar,csw_u,csw,alaw_u,alaw, &
405 ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
406 strd_u,strd,drst_u,drst,ss_u,ss_urb,ss,pb_u, &
407 pb_urb,pb,lp_urb2d(ix,iy), &
408 lb_urb2d(ix,iy),hgt_urb2d(ix,iy),FRC_URB2D(ix,iy))
409 !
410 !We compute the view factors in the icBEP_XY routine
411 !
413 call icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u,fws_u,fsg_u, &
414 nd_u(iurb),strd,ws,nzurban(iurb),z_u)
416 do iz= kts,kte
417 ua1D(iz)=u_phy(ix,iz,iy)
418 va1D(iz)=v_phy(ix,iz,iy)
419 pt1D(iz)=th_phy(ix,iz,iy)
420 da1D(iz)=rho(ix,iz,iy)
421 pr1D(iz)=p_phy(ix,iz,iy)
422 ! pt01D(iz)=th_phy(ix,iz,iy)
423 pt01D(iz)=300.
424 z1D(iz)=z(ix,iz,iy)
425 a_u1D(iz)=0.
426 a_v1D(iz)=0.
427 a_t1D(iz)=0.
428 a_e1D(iz)=0.
429 b_u1D(iz)=0.
430 b_v1D(iz)=0.
431 b_t1D(iz)=0.
432 b_e1D(iz)=0.
433 enddo
434 z1D(kte+1)=z(ix,kte+1,iy)
436 do id=1,ndm
437 do iz_u=1,nz_um
438 do iw=1,nwr_u
439 ! tw1D(2*id-1,iz_u,iw)=tw1_u(ix,iy,ind_zwd(iz_u,iw,id))
440 ! tw1D(2*id,iz_u,iw)=tw2_u(ix,iy,ind_zwd(iz_u,iw,id))
441 if(ind_zwd(iz_u,iw,id).gt.urban_map_zwd)write(*,*)'ind_zwd too big w',ind_zwd(iz_u,iw,id)
442 tw1D(2*id-1,iz_u,iw)=tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
443 tw1D(2*id,iz_u,iw)=tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
444 enddo
445 enddo
446 enddo
448 do id=1,ndm
449 do ig=1,ng_u
450 ! tg1D(id,ig)=tg_u(ix,iy,ind_gd(ig,id))
451 tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
452 enddo
453 do iz_u=1,nz_um
454 do ir=1,nwr_u
455 ! tr1D(id,iz_u,ir)=tr_u(ix,iy,ind_zwd(iz_u,ir,id))
456 if(ind_zwd(iz_u,ir,id).gt.urban_map_zwd)write(*,*)'ind_zwd too big r',ind_zwd(iz_u,ir,id)
457 tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)
458 enddo
459 enddo
460 enddo
462 do id=1,ndm
463 do iz=1,nz_um
464 ! sfw1D(2*id-1,iz)=sfw1(ix,iy,ind_zd(iz,id))
465 ! sfw1D(2*id,iz)=sfw2(ix,iy,ind_zd(iz,id))
466 sfw1D(2*id-1,iz)=sfw1_urb3d(ix,ind_zd(iz,id),iy)
467 sfw1D(2*id,iz)=sfw2_urb3d(ix,ind_zd(iz,id),iy)
468 enddo
469 enddo
471 do id=1,ndm
472 ! sfg1D(id)=sfg(ix,iy,id)
473 sfg1D(id)=sfg_urb3d(ix,id,iy)
474 enddo
476 do id=1,ndm
477 do iz=1,nz_um
478 ! sfr1D(id,iz)=sfr(ix,iy,ind_zd(iz,id))
479 sfr1D(id,iz)=sfr_urb3d(ix,ind_zd(iz,id),iy)
480 enddo
481 enddo
484 rs1D=swdown(ix,iy)
485 rld1D=glw(ix,iy)
486 time_h=(itimestep*dt)/3600.+gmt
488 zr1D=acos(COSZ_URB2D(ix,iy))
489 deltar1D=DECLIN_URB
490 ah1D=OMG_URB2D(ix,iy)
491 ! call angle(xlong(ix,iy),xlat(ix,iy),julday,time_h,zr1D,deltar1D,ah1D)
492 ! write(*,*) 'entro en BEP1D'
493 call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
494 zr1D,deltar1D,ah1D,rs1D,rld1D, &
495 alag,alaw,alar,csg,csw,csr, &
496 albg_u(iurb),albw_u(iurb),albr_u(iurb), &
497 emg_u(iurb),emw_u(iurb),emr_u(iurb), &
498 fww_u,fwg_u,fgw_u,fsw_u, &
499 fws_u,fsg_u,z0, &
500 nd_u(iurb),strd,drst,ws,bs,ss,pb, &
501 nzurban(iurb),z_u, &
502 tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D, &
503 a_u1D,a_v1D,a_t1D,a_e1D, &
504 b_u1D,b_v1D,b_t1D,b_e1D, &
505 dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy), &
506 rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy))
507 ! write(*,*) 'salgo de BEP1D'
508 do id=1,ndm
509 do iz=1,nz_um
510 sfw1_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id-1,iz)
511 sfw2_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id,iz)
512 enddo
513 enddo
515 do id=1,ndm
516 sfg_urb3d(ix,id,iy)=sfg1D(id)
517 enddo
519 do id=1,ndm
520 do iz=1,nz_um
521 sfr_urb3d(ix,ind_zd(iz,id),iy)=sfr1D(id,iz)
522 enddo
523 enddo
524 !
525 do id=1,ndm
526 do iz_u=1,nz_um
527 do iw=1,nwr_u
528 tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw)
529 tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw)
530 enddo
531 enddo
532 enddo
534 do id=1,ndm
535 do ig=1,ng_u
536 tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
537 enddo
538 do iz_u=1,nz_um
539 do ir=1,nwr_u
540 trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
541 enddo
542 enddo
543 enddo
545 sf(ix,kts:kte,iy)=0.
546 vl(ix,kts:kte,iy)=0.
547 a_u(ix,kts:kte,iy)=0.
548 a_v(ix,kts:kte,iy)=0.
549 a_t(ix,kts:kte,iy)=0.
550 a_e(ix,kts:kte,iy)=0.
551 b_u(ix,kts:kte,iy)=0.
552 b_v(ix,kts:kte,iy)=0.
553 b_t(ix,kts:kte,iy)=0.
554 b_e(ix,kts:kte,iy)=0.
555 b_q(ix,kts:kte,iy)=0.
556 dlg(ix,kts:kte,iy)=0.
557 dl_u(ix,kts:kte,iy)=0.
559 do iz= kts,kte
560 sf(ix,iz,iy)=sf1D(iz)
561 vl(ix,iz,iy)=vl1D(iz)
562 a_u(ix,iz,iy)=a_u1D(iz)
563 a_v(ix,iz,iy)=a_v1D(iz)
564 a_t(ix,iz,iy)=a_t1D(iz)
565 a_e(ix,iz,iy)=a_e1D(iz)
566 b_u(ix,iz,iy)=b_u1D(iz)
567 b_v(ix,iz,iy)=b_v1D(iz)
568 b_t(ix,iz,iy)=b_t1D(iz)
569 b_e(ix,iz,iy)=b_e1D(iz)
570 dlg(ix,iz,iy)=dlg1D(iz)
571 dl_u(ix,iz,iy)=dl_u1D(iz)
572 enddo
573 sf(ix,kte+1,iy)=sf1D(kte+1)
574 !
575 endif ! FRC_URB2D
577 enddo ! iy
578 enddo ! ix
581 time_bep=time_bep+dt
584 return
585 end subroutine BEP
587 ! ===6=8===============================================================72
589 subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, &
590 zr,deltar,ah,rs,rld, &
591 alag,alaw,alar,csg,csw,csr, &
592 albg,albw,albr,emg,emw,emr, &
593 fww,fwg,fgw,fsw,fws,fsg,z0, &
594 ndu,strd,drst,ws,bs,ss,pb, &
595 nzu,z_u, &
596 tw,tg,tr,sfw,sfg,sfr, &
597 a_u,a_v,a_t,a_e, &
598 b_u,b_v,b_t,b_e, &
599 dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb)
601 ! ----------------------------------------------------------------------
602 ! This routine computes the effects of buildings on momentum, heat and
603 ! TKE (turbulent kinetic energy) sources or sinks and on the mixing length.
604 ! It provides momentum, heat and TKE sources or sinks at different levels of a
605 ! mesoscale grid defined by the altitude of its cell interfaces "z" and
606 ! its number of levels "nz".
607 ! The meteorological input parameters (wind, temperature, solar radiation)
608 ! are specified on the "mesoscale grid".
609 ! The inputs concerning the building and street charateristics are defined
610 ! on a "urban grid". The "urban grid" is defined with its number of levels
611 ! "nz_u" and its space step "dz_u".
612 ! The input parameters are interpolated on the "urban grid". The sources or sinks
613 ! are calculated on the "urban grid". Finally the sources or sinks are
614 ! interpolated on the "mesoscale grid".
617 ! Mesoscale grid Urban grid Mesoscale grid
618 !
619 ! z(4) --- ---
620 ! | |
621 ! | |
622 ! | Interpolation Interpolation |
623 ! | Sources or sinks calculation |
624 ! z(3) --- ---
625 ! | ua ua_u --- uv_a a_u |
626 ! | va va_u | uv_b b_u |
627 ! | pt pt_u --- uh_b a_v |
628 ! z(2) --- | etc... etc...---
629 ! | z_u(1) --- |
630 ! | | |
631 ! z(1) ------------------------------------------------------------
633 !
634 ! Reference:
635 ! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
636 ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
637 ! 261-304
639 ! ----------------------------------------------------------------------
641 implicit none
645 ! ----------------------------------------------------------------------
646 ! INPUT:
647 ! ----------------------------------------------------------------------
649 ! Data relative to the "mesoscale grid"
651 ! integer nz ! Number of vertical levels
652 integer kms,kme,kts,kte
653 real z(kms:kme) ! Altitude above the ground of the cell interfaces.
654 real ua(kms:kme) ! Wind speed in the x direction
655 real va(kms:kme) ! Wind speed in the y direction
656 real pt(kms:kme) ! Potential temperature
657 real da(kms:kme) ! Air density
658 real pr(kms:kme) ! Air pressure
659 real pt0(kms:kme) ! Reference potential temperature (could be equal to "pt")
660 real dt ! Time step
661 real zr ! Zenith angle
662 real deltar ! Declination of the sun
663 real ah ! Hour angle
664 real rs ! Solar radiation
665 real rld ! Downward flux of the longwave radiation
667 ! Data relative to the "urban grid"
669 integer iurb ! Current urban class
671 ! Building parameters
672 real alag(ng_u) ! Ground thermal diffusivity [m^2 s^-1]
673 real alaw(nwr_u) ! Wall thermal diffusivity [m^2 s^-1]
674 real alar(nwr_u) ! Roof thermal diffusivity [m^2 s^-1]
675 real csg(ng_u) ! Specific heat of the ground material [J m^3 K^-1]
676 real csw(nwr_u) ! Specific heat of the wall material [J m^3 K^-1]
677 real csr(nwr_u) ! Specific heat of the roof material [J m^3 K^-1]
679 ! Radiation parameters
680 real albg ! Albedo of the ground
681 real albw ! Albedo of the wall
682 real albr ! Albedo of the roof
683 real emg ! Emissivity of ground
684 real emw ! Emissivity of wall
685 real emr ! Emissivity of roof
687 ! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and
688 ! short wave radation.
689 ! The calculation of these factor is explained in the Appendix A of the BLM paper
690 real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
691 real fwg(nz_um,ndm,nurbm) ! from wall to ground
692 real fgw(nz_um,ndm,nurbm) ! from ground to wall
693 real fsw(nz_um,ndm,nurbm) ! from sky to wall
694 real fws(nz_um,ndm,nurbm) ! from wall to sky
695 real fsg(ndm,nurbm) ! from sky to ground
697 ! Roughness parameters
698 real z0(ndm,nz_um) ! Roughness lengths "profiles"
700 ! Street parameters
701 integer ndu ! Number of street direction for each urban class
702 real strd(ndm) ! Street length (set to a greater value then the horizontal length of the cells)
703 real drst(ndm) ! Street direction
704 real ws(ndm) ! Street width
705 real bs(ndm) ! Building width
706 real ss(nz_um) ! The probability that a building has an height equal to "z"
707 real pb(nz_um) ! The probability that a building has an height greater or equal to "z"
709 ! Grid parameters
710 integer nzu ! Number of layer in the urban grid
711 real z_u(nz_um) ! Height of the urban grid levels
714 ! ----------------------------------------------------------------------
715 ! INPUT-OUTPUT
716 ! ----------------------------------------------------------------------
718 ! Data relative to the "urban grid" which should be stored from the current time step to the next one
720 real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
721 real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
722 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
723 real sfw(2*ndm,nz_um) ! Sensible heat flux from walls
724 real sfg(ndm) ! Sensible heat flux from ground (road)
725 real sfr(ndm,nz_um) ! Sensible heat flux from roofs
726 real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) towards the interior
727 real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof towards the interior
728 real gfw(2*ndm,nz_um) ! Heat flux transfered from the surface of the walls towards the interior
729 ! ----------------------------------------------------------------------
730 ! OUTPUT:
731 ! ----------------------------------------------------------------------
733 ! Data relative to the "mesoscale grid"
735 real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
736 real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
738 ! Implicit and explicit components of the source and sink terms at each levels,
739 ! the fluxes can be computed as follow: FX = A*X + B example: Heat fluxes = a_t * pt + b_t
740 real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
741 real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
742 real a_t(kms:kme) ! Implicit component of the heat sources or sinks
743 real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
744 real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
745 real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
746 real b_t(kms:kme) ! Explicit component of the heat sources or sinks
747 real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
748 real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
749 real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
751 ! ----------------------------------------------------------------------
752 ! LOCAL:
753 ! ----------------------------------------------------------------------
755 real dz(kms:kme) ! vertical space steps of the "mesoscale grid"
757 ! Data interpolated from the "mesoscale grid" to the "urban grid"
759 real ua_u(nz_um) ! Wind speed in the x direction
760 real va_u(nz_um) ! Wind speed in the y direction
761 real pt_u(nz_um) ! Potential temperature
762 real da_u(nz_um) ! Air density
763 real pt0_u(nz_um) ! Reference potential temperature
764 real pr_u(nz_um) ! Air pressure
766 ! Solar radiation at each level of the "urban grid"
768 real rsg(ndm) ! Short wave radiation from the ground
769 real rsw(2*ndm,nz_um) ! Short wave radiation from the walls
770 real rlg(ndm) ! Long wave radiation from the ground
771 real rlw(2*ndm,nz_um) ! Long wave radiation from the walls
773 ! Potential temperature of the surfaces at each level of the "urban grid"
775 real ptg(ndm) ! Ground potential temperatures
776 real ptr(ndm,nz_um) ! Roof potential temperatures
777 real ptw(2*ndm,nz_um) ! Walls potential temperatures
780 ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
781 ! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
782 ! The fluxes can be computed as follow: Fluxes of X = A*X + B
783 ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
785 real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
786 real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
787 real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
788 real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
789 real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
790 real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
791 real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
792 real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
793 real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
794 real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
795 real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
797 !
798 real rs_abs ! solar radiation absorbed by urban surfaces
799 real rl_up ! longwave radiation emitted by urban surface to the atmosphere
800 real emiss ! mean emissivity of the urban surface
801 real grdflx_urb ! ground heat flux
802 integer iz,id
803 integer iw,ix,iy
805 ! ----------------------------------------------------------------------
806 ! END VARIABLES DEFINITIONS
807 ! ----------------------------------------------------------------------
809 ! Fix some usefull parameters for the computation of the sources or sinks
811 do iz=kts,kte
812 dz(iz)=z(iz+1)-z(iz)
813 end do
815 ! Interpolation on the "urban grid"
816 call interpol(kms,kme,kts,kte,nzu,z,z_u,ua,ua_u)
817 call interpol(kms,kme,kts,kte,nzu,z,z_u,va,va_u)
818 call interpol(kms,kme,kts,kte,nzu,z,z_u,pt,pt_u)
819 call interpol(kms,kme,kts,kte,nzu,z,z_u,pt0,pt0_u)
820 call interpol(kms,kme,kts,kte,nzu,z,z_u,pr,pr_u)
821 call interpol(kms,kme,kts,kte,nzu,z,z_u,da,da_u)
824 ! Compute the modification of the radiation due to the buildings
826 call modif_rad(iurb,ndu,nzu,z_u,ws, &
827 drst,strd,ss,pb, &
828 tw,tg,albg,albw,emw,emg, &
829 fww,fwg,fgw,fsw,fsg, &
830 zr,deltar,ah, &
831 rs,rld,rsw,rsg,rlw,rlg)
833 ! calculation of the urban albedo and the upward long wave radiation
835 call upward_rad(ndu,nzu,ws,bs, &
836 sigma,pb,ss, &
837 tg,emg,albg,rlg,rsg,sfg, &
838 tw,emw,albw,rlw,rsw,sfw, &
839 tr,emr,albr,rld,rs,sfr, &
840 rs_abs,rl_up,emiss,grdflx_urb)
842 ! Compute the surface temperatures
845 call surf_temp(nzu,ndu,pr_u,dt,ss, &
846 rs,rld,rsg,rlg,rsw,rlw, &
847 tg,alag,csg,emg,albg,ptg,sfg,gfg, &
848 tr,alar,csr,emr,albr,ptr,sfr,gfr, &
849 tw,alaw,csw,emw,albw,ptw,sfw,gfw)
852 ! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
854 call buildings(ndu,nzu,z0,ua_u,va_u, &
855 pt_u,pt0_u,ptg,ptr,da_u,ptw,drst, &
856 uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
857 uhb_u,vhb_u,thb_u,ehb_u,ss,dt)
860 ! Calculation of the sensible heat fluxes for the ground, the wall and roof
861 ! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
862 ! where A and B are the implicit and explicit components of the heat sources or sinks.
863 !
864 !
866 do id=1,ndu
867 sfg(id)=-da_u(1)*cp_u*thb_u(id,1)
868 do iz=2,nzu
869 sfr(id,iz)=-da_u(iz)*cp_u*thb_u(id,iz)
870 enddo
872 do iz=1,nzu
873 sfw(2*id-1,iz)=-da_u(iz)*cp_u*(tvb_u(2*id-1,iz)+ &
874 tva_u(2*id-1,iz)*pt_u(iz))
875 sfw(2*id,iz)=-da_u(iz)*cp_u*(tvb_u(2*id,iz)+ &
876 tva_u(2*id,iz)*pt_u(iz))
877 enddo
878 enddo
880 ! calculation of the urban albedo and the upward long wave radiation
882 !! call upward_rad(ndu,nzu,ws,bs, &
883 !! sigma,pb,ss, &
884 !! tg,emg,albg,rlg,rsg,sfg, &
885 !! tw,emw,albw,rlw,rsw,sfw, &
886 !! tr,emr,albr,rld,rs,sfr, &
887 !! rs_abs,rl_up,emiss,grdflx_urb)
889 ! Interpolation on the "mesoscale grid"
891 call urban_meso(ndu,kms,kme,kts,kte,nzu,z,dz,z_u,pb,ss,bs,ws,sf, &
892 vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
893 uhb_u,vhb_u,thb_u,ehb_u, &
894 a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)
897 ! computation of the mean road temperature tsk (this value could be used
898 ! to replace the surface temperature in the radiation routines, if needed).
900 ! Calculation of the length scale taking into account the buildings effects
902 call interp_length(ndu,kms,kme,kts,kte,nzu,z_u,z,ss,ws,bs,dlg,dl_u)
904 return
905 end subroutine BEP1D
907 ! ===6=8===============================================================72
908 ! ===6=8===============================================================72
910 subroutine param(iurb,nzu,nzurb,nzurban,ndu, &
911 csg_u,csg,alag_u,alag,csr_u,csr, &
912 alar_u,alar,csw_u,csw,alaw_u,alaw, &
913 ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
914 strd_u,strd,drst_u,drst,ss_u,ss_urb,ss,pb_u, &
915 pb_urb,pb,lp_urb,lb_urb,hgt_urb,frc_urb)
917 ! ----------------------------------------------------------------------
918 ! This routine prepare some usefull parameters
919 ! ----------------------------------------------------------------------
921 implicit none
924 ! ----------------------------------------------------------------------
925 ! INPUT:
926 ! ----------------------------------------------------------------------
927 integer iurb ! Current urban class
928 integer nzu ! Number of vertical urban levels in the current class
929 integer nzurb ! Number of vertical urban levels in the current class
930 integer ndu ! Number of street direction for the current urban class
931 real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
932 real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
933 real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
934 real bs_u(ndm,nurbm) ! Building width
935 real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
936 real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
937 real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
938 real drst_u(ndm,nurbm) ! Street direction
939 real strd_u(ndm,nurbm) ! Street length
940 real ws_u(ndm,nurbm) ! Street width
941 real z0g_u(nurbm) ! The ground's roughness length
942 real z0r_u(nurbm) ! The roof's roughness length
943 real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
944 real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
945 real ss_urb(nz_um) ! The probability that a building has an height equal to "z"
946 real pb_urb(nz_um) ! The probability that a building has an height greater or equal to "z"
947 real lp_urb ! Building plan area density
948 real lb_urb ! Building surface area to plan area ratio
949 real hgt_urb ! Average building height weighted by building plan area [m]
950 real frc_urb ! Urban fraction
952 ! ----------------------------------------------------------------------
953 ! OUTPUT:
954 ! ----------------------------------------------------------------------
955 real alag(ng_u) ! Ground thermal diffusivity at each ground levels
956 real alar(nwr_u) ! Roof thermal diffusivity at each roof levels
957 real alaw(nwr_u) ! Wall thermal diffusivity at each wall levels
958 real csg(ng_u) ! Specific heat of the ground material at each ground levels
959 real csr(nwr_u) ! Specific heat of the roof material at each roof levels
960 real csw(nwr_u) ! Specific heat of the wall material at each wall levels
961 real bs(ndm) ! Building width for the current urban class
962 real drst(ndm) ! street directions for the current urban class
963 real strd(ndm) ! Street lengths for the current urban class
964 real ws(ndm) ! Street widths of the current urban class
965 real z0(ndm,nz_um) ! Roughness lengths "profiles"
966 real ss(nz_um) ! Probability to have a building with height h
967 real pb(nz_um) ! Probability to have a building with an height equal
968 integer nzurban
970 ! ----------------------------------------------------------------------
971 ! LOCAL:
972 ! ----------------------------------------------------------------------
973 integer id,ig,ir,iw,iz,ihu
975 ! ----------------------------------------------------------------------
976 ! END VARIABLES DEFINITIONS
977 ! ----------------------------------------------------------------------
978 !
979 !Initialize
980 !
981 ss=0.
982 pb=0.
983 csg=0.
984 alag=0.
985 csr=0.
986 alar=0.
987 csw=0.
988 alaw=0.
989 z0=0.
990 ws=0.
991 bs=0.
992 strd=0.
993 drst=0.
994 nzurban=0
996 ihu=0
998 do iz=1,nz_um
999 if (ss_urb(iz)/=0.) then
1000 ihu=1
1001 exit
1002 else
1003 continue
1004 endif
1005 enddo
1007 if (ihu==1) then
1008 do iz=1,nzurb+1
1009 ss(iz)=ss_urb(iz)
1010 pb(iz)=pb_urb(iz)
1011 enddo
1012 nzurban=nzurb
1013 else
1014 do iz=1,nzu+1
1015 ss(iz)=ss_u(iz,iurb)
1016 pb(iz)=pb_u(iz,iurb)
1017 end do
1018 nzurban=nzu
1019 endif
1021 do id=1,ndu
1022 z0(id,1)=z0g_u(iurb)
1023 do iz=2,nzurban+1
1024 z0(id,iz)=z0r_u(iurb)
1025 enddo
1026 enddo
1028 do ig=1,ng_u
1029 csg(ig)=csg_u(iurb)
1030 alag(ig)=alag_u(iurb)
1031 enddo
1033 do ir=1,nwr_u
1034 csr(ir)=csr_u(iurb)
1035 alar(ir)=alar_u(iurb)
1036 enddo
1038 do iw=1,nwr_u
1039 csw(iw)=csw_u(iurb)
1040 alaw(iw)=alaw_u(iurb)
1041 enddo
1043 do id=1,ndu
1044 strd(id)=strd_u(id,iurb)
1045 drst(id)=drst_u(id,iurb)
1046 enddo
1048 do id=1,ndu
1049 if ((hgt_urb<=0.).OR.(lp_urb<=0.).OR.(lb_urb<=0.)) then
1050 ws(id)=ws_u(id,iurb)
1051 bs(id)=bs_u(id,iurb)
1052 else if ((lp_urb/frc_urb<1.).and.(lp_urb<lb_urb)) then
1053 bs(id)=2.*hgt_urb*lp_urb/(lb_urb-lp_urb)
1054 ws(id)=2.*hgt_urb*lp_urb*((frc_urb/lp_urb)-1.)/(lb_urb-lp_urb)
1055 else
1056 ws(id)=ws_u(id,iurb)
1057 bs(id)=bs_u(id,iurb)
1058 endif
1059 enddo
1060 do id=1,ndu
1061 if ((bs(id)<=1.).OR.(bs(id)>=150.)) then
1062 ! write(*,*) 'WARNING, WIDTH OF THE BUILDING WRONG',id,bs(id)
1063 ! write(*,*) 'WIDTH OF THE STREET',id,ws(id)
1064 bs(id)=bs_u(id,iurb)
1065 ws(id)=ws_u(id,iurb)
1066 endif
1067 if ((ws(id)<=1.).OR.(ws(id)>=150.)) then
1068 ! write(*,*) 'WARNING, WIDTH OF THE STREET WRONG',id,ws(id)
1069 ! write(*,*) 'WIDTH OF THE BUILDING',id,bs(id)
1070 bs(id)=bs_u(id,iurb)
1071 ws(id)=ws_u(id,iurb)
1072 endif
1073 enddo
1074 return
1075 end subroutine param
1077 ! ===6=8===============================================================72
1078 ! ===6=8===============================================================72
1080 subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
1082 ! ----------------------------------------------------------------------
1083 ! This routine interpolate para
1084 ! meters from the "mesoscale grid" to
1085 ! the "urban grid".
1086 ! See p300 Appendix B.1 of the BLM paper.
1087 ! ----------------------------------------------------------------------
1089 implicit none
1091 ! ----------------------------------------------------------------------
1092 ! INPUT:
1093 ! ----------------------------------------------------------------------
1094 ! Data relative to the "mesoscale grid"
1095 integer kts,kte,kms,kme
1096 real z(kms:kme) ! Altitude of the cell interface
1097 real c(kms:kme) ! Parameter which has to be interpolated
1098 ! Data relative to the "urban grid"
1099 integer nz_u ! Number of levels
1100 !! real z_u(nz_u+1) ! Altitude of the cell interface
1101 real z_u(nz_um) ! Altitude of the cell interface
1102 ! ----------------------------------------------------------------------
1103 ! OUTPUT:
1104 ! ----------------------------------------------------------------------
1105 !! real c_u(nz_u) ! Interpolated paramters in the "urban grid"
1106 real c_u(nz_um) ! Interpolated paramters in the "urban grid"
1109 ! LOCAL:
1110 ! ----------------------------------------------------------------------
1111 integer iz_u,iz
1112 real ctot,dz
1114 ! ----------------------------------------------------------------------
1115 ! END VARIABLES DEFINITIONS
1116 ! ----------------------------------------------------------------------
1118 do iz_u=1,nz_u
1119 ctot=0.
1120 do iz=kts,kte
1121 dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
1122 ctot=ctot+c(iz)*dz
1123 enddo
1124 c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
1125 enddo
1127 return
1128 end subroutine interpol
1130 ! ===6=8===============================================================72
1131 ! ===6=8===============================================================72
1133 subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, &
1134 tw,tg,albg,albw,emw,emg, &
1135 fww,fwg,fgw,fsw,fsg, &
1136 zr,deltar,ah, &
1137 rs,rl,rsw,rsg,rlw,rlg)
1139 ! ----------------------------------------------------------------------
1140 ! This routine computes the modification of the short wave and
1141 ! long wave radiation due to the buildings.
1142 ! ----------------------------------------------------------------------
1144 implicit none
1147 ! ----------------------------------------------------------------------
1148 ! INPUT:
1149 ! ----------------------------------------------------------------------
1150 integer iurb ! current urban class
1151 integer nd ! Number of street direction for the current urban class
1152 integer nz_u ! Number of layer in the urban grid
1153 real z(nz_um) ! Height of the urban grid levels
1154 real ws(ndm) ! Street widths of the current urban class
1155 real drst(ndm) ! street directions for the current urban class
1156 real strd(ndm) ! Street lengths for the current urban class
1157 real ss(nz_um) ! probability to have a building with height h
1158 real pb(nz_um) ! probability to have a building with an height equal
1159 real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
1160 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1161 real albg ! Albedo of the ground for the current urban class
1162 real albw ! Albedo of the wall for the current urban class
1163 real emg ! Emissivity of ground for the current urban class
1164 real emw ! Emissivity of wall for the current urban class
1165 real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
1166 real fsg(ndm,nurbm) ! View factors from sky to ground
1167 real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
1168 real fws(nz_um,ndm,nurbm) ! View factors from wall to sky
1169 real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
1170 real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
1171 real ah ! Hour angle (it should come from the radiation routine)
1172 real zr ! zenith angle
1173 real deltar ! Declination of the sun
1174 real rs ! solar radiation
1175 real rl ! downward flux of the longwave radiation
1176 ! ----------------------------------------------------------------------
1177 ! OUTPUT:
1178 ! ----------------------------------------------------------------------
1179 real rlg(ndm) ! Long wave radiation at the ground
1180 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1181 real rsg(ndm) ! Short wave radiation at the ground
1182 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1184 ! ----------------------------------------------------------------------
1185 ! LOCAL:
1186 ! ----------------------------------------------------------------------
1188 integer id,iz
1190 ! Calculation of the shadow effects
1192 call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
1193 rs,rsw,rsg)
1195 ! Calculation of the reflection effects
1196 do id=1,nd
1197 call long_rad(iurb,nz_u,id,emw,emg, &
1198 fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
1200 call short_rad(iurb,nz_u,id,albw,albg,fwg,fww,fgw,rsg,rsw,pb)
1202 enddo
1204 return
1205 end subroutine modif_rad
1208 ! ===6=8===============================================================72
1209 ! ===6=8===============================================================72
1211 subroutine surf_temp(nz_u,nd,pr,dt,ss,rs,rl,rsg,rlg,rsw,rlw, &
1212 tg,alag,csg,emg,albg,ptg,sfg,gfg, &
1213 tr,alar,csr,emr,albr,ptr,sfr,gfr, &
1214 tw,alaw,csw,emw,albw,ptw,sfw,gfw)
1217 ! ----------------------------------------------------------------------
1218 ! Computation of the surface temperatures for walls, ground and roofs
1219 ! ----------------------------------------------------------------------
1221 implicit none
1225 ! ----------------------------------------------------------------------
1226 ! INPUT:
1227 ! ----------------------------------------------------------------------
1228 integer nz_u ! Number of vertical layers defined in the urban grid
1229 integer nd ! Number of street direction for the current urban class
1230 real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
1231 real alar(nwr_u) ! Roof thermal diffusivity for the current urban class [m^2 s^-1]
1232 real alaw(nwr_u) ! Wall thermal diffusivity for the current urban class [m^2 s^-1]
1233 real albg ! Albedo of the ground for the current urban class
1234 real albr ! Albedo of the roof for the current urban class
1235 real albw ! Albedo of the wall for the current urban class
1236 real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
1237 real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
1238 real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
1239 real dt ! Time step
1240 real emg ! Emissivity of ground for the current urban class
1241 real emr ! Emissivity of roof for the current urban class
1242 real emw ! Emissivity of wall for the current urban class
1243 real pr(nz_um) ! Air pressure
1244 real rs ! Solar radiation
1245 real rl ! Downward flux of the longwave radiation
1246 real rlg(ndm) ! Long wave radiation at the ground
1247 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1248 real rsg(ndm) ! Short wave radiation at the ground
1249 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1250 real sfg(ndm) ! Sensible heat flux from ground (road)
1251 real sfr(ndm,nz_um) ! Sensible heat flux from roofs
1252 real sfw(2*ndm,nz_um) ! Sensible heat flux from walls
1253 real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) toward the interior
1254 real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof toward the interior
1255 real gfw(2*ndm,nz_um) ! Heat flux transfered from the surface of the walls toward the interior
1256 real ss(nz_um) ! Probability to have a building with height h
1257 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1258 real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
1259 real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
1262 ! ----------------------------------------------------------------------
1263 ! OUTPUT:
1264 ! ----------------------------------------------------------------------
1265 real ptg(ndm) ! Ground potential temperatures
1266 real ptr(ndm,nz_um) ! Roof potential temperatures
1267 real ptw(2*ndm,nz_um) ! Walls potential temperatures
1269 ! ----------------------------------------------------------------------
1270 ! LOCAL:
1271 ! ----------------------------------------------------------------------
1272 integer id,ig,ir,iw,iz
1274 real rtg(ndm) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
1275 real rtr(ndm,nz_um) ! Total radiation at roof surface (solar+incoming long+outgoing long)
1276 real rtw(2*ndm,nz_um) ! Radiation at walls surface (solar+incoming long+outgoing long)
1277 real tg_tmp(ng_u)
1278 real tr_tmp(nwr_u)
1279 real tw_tmp(nwr_u)
1281 real dzg_u(ng_u) ! Layer sizes in the ground
1283 real dzr_u(nwr_u) ! Layers sizes in the roof
1285 real dzw_u(nwr_u) ! Layer sizes in the wall
1288 data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
1289 data dzr_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
1290 data dzw_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
1291 ! ----------------------------------------------------------------------
1292 ! END VARIABLES DEFINITIONS
1293 ! ----------------------------------------------------------------------
1297 do id=1,nd
1299 ! Calculation for the ground surfaces
1300 do ig=1,ng_u
1301 tg_tmp(ig)=tg(id,ig)
1302 end do
1303 !
1304 call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg, &
1305 rsg(id),rlg(id),pr(1), &
1306 dt,emg,albg, &
1307 rtg(id),sfg(id),gfg(id))
1308 do ig=1,ng_u
1309 tg(id,ig)=tg_tmp(ig)
1310 end do
1312 ! Calculation for the roofs surfaces
1314 do iz=2,nz_u
1316 if(ss(iz).gt.0.)then
1317 do ir=1,nwr_u
1318 tr_tmp(ir)=tr(id,iz,ir)
1319 end do
1321 call soil_temp(nwr_u,dzr_u,tr_tmp,ptr(id,iz), &
1322 alar,csr,rs,rl,pr(iz),dt,emr,albr, &
1323 rtr(id,iz),sfr(id,iz),gfr(id,iz))
1324 do ir=1,nwr_u
1325 tr(id,iz,ir)=tr_tmp(ir)
1326 end do
1328 end if
1330 end do !iz
1332 ! Calculation for the walls surfaces
1334 do iz=1,nz_u
1336 do iw=1,nwr_u
1337 tw_tmp(iw)=tw(2*id-1,iz,iw)
1338 end do
1339 call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id-1,iz),alaw, &
1340 csw, &
1341 rsw(2*id-1,iz),rlw(2*id-1,iz), &
1342 pr(iz),dt,emw, &
1343 albw,rtw(2*id-1,iz),sfw(2*id-1,iz),gfw(2*id-1,iz))
1345 do iw=1,nwr_u
1346 tw(2*id-1,iz,iw)=tw_tmp(iw)
1347 end do
1349 do iw=1,nwr_u
1350 tw_tmp(iw)=tw(2*id,iz,iw)
1351 end do
1353 call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id,iz),alaw, &
1354 csw, &
1355 rsw(2*id,iz),rlw(2*id,iz), &
1356 pr(iz),dt,emw, &
1357 albw,rtw(2*id,iz),sfw(2*id,iz),gfw(2*id,iz))
1358 do iw=1,nwr_u
1359 tw(2*id,iz,iw)=tw_tmp(iw)
1360 end do
1362 end do !iz
1364 end do !id
1366 return
1367 end subroutine surf_temp
1369 ! ===6=8===============================================================72
1370 ! ===6=8===============================================================72
1372 subroutine buildings(nd,nz,z0,ua_u,va_u,pt_u,pt0_u, &
1373 ptg,ptr,da_u,ptw, &
1374 drst,uva_u,vva_u,uvb_u,vvb_u, &
1375 tva_u,tvb_u,evb_u, &
1376 uhb_u,vhb_u,thb_u,ehb_u,ss,dt)
1378 ! ----------------------------------------------------------------------
1379 ! This routine computes the sources or sinks of the different quantities
1380 ! on the urban grid. The actual calculation is done in the subroutines
1381 ! called flux_wall, and flux_flat.
1382 ! ----------------------------------------------------------------------
1384 implicit none
1387 ! ----------------------------------------------------------------------
1388 ! INPUT:
1389 ! ----------------------------------------------------------------------
1390 integer nd ! Number of street direction for the current urban class
1391 integer nz ! number of vertical space steps
1392 real ua_u(nz_um) ! Wind speed in the x direction on the urban grid
1393 real va_u(nz_um) ! Wind speed in the y direction on the urban grid
1394 real da_u(nz_um) ! air density on the urban grid
1395 real drst(ndm) ! Street directions for the current urban class
1396 real dz
1397 real pt_u(nz_um) ! Potential temperature on the urban grid
1398 real pt0_u(nz_um) ! reference potential temperature on the urban grid
1399 real ptg(ndm) ! Ground potential temperatures
1400 real ptr(ndm,nz_um) ! Roof potential temperatures
1401 real ptw(2*ndm,nz_um) ! Walls potential temperatures
1402 real ss(nz_um) ! probability to have a building with height h
1403 real z0(ndm,nz_um) ! Roughness lengths "profiles"
1404 real dt ! time step
1407 ! ----------------------------------------------------------------------
1408 ! OUTPUT:
1409 ! ----------------------------------------------------------------------
1410 ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
1411 ! vertical surfaces (walls) and horizontal surfaces (roofs and street)
1412 ! The fluxes can be computed as follow: Fluxes of X = A*X + B
1413 ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
1415 real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
1416 real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
1417 real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
1418 real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
1419 real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
1420 real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
1421 real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1422 real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1423 real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1424 real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1425 real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1427 ! ----------------------------------------------------------------------
1428 ! LOCAL:
1429 ! ----------------------------------------------------------------------
1430 integer id,iz
1432 ! ----------------------------------------------------------------------
1433 ! END VARIABLES DEFINITIONS
1434 ! ----------------------------------------------------------------------
1435 dz=dz_u
1437 do id=1,nd
1439 ! Calculation at the ground surfaces
1440 call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1), &
1441 ptg(id),uhb_u(id,1), &
1442 vhb_u(id,1),thb_u(id,1),ehb_u(id,1))
1444 ! Calculation at the roof surfaces
1445 do iz=2,nz
1446 if(ss(iz).gt.0)then
1447 call flux_flat(dz,z0(id,iz),ua_u(iz), &
1448 va_u(iz),pt_u(iz),pt0_u(iz), &
1449 ptr(id,iz),uhb_u(id,iz), &
1450 vhb_u(id,iz),thb_u(id,iz),ehb_u(id,iz))
1451 else
1452 uhb_u(id,iz) = 0.0
1453 vhb_u(id,iz) = 0.0
1454 thb_u(id,iz) = 0.0
1455 ehb_u(id,iz) = 0.0
1456 endif
1457 end do
1459 ! Calculation at the wall surfaces
1460 do iz=1,nz
1461 call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
1462 ptw(2*id-1,iz), &
1463 uva_u(2*id-1,iz),vva_u(2*id-1,iz), &
1464 uvb_u(2*id-1,iz),vvb_u(2*id-1,iz), &
1465 tva_u(2*id-1,iz),tvb_u(2*id-1,iz), &
1466 evb_u(2*id-1,iz),drst(id),dt)
1468 call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
1469 ptw(2*id,iz), &
1470 uva_u(2*id,iz),vva_u(2*id,iz), &
1471 uvb_u(2*id,iz),vvb_u(2*id,iz), &
1472 tva_u(2*id,iz),tvb_u(2*id,iz), &
1473 evb_u(2*id,iz),drst(id),dt)
1474 !
1476 end do
1478 end do
1480 return
1481 end subroutine buildings
1484 ! ===6=8===============================================================72
1485 ! ===6=8===============================================================72
1487 subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl, &
1488 uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
1489 uhb_u,vhb_u,thb_u,ehb_u, &
1490 a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)
1492 ! ----------------------------------------------------------------------
1493 ! This routine interpolates the parameters from the "urban grid" to the
1494 ! "mesoscale grid".
1495 ! See p300-301 Appendix B.2 of the BLM paper.
1496 ! ----------------------------------------------------------------------
1498 implicit none
1500 ! ----------------------------------------------------------------------
1501 ! INPUT:
1502 ! ----------------------------------------------------------------------
1503 ! Data relative to the "mesoscale grid"
1504 integer kms,kme,kts,kte
1505 real z(kms:kme) ! Altitude above the ground of the cell interface
1506 real dz(kms:kme) ! Vertical space steps
1508 ! Data relative to the "uban grid"
1509 integer nz_u ! Number of layer in the urban grid
1510 integer nd ! Number of street direction for the current urban class
1511 real bs(ndm) ! Building widths of the current urban class
1512 real ws(ndm) ! Street widths of the current urban class
1513 real z_u(nz_um) ! Height of the urban grid levels
1514 real pb(nz_um) ! Probability to have a building with an height equal
1515 real ss(nz_um) ! Probability to have a building with height h
1516 real uhb_u(ndm,nz_um) ! U (x-wind component) Horizontal surfaces, B (explicit) term
1517 real uva_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, A (implicit) term
1518 real uvb_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, B (explicit) term
1519 real vhb_u(ndm,nz_um) ! V (y-wind component) Horizontal surfaces, B (explicit) term
1520 real vva_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, A (implicit) term
1521 real vvb_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, B (explicit) term
1522 real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1523 real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1524 real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1525 real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1526 real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1528 ! ----------------------------------------------------------------------
1529 ! OUTPUT:
1530 ! ----------------------------------------------------------------------
1531 ! Data relative to the "mesoscale grid"
1532 real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
1533 real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
1534 real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
1535 real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
1536 real a_t(kms:kme) ! Implicit component of the heat sources or sinks
1537 real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
1538 real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
1539 real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
1540 real b_t(kms:kme) ! Explicit component of the heat sources or sinks
1541 real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
1543 ! ----------------------------------------------------------------------
1544 ! LOCAL:
1545 ! ----------------------------------------------------------------------
1546 real dzz
1547 real fact
1548 integer id,iz,iz_u
1549 real se,sr,st,su,sv
1550 real uet(kms:kme) ! Contribution to TKE due to walls
1551 real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb
1554 ! ----------------------------------------------------------------------
1555 ! END VARIABLES DEFINITIONS
1556 ! ----------------------------------------------------------------------
1558 ! initialisation
1560 do iz=kts,kte
1561 a_u(iz)=0.
1562 a_v(iz)=0.
1563 a_t(iz)=0.
1564 a_e(iz)=0.
1565 b_u(iz)=0.
1566 b_v(iz)=0.
1567 b_e(iz)=0.
1568 b_t(iz)=0.
1569 uet(iz)=0.
1570 end do
1572 ! horizontal surfaces
1573 do iz=kts,kte
1574 sf(iz)=0.
1575 vl(iz)=0.
1576 enddo
1577 sf(kte+1)=0.
1579 do id=1,nd
1580 do iz=kts+1,kte+1
1581 sr=0.
1582 do iz_u=2,nz_u
1583 if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
1584 sr=pb(iz_u)
1585 endif
1586 enddo
1587 sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
1588 enddo
1589 enddo
1591 ! volume
1592 do iz=kts,kte
1593 do id=1,nd
1594 vtot=0.
1595 do iz_u=1,nz_u
1596 dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
1597 vtot=vtot+pb(iz_u+1)*dzz
1598 enddo
1599 vtot=vtot/(z(iz+1)-z(iz))
1600 vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
1601 enddo
1602 enddo
1604 ! horizontal surface impact
1606 do id=1,nd
1608 fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
1609 b_t(kts)=b_t(kts)+thb_u(id,1)*fact
1610 b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
1611 b_v(kts)=b_v(kts)+vhb_u(id,1)*fact
1612 b_e(kts)=b_e(kts)+ehb_u(id,1)*fact*(z_u(2)-z_u(1))
1614 do iz=kts,kte
1615 st=0.
1616 su=0.
1617 sv=0.
1618 se=0.
1619 do iz_u=2,nz_u
1620 if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
1621 st=st+ss(iz_u)*thb_u(id,iz_u)
1622 su=su+ss(iz_u)*uhb_u(id,iz_u)
1623 sv=sv+ss(iz_u)*vhb_u(id,iz_u)
1624 se=se+ss(iz_u)*ehb_u(id,iz_u)*(z_u(iz_u+1)-z_u(iz_u))
1625 endif
1626 enddo
1628 fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
1629 b_t(iz)=b_t(iz)+st*fact
1630 b_u(iz)=b_u(iz)+su*fact
1631 b_v(iz)=b_v(iz)+sv*fact
1632 b_e(iz)=b_e(iz)+se*fact
1633 enddo
1634 enddo
1636 ! vertical surface impact
1638 do iz=kts,kte
1639 uet(iz)=0.
1640 do id=1,nd
1641 vtb=0.
1642 vta=0.
1643 vua=0.
1644 vub=0.
1645 vva=0.
1646 vvb=0.
1647 veb=0.
1648 vte=0.
1649 do iz_u=1,nz_u
1650 dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
1651 fact=dzz/(ws(id)+bs(id))
1652 vtb=vtb+pb(iz_u+1)* &
1653 (tvb_u(2*id-1,iz_u)+tvb_u(2*id,iz_u))*fact
1654 vta=vta+pb(iz_u+1)* &
1655 (tva_u(2*id-1,iz_u)+tva_u(2*id,iz_u))*fact
1656 vua=vua+pb(iz_u+1)* &
1657 (uva_u(2*id-1,iz_u)+uva_u(2*id,iz_u))*fact
1658 vva=vva+pb(iz_u+1)* &
1659 (vva_u(2*id-1,iz_u)+vva_u(2*id,iz_u))*fact
1660 vub=vub+pb(iz_u+1)* &
1661 (uvb_u(2*id-1,iz_u)+uvb_u(2*id,iz_u))*fact
1662 vvb=vvb+pb(iz_u+1)* &
1663 (vvb_u(2*id-1,iz_u)+vvb_u(2*id,iz_u))*fact
1664 veb=veb+pb(iz_u+1)* &
1665 (evb_u(2*id-1,iz_u)+evb_u(2*id,iz_u))*fact
1666 enddo
1668 fact=1./vl(iz)/dz(iz)/nd
1669 b_t(iz)=b_t(iz)+vtb*fact
1670 a_t(iz)=a_t(iz)+vta*fact
1671 a_u(iz)=a_u(iz)+vua*fact
1672 a_v(iz)=a_v(iz)+vva*fact
1673 b_u(iz)=b_u(iz)+vub*fact
1674 b_v(iz)=b_v(iz)+vvb*fact
1675 b_e(iz)=b_e(iz)+veb*fact
1676 uet(iz)=uet(iz)+vte*fact
1677 enddo
1678 enddo
1682 return
1683 end subroutine urban_meso
1685 ! ===6=8===============================================================72
1687 ! ===6=8===============================================================72
1689 subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
1690 dlg,dl_u)
1692 ! ----------------------------------------------------------------------
1693 ! Calculation of the length scales
1694 ! See p272-274 formula (22) and (24) of the BLM paper
1695 ! ----------------------------------------------------------------------
1697 implicit none
1700 ! ----------------------------------------------------------------------
1701 ! INPUT:
1702 ! ----------------------------------------------------------------------
1703 integer kms,kme,kts,kte
1704 real z(kms:kme) ! Altitude above the ground of the cell interface
1705 integer nd ! Number of street direction for the current urban class
1706 integer nz_u ! Number of levels in the "urban grid"
1707 real z_u(nz_um) ! Height of the urban grid levels
1708 real bs(ndm) ! Building widths of the current urban class
1709 real ss(nz_um) ! Probability to have a building with height h
1710 real ws(ndm) ! Street widths of the current urban class
1713 ! ----------------------------------------------------------------------
1714 ! OUTPUT:
1715 ! ----------------------------------------------------------------------
1716 real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
1717 real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
1719 ! ----------------------------------------------------------------------
1720 ! LOCAL:
1721 ! ----------------------------------------------------------------------
1722 real dlgtmp
1723 integer id,iz,iz_u
1724 real sftot
1725 real ulu,ssl
1727 ! ----------------------------------------------------------------------
1728 ! END VARIABLES DEFINITIONS
1729 ! ----------------------------------------------------------------------
1731 do iz=kts,kte
1732 ulu=0.
1733 ssl=0.
1734 do id=1,nd
1735 do iz_u=2,nz_u
1736 if(z_u(iz_u).gt.z(iz))then
1737 ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
1738 ssl=ssl+ss(iz_u)/nd
1739 endif
1740 enddo
1741 enddo
1743 if(ulu.ne.0)then
1744 dl_u(iz)=ssl/ulu
1745 else
1746 dl_u(iz)=0.
1747 endif
1748 enddo
1751 do iz=kts,kte
1752 dlg(iz)=0.
1753 do id=1,nd
1754 sftot=ws(id)
1755 dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
1756 do iz_u=1,nz_u
1757 if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
1758 dlgtmp=dlgtmp+ss(iz_u)*bs(id)/ &
1759 ((z(iz)+z(iz+1))/2.-z_u(iz_u))
1760 sftot=sftot+ss(iz_u)*bs(id)
1761 endif
1762 enddo
1763 dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
1764 enddo
1765 dlg(iz)=1./dlg(iz)
1766 enddo
1768 return
1769 end subroutine interp_length
1771 ! ===6=8===============================================================72
1772 ! ===6=8===============================================================72
1774 subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
1775 rs,rsw,rsg)
1777 ! ----------------------------------------------------------------------
1778 ! Modification of short wave radiation to take into account
1779 ! the shadow produced by the buildings
1780 ! ----------------------------------------------------------------------
1782 implicit none
1784 ! ----------------------------------------------------------------------
1785 ! INPUT:
1786 ! ----------------------------------------------------------------------
1787 integer nd ! Number of street direction for the current urban class
1788 integer nz_u ! number of vertical layers defined in the urban grid
1789 real ah ! Hour angle (it should come from the radiation routine)
1790 real deltar ! Declination of the sun
1791 real drst(ndm) ! street directions for the current urban class
1792 real rs ! solar radiation
1793 real ss(nz_um) ! probability to have a building with height h
1794 real pb(nz_um) ! Probability that a building has an height greater or equal to h
1795 real ws(ndm) ! Street width of the current urban class
1796 real z(nz_um) ! Height of the urban grid levels
1797 real zr ! zenith angle
1799 ! ----------------------------------------------------------------------
1800 ! OUTPUT:
1801 ! ----------------------------------------------------------------------
1802 real rsg(ndm) ! Short wave radiation at the ground
1803 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1805 ! ----------------------------------------------------------------------
1806 ! LOCAL:
1807 ! ----------------------------------------------------------------------
1808 integer id,iz,jz
1809 real aae,aaw,bbb,phix,rd,rtot,wsd
1811 ! ----------------------------------------------------------------------
1812 ! END VARIABLES DEFINITIONS
1813 ! ----------------------------------------------------------------------
1815 if(rs.eq.0.or.sin(zr).eq.1)then
1816 do id=1,nd
1817 rsg(id)=0.
1818 do iz=1,nz_u
1819 rsw(2*id-1,iz)=0.
1820 rsw(2*id,iz)=0.
1821 enddo
1822 enddo
1823 else
1824 !test
1825 if(abs(sin(zr)).gt.1.e-10)then
1826 if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
1827 bbb=pi/2.
1828 elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
1829 bbb=-pi/2.
1830 else
1831 bbb=asin(cos(deltar)*sin(ah)/sin(zr))
1832 endif
1833 else
1834 if(cos(deltar)*sin(ah).ge.0)then
1835 bbb=pi/2.
1836 elseif(cos(deltar)*sin(ah).lt.0)then
1837 bbb=-pi/2.
1838 endif
1839 endif
1841 phix=zr
1843 do id=1,nd
1845 rsg(id)=0.
1847 aae=bbb-drst(id)
1848 aaw=bbb-drst(id)+pi
1850 do iz=1,nz_u
1851 rsw(2*id-1,iz)=0.
1852 rsw(2*id,iz)=0.
1853 if(pb(iz+1).gt.0.)then
1854 do jz=1,nz_u
1855 if(abs(sin(aae)).gt.1.e-10)then
1856 call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae, &
1857 ws(id),rd)
1858 rsw(2*id-1,iz)=rsw(2*id-1,iz)+rs*rd*ss(jz+1)/pb(iz+1)
1859 endif
1861 if(abs(sin(aaw)).gt.1.e-10)then
1862 call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw, &
1863 ws(id),rd)
1864 rsw(2*id,iz)=rsw(2*id,iz)+rs*rd*ss(jz+1)/pb(iz+1)
1865 endif
1866 enddo
1867 endif
1868 enddo
1869 if(abs(sin(aae)).gt.1.e-10)then
1870 wsd=abs(ws(id)/sin(aae))
1872 do jz=1,nz_u
1873 rd=max(0.,wsd-z(jz+1)*tan(phix))
1874 rsg(id)=rsg(id)+rs*rd*ss(jz+1)/wsd
1875 enddo
1876 rtot=0.
1878 do iz=1,nz_u
1879 rtot=rtot+(rsw(2*id,iz)+rsw(2*id-1,iz))* &
1880 (z(iz+1)-z(iz))
1881 enddo
1882 rtot=rtot+rsg(id)*ws(id)
1883 else
1884 rsg(id)=rs
1885 endif
1887 enddo
1888 endif
1890 return
1891 end subroutine shadow_mas
1893 ! ===6=8===============================================================72
1894 ! ===6=8===============================================================72
1896 subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
1898 ! ----------------------------------------------------------------------
1899 ! This routine computes the effects of a shadow induced by a building of
1900 ! height hu, on a portion of wall between z1 and z2. See equation A10,
1901 ! and correction described below formula A11, and figure A1. Basically rd
1902 ! is the ratio between the horizontal surface illuminated and the portion
1903 ! of wall. Referring to figure A1, multiplying radiation flux density on
1904 ! a horizontal surface (rs) by x1-x2 we have the radiation energy per
1905 ! unit time. Dividing this by z2-z1, we obtain the radiation flux
1906 ! density reaching the portion of the wall between z2 and z1
1907 ! (everything is assumed in 2D)
1908 ! ----------------------------------------------------------------------
1910 implicit none
1912 ! ----------------------------------------------------------------------
1913 ! INPUT:
1914 ! ----------------------------------------------------------------------
1915 real aa ! Angle between the sun direction and the face of the wall (A12)
1916 real hu ! Height of the building that generates the shadow
1917 real phix ! Solar zenith angle
1918 real ws ! Width of the street
1919 real z1 ! Height of the level z(iz)
1920 real z2 ! Height of the level z(iz+1)
1922 ! ----------------------------------------------------------------------
1923 ! OUTPUT:
1924 ! ----------------------------------------------------------------------
1925 real rd ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A.
1926 ! Multiplying rd by rs (radiation flux
1927 ! density on a horizontal surface) gives
1928 ! the radiation flux density on the
1929 ! portion of wall between z1 and z2.
1930 ! ----------------------------------------------------------------------
1931 ! LOCAL:
1932 ! ----------------------------------------------------------------------
1933 real x1,x2 ! x1,x2 see Fig. A1.
1935 ! ----------------------------------------------------------------------
1936 ! END VARIABLES DEFINITIONS
1937 ! ----------------------------------------------------------------------
1939 x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
1941 x2=max((hu-z2)*tan(phix),0.)
1943 rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
1945 return
1946 end subroutine shade_wall
1948 ! ===6=8===============================================================72
1949 ! ===6=8===============================================================72
1951 subroutine long_rad(iurb,nz_u,id,emw,emg, &
1952 fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
1954 ! ----------------------------------------------------------------------
1955 ! This routine computes the effects of the reflections of long-wave
1956 ! radiation in the street canyon by solving the system
1957 ! of 2*nz_u+1 eqn. in 2*nz_u+1
1958 ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
1959 ! The system is solved by solving A X= B,
1960 ! with A matrix B vector, and X solution.
1961 ! ----------------------------------------------------------------------
1963 implicit none
1967 ! ----------------------------------------------------------------------
1968 ! INPUT:
1969 ! ----------------------------------------------------------------------
1970 real emg ! Emissivity of ground for the current urban class
1971 real emw ! Emissivity of wall for the current urban class
1972 real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
1973 real fsg(ndm,nurbm) ! View factors from sky to ground
1974 real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
1975 real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
1976 real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
1977 integer id ! Current street direction
1978 integer iurb ! Current urban class
1979 integer nz_u ! Number of layer in the urban grid
1980 real pb(nz_um) ! Probability to have a building with an height equal
1981 real rl ! Downward flux of the longwave radiation
1982 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1983 real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
1986 ! ----------------------------------------------------------------------
1987 ! OUTPUT:
1988 ! ----------------------------------------------------------------------
1989 real rlg(ndm) ! Long wave radiation at the ground
1990 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1992 ! ----------------------------------------------------------------------
1993 ! LOCAL:
1994 ! ----------------------------------------------------------------------
1995 integer i,j
1996 real aaa(2*nz_um+1,2*nz_um+1) ! terms of the matrix
1997 real bbb(2*nz_um+1) ! terms of the vector
1999 ! ----------------------------------------------------------------------
2000 ! END VARIABLES DEFINITIONS
2001 ! ----------------------------------------------------------------------
2004 ! west wall
2006 do i=1,nz_u
2008 do j=1,nz_u
2009 aaa(i,j)=0.
2010 enddo
2012 aaa(i,i)=1.
2014 do j=nz_u+1,2*nz_u
2015 aaa(i,j)=-(1.-emw)*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
2016 enddo
2018 !! aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)*pb(i+1)
2019 aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
2021 bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg(id,ng_u)**4
2022 do j=1,nz_u
2023 bbb(i)=bbb(i)+pb(j+1)*emw*sigma*fww(j,i,id,iurb)* &
2024 tw(2*id,j,nwr_u)**4+ &
2025 fww(j,i,id,iurb)*rl*(1.-pb(j+1))
2026 enddo
2028 enddo
2030 ! east wall
2032 do i=1+nz_u,2*nz_u
2034 do j=1,nz_u
2035 aaa(i,j)=-(1.-emw)*fww(j,i-nz_u,id,iurb)*pb(j+1)
2036 enddo
2038 do j=1+nz_u,2*nz_u
2039 aaa(i,j)=0.
2040 enddo
2042 aaa(i,i)=1.
2044 !! aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)*pb(i-nz_u+1)
2045 aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
2047 bbb(i)=fsw(i-nz_u,id,iurb)*rl+ &
2048 emg*fgw(i-nz_u,id,iurb)*sigma*tg(id,ng_u)**4
2050 do j=1,nz_u
2051 bbb(i)=bbb(i)+pb(j+1)*emw*sigma*fww(j,i-nz_u,id,iurb)* &
2052 tw(2*id-1,j,nwr_u)**4+ &
2053 fww(j,i-nz_u,id,iurb)*rl*(1.-pb(j+1))
2054 enddo
2056 enddo
2058 ! ground
2059 do j=1,nz_u
2060 aaa(2*nz_u+1,j)=-(1.-emw)*fwg(j,id,iurb)*pb(j+1)
2061 enddo
2063 do j=nz_u+1,2*nz_u
2064 aaa(2*nz_u+1,j)=-(1.-emw)*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
2065 enddo
2067 aaa(2*nz_u+1,2*nz_u+1)=1.
2069 bbb(2*nz_u+1)=fsg(id,iurb)*rl
2071 do i=1,nz_u
2072 bbb(2*nz_u+1)=bbb(2*nz_u+1)+emw*sigma*fwg(i,id,iurb)*pb(i+1)* &
2073 (tw(2*id-1,i,nwr_u)**4+tw(2*id,i,nwr_u)**4)+ &
2074 2.*fwg(i,id,iurb)*(1.-pb(i+1))*rl
2075 enddo
2079 call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
2081 do i=1,nz_u
2082 rlw(2*id-1,i)=bbb(i)
2083 enddo
2085 do i=nz_u+1,2*nz_u
2086 rlw(2*id,i-nz_u)=bbb(i)
2087 enddo
2089 rlg(id)=bbb(2*nz_u+1)
2091 return
2092 end subroutine long_rad
2094 ! ===6=8===============================================================72
2095 ! ===6=8===============================================================72
2097 subroutine short_rad(iurb,nz_u,id,albw, &
2098 albg,fwg,fww,fgw,rsg,rsw,pb)
2100 ! ----------------------------------------------------------------------
2101 ! This routine computes the effects of the reflections of short-wave
2102 ! (solar) radiation in the street canyon by solving the system
2103 ! of 2*nz_u+1 eqn. in 2*nz_u+1
2104 ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
2105 ! The system is solved by solving A X= B,
2106 ! with A matrix B vector, and X solution.
2107 ! ----------------------------------------------------------------------
2109 implicit none
2113 ! ----------------------------------------------------------------------
2114 ! INPUT:
2115 ! ----------------------------------------------------------------------
2116 real albg ! Albedo of the ground for the current urban class
2117 real albw ! Albedo of the wall for the current urban class
2118 real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
2119 real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
2120 real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
2121 integer id ! current street direction
2122 integer iurb ! current urban class
2123 integer nz_u ! Number of layer in the urban grid
2124 real pb(nz_um) ! probability to have a building with an height equal
2126 ! ----------------------------------------------------------------------
2127 ! OUTPUT:
2128 ! ----------------------------------------------------------------------
2129 real rsg(ndm) ! Short wave radiation at the ground
2130 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
2132 ! ----------------------------------------------------------------------
2133 ! LOCAL:
2134 ! ----------------------------------------------------------------------
2135 integer i,j
2136 real aaa(2*nz_um+1,2*nz_um+1) ! terms of the matrix
2137 real bbb(2*nz_um+1) ! terms of the vector
2139 ! ----------------------------------------------------------------------
2140 ! END VARIABLES DEFINITIONS
2141 ! ----------------------------------------------------------------------
2144 ! west wall
2146 do i=1,nz_u
2147 do j=1,nz_u
2148 aaa(i,j)=0.
2149 enddo
2151 aaa(i,i)=1.
2153 do j=nz_u+1,2*nz_u
2154 aaa(i,j)=-albw*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
2155 enddo
2157 aaa(i,2*nz_u+1)=-albg*fgw(i,id,iurb)
2158 bbb(i)=rsw(2*id-1,i)
2160 enddo
2162 ! east wall
2164 do i=1+nz_u,2*nz_u
2165 do j=1,nz_u
2166 aaa(i,j)=-albw*fww(j,i-nz_u,id,iurb)*pb(j+1)
2167 enddo
2169 do j=1+nz_u,2*nz_u
2170 aaa(i,j)=0.
2171 enddo
2173 aaa(i,i)=1.
2174 aaa(i,2*nz_u+1)=-albg*fgw(i-nz_u,id,iurb)
2175 bbb(i)=rsw(2*id,i-nz_u)
2177 enddo
2179 ! ground
2181 do j=1,nz_u
2182 aaa(2*nz_u+1,j)=-albw*fwg(j,id,iurb)*pb(j+1)
2183 enddo
2185 do j=nz_u+1,2*nz_u
2186 aaa(2*nz_u+1,j)=-albw*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
2187 enddo
2189 aaa(2*nz_u+1,2*nz_u+1)=1.
2190 bbb(2*nz_u+1)=rsg(id)
2192 call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
2194 do i=1,nz_u
2195 rsw(2*id-1,i)=bbb(i)
2196 enddo
2198 do i=nz_u+1,2*nz_u
2199 rsw(2*id,i-nz_u)=bbb(i)
2200 enddo
2202 rsg(id)=bbb(2*nz_u+1)
2204 return
2205 end subroutine short_rad
2208 ! ===6=8===============================================================72
2209 ! ===6=8===============================================================72
2211 subroutine gaussj(a,n,b,np)
2213 ! ----------------------------------------------------------------------
2214 ! This routine solve a linear system of n equations of the form
2215 ! A X = B
2216 ! where A is a matrix a(i,j)
2217 ! B a vector and X the solution
2218 ! In output b is replaced by the solution
2219 ! ----------------------------------------------------------------------
2221 implicit none
2223 ! ----------------------------------------------------------------------
2224 ! INPUT:
2225 ! ----------------------------------------------------------------------
2226 integer np
2227 real a(np,np)
2229 ! ----------------------------------------------------------------------
2230 ! OUTPUT:
2231 ! ----------------------------------------------------------------------
2232 real b(np)
2234 ! ----------------------------------------------------------------------
2235 ! LOCAL:
2236 ! ----------------------------------------------------------------------
2237 integer nmax
2238 parameter (nmax=150)
2240 real big,dum
2241 integer i,icol,irow
2242 integer j,k,l,ll,n
2243 integer ipiv(nmax)
2244 real pivinv
2246 ! ----------------------------------------------------------------------
2247 ! END VARIABLES DEFINITIONS
2248 ! ----------------------------------------------------------------------
2250 do j=1,n
2251 ipiv(j)=0.
2252 enddo
2254 do i=1,n
2255 big=0.
2256 do j=1,n
2257 if(ipiv(j).ne.1)then
2258 do k=1,n
2259 if(ipiv(k).eq.0)then
2260 if(abs(a(j,k)).ge.big)then
2261 big=abs(a(j,k))
2262 irow=j
2263 icol=k
2264 endif
2265 elseif(ipiv(k).gt.1)then
2266 FATAL_ERROR('singular matrix in gaussj')
2267 endif
2268 enddo
2269 endif
2270 enddo
2272 ipiv(icol)=ipiv(icol)+1
2274 if(irow.ne.icol)then
2275 do l=1,n
2276 dum=a(irow,l)
2277 a(irow,l)=a(icol,l)
2278 a(icol,l)=dum
2279 enddo
2281 dum=b(irow)
2282 b(irow)=b(icol)
2283 b(icol)=dum
2285 endif
2287 if(a(icol,icol).eq.0) FATAL_ERROR('singular matrix in gaussj')
2289 pivinv=1./a(icol,icol)
2290 a(icol,icol)=1
2292 do l=1,n
2293 a(icol,l)=a(icol,l)*pivinv
2294 enddo
2296 b(icol)=b(icol)*pivinv
2298 do ll=1,n
2299 if(ll.ne.icol)then
2300 dum=a(ll,icol)
2301 a(ll,icol)=0.
2302 do l=1,n
2303 a(ll,l)=a(ll,l)-a(icol,l)*dum
2304 enddo
2306 b(ll)=b(ll)-b(icol)*dum
2308 endif
2309 enddo
2310 enddo
2312 return
2313 end subroutine gaussj
2315 ! ===6=8===============================================================72
2316 ! ===6=8===============================================================72
2318 subroutine soil_temp(nz,dz,temp,pt,ala,cs, &
2319 rs,rl,press,dt,em,alb,rt,sf,gf)
2321 ! ----------------------------------------------------------------------
2322 ! This routine solves the Fourier diffusion equation for heat in
2323 ! the material (wall, roof, or ground). Resolution is done implicitely.
2324 ! Boundary conditions are:
2325 ! - fixed temperature at the interior
2326 ! - energy budget at the surface
2327 ! ----------------------------------------------------------------------
2329 implicit none
2333 ! ----------------------------------------------------------------------
2334 ! INPUT:
2335 ! ----------------------------------------------------------------------
2336 integer nz ! Number of layers
2337 real ala(nz) ! Thermal diffusivity in each layers [m^2 s^-1]
2338 real alb ! Albedo of the surface
2339 real cs(nz) ! Specific heat of the material [J m^3 K^-1]
2340 real dt ! Time step
2341 real em ! Emissivity of the surface
2342 real press ! Pressure at ground level
2343 real rl ! Downward flux of the longwave radiation
2344 real rs ! Solar radiation
2345 real sf ! Sensible heat flux at the surface
2346 real temp(nz) ! Temperature in each layer [K]
2347 real dz(nz) ! Layer sizes [m]
2350 ! ----------------------------------------------------------------------
2351 ! OUTPUT:
2352 ! ----------------------------------------------------------------------
2353 real gf ! Heat flux transferred from the surface toward the interior
2354 real pt ! Potential temperature at the surface
2355 real rt ! Total radiation at the surface (solar+incoming long+outgoing long)
2357 ! ----------------------------------------------------------------------
2358 ! LOCAL:
2359 ! ----------------------------------------------------------------------
2360 integer iz
2361 real a(nz,3)
2362 real alpha
2363 real c(nz)
2364 real cddz(nz+2)
2365 real tsig
2367 ! ----------------------------------------------------------------------
2368 ! END VARIABLES DEFINITIONS
2369 ! ----------------------------------------------------------------------
2371 tsig=temp(nz)
2372 alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
2373 ! Compute cddz=2*cd/dz
2375 cddz(1)=ala(1)/dz(1)
2376 do iz=2,nz
2377 cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
2378 enddo
2379 ! cddz(nz+1)=ala(nz+1)/dz(nz)
2381 a(1,1)=0.
2382 a(1,2)=1.
2383 a(1,3)=0.
2384 c(1)=temp(1)
2386 do iz=2,nz-1
2387 a(iz,1)=-cddz(iz)*dt/dz(iz)
2388 a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
2389 a(iz,3)=-cddz(iz+1)*dt/dz(iz)
2390 c(iz)=temp(iz)
2391 enddo
2393 a(nz,1)=-dt*cddz(nz)/dz(nz)
2394 a(nz,2)=1.+dt*cddz(nz)/dz(nz)
2395 a(nz,3)=0.
2396 c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
2399 call invert(nz,a,c,temp)
2402 pt=temp(nz)*(press/1.e+5)**(-rcp_u)
2404 rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)
2406 ! gf=-cddz(nz)*(temp(nz)-temp(nz-1))*cs(nz)
2407 gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
2408 return
2409 end subroutine soil_temp
2411 ! ===6=8===============================================================72
2412 ! ===6=8===============================================================72
2414 subroutine invert(n,a,c,x)
2416 ! ----------------------------------------------------------------------
2417 ! Inversion and resolution of a tridiagonal matrix
2418 ! A X = C
2419 ! ----------------------------------------------------------------------
2421 implicit none
2423 ! ----------------------------------------------------------------------
2424 ! INPUT:
2425 ! ----------------------------------------------------------------------
2426 integer n
2427 real a(n,3) ! a(*,1) lower diagonal (Ai,i-1)
2428 ! a(*,2) principal diagonal (Ai,i)
2429 ! a(*,3) upper diagonal (Ai,i+1)
2430 real c(n)
2432 ! ----------------------------------------------------------------------
2433 ! OUTPUT:
2434 ! ----------------------------------------------------------------------
2435 real x(n)
2437 ! ----------------------------------------------------------------------
2438 ! LOCAL:
2439 ! ----------------------------------------------------------------------
2440 integer i
2442 ! ----------------------------------------------------------------------
2443 ! END VARIABLES DEFINITIONS
2444 ! ----------------------------------------------------------------------
2446 do i=n-1,1,-1
2447 c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
2448 a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
2449 enddo
2451 do i=2,n
2452 c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
2453 enddo
2455 do i=1,n
2456 x(i)=c(i)/a(i,2)
2457 enddo
2459 return
2460 end subroutine invert
2463 ! ===6=8===============================================================72
2464 ! ===6=8===============================================================72
2466 subroutine flux_wall(ua,va,pt,da,ptw,uva,vva,uvb,vvb, &
2467 tva,tvb,evb,drst,dt)
2469 ! ----------------------------------------------------------------------
2470 ! This routine computes the surface sources or sinks of momentum, tke,
2471 ! and heat from vertical surfaces (walls).
2472 ! ----------------------------------------------------------------------
2474 implicit none
2478 ! INPUT:
2479 ! -----
2480 real drst ! street directions for the current urban class
2481 real da ! air density
2482 real pt ! potential temperature
2483 real ptw ! Walls potential temperatures
2484 real ua ! wind speed
2485 real va ! wind speed
2487 real dt !time step
2488 ! OUTPUT:
2489 ! ------
2490 ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
2491 ! vertical surfaces (walls).
2492 ! The fluxes can be computed as follow: Fluxes of X = A*X + B
2493 ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
2494 real uva ! U (wind component) Vertical surfaces, A (implicit) term
2495 real uvb ! U (wind component) Vertical surfaces, B (explicit) term
2496 real vva ! V (wind component) Vertical surfaces, A (implicit) term
2497 real vvb ! V (wind component) Vertical surfaces, B (explicit) term
2498 real tva ! Temperature Vertical surfaces, A (implicit) term
2499 real tvb ! Temperature Vertical surfaces, B (explicit) term
2500 real evb ! Energy (TKE) Vertical surfaces, B (explicit) term
2502 ! LOCAL:
2503 ! -----
2504 real hc
2505 real u_ort
2506 real vett
2508 ! -------------------------
2509 ! END VARIABLES DEFINITIONS
2510 ! -------------------------
2512 vett=(ua**2+va**2)**.5
2514 u_ort=abs((cos(drst)*ua-sin(drst)*va))
2516 uva=-cdrag*u_ort/2.*cos(drst)*cos(drst)
2517 vva=-cdrag*u_ort/2.*sin(drst)*sin(drst)
2519 uvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*va
2520 vvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*ua
2522 hc=5.678*(1.09+0.23*(vett/0.3048))
2524 if(hc.gt.da*cp_u/dt)then
2525 hc=da*cp_u/dt
2526 endif
2528 ! tvb=hc*ptw/da/cp_u
2529 ! tva=-hc/da/cp_u
2530 !!!!!!!!!!!!!!!!!!!!
2531 ! explicit
2532 tvb=hc*ptw/da/cp_u-hc/da/cp_u*pt !c
2533 tva = 0. !c
2535 evb=cdrag*(abs(u_ort)**3.)/2.
2537 return
2538 end subroutine flux_wall
2540 ! ===6=8===============================================================72
2542 ! ===6=8===============================================================72
2544 subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg, &
2545 uhb,vhb,thb,ehb)
2547 ! ----------------------------------------------------------------------
2548 ! Calculation of the flux at the ground
2549 ! Formulation of Louis (Louis, 1979)
2550 ! ----------------------------------------------------------------------
2552 implicit none
2556 ! ----------------------------------------------------------------------
2557 ! INPUT:
2558 ! ----------------------------------------------------------------------
2559 real dz ! first vertical level
2560 real pt ! potential temperature
2561 real pt0 ! reference potential temperature
2562 real ptg ! ground potential temperature
2563 real ua ! wind speed
2564 real va ! wind speed
2565 real z0 ! Roughness length
2567 ! ----------------------------------------------------------------------
2568 ! OUTPUT:
2569 ! ----------------------------------------------------------------------
2570 ! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
2571 ! surfaces (roofs and street)
2572 ! The fluxes can be computed as follow: Fluxes of X = B
2573 ! Example: Momentum fluxes on horizontal surfaces = uhb_u
2574 real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
2575 real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
2576 real thb ! Temperature Horizontal surfaces, B (explicit) term
2577 real tva ! Temperature Vertical surfaces, A (implicit) term
2578 real tvb ! Temperature Vertical surfaces, B (explicit) term
2579 real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
2582 ! ----------------------------------------------------------------------
2583 ! LOCAL:
2584 ! ----------------------------------------------------------------------
2585 real aa
2586 real al
2587 real buu
2588 real c
2589 real fbuw
2590 real fbpt
2591 real fh
2592 real fm
2593 real ric
2594 real tstar
2595 real ustar
2596 real utot
2597 real wstar
2598 real zz
2600 real b,cm,ch,rr,tol
2601 parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
2603 ! ----------------------------------------------------------------------
2604 ! END VARIABLES DEFINITIONS
2605 ! ----------------------------------------------------------------------
2608 ! computation of the ground temperature
2610 utot=(ua**2+va**2)**.5
2613 !!!! Louis formulation
2614 !
2615 ! compute the bulk Richardson Number
2617 zz=dz/2.
2619 ! if(tstar.lt.0.)then
2620 ! wstar=(-ustar*tstar*g*hii/pt)**(1./3.)
2621 ! else
2622 ! wstar=0.
2623 ! endif
2624 !
2625 ! if (utot.le.0.7*wstar) utot=max(0.7*wstar,0.00001)
2627 utot=max(utot,0.01)
2629 ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
2631 aa=vk/log(zz/z0)
2633 ! determine the parameters fm and fh for stable, neutral and unstable conditions
2635 if(ric.gt.0)then
2636 fm=1/(1+0.5*b*ric)**2
2637 fh=fm
2638 else
2639 c=b*cm*aa*aa*(zz/z0)**.5
2640 fm=1-b*ric/(1+c*(-ric)**.5)
2641 c=c*ch/cm
2642 fh=1-b*ric/(1+c*(-ric)**.5)
2643 endif
2645 fbuw=-aa*aa*utot*utot*fm
2646 fbpt=-aa*aa*utot*(pt-ptg)*fh/rr
2648 ustar=(-fbuw)**.5
2649 tstar=-fbpt/ustar
2651 al=(vk*g_u*tstar)/(pt*ustar*ustar)
2653 buu=-g_u/pt0*ustar*tstar
2655 uhb=-ustar*ustar*ua/utot
2656 vhb=-ustar*ustar*va/utot
2657 thb=-ustar*tstar
2658 ! thb= 0.
2659 ehb=buu
2660 !!!!!!!!!!!!!!!
2662 return
2663 end subroutine flux_flat
2665 ! ===6=8===============================================================72
2666 ! ===6=8===============================================================72
2668 subroutine icBEP (nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u)
2670 implicit none
2673 ! Street parameters
2674 integer nd_u(nurbm) ! Number of street direction for each urban class
2675 real h_b(nz_um,nurbm) ! Bulding's heights [m]
2676 real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
2677 ! -----------------------------------------------------------------------
2678 ! Output
2679 !------------------------------------------------------------------------
2681 real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to z
2682 real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to z
2684 ! Grid parameters
2685 integer nz_u(nurbm) ! Number of layer in the urban grid
2686 real z_u(nz_um) ! Height of the urban grid levels
2689 ! -----------------------------------------------------------------------
2690 ! Local
2691 !------------------------------------------------------------------------
2693 integer iz_u,id,ilu,iurb
2695 real dtot
2696 real hbmax
2698 !------------------------------------------------------------------------
2701 ! -----------------------------------------------------------------------
2702 ! This routine initialise the urban paramters for the BEP module
2703 !------------------------------------------------------------------------
2704 !
2705 !Initialize variables
2706 !
2707 z_u=0.
2708 nz_u=0
2709 ss_u=0.
2710 pb_u=0.
2712 ! Computation of the urban levels height
2714 z_u(1)=0.
2716 do iz_u=1,nz_um-1
2717 z_u(iz_u+1)=z_u(iz_u)+dz_u
2718 enddo
2720 ! Normalisation of the building density
2722 do iurb=1,nurbm
2723 dtot=0.
2724 do ilu=1,nz_um
2725 dtot=dtot+d_b(ilu,iurb)
2726 enddo
2727 do ilu=1,nz_um
2728 d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
2729 enddo
2730 enddo
2732 ! Compute the view factors, pb and ss
2734 do iurb=1,nurbm
2735 hbmax=0.
2736 nz_u(iurb)=0
2737 do ilu=1,nz_um
2738 if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
2739 enddo
2741 do iz_u=1,nz_um-1
2742 if(z_u(iz_u+1).gt.hbmax)go to 10
2743 enddo
2745 10 continue
2746 nz_u(iurb)=iz_u+1
2748 do id=1,nd_u(iurb)
2750 do iz_u=1,nz_u(iurb)
2751 ss_u(iz_u,iurb)=0.
2752 do ilu=1,nz_um
2753 if(z_u(iz_u).le.h_b(ilu,iurb) &
2754 .and.z_u(iz_u+1).gt.h_b(ilu,iurb))then
2755 ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
2756 endif
2757 enddo
2758 enddo
2760 pb_u(1,iurb)=1.
2761 do iz_u=1,nz_u(iurb)
2762 pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
2763 enddo
2765 enddo
2766 end do
2769 return
2770 end subroutine icBEP
2772 ! ===6=8===============================================================72
2773 ! ===6=8===============================================================72
2775 subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws)
2777 implicit none
2781 ! -----------------------------------------------------------------------
2782 ! Input
2783 !------------------------------------------------------------------------
2785 integer iurb ! Number of the urban class
2786 integer nz_u ! Number of levels in the urban grid
2787 integer id ! Street direction number
2788 real ws ! Street width
2789 real z(nz_um) ! Height of the urban grid levels
2790 real dxy ! Street lenght
2793 ! -----------------------------------------------------------------------
2794 ! Output
2795 !------------------------------------------------------------------------
2797 ! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
2798 ! and the short wave radation. They are the part of radiation from a surface
2799 ! or from the sky to another surface.
2801 real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
2802 real fwg(nz_um,ndm,nurbm) ! from wall to ground
2803 real fgw(nz_um,ndm,nurbm) ! from ground to wall
2804 real fsw(nz_um,ndm,nurbm) ! from sky to wall
2805 real fws(nz_um,ndm,nurbm) ! from wall to sky
2806 real fsg(ndm,nurbm) ! from sky to ground
2809 ! -----------------------------------------------------------------------
2810 ! Local
2811 !------------------------------------------------------------------------
2813 integer jz,iz
2815 real hut
2816 real f1,f2,f12,f23,f123,ftot
2817 real fprl,fnrm
2818 real a1,a2,a3,a4,a12,a23,a123
2820 ! -----------------------------------------------------------------------
2821 ! This routine calculates the view factors
2822 !------------------------------------------------------------------------
2824 hut=z(nz_u+1)
2826 do jz=1,nz_u
2828 ! radiation from wall to wall
2830 do iz=1,nz_u
2832 call fprls (fprl,dxy,abs(z(jz+1)-z(iz )),ws)
2833 f123=fprl
2834 call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
2835 f23=fprl
2836 call fprls (fprl,dxy,abs(z(jz )-z(iz )),ws)
2837 f12=fprl
2838 call fprls (fprl,dxy,abs(z(jz )-z(iz+1)),ws)
2839 f2 = fprl
2841 a123=dxy*(abs(z(jz+1)-z(iz )))
2842 a12 =dxy*(abs(z(jz )-z(iz )))
2843 a23 =dxy*(abs(z(jz+1)-z(iz+1)))
2844 a1 =dxy*(abs(z(iz+1)-z(iz )))
2845 a2 =dxy*(abs(z(jz )-z(iz+1)))
2846 a3 =dxy*(abs(z(jz+1)-z(jz )))
2848 ftot=0.5*(a123*f123-a23*f23-a12*f12+a2*f2)/a1
2850 fww(iz,jz,id,iurb)=ftot*a1/a3
2852 enddo
2854 ! radiation from ground to wall
2856 call fnrms (fnrm,z(jz+1),dxy,ws)
2857 f12=fnrm
2858 call fnrms (fnrm,z(jz) ,dxy,ws)
2859 f1=fnrm
2861 a1 = ws*dxy
2863 a12= ws*dxy
2865 a4=(z(jz+1)-z(jz))*dxy
2867 ftot=(a12*f12-a12*f1)/a1
2869 fgw(jz,id,iurb)=ftot*a1/a4
2871 ! radiation from sky to wall
2873 call fnrms(fnrm,hut-z(jz) ,dxy,ws)
2874 f12 = fnrm
2875 call fnrms (fnrm,hut-z(jz+1),dxy,ws)
2876 f1 =fnrm
2878 a1 = ws*dxy
2880 a12= ws*dxy
2882 a4 = (z(jz+1)-z(jz))*dxy
2884 ftot=(a12*f12-a12*f1)/a1
2886 fsw(jz,id,iurb)=ftot*a1/a4
2888 enddo
2890 ! radiation from wall to sky
2891 do iz=1,nz_u
2892 call fnrms(fnrm,ws,dxy,hut-z(iz))
2893 f12=fnrm
2894 call fnrms(fnrm,ws,dxy,hut-z(iz+1))
2895 f1=fnrm
2896 a1 = (z(iz+1)-z(iz))*dxy
2897 a2 = (hut-z(iz+1))*dxy
2898 a12= (hut-z(iz))*dxy
2899 a4 = ws*dxy
2900 ftot=(a12*f12-a2*f1)/a1
2901 fws(iz,id,iurb)=ftot*a1/a4
2903 enddo
2904 !!!!!!!!!!!!!
2907 do iz=1,nz_u
2909 ! radiation from wall to ground
2911 call fnrms (fnrm,ws,dxy,z(iz+1))
2912 f12=fnrm
2913 call fnrms (fnrm,ws,dxy,z(iz ))
2914 f1 =fnrm
2916 a1= (z(iz+1)-z(iz) )*dxy
2918 a2 = z(iz)*dxy
2919 a12= z(iz+1)*dxy
2920 a4 = ws*dxy
2922 ftot=(a12*f12-a2*f1)/a1
2924 fwg(iz,id,iurb)=ftot*a1/a4
2926 enddo
2928 ! radiation from sky to ground
2930 call fprls (fprl,dxy,ws,hut)
2931 fsg(id,iurb)=fprl
2933 return
2934 end subroutine view_factors
2936 ! ===6=8===============================================================72
2937 ! ===6=8===============================================================72
2939 SUBROUTINE fprls (fprl,a,b,c)
2941 implicit none
2945 real a,b,c
2946 real x,y
2947 real fprl
2950 x=a/c
2951 y=b/c
2953 if(a.eq.0.or.b.eq.0.)then
2954 fprl=0.
2955 else
2956 fprl=log( ( (1.+x**2)*(1.+y**2)/(1.+x**2+y**2) )**.5)+ &
2957 y*((1.+x**2)**.5)*atan(y/((1.+x**2)**.5))+ &
2958 x*((1.+y**2)**.5)*atan(x/((1.+y**2)**.5))- &
2959 y*atan(y)-x*atan(x)
2960 fprl=fprl*2./(pi*x*y)
2961 endif
2963 return
2964 end subroutine fprls
2966 ! ===6=8===============================================================72
2967 ! ===6=8===============================================================72
2969 SUBROUTINE fnrms (fnrm,a,b,c)
2971 implicit none
2975 real a,b,c
2976 real x,y,z,a1,a2,a3,a4,a5,a6
2977 real fnrm
2979 x=a/b
2980 y=c/b
2981 z=x**2+y**2
2983 if(y.eq.0.or.x.eq.0)then
2984 fnrm=0.
2985 else
2986 a1=log( (1.+x*x)*(1.+y*y)/(1.+z) )
2987 a2=y*y*log(y*y*(1.+z)/z/(1.+y*y) )
2988 a3=x*x*log(x*x*(1.+z)/z/(1.+x*x) )
2989 a4=y*atan(1./y)
2990 a5=x*atan(1./x)
2991 a6=sqrt(z)*atan(1./sqrt(z))
2992 fnrm=0.25*(a1+a2+a3)+a4+a5-a6
2993 fnrm=fnrm/(pi*y)
2994 endif
2996 return
2997 end subroutine fnrms
2998 ! ===6=8===============================================================72
3000 SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
3001 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
3002 emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
3004 ! initialization routine, where the variables from the table are read
3006 implicit none
3008 integer iurb ! urban class number
3009 ! Building parameters
3010 real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
3011 real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
3012 real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
3013 real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
3014 real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
3015 real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
3016 real twini_u(nurbm) ! Temperature inside the buildings behind the wall [K]
3017 real trini_u(nurbm) ! Temperature inside the buildings behind the roof [K]
3018 real tgini_u(nurbm) ! Initial road temperature
3020 ! Radiation parameters
3021 real albg_u(nurbm) ! Albedo of the ground
3022 real albw_u(nurbm) ! Albedo of the wall
3023 real albr_u(nurbm) ! Albedo of the roof
3024 real emg_u(nurbm) ! Emissivity of ground
3025 real emw_u(nurbm) ! Emissivity of wall
3026 real emr_u(nurbm) ! Emissivity of roof
3028 ! Roughness parameters
3029 real z0g_u(nurbm) ! The ground's roughness length
3030 real z0r_u(nurbm) ! The roof's roughness length
3032 ! Street parameters
3033 integer nd_u(nurbm) ! Number of street direction for each urban class
3035 real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
3036 real drst_u(ndm,nurbm) ! Street direction [degree]
3037 real ws_u(ndm,nurbm) ! Street width [m]
3038 real bs_u(ndm,nurbm) ! Building width [m]
3039 real h_b(nz_um,nurbm) ! Bulding's heights [m]
3040 real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
3042 integer i,iu
3043 integer nurb ! number of urban classes used
3045 !
3046 !Initialize some variables
3047 !
3049 h_b=0.
3050 d_b=0.
3052 nurb=ICATE
3053 do iu=1,nurb
3054 nd_u(iu)=0
3055 enddo
3057 csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
3058 csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
3059 csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
3060 do i=1,icate
3061 alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3062 alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3063 alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3064 enddo
3065 twini_u=TBLEND_TBL
3066 trini_u=TRLEND_TBL
3067 tgini_u=TGLEND_TBL
3068 albw_u=ALBB_TBL
3069 albr_u=ALBR_TBL
3070 albg_u=ALBG_TBL
3071 emw_u=EPSB_TBL
3072 emr_u=EPSR_TBL
3073 emg_u=EPSG_TBL
3074 z0r_u=Z0R_TBL
3075 z0g_u=Z0G_TBL
3076 nd_u=NUMDIR_TBL
3077 do iu=1,icate
3078 if(ndm.lt.nd_u(iu))then
3079 write(*,*)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
3080 write(*,*)'remember also that urban_map_zrd should be equal or greater than nz_um*ndm*nwr-u!'
3081 stop
3082 endif
3083 do i=1,nd_u(iu)
3084 drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
3085 ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
3086 bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
3087 enddo
3088 enddo
3089 do iu=1,ICATE
3090 if(nz_um.lt.numhgt_tbl(iu)+3)then
3091 write(*,*)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
3092 write(*,*)'remember also that urban_map_zrd should be equal or greater than nz_um*ndm*nwr-u!'
3093 stop
3094 endif
3095 do i=1,NUMHGT_TBL(iu)
3096 h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
3097 d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
3098 enddo
3099 enddo
3101 do i=1,ndm
3102 do iu=1,nurbm
3103 strd_u(i,iu)=100000.
3104 enddo
3105 enddo
3107 return
3108 END SUBROUTINE init_para
3109 !==============================================================
3111 !==============================================================
3112 subroutine angle(along,alat,day,realt,zr,deltar,ah)
3113 ! ----------------
3114 !
3115 ! Computation of the solar angles
3116 ! schayes (1982,atm. env. , p1407)
3117 ! Inputs
3118 !========================
3119 ! along=longitud
3120 ! alat=latitude
3121 ! day=julian day (from the beginning of the year)
3122 ! realt= time GMT in hours
3123 ! Outputs
3124 !============================
3125 ! zr=solar zenith angle
3126 ! deltar=declination angle
3127 ! ah=hour angle
3128 !===============================
3130 implicit none
3131 real along,alat, realt, zr, deltar, ah, arg
3132 real rad,om,radh,initt, pii, drad, alongt, cphi, sphi
3133 real c1, c2, c3, s1, s2, s3, delta, rmsr2, cd, sid
3134 real et, ahor, chor, coznt
3135 integer day
3138 data rad,om,radh,initt/0.0174533,0.0172142,0.26179939,0/
3140 zr=0.
3141 deltar=0.
3142 ah=0.
3144 pii = 3.14159265358979312
3145 drad = pii/180.
3147 alongt=along/15.
3148 cphi=cos(alat*drad)
3149 sphi=sin(alat*drad)
3150 !
3151 ! declination
3152 !
3153 arg=om*day
3154 c1=cos(arg)
3155 c2=cos(2.*arg)
3156 c3=cos(3.*arg)
3157 s1=sin(arg)
3158 s2=sin(2.*arg)
3159 s3=sin(3.*arg)
3160 delta=0.33281-22.984*c1-0.3499*c2-0.1398*c3+3.7872*s1+0.03205*s2+0.07187*s3
3161 rmsr2=(1./(1.-0.01673*c1))**2
3162 deltar=delta*rad
3163 cd=cos(deltar)
3164 sid=sin(deltar)
3165 !
3166 ! time equation in hours
3167 !
3168 et=0.0072*c1-0.0528*c2-0.0012*c3-0.1229*s1-0.1565*s2-0.0041*s3
3169 !
3170 !
3171 ! hour angle
3172 !
3174 ! ifh=0
3176 ! ahor=realt-12.+ifh+et+alongt
3177 ahor=realt-12.+et+alongt
3178 ah=ahor*radh
3179 chor=cos(ah)
3180 !
3181 ! zenith angle
3182 !
3183 coznt=sphi*sid+cphi*cd*chor
3185 zr=acos(coznt)
3187 return
3189 END SUBROUTINE angle
3190 !
3191 !====6=8===============================================================72
3192 !====6=8===============================================================72
3194 subroutine upward_rad(ndu,nzu,ws,bs,sigma,pb,ss, &
3195 tg,emg_u,albg_u,rlg,rsg,sfg, &
3196 tw,emw_u,albw_u,rlw,rsw,sfw, &
3197 tr,emr_u,albr_u,rld,rs, sfr, &
3198 rs_abs,rl_up,emiss,grdflx_urb)
3199 !
3200 ! IN this surboutine we compute the upward longwave flux, and the albedo
3201 ! needed for the radiation scheme
3202 !
3203 implicit none
3205 !
3206 !INPUT VARIABLES
3207 !
3208 real rsw(2*ndm,nz_um) ! Short wave radiation at the wall for a given canyon direction [W/m2]
3209 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls for a given canyon direction [W/m2]
3210 real rsg(ndm) ! Short wave radiation at the canyon for a given canyon direction [W/m2]
3211 real rlg(ndm) ! Long wave radiation at the ground for a given canyon direction [W/m2]
3212 real rs ! Short wave radiation at the horizontal surface from the sun [W/m2]
3213 real sfw(2*ndm,nz_um) ! Sensible heat flux from walls [W/m2]
3214 real sfg(ndm) ! Sensible heat flux from ground (road) [W/m2]
3215 real sfr(ndm,nz_um) ! Sensible heat flux from roofs [W/m2]
3216 real rld ! Long wave radiation from the sky [W/m2]
3217 real albg_u ! albedo of the ground/street
3218 real albw_u ! albedo of the walls
3219 real albr_u ! albedo of the roof
3220 real ws(ndm) ! width of the street
3221 real bs(ndm)
3222 ! building size
3223 real pb(nz_um) ! Probability to have a building with an height equal or higher
3224 integer nzu
3225 real ss(nz_um) ! Probability to have a building of a given height
3226 real sigma
3227 real emg_u ! emissivity of the street
3228 real emw_u ! emissivity of the wall
3229 real emr_u ! emissivity of the roof
3230 real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
3231 real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
3232 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
3233 integer id ! street direction
3234 integer ndu ! number of street directions
3235 !OUTPUT/INPUT
3236 real rs_abs ! absrobed solar radiationfor this street direction
3237 real rl_up ! upward longwave radiation for this street direction
3238 real emiss ! mean emissivity
3239 real grdflx_urb ! ground heat flux
3240 !LOCAL
3241 integer iz,iw
3242 real rl_inc,rl_emit
3243 real gfl
3244 integer ix,iy,iwrong
3246 iwrong=1
3247 do iz=1,nzu+1
3248 do id=1,ndu
3249 do iw=1,nwr_u
3250 if(tr(id,iz,iw).lt.100.)then
3251 write(*,*)'in upward_rad ',iz,id,iw,tr(id,iz,iw)
3252 iwrong=0
3253 endif
3254 if(tw(2*id-1,iz,iw).lt.100.) then
3255 write(*,*)'in upward_rad ',iz,id,iw,tw(2*id-1,iz,iw)
3256 iwrong=0
3257 endif
3258 if(tw(2*id,iz,iw).lt.100.) then
3259 write(*,*)'in upward_rad ',iz,id,iw,tw(2*id,iz,iw)
3260 iwrong=0
3261 endif
3262 enddo
3263 enddo
3264 enddo
3265 do id=1,ndu
3266 do iw=1,ng_u
3267 if(tg(id,iw).lt.100.) then
3268 write(*,*)'in upward_rad ',id,iw,tg(id,iw)
3269 iwrong=0
3270 endif
3271 enddo
3272 enddo
3273 if(iwrong.eq.0)stop
3275 rl_up=0.
3277 rs_abs=0.
3278 rl_inc=0.
3279 emiss=0.
3280 rl_emit=0.
3281 grdflx_urb=0.
3282 do id=1,ndu
3283 rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/ndu
3284 rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/ndu
3285 rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/ndu
3286 gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg(id,ng_u)**4.)+sfg(id)
3287 grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/ndu
3289 do iz=2,nzu
3290 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))/ndu
3291 rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
3292 rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
3293 gfl=(1.-albr_u)*rs+emr_u*rld-emr_u*sigma*(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
3294 grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
3295 enddo
3297 do iz=1,nzu
3298 rl_emit=rl_emit-(emw_u*sigma*( tw(2*id-1,iz,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+ &
3299 (1-emw_u)*( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
3300 rl_inc=rl_inc+(( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
3301 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))/ndu
3302 gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) ) &
3303 -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))
3304 grdflx_urb=grdflx_urb-gfl*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
3305 enddo
3307 enddo
3308 emiss=(emg_u+emw_u+emr_u)/3.
3309 rl_up=(rl_inc+rl_emit)-rld
3312 return
3314 END SUBROUTINE upward_rad
3316 !====6=8===============================================================72
3317 !====6=8===============================================================72
3318 ! ===6=8===============================================================72
3319 ! ===6=8===============================================================72
3321 subroutine icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u, &
3322 fws_u,fsg_u,ndu,strd,ws,nzu,z_u)
3324 implicit none
3326 ! Street parameters
3327 integer ndu ! Number of street direction for each urban class
3328 integer iurb
3330 real strd(ndm) ! Street length (fix to greater value to the horizontal length of the cells)
3331 real ws(ndm) ! Street width [m]
3333 ! Grid parameters
3334 integer nzu ! Number of layer in the urban grid
3335 real z_u(nz_um) ! Height of the urban grid levels
3336 ! -----------------------------------------------------------------------
3337 ! Output
3338 !------------------------------------------------------------------------
3340 ! fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
3341 ! and the short wave radation. They are the part of radiation from a surface
3342 ! or from the sky to another surface.
3344 real fww_u(nz_um,nz_um,ndm,nurbm) ! from wall to wall
3345 real fwg_u(nz_um,ndm,nurbm) ! from wall to ground
3346 real fgw_u(nz_um,ndm,nurbm) ! from ground to wall
3347 real fsw_u(nz_um,ndm,nurbm) ! from sky to wall
3348 real fws_u(nz_um,ndm,nurbm) ! from sky to wall
3349 real fsg_u(ndm,nurbm) ! from sky to ground
3351 ! -----------------------------------------------------------------------
3352 ! Local
3353 !------------------------------------------------------------------------
3355 integer id
3357 ! -----------------------------------------------------------------------
3358 ! This routine compute the view factors
3359 !------------------------------------------------------------------------
3360 !
3361 !Initialize
3362 !
3363 fww_u=0.
3364 fwg_u=0.
3365 fgw_u=0.
3366 fsw_u=0.
3367 fws_u=0.
3368 fsg_u=0.
3370 do id=1,ndu
3372 call view_factors(iurb,nzu,id,strd(id),z_u,ws(id), &
3373 fww_u,fwg_u,fgw_u,fsg_u,fsw_u,fws_u)
3375 enddo
3376 return
3377 end subroutine icBEP_XY
3378 ! ===6=8===============================================================72
3379 ! ===6=8===============================================================72
3381 subroutine icBEPHI_XY(hb_u,hi_urb1D,ss_u,pb_u,nzu,z_u)
3383 implicit none
3384 !-----------------------------------------------------------------------
3385 ! Inputs
3386 !-----------------------------------------------------------------------
3387 ! Street parameters
3388 !
3389 real hi_urb1D(nz_um) ! The probability that a building has an height h_b
3390 !
3391 ! Grid parameters
3392 !
3393 real z_u(nz_um) ! Height of the urban grid levels
3394 ! -----------------------------------------------------------------------
3395 ! Output
3396 !------------------------------------------------------------------------
3398 real ss_u(nz_um) ! The probability that a building has an height equal to z
3399 real pb_u(nz_um) ! The probability that a building has an height greater or equal to z
3400 !
3401 ! Grid parameters
3402 !
3403 integer nzu ! Number of layer in the urban grid
3405 ! -----------------------------------------------------------------------
3406 ! Local
3407 !------------------------------------------------------------------------
3408 real hb_u(nz_um) ! Bulding's heights [m]
3409 integer iz_u,id,ilu
3411 real dtot
3412 real hbmax
3414 !------------------------------------------------------------------------
3416 !Initialize variables
3417 !
3419 nzu=0
3420 ss_u=0.
3421 pb_u=0.
3423 ! Normalisation of the building density
3425 dtot=0.
3426 hb_u=0.
3428 do ilu=1,nz_um
3429 dtot=dtot+hi_urb1D(ilu)
3430 enddo
3432 do ilu=1,nz_um
3433 if (hi_urb1D(ilu)<0.) then
3434 ! write(*,*) 'WARNING, HI_URB1D(ilu) < 0 IN BEP'
3435 go to 20
3436 endif
3437 enddo
3439 if (dtot.gt.0.) then
3440 continue
3441 else
3442 ! write(*,*) 'WARNING, HI_URB1D <= 0 IN BEP'
3443 go to 20
3444 endif
3446 do ilu=1,nz_um
3447 hi_urb1D(ilu)=hi_urb1D(ilu)/dtot
3448 enddo
3450 hb_u(1)=dz_u
3451 do ilu=2,nz_um
3452 hb_u(ilu)=dz_u+hb_u(ilu-1)
3453 enddo
3456 ! Compute pb and ss
3459 hbmax=0.
3461 do ilu=1,nz_um
3462 if (hi_urb1D(ilu)>0.and.hi_urb1D(ilu)<=1.) then
3463 hbmax=hb_u(ilu)
3464 endif
3465 enddo
3467 do iz_u=1,nz_um-1
3468 if(z_u(iz_u+1).gt.hbmax)go to 10
3469 enddo
3471 10 continue
3473 nzu=iz_u+1
3475 if ((nzu+1).gt.nz_um) then
3476 write(*,*) 'error, nz_um has to be increased to at least',nzu+1
3477 stop
3478 endif
3480 do iz_u=1,nzu
3481 ss_u(iz_u)=0.
3482 do ilu=1,nz_um
3483 if(z_u(iz_u).le.hb_u(ilu) &
3484 .and.z_u(iz_u+1).gt.hb_u(ilu))then
3485 ss_u(iz_u)=ss_u(iz_u)+hi_urb1D(ilu)
3486 endif
3487 enddo
3488 enddo
3490 pb_u(1)=1.
3491 do iz_u=1,nzu
3492 pb_u(iz_u+1)=max(0.,pb_u(iz_u)-ss_u(iz_u))
3493 enddo
3495 20 continue
3496 return
3497 end subroutine icBEPHI_XY
3498 ! ===6=8===============================================================72
3499 ! ===6=8===============================================================72
3500 END MODULE module_sf_bep
3502 FUNCTION bep_nurbm () RESULT (bep_val_nurbm)
3503 USE module_sf_bep
3504 IMPLICIT NONE
3505 INTEGER :: bep_val_nurbm
3506 bep_val_nurbm = nurbm
3507 END FUNCTION bep_nurbm
3509 FUNCTION bep_ndm () RESULT (bep_val_ndm)
3510 USE module_sf_bep
3511 IMPLICIT NONE
3512 INTEGER :: bep_val_ndm
3513 bep_val_ndm = ndm
3514 END FUNCTION bep_ndm
3516 FUNCTION bep_nz_um () RESULT (bep_val_nz_um)
3517 USE module_sf_bep
3518 IMPLICIT NONE
3519 INTEGER :: bep_val_nz_um
3520 bep_val_nz_um = nz_um
3521 END FUNCTION bep_nz_um
3523 FUNCTION bep_ng_u () RESULT (bep_val_ng_u)
3524 USE module_sf_bep
3525 IMPLICIT NONE
3526 INTEGER :: bep_val_ng_u
3527 bep_val_ng_u = ng_u
3528 END FUNCTION bep_ng_u
3530 FUNCTION bep_nwr_u () RESULT (bep_val_nwr_u)
3531 USE module_sf_bep
3532 IMPLICIT NONE
3533 INTEGER :: bep_val_nwr_u
3534 bep_val_nwr_u = nwr_u
3535 END FUNCTION bep_nwr_u