SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
ehc_module Module Reference

Functions/Subroutines

subroutine ehc (tstep, nlayer, qg_surf, qg_roof, qg_wall, tsfc_roof, tin_roof, temp_in_roof, k_roof, cp_roof, dz_roof, sfr_roof, tsfc_wall, tin_wall, temp_in_wall, k_wall, cp_wall, dz_wall, sfr_wall, tsfc_surf, tin_surf, temp_in_surf, k_surf, cp_surf, dz_surf, sfr_surf, temp_out_roof, qs_roof, temp_out_wall, qs_wall, temp_out_surf, qs_surf, qs)
 

Function/Subroutine Documentation

◆ ehc()

subroutine ehc_module::ehc ( integer, intent(in) tstep,
integer, intent(in) nlayer,
real(kind(1d0)), dimension(nsurf), intent(in) qg_surf,
real(kind(1d0)), dimension(nlayer), intent(in) qg_roof,
real(kind(1d0)), dimension(nlayer), intent(in) qg_wall,
real(kind(1d0)), dimension(nlayer), intent(in) tsfc_roof,
real(kind(1d0)), dimension(nlayer), intent(in) tin_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) temp_in_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) k_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) cp_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) dz_roof,
real(kind(1d0)), dimension(nlayer), intent(in) sfr_roof,
real(kind(1d0)), dimension(nlayer), intent(in) tsfc_wall,
real(kind(1d0)), dimension(nlayer), intent(in) tin_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) temp_in_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) k_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) cp_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in) dz_wall,
real(kind(1d0)), dimension(nlayer), intent(in) sfr_wall,
real(kind(1d0)), dimension(nsurf), intent(in) tsfc_surf,
real(kind(1d0)), dimension(nsurf), intent(in) tin_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in) temp_in_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in) k_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in) cp_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in) dz_surf,
real(kind(1d0)), dimension(nsurf), intent(in) sfr_surf,
real(kind(1d0)), dimension(nlayer, ndepth), intent(out) temp_out_roof,
real(kind(1d0)), dimension(nlayer), intent(out) qs_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(out) temp_out_wall,
real(kind(1d0)), dimension(nlayer), intent(out) qs_wall,
real(kind(1d0)), dimension(nsurf, ndepth), intent(out) temp_out_surf,
real(kind(1d0)), dimension(nsurf), intent(out) qs_surf,
real(kind(1d0)), intent(out) qs )

Definition at line 104 of file suews_phys_ehc.f95.

115 USE allocatearray, ONLY: &
116 nsurf, ndepth, &
118 USE heatflux, ONLY: heatcond1d_vstep
119
120 IMPLICIT NONE
121 INTEGER, INTENT(in) :: tstep
122 INTEGER, INTENT(in) :: nlayer ! number of vertical levels in urban canopy
123
124 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: QG_surf ! ground heat flux
125 ! extended for ESTM_ehc
126
127 ! keys:
128 ! tsfc: surface temperature
129 ! tin: indoor/deep bottom temperature
130 ! temp_in: temperature at inner interfaces
131 ! k: thermal conductivity
132 ! cp: heat capacity
133 ! dz: thickness of each layer
134 ! roof/wall/surf: roof/wall/ground surface types
135
136 ! input arrays: roof facets
137 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: qg_roof
138 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tsfc_roof
139 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tin_roof
140 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: sfr_roof
141 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: temp_in_roof
142 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: k_roof
143 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: cp_roof
144 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: dz_roof
145 ! input arrays: wall facets
146 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: qg_wall
147 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tsfc_wall
148 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tin_wall
149 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: sfr_wall
150 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: temp_in_wall
151 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: k_wall
152 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: cp_wall
153 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: dz_wall
154 ! input arrays: standard suews surfaces
155 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: tsfc_surf
156 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: tin_surf
157 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
158 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: temp_in_surf
159 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: k_surf
160 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: cp_surf
161 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: dz_surf
162
163 ! output arrays
164 ! roof facets
165 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: QS_roof
166 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(out) :: temp_out_roof !interface temperature between depth layers
167 ! wall facets
168 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: QS_wall
169 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(out) :: temp_out_wall !interface temperature between depth layers
170 ! standard suews surfaces
171 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: QS_surf
172 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(out) :: temp_out_surf !interface temperature between depth layers
173
174 ! grid aggregated results
175 REAL(KIND(1D0)), INTENT(out) :: QS
176
177 ! internal use
178 ! temporary arrays in calculations
179 ! note: nsurf is used as the maximum number of surfaces
180 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tsfc_cal
181 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tin_cal
182 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: temp_cal
183 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: k_cal
184 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: cp_cal
185 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: dz_cal
186 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: qs_cal
187
188 ! REAL(KIND(1D0)), DIMENSION(ndepth) :: temp_all_cal
189 ! surface temperatures at innermost depth layer
190 ! REAL(KIND(1D0)) :: temp_IBC
191 INTEGER :: i_facet, i_group, nfacet, i_depth
192
193 ! settings for boundary conditions
194 REAL(KIND(1D0)), DIMENSION(2) :: bc
195
196 ! use temperature as boundary conditions
197 LOGICAL, DIMENSION(2) :: bctype
198 LOGICAL :: debug
199
200 ! use finite depth heat conduction solver
201 LOGICAL :: use_heatcond1d, use_heatcond1d_water
202
203 ! normalised surface fractions
204 REAL(KIND(1D0)), DIMENSION(nlayer) :: sfr_roof_n
205 REAL(KIND(1D0)), DIMENSION(nlayer) :: sfr_wall_n
206
207 ! initialise solver flags
208 use_heatcond1d = .true.
209 use_heatcond1d_water = .false.
210 debug = .false.
211
212 ! normalised surface fractions
213 ! NB: the sum of sfr_roof (sfr_wall) is NOT unity; so need to be normalised by the sum
214 sfr_roof_n = sfr_roof/sum(sfr_roof)
215 sfr_wall_n = sfr_wall/sum(sfr_wall)
216 ! sub-facets of buildings: e.g. walls, roofs, etc.
217 DO i_group = 1, 3
218
219 ! allocate arrays
220 IF (i_group == 1) THEN
221 ! PRINT *, 'group: ', 'roof'
222 nfacet = nlayer
223 ELSE IF (i_group == 2) THEN
224 ! PRINT *, 'group: ', 'wall'
225 nfacet = nlayer
226 ELSE IF (i_group == 3) THEN
227 ! PRINT *, 'group: ', 'surf'
228 nfacet = nsurf
229 END IF
230 ! PRINT *, 'nfacet here: ', nfacet
231 ALLOCATE (tsfc_cal(nfacet))
232 ALLOCATE (tin_cal(nfacet))
233 ALLOCATE (qs_cal(nfacet))
234 ALLOCATE (temp_cal(nfacet, ndepth))
235 ALLOCATE (k_cal(nfacet, ndepth))
236 ALLOCATE (cp_cal(nfacet, ndepth))
237 ALLOCATE (dz_cal(nfacet, ndepth))
238 ! PRINT *, 'allocation done! '
239
240 ! translate input arrays of facet groups to internal use arrays
241 IF (i_group == 1) THEN
242 ! PRINT *, 'translation for roof! '
243 ! TODO: to update with actual values from input files
244 tsfc_cal(1:nfacet) = tsfc_roof(1:nfacet)
245 ! PRINT *, 'tsfc_cal for roof! ', tsfc_cal
246 tin_cal(1:nfacet) = tin_roof(1:nfacet)
247 ! PRINT *, 'tin_cal for roof! ', tin_cal
248 temp_cal(1:nfacet, 1:ndepth) = temp_in_roof(1:nfacet, 1:ndepth)
249 ! PRINT *, 'temp_cal for roof! ',temp_cal
250 k_cal(1:nfacet, 1:ndepth) = k_roof(1:nfacet, 1:ndepth)
251 ! PRINT *, 'k_roof for roof! ',k_roof(1,:)
252 cp_cal(1:nfacet, 1:ndepth) = cp_roof(1:nfacet, 1:ndepth)
253 dz_cal(1:nfacet, 1:ndepth) = dz_roof(1:nfacet, 1:ndepth)
254 ! qs_cal(1:nfacet) = qs_roof(1:nfacet)
255 ELSE IF (i_group == 2) THEN
256 ! PRINT *, 'translation for wall! '
257 ! TODO: to update with actual values from input files
258 tsfc_cal(1:nfacet) = tsfc_wall(1:nfacet)
259 tin_cal(1:nfacet) = tin_wall(1:nfacet)
260 temp_cal(1:nfacet, 1:ndepth) = temp_in_wall(1:nfacet, 1:ndepth)
261 ! TODO: temporarily set this for testing
262 k_cal(1:nfacet, 1:ndepth) = k_wall(1:nfacet, 1:ndepth)
263 cp_cal(1:nfacet, 1:ndepth) = cp_wall(1:nfacet, 1:ndepth)
264 dz_cal(1:nfacet, 1:ndepth) = dz_wall(1:nfacet, 1:ndepth)
265 ! qs_cal(1:nfacet) = qs_wall(1:nfacet)
266
267 ELSE IF (i_group == 3) THEN
268 ! PRINT *, 'translation for surf! '
269 ! nfacet = nsurf
270 tsfc_cal(1:nfacet) = tsfc_surf(1:nfacet)
271 tin_cal(1:nfacet) = tin_surf(1:nfacet)
272 temp_cal(1:nfacet, 1:ndepth) = temp_in_surf(1:nfacet, 1:ndepth)
273 ! TODO: to update with actual values from input files
274 k_cal(1:nfacet, 1:ndepth) = k_surf(1:nfacet, 1:ndepth)
275 cp_cal(1:nfacet, 1:ndepth) = cp_surf(1:nfacet, 1:ndepth)
276 dz_cal(1:nfacet, 1:ndepth) = dz_surf(1:nfacet, 1:ndepth)
277 ! k_cal(1:nfacet, 1:ndepth) = 1.2
278 ! cp_cal(1:nfacet, 1:ndepth) = 2E6
279 ! dz_cal(1:nfacet, 1:ndepth) = 0.1
280 ! qs_cal(1:nfacet) = QS_surf(1:nfacet)
281 END IF
282 ! PRINT *, 'translation done! '
283 ! TODO: temporary setting
284 ! k_cal(1:nfacet, 1:ndepth) = 1.2
285 ! cp_cal(1:nfacet, 1:ndepth) = 2E6
286 ! dz_cal(1:nfacet, 1:ndepth) = 0.2
287
288 ! PRINT *, 'nfacet: ', nfacet
289 DO i_facet = 1, nfacet
290 ! PRINT *, 'i_facet: ', i_facet
291 ! ASSOCIATE (v => dz_cal(i_facet, 1:ndepth))
292 ! PRINT *, 'dz_cal in ESTM_ehc', v, SIZE(v)
293 ! END ASSOCIATE
294
295 ! determine the calculation method
296 IF (i_group == 3) THEN
297 use_heatcond1d = .true.
298 IF (i_facet == bldgsurf) THEN
299 ! building surface needs a different treatment
300 use_heatcond1d = .false.
301 ELSE IF (i_facet == watersurf) THEN
302 ! water surface needs a different treatment
303 use_heatcond1d = .false.
304 use_heatcond1d_water = .true.
305 END IF
306 ELSE
307 use_heatcond1d = .true.
308 END IF
309
310 ! actual heat conduction calculations
311 IF (dz_cal(i_facet, 1) /= -999.0 .AND. use_heatcond1d) THEN
312
313 ! surface heat flux
314 IF (i_group == 1) THEN
315 bc(1) = qg_roof(i_facet)
316 ELSE IF (i_group == 2) THEN
317 bc(1) = qg_wall(i_facet)
318 ELSE IF (i_group == 3) THEN
319 bc(1) = qg_surf(i_facet)
320 END IF
321 ! bctype(1) = .TRUE.
322 bc(1) = tsfc_cal(i_facet)
323 bctype(1) = .false.
324
325 !TODO: should be a prescribed temperature of the innermost boundary
326 bc(2) = tin_cal(i_facet)
327 bctype(2) = .false.
328
329 ! IF (i_group == 3 .AND. i_facet == 3) THEN
330 ! PRINT *, 'temp_cal before: ', temp_cal(i_facet, :)
331 ! PRINT *, 'k_cal: ', k_cal(i_facet, 1:ndepth)
332 ! PRINT *, 'cp_cal: ', cp_cal(i_facet, 1:ndepth)
333 ! PRINT *, 'dz_cal: ', dz_cal(i_facet, 1:ndepth)
334 ! PRINT *, 'bc: ', bc
335
336 ! END IF
337
338 ! 1D heat conduction for finite depth layers
339
340 ! IF ((i_group == 3) .AND. (i_facet == 1)) THEN
341 ! debug = .FALSE.
342 ! ELSE
343 ! debug = .FALSE.
344 ! END IF
345 ! CALL heatcond1d_ext( &
346 CALL heatcond1d_vstep( &
347 temp_cal(i_facet, :), &
348 qs_cal(i_facet), &
349 tsfc_cal(i_facet), &
350 dz_cal(i_facet, 1:ndepth), &
351 tstep*1.d0, &
352 k_cal(i_facet, 1:ndepth), &
353 cp_cal(i_facet, 1:ndepth), &
354 bc, &
355 bctype, debug)
356
357 ! update temperature at all inner interfaces
358 ! tin_cal(i_facet, :) = temp_all_cal
359 ! IF (i_group == 3 .AND. i_facet == 3) THEN
360 ! PRINT *, 'temp_cal after: ', temp_cal(i_facet, :)
361 ! PRINT *, 'QS_cal after: ', QS_cal(i_facet)
362 ! PRINT *, 'k_cal: ', k_cal(i_facet, 1:ndepth)
363 ! PRINT *, 'cp_cal: ', cp_cal(i_facet, 1:ndepth)
364 ! PRINT *, 'dz_cal: ', dz_cal(i_facet, 1:ndepth)
365 ! PRINT *, 'bc: ', bc
366 ! PRINT *, ''
367
368 ! END IF
369 END IF
370
371 IF (dz_cal(i_facet, 1) /= -999.0 .AND. use_heatcond1d_water) THEN
372 ! temperatures at all interfaces, including the outmost surface
373 ! temp_all_cal = temp_cal(i_facet, :)
374
375 ! outermost surface temperature
376 ! bc(1) = tsfc_cal(i_facet)
377 bc(1) = tsfc_cal(i_facet)
378 bctype(1) = .false.
379
380 !TODO: should be a prescribed temperature of the innermost boundary
381 bc(2) = tin_cal(i_facet)
382 bctype(2) = .false.
383
384 ! 1D heat conduction for finite depth layers
385 ! TODO: this should be a water specific heat conduction solver: to implement
386 ! CALL heatcond1d_ext( &
387 CALL heatcond1d_vstep( &
388 temp_cal(i_facet, :), &
389 qs_cal(i_facet), &
390 tsfc_cal(i_facet), &
391 dz_cal(i_facet, 1:ndepth), &
392 tstep*1.d0, &
393 k_cal(i_facet, 1:ndepth), &
394 cp_cal(i_facet, 1:ndepth), &
395 bc, &
396 bctype, debug)
397
398 ! ! update temperature at all inner interfaces
399 ! temp_cal(i_facet, :) = temp_all_cal
400 END IF
401
402 END DO ! end of i_facet loop
403
404 ! translate results to back to the output arrays of facet groups
405 IF (i_group == 1) THEN
406 qs_roof = qs_cal(1:nfacet)
407 ! tsfc_roof = tsfc_cal(1:nfacet)
408 temp_out_roof = temp_cal(:nfacet, :)
409 ELSE IF (i_group == 2) THEN
410 qs_wall = qs_cal(1:nfacet)
411 ! tsfc_wall = tsfc_cal(1:nfacet)
412 temp_out_wall = temp_cal(:nfacet, :)
413 ELSE IF (i_group == 3) THEN
414 qs_surf = qs_cal(1:nfacet)
415 ! tsfc_surf = tsfc_cal(1:nfacet)
416 temp_out_surf = temp_cal(:nfacet, :)
417 END IF
418
419 ! deallocate memory
420 DEALLOCATE (tsfc_cal)
421 DEALLOCATE (tin_cal)
422 DEALLOCATE (qs_cal)
423 DEALLOCATE (temp_cal)
424 DEALLOCATE (k_cal)
425 DEALLOCATE (cp_cal)
426 DEALLOCATE (dz_cal)
427
428 END DO ! end do i_group
429
430 ! aggregated results
431 ! building surface
432 ! PRINT *, 'QS_roof in ESTM_ehc', DOT_PRODUCT(QS_roof, sfr_roof), 'for', sfr_roof
433 ! PRINT *, 'QS_wall in ESTM_ehc', DOT_PRODUCT(QS_wall, sfr_wall), 'for', sfr_wall
434 IF (sfr_surf(bldgsurf) < 1.0e-8) THEN
435 qs_surf(bldgsurf) = 0.0
436 ELSE
437 qs_surf(bldgsurf) = (dot_product(qs_roof, sfr_roof) + dot_product(qs_wall, sfr_wall))/sfr_surf(bldgsurf)
438 END IF
439
440 DO i_depth = 1, ndepth
441 temp_out_surf(bldgsurf, i_depth) = &
442 (dot_product(temp_out_roof(:, i_depth), sfr_roof) &
443 + dot_product(temp_out_wall(:, i_depth), sfr_wall)) &
444 /(sum(sfr_roof) + sum(sfr_wall))
445 END DO
446
447 ! all standard suews surfaces
448 qs = dot_product(qs_surf, sfr_surf)
449
integer, parameter bldgsurf
real(kind(1d0)), dimension(:, :), allocatable dz_wall
integer, parameter conifsurf
real(kind(1d0)), dimension(:), allocatable tsfc_roof
real(kind(1d0)), dimension(:, :), allocatable dz_surf
real(kind(1d0)), dimension(nsurf) sfr_surf
real(kind(1d0)), dimension(:), allocatable sfr_roof
real(kind(1d0)), dimension(:, :), allocatable cp_wall
real(kind(1d0)), dimension(:), allocatable tsfc_surf
real(kind(1d0)), dimension(:, :), allocatable cp_surf
real(kind(1d0)), dimension(:, :), allocatable k_roof
real(kind(1d0)), dimension(:, :), allocatable k_surf
real(kind(1d0)), dimension(:), allocatable tin_surf
integer, parameter watersurf
real(kind(1d0)), dimension(:), allocatable tin_wall
integer, parameter grasssurf
real(kind(1d0)), dimension(:, :), allocatable dz_roof
integer, parameter bsoilsurf
integer, parameter nsurf
integer, parameter pavsurf
real(kind(1d0)), dimension(:, :), allocatable cp_roof
integer, parameter ndepth
real(kind(1d0)), dimension(:, :), allocatable k_wall
real(kind(1d0)), dimension(:), allocatable sfr_wall
real(kind(1d0)), dimension(:), allocatable tin_roof
integer, parameter decidsurf
real(kind(1d0)), dimension(:), allocatable tsfc_wall
subroutine heatcond1d_vstep(t, qs, tsfc, dx, dt, k, rhocp, bc, bctype, debug)

References allocatearray::bldgsurf, allocatearray::bsoilsurf, allocatearray::conifsurf, allocatearray::decidsurf, allocatearray::grasssurf, heatflux::heatcond1d_vstep(), allocatearray::ndepth, allocatearray::nsurf, allocatearray::pavsurf, and allocatearray::watersurf.

Referenced by suews_driver::suews_cal_qs(), and suews_driver::suews_cal_qs_dts().

Here is the call graph for this function:
Here is the caller graph for this function: