Changeset 2696 for palm/trunk/SOURCE/urban_surface_mod.f90
- Timestamp:
- Dec 14, 2017 5:12:51 PM (5 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
-
palm/trunk/SOURCE
-
palm/trunk/SOURCE/urban_surface_mod.f90
r2583 r2696 1 1 !> @file urban_surface_mod.f90 2 !------------------------------------------------------------------------------ --!3 ! This file is part of PALM.2 !------------------------------------------------------------------------------! 3 ! This file is part of the PALM model system. 4 4 ! 5 5 ! PALM is free software: you can redistribute it and/or modify it under the … … 26 26 ! ----------------- 27 27 ! $Id$ 28 ! - Bugfix in calculation of pt_surface and related fluxes. (BM) 29 ! - Do not write surface temperatures onto pt array as this might cause 30 ! problems with nesting. (MS) 31 ! - Revised calculation of pt1 (now done in surface_layer_fluxes). 32 ! Bugfix, f_shf_window and f_shf_green were not set at vertical surface 33 ! elements. (MS) 34 ! - merged with branch ebsolver 35 ! green building surfaces do not evaporate yet 36 ! properties of green wall layers and window layers are taken from wall layers 37 ! this input data is missing. (RvT) 38 ! - Merged with branch radiation (developed by Mohamed Salim) 39 ! - Revised initialization. (MS) 40 ! - Rename emiss_surf into emissivity, roughness_wall into z0, albedo_surf into 41 ! albedo. (MS) 42 ! - Move first call of usm_radiatin from usm_init to init_3d_model 43 ! - fixed problem with near surface temperature 44 ! - added near surface temperature t_surf_10cm_h(m), t_surf_10cm_v(l)%t(m) 45 ! - does not work with temp profile including stability, ol 46 ! t_surf_10cm = pt1 now 47 ! - merged with 2357 bugfix, error message for nopointer version 48 ! - added indoor model coupling with wall heat flux 49 ! - added green substrate/ dry vegetation layer for buildings 50 ! - merged with 2232 new surface-type structure 51 ! - added transmissivity of window tiles 52 ! - added MOSAIK tile approach for 3 different surfaces (RvT) 53 ! 54 ! 2583 2017-10-26 13:58:38Z knoop 28 55 ! Date and time quantities are now read from date_and_time_mod. Solar constant is 29 56 ! read from radiation_model_mod … … 142 169 !> 1.2 Km/s, which seem to be not realistic. 143 170 !> 171 !> @todo Output of _av variables in case of restarts 144 172 !> @todo Revise flux conversion in energy-balance solver 145 173 !> @todo Bugfixing in nopointer branch … … 148 176 !> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog 149 177 !> factor 3 between min and max time 178 !> @todo Move setting of flag indoor_model to indoor_model_mod once available 179 !> @todo Check divisions in wtend (etc.) calculations for possible division 180 !> by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?! 181 !> @todo Use unit 90 for OPEN/CLOSE of input files (FK) 150 182 !------------------------------------------------------------------------------! 151 183 MODULE urban_surface_mod … … 170 202 g, pt_surface, large_scale_forcing, lsf_surf, spinup, & 171 203 spinup_pt_mean, spinup_time, time_do3d, dt_do3d, & 172 average_count_3d, varnamelength, urban_surface 204 average_count_3d, varnamelength, urban_surface, kappa 173 205 174 206 USE cpulog, & 175 207 ONLY: cpu_log, log_point, log_point_s 176 208 177 209 USE date_and_time_mod, & 178 210 ONLY: d_seconds_year, day_of_year_init, time_utc_init 179 211 180 212 USE grid_variables, & 181 213 ONLY: dx, dy, ddx, ddy, ddx2, ddy2 182 214 183 215 USE indices, & 184 216 ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, & … … 192 224 193 225 USE plant_canopy_model_mod, & 194 ONLY: plant_canopy, pch_index, & 195 pc_heating_rate, lad_s 226 ONLY: pc_heating_rate, plant_canopy, usm_lad_rma 196 227 197 228 USE radiation_model_mod, & 198 ONLY: radiation, calc_zenith, zenith,&199 rad_ net, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,&229 ONLY: albedo_type, radiation, calc_zenith, zenith, & 230 rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, & 200 231 sigma_sb, solar_constant, sun_direction, sun_dir_lat, & 201 sun_dir_lon, force_radiation_call 232 sun_dir_lon, & 233 force_radiation_call, surfinsw, surfinlw, surfinswdir, & 234 surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, & 235 surfinl, surfinlwdif, energy_balance_surf_h, & 236 energy_balance_surf_v, rad_sw_in_dir, rad_sw_in_diff, & 237 rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, & 238 surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, pcbl, npcbl, startenergy, & 239 endenergy, nenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, & 240 inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,id, & 241 iz, iy, ix, idir, jdir, kdir, startborder, endborder, nsurf_type, nzub, nzut, & 242 isky, inorth_b,idown_a, isouth_b, ieast_b, iwest_b, nzu, pch, nsurf, & 243 iup_a, inorth_a, isouth_a, ieast_a, iwest_a, idsvf, ndsvf, & 244 idcsf, ndcsf, kdcsf, pct, startland, endland, startwall, endwall 202 245 203 246 USE statistics, & … … 206 249 USE surface_mod 207 250 208 209 251 210 252 IMPLICIT NONE 211 253 254 212 255 !-- configuration parameters (they can be setup in PALM config) 213 LOGICAL :: split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation214 !< (.F. in case the radiation model already does it)215 LOGICAL :: usm_energy_balance_land = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs216 LOGICAL :: usm_energy_balance_wall = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs217 256 LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used 218 257 LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used 219 258 LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls 220 LOGICAL :: mrt_factors = .FALSE. !< whether to generate MRT factor files during init 221 LOGICAL :: write_svf_on_init = .FALSE. 222 LOGICAL :: read_svf_on_init = .FALSE. 223 LOGICAL :: usm_lad_rma = .TRUE. !< use MPI RMA to access LAD for raytracing (instead of global array) 224 225 INTEGER(iwp) :: nrefsteps = 0 !< number of reflection steps to perform 226 259 LOGICAL :: indoor_model = .FALSE. !< whether to use the indoor model 260 261 262 INTEGER(iwp) :: building_type = 1 !< default building type (preleminary setting) 227 263 INTEGER(iwp) :: land_category = 2 !< default category for land surface 228 264 INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone 229 265 INTEGER(iwp) :: pedestrant_category = 2 !< default category for wall surface in pedestrian zone 230 266 INTEGER(iwp) :: roof_category = 2 !< default category for root surface 231 REAL(wp) :: roof_height_limit = 4._wp !< height for distinguish between land surfaces and roofs 232 233 REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha) 234 REAL(wp) :: ra_horiz_coef = 5.0_wp !< mysterious coefficient for correction of overestimation 267 268 ! 269 !-- Indices of input attributes for (above) ground floor level 270 INTEGER(iwp) :: ind_alb_wall = 38 !< index in input list for albedo_type of wall fraction 271 INTEGER(iwp) :: ind_alb_green = 39 !< index in input list for albedo_type of green fraction 272 INTEGER(iwp) :: ind_alb_win = 40 !< index in input list for albedo_type of window fraction 273 INTEGER(iwp) :: ind_emis_wall_agfl = 14 !< index in input list for wall emissivity, above ground floor level 274 INTEGER(iwp) :: ind_emis_wall_gfl = 32 !< index in input list for wall emissivity, ground floor level 275 INTEGER(iwp) :: ind_emis_green_agfl = 15 !< index in input list for green emissivity, above ground floor level 276 INTEGER(iwp) :: ind_emis_green_gfl = 33 !< index in input list for green emissivity, ground floor level 277 INTEGER(iwp) :: ind_emis_win_agfl = 16 !< index in input list for window emissivity, above ground floor level 278 INTEGER(iwp) :: ind_emis_win_gfl = 34 !< index in input list for window emissivity, ground floor level 279 INTEGER(iwp) :: ind_green_frac_w_agfl = 2 !< index in input list for green fraction on wall, above ground floor level 280 INTEGER(iwp) :: ind_green_frac_w_gfl = 23 !< index in input list for green fraction on wall, ground floor level 281 INTEGER(iwp) :: ind_green_frac_r_agfl = 3 !< index in input list for green fraction on roof, above ground floor level 282 INTEGER(iwp) :: ind_green_frac_r_gfl = 24 !< index in input list for green fraction on roof, ground floor level 283 INTEGER(iwp) :: ind_hc1_agfl = 6 !< index in input list for heat capacity at first wall layer, above ground floor level 284 INTEGER(iwp) :: ind_hc1_gfl = 26 !< index in input list for heat capacity at first wall layer, ground floor level 285 INTEGER(iwp) :: ind_hc2_agfl = 7 !< index in input list for heat capacity at second wall layer, above ground floor level 286 INTEGER(iwp) :: ind_hc2_gfl = 27 !< index in input list for heat capacity at second wall layer, ground floor level 287 INTEGER(iwp) :: ind_hc3_agfl = 8 !< index in input list for heat capacity at third wall layer, above ground floor level 288 INTEGER(iwp) :: ind_hc3_gfl = 28 !< index in input list for heat capacity at third wall layer, ground floor level 289 INTEGER(iwp) :: ind_gflh = 20 !< index in input list for ground floor level height 290 INTEGER(iwp) :: ind_lai_r_agfl = 4 !< index in input list for LAI on roof, above ground floor level 291 INTEGER(iwp) :: ind_lai_r_gfl = 4 !< index in input list for LAI on roof, ground floor level 292 INTEGER(iwp) :: ind_lai_w_agfl = 5 !< index in input list for LAI on wall, above ground floor level 293 INTEGER(iwp) :: ind_lai_w_gfl = 25 !< index in input list for LAI on wall, ground floor level 294 INTEGER(iwp) :: ind_tc1_agfl = 9 !< index in input list for thermal conductivity at first wall layer, above ground floor level 295 INTEGER(iwp) :: ind_tc1_gfl = 29 !< index in input list for thermal conductivity at first wall layer, ground floor level 296 INTEGER(iwp) :: ind_tc2_agfl = 10 !< index in input list for thermal conductivity at second wall layer, above ground floor level 297 INTEGER(iwp) :: ind_tc2_gfl = 30 !< index in input list for thermal conductivity at second wall layer, ground floor level 298 INTEGER(iwp) :: ind_tc3_agfl = 11 !< index in input list for thermal conductivity at third wall layer, above ground floor level 299 INTEGER(iwp) :: ind_tc3_gfl = 31 !< index in input list for thermal conductivity at third wall layer, ground floor level 300 INTEGER(iwp) :: ind_thick_1 = 41 !< index for wall layer thickness - 1st layer 301 INTEGER(iwp) :: ind_thick_2 = 42 !< index for wall layer thickness - 2nd layer 302 INTEGER(iwp) :: ind_thick_3 = 43 !< index for wall layer thickness - 3rd layer 303 INTEGER(iwp) :: ind_thick_4 = 44 !< index for wall layer thickness - 4th layer 304 INTEGER(iwp) :: ind_trans_agfl = 17 !< index in input list for window transmissivity, above ground floor level 305 INTEGER(iwp) :: ind_trans_gfl = 35 !< index in input list for window transmissivity, ground floor level 306 INTEGER(iwp) :: ind_wall_frac_agfl = 0 !< index in input list for wall fraction, above ground floor level 307 INTEGER(iwp) :: ind_wall_frac_gfl = 21 !< index in input list for wall fraction, ground floor level 308 INTEGER(iwp) :: ind_win_frac_agfl = 1 !< index in input list for window fraction, above ground floor level 309 INTEGER(iwp) :: ind_win_frac_gfl = 22 !< index in input list for window fraction, ground floor level 310 INTEGER(iwp) :: ind_z0_agfl = 18 !< index in input list for z0, above ground floor level 311 INTEGER(iwp) :: ind_z0_gfl = 36 !< index in input list for z0, ground floor level 312 INTEGER(iwp) :: ind_z0qh_agfl = 19 !< index in input list for z0h / z0q, above ground floor level 313 INTEGER(iwp) :: ind_z0qh_gfl = 37 !< index in input list for z0h / z0q, ground floor level 314 315 316 REAL(wp) :: roof_height_limit = 4._wp !< height for distinguish between land surfaces and roofs 317 REAL(wp) :: ground_floor_level = 4.0_wp !< default ground floor level 318 REAL(wp) :: ra_horiz_coef = 5.0_wp !< mysterious coefficient for correction of overestimation 235 319 !< of r_a for horizontal surfaces -> TODO 236 237 !-- parameters of urban surface model 238 INTEGER(iwp), PARAMETER :: usm_version_len = 10 !< length of identification string of usm version 239 CHARACTER(usm_version_len), PARAMETER :: usm_version = 'USM v. 1.0' !< identification of version of binary svf and restart files 240 INTEGER(iwp), PARAMETER :: svf_code_len = 15 !< length of code for verification of the end of svf file 241 CHARACTER(svf_code_len), PARAMETER :: svf_code = '*** end svf ***' !< code for verification of the end of svf file 242 INTEGER(iwp) :: nzu !< number of layers of urban surface (will be calculated) 243 INTEGER(iwp) :: nzub,nzut !< bottom and top layer of urban surface (will be calculated) 244 INTEGER(iwp), PARAMETER :: nzut_free = 3 !< number of free layers in urban surface layer above top of buildings 245 INTEGER(iwp), PARAMETER :: ndsvf = 2 !< number of dimensions of real values in SVF 246 INTEGER(iwp), PARAMETER :: idsvf = 2 !< number of dimensions of integer values in SVF 247 INTEGER(iwp), PARAMETER :: ndcsf = 2 !< number of dimensions of real values in CSF 248 INTEGER(iwp), PARAMETER :: idcsf = 2 !< number of dimensions of integer values in CSF 249 INTEGER(iwp), PARAMETER :: kdcsf = 4 !< number of dimensions of integer values in CSF calculation array 250 INTEGER(iwp), PARAMETER :: id = 1 !< position of d-index in surfl and surf 251 INTEGER(iwp), PARAMETER :: iz = 2 !< position of k-index in surfl and surf 252 INTEGER(iwp), PARAMETER :: iy = 3 !< position of j-index in surfl and surf 253 INTEGER(iwp), PARAMETER :: ix = 4 !< position of i-index in surfl and surf 254 INTEGER(iwp), PARAMETER :: iroof = 0 !< 0 - index of ground or roof 255 INTEGER(iwp), PARAMETER :: isouth = 1 !< 1 - index of south facing wall 256 INTEGER(iwp), PARAMETER :: inorth = 2 !< 2 - index of north facing wall 257 INTEGER(iwp), PARAMETER :: iwest = 3 !< 3 - index of west facing wall 258 INTEGER(iwp), PARAMETER :: ieast = 4 !< 4 - index of east facing wall 259 INTEGER(iwp), PARAMETER :: isky = 5 !< 5 - index of top border of the urban surface layer ("urban sky") 260 INTEGER(iwp), PARAMETER :: inorthb = 6 !< 6 - index of free north border of the domain (south facing) 261 INTEGER(iwp), PARAMETER :: isouthb = 7 !< 7 - index of north south border of the domain (north facing) 262 INTEGER(iwp), PARAMETER :: ieastb = 8 !< 8 - index of east border of the domain (west facing) 263 INTEGER(iwp), PARAMETER :: iwestb = 9 !< 9 - index of wast border of the domain (east facing) 264 INTEGER(iwp), DIMENSION(0:9), PARAMETER :: idir = (/0,0,0,-1,1,0,0,0,-1,1/) !< surface normal direction x indices 265 INTEGER(iwp), DIMENSION(0:9), PARAMETER :: jdir = (/0,-1,1,0,0,0,-1,1,0,0/) !< surface normal direction y indices 266 INTEGER(iwp), DIMENSION(0:9), PARAMETER :: kdir = (/1,0,0,0,0,-1,0,0,0,0/) !< surface normal direction z indices 267 REAL(wp), DIMENSION(1:4) :: ddxy2 !< 1/dx^2 or 1/dy^2 (in surface normal direction) 268 INTEGER(iwp), DIMENSION(1:4,6:9) :: ijdb !< start and end of the local domain border coordinates (set in code) 269 LOGICAL, DIMENSION(6:9) :: isborder !< is PE on the border of the domain in four corresponding directions 270 !< parameter but set in the code 271 272 !-- indices and sizes of urban surface model 273 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surfl !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x] 274 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surf !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x] 275 INTEGER(iwp) :: nsurfl !< number of all surfaces in local processor 276 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nsurfs !< array of number of all surfaces in individual processors 277 INTEGER(iwp) :: startsky !< start index of block of sky 278 INTEGER(iwp) :: endsky !< end index of block of sky 279 INTEGER(iwp) :: nskys !< number of sky surfaces in local processor 280 INTEGER(iwp) :: startland !< start index of block of land and roof surfaces 281 INTEGER(iwp) :: endland !< end index of block of land and roof surfaces 282 INTEGER(iwp) :: nlands !< number of land and roof surfaces in local processor 283 INTEGER(iwp) :: startwall !< start index of block of wall surfaces 284 INTEGER(iwp) :: endwall !< end index of block of wall surfaces 285 INTEGER(iwp) :: nwalls !< number of wall surfaces in local processor 286 INTEGER(iwp) :: startenergy !< start index of block of real surfaces (land, walls and roofs) 287 INTEGER(iwp) :: endenergy !< end index of block of real surfaces (land, walls and roofs) 288 INTEGER(iwp) :: nenergy !< number of real surfaces in local processor 289 INTEGER(iwp) :: nsurf !< global number of surfaces in index array of surfaces (nsurf = Σproc nsurfs) 290 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surfstart !< starts of blocks of surfaces for individual processors in array surf 291 !< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1] 292 INTEGER(iwp) :: nsvfl !< number of svf for local processor 293 INTEGER(iwp) :: ncsfl !< no. of csf in local processor 294 !< needed only during calc_svf but must be here because it is 295 !< shared between subroutines usm_calc_svf and usm_raytrace 296 297 !-- type for calculation of svf 298 TYPE t_svf 299 INTEGER(iwp) :: isurflt !< 300 INTEGER(iwp) :: isurfs !< 301 REAL(wp) :: rsvf !< 302 REAL(wp) :: rtransp !< 303 END TYPE 304 305 !-- type for calculation of csf 306 TYPE t_csf 307 INTEGER(iwp) :: ip !< 308 INTEGER(iwp) :: itx !< 309 INTEGER(iwp) :: ity !< 310 INTEGER(iwp) :: itz !< 311 INTEGER(iwp) :: isurfs !< 312 REAL(wp) :: rsvf !< 313 REAL(wp) :: rtransp !< 314 END TYPE 320 321 CHARACTER(37), DIMENSION(0:6), PARAMETER :: building_type_name = (/ & 322 'user-defined ', & ! 0 323 'residential - 1950 ', & ! 1 324 'residential 1951 - 2000 ', & ! 2 325 'residential 2001 - ', & ! 3 326 'office - 1950 ', & ! 4 327 'office 1951 - 2000 ', & ! 5 328 'office 2001 - ' & ! 6 329 /) 330 ! 331 !-- building parameters, 4 different types 332 !-- 0 - wall fraction, 1- window fraction, 2 - green fraction on wall, 3- green fraction 333 !-- at roof, 4 - lai of green fraction at roof, 5 - lai of green fraction at wall, 334 !-- 6 - heat capacity of wall layer 1, 7 - heat capacity of wall layer 2, 335 !-- 8 - heat capacity of wall layer 3, 9 - thermal conductivity of wall layer 1, 336 !-- 10 - thermal conductivity of wall layer 2, 11 - thermal conductivity of wall layer 3, 337 !-- 12 - indoor target summer temperature ( K ), 13 - indoor target winter temperature (K), 338 !-- 14 - emissivity of wall fraction, 15 - emissivity of green fraction, 16 - emissivity of window fraction, 339 !-- 17 - transmissivity of window fraction, 18 - z0, 19 - z0h/z0q, 20 - ground floor height, 340 !-- 21 - ground floor wall fraction, 22 - ground floor window fraction, 23 ground floor green fraction, 341 !-- 24 - ground floor green fraction on roof, 25 - ground floor lai of green fraction, 342 !-- 26 - ground floor heat capacity of wall layer 1, 27 - ground floor heat capacity of wall layer 1, 343 !-- 28 - ground floor heat capacity of wall layer 3, 29 - ground floor thermal conductivity of wall layer 1, 344 !-- 30 - ground floor thermal conductivity of wall layer 2, 31 - ground floor thermal conductivity of wall layer 3, 345 !-- 32 - ground floor emissivity of wall fraction, 33 - ground floor emissivity of green fraction, 346 !-- 34 - ground floor emissivity of window fraction, 35 - ground floor transmissivity of window fraction, 347 !-- 36 - ground floor z0, 37 - ground floor z0h/z0q, 38 - albedo type wall fraction 348 !-- 39 - albedo type green fraction, 40 - albedo type window fraction 349 !-- 41 - wall layer thickness - 1st layer, 42 - wall layer thickness - 2nd layer, 350 !-- 43 - wall layer thickness - 3rd layer, 44 - wall layer thickness - 4th layer, 351 !-- 45 - heat capacity of the wall surface, 46 - heat conductivity 352 !-- Please note, only preleminary dummy values so far! 353 REAL(wp), DIMENSION(0:46,1:6), PARAMETER :: building_pars = RESHAPE( (/ & 354 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5 355 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11 356 296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17 357 0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20 358 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25 359 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28 360 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31 361 0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37 362 24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40 363 0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44 364 20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 1 365 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5 366 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11 367 296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17 368 0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20 369 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25 370 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28 371 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31 372 0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37 373 24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40 374 0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44 375 20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 2 376 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5 377 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11 378 296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17 379 0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20 380 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25 381 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28 382 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31 383 0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37 384 24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40 385 0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44 386 20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 3 387 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5 388 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11 389 296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17 390 0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20 391 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25 392 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28 393 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31 394 0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37 395 24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40 396 0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44 397 20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 4 398 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5 399 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11 400 296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17 401 0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20 402 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25 403 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28 404 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31 405 0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37 406 24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40 407 0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44 408 20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 5 409 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5 410 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11 411 296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17 412 0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20 413 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25 414 1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28 415 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31 416 0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37 417 24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40 418 0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44 419 20000.0_wp, 10.0_wp & !parameter 45-46 - end of type 6 420 /), & 421 (/47, 6/) ) 422 315 423 ! 316 424 !-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls. … … 324 432 END TYPE t_wall_vertical 325 433 326 !-- arrays for calculation of svf and csf327 TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array328 TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array329 TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: asvf1, asvf2 !< realizations of svf array330 TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET :: acsf1, acsf2 !< realizations of csf array331 INTEGER(iwp) :: nsvfla !< dimmension of array allocated for storage of svf in local processor332 INTEGER(iwp) :: ncsfla !< dimmension of array allocated for storage of csf in local processor333 INTEGER(iwp) :: msvf, mcsf !< mod for swapping the growing array334 INTEGER(iwp), PARAMETER :: gasize = 10000 !< initial size of growing arrays335 !-- temporary arrays for calculation of csf in raytracing336 INTEGER(iwp) :: maxboxesg !< max number of boxes ray can cross in the domain337 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: boxes !< coordinates of gridboxes being crossed by ray338 REAL(wp), DIMENSION(:), ALLOCATABLE :: crlens !< array of crossing lengths of ray for particular grid boxes339 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored340 #if defined( __parallel )341 INTEGER(kind=MPI_ADDRESS_KIND), &342 DIMENSION(:), ALLOCATABLE :: lad_disp !< array of displaycements of lad in local array of proc lad_ip343 #endif344 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad_s_ray !< array of received lad_s for appropriate gridboxes crossed by ray345 346 !-- arrays storing the values of USM347 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of source and target surface for svf[isvf]348 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: svf !< array of shape view factors+direct irradiation factors for local surfaces349 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins !< array of sw radiation falling to local surface after i-th reflection350 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl !< array of lw radiation for local surface after i-th reflection351 352 !< Inward radiation is also valid for virtual surfaces (radiation leaving domain)353 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< array of sw radiation falling to local surface including radiation from reflections354 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< array of lw radiation falling to local surface including radiation from reflections355 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir !< array of direct sw radiation falling to local surface356 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif !< array of diffuse sw radiation from sky and model boundary falling to local surface357 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif !< array of diffuse lw radiation from sky and model boundary falling to local surface358 359 !< Outward radiation is only valid for nonvirtual surfaces360 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< array of reflected sw radiation for local surface in i-th reflection361 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< array of reflected + emitted lw radiation for local surface in i-th reflection362 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection363 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection364 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection365 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection366 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf !< array of total radiation flux incoming to minus outgoing from local surface367 REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< local copy of rad_net (net radiation at surface)368 434 369 435 !-- arrays for time averages 436 !-- Attention: the variable rad_net_av is also used in the 3d field variable in radiation_model_mod.f90. It may be better to rename it 370 437 REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l 371 438 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections … … 381 448 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection 382 449 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface 450 REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_av !< average of wghf_eb 451 REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb_av !< average of wshf_eb 452 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_wall_av !< Average of t_wall 453 REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_green_av !< average of wghf_eb_green 454 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_green_av !< Average of t_green 455 REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_window_av !< average of wghf_eb_window 456 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_window_av !< Average of t_window 383 457 384 !-- block variables needed for calculation of the plant canopy model inside the urban surface model385 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: csfsurf !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]386 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csf !< array of plant canopy sink fators + direct irradiation factors (transparency)387 !< for local surfaces388 INTEGER(wp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]389 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< index of local pcb[k,j,i]390 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw !< array of absorbed sw radiation for local plant canopy box391 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw !< array of absorbed lw radiation for local plant canopy box392 INTEGER(iwp) :: npcbl !< number of the plant canopy gridboxes in local processor393 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy394 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy395 REAL(wp), DIMENSION(:,:,:), POINTER :: usm_lad !< subset of lad_s within urban surface, transformed to plain Z coordinate396 REAL(wp), DIMENSION(:), POINTER :: usm_lad_g !< usm_lad globalized (used to avoid MPI RMA calls in raytracing)397 REAL(wp) :: prototype_lad !< prototype leaf area density for computing effective optical depth398 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing399 400 !-- radiation related arrays (it should be better in interface of radiation module of PALM401 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_dir !< direct sw radiation402 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_diff !< diffusion sw radiation403 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_diff !< diffusion lw radiation404 458 405 459 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 419 473 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /) 420 474 !< normalized soil, wall and roof layer depths (m/m) 475 ! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.33_wp, 0.66_wp, 1.0_wp /) 476 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /) 477 ! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.33_wp, 0.66_wp, 1.0_wp /) 478 ! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /) 479 !< normalized window layer depths (m/m) 480 ! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /) 481 !< normalized green layer depths (m/m) 482 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /) 483 ! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /) 484 421 485 422 486 REAL(wp) :: wall_inner_temperature = 296.0_wp !< temperature of the inner wall surface (~23 degrees C) (K) 423 487 REAL(wp) :: roof_inner_temperature = 296.0_wp !< temperature of the inner roof surface (~23 degrees C) (K) 424 488 REAL(wp) :: soil_inner_temperature = 283.0_wp !< temperature of the deep soil (~10 degrees C) (K) 489 REAL(wp) :: window_inner_temperature = 296.0_wp !< temperature of the inner window surface (~23 degrees C) (K) 425 490 426 491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 428 493 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 429 494 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m) 495 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_window !< normalized window layer depths (m) 496 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_green !< normalized green layer depths (m) 430 497 431 498 #if defined( __nopointer ) 432 499 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h !< wall surface temperature (K) at horizontal walls 433 500 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_p !< progn. wall surface temperature (K) at horizontal walls 434 501 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h !< window surface temperature (K) at horizontal walls 502 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_p !< progn. window surface temperature (K) at horizontal walls 503 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h !< green surface temperature (K) at horizontal walls 504 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_p !< progn. green surface temperature (K) at horizontal walls 505 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_whole_h !< whole surface temperature (K) at horizontal walls 506 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_whole_h_p !< progn. whole surface temperature (K) at horizontal walls 507 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h !< near surface temperature (10cm) (K) at horizontal walls 508 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h_p !< progn. near surface temperature (10cm) (K) at horizontal walls 435 509 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v 436 510 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_p 511 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v 512 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_p 513 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v 514 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_p 515 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_whole_v 516 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_whole_v_p 517 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v 518 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v_p 437 519 #else 438 520 REAL(wp), DIMENSION(:), POINTER :: t_surf_h 439 521 REAL(wp), DIMENSION(:), POINTER :: t_surf_h_p 522 REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h 523 REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h_p 524 REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h 525 REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h_p 526 REAL(wp), DIMENSION(:), POINTER :: t_surf_whole_h 527 REAL(wp), DIMENSION(:), POINTER :: t_surf_whole_h_p 528 REAL(wp), DIMENSION(:), POINTER :: t_surf_10cm_h 529 REAL(wp), DIMENSION(:), POINTER :: t_surf_10cm_h_p 440 530 441 531 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_1 442 532 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_2 533 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_1 534 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_2 535 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_1 536 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_2 537 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_whole_h_1 538 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_whole_h_2 539 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h_1 540 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h_2 443 541 444 542 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v 445 543 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v_p 544 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v 545 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v_p 546 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v 547 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v_p 548 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_whole_v 549 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_whole_v_p 550 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_10cm_v 551 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_10cm_v_p 446 552 447 553 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_1 448 554 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_2 555 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_1 556 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_2 557 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_1 558 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_2 559 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_whole_v_1 560 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_whole_v_2 561 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v_1 562 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v_2 563 449 564 #endif 450 565 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_av !< average of wall surface temperature (K) 566 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_av !< average of window surface temperature (K) 567 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_av !< average of green wall surface temperature (K) 568 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_whole_av !< average of whole wall surface temperature (K) 569 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_av !< average of whole wall surface temperature (K) 451 570 452 571 !-- Temporal tendencies for time stepping 453 REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K) 572 REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency of wall (K) 573 REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_window_m !< surface temperature tendency of window (K) 574 REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_green_m !< surface temperature tendency of green wall (K) 454 575 455 576 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 457 578 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 458 579 !-- parameters of the land, roof and wall surfaces 459 REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface460 !-- parameters of the wall surfaces461 REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface462 580 463 581 #if defined( __nopointer ) … … 465 583 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av !< Average of t_wall 466 584 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_p !< Prog. wall temperature (K) 585 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h !< Window temperature (K) 586 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_av !< Average of t_window 587 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_p !< Prog. window temperature (K) 588 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h !< Green temperature (K) 589 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_av !< Average of t_green 590 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_p !< Prog. green temperature (K) 467 591 468 592 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v !< Wall temperature (K) 469 593 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av !< Average of t_wall 470 594 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_p !< Prog. wall temperature (K) 595 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v !< Window temperature (K) 596 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_av !< Average of t_window 597 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_p !< Prog. window temperature (K) 598 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v !< Green temperature (K) 599 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_av !< Average of t_green 600 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_p !< Prog. green temperature (K) 471 601 #else 472 602 REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p 473 603 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av, t_wall_h_1, t_wall_h_2 604 REAL(wp), DIMENSION(:,:), POINTER :: t_window_h, t_window_h_p 605 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_av, t_window_h_1, t_window_h_2 606 REAL(wp), DIMENSION(:,:), POINTER :: t_green_h, t_green_h_p 607 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_av, t_green_h_1, t_green_h_2 474 608 475 609 TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p 476 610 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av, t_wall_v_1, t_wall_v_2 611 TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_window_v, t_window_v_p 612 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_av, t_window_v_1, t_window_v_2 613 TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_green_v, t_green_v_p 614 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_av, t_green_v_1, t_green_v_2 477 615 #endif 478 616 479 617 !-- Wall temporal tendencies for time stepping 480 618 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array 619 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_window_m !< t_window prognostic array 620 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_green_m !< t_green prognostic array 481 621 482 622 !-- Surface and material parameters classes (surface_type) … … 495 635 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types 496 636 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: surface_params !< parameters of wall types 497 498 CHARACTER(len=*), PARAMETER :: svf_file_name='usm_svf' 637 499 638 500 639 !-- interfaces of subroutines accessed from outside of this module 640 INTERFACE usm_boundary_condition 641 MODULE PROCEDURE usm_boundary_condition 642 END INTERFACE usm_boundary_condition 643 501 644 INTERFACE usm_check_data_output 502 645 MODULE PROCEDURE usm_check_data_output … … 523 666 END INTERFACE usm_material_heat_model 524 667 668 INTERFACE usm_green_heat_model 669 MODULE PROCEDURE usm_green_heat_model 670 END INTERFACE usm_green_heat_model 671 525 672 INTERFACE usm_parin 526 673 MODULE PROCEDURE usm_parin 527 674 END INTERFACE usm_parin 528 529 INTERFACE usm_radiation530 MODULE PROCEDURE usm_radiation531 END INTERFACE usm_radiation532 675 676 INTERFACE usm_temperature_near_surface 677 MODULE PROCEDURE usm_temperature_near_surface 678 END INTERFACE usm_temperature_near_surface 679 533 680 INTERFACE usm_read_restart_data 534 681 MODULE PROCEDURE usm_read_restart_data … … 546 693 MODULE PROCEDURE usm_write_restart_data 547 694 END INTERFACE usm_write_restart_data 695 696 INTERFACE usm_allocate_surface 697 MODULE PROCEDURE usm_allocate_surface 698 END INTERFACE usm_allocate_surface 699 700 INTERFACE usm_average_3d_data 701 MODULE PROCEDURE usm_average_3d_data 702 END INTERFACE usm_average_3d_data 703 548 704 549 705 SAVE … … 551 707 PRIVATE 552 708 709 !-- Public functions 710 PUBLIC usm_boundary_condition, usm_check_parameters, usm_init_urban_surface,& 711 usm_read_restart_data, & 712 usm_surface_energy_balance, usm_material_heat_model, & 713 usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, & 714 usm_data_output_3d, usm_define_netcdf_grid, usm_parin, & 715 usm_write_restart_data, usm_allocate_surface 716 553 717 !-- Public parameters, constants and initial values 554 PUBLIC split_diffusion_radiation, & 555 usm_anthropogenic_heat, usm_material_model, mrt_factors, & 556 usm_check_parameters, & 557 usm_energy_balance_land, usm_energy_balance_wall, nrefsteps, & 558 usm_init_urban_surface, usm_radiation, usm_read_restart_data, & 559 usm_surface_energy_balance, usm_material_heat_model, & 560 usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, & 561 usm_data_output_3d, usm_define_netcdf_grid, usm_parin, & 562 usm_write_restart_data, & 563 nzub, nzut, ra_horiz_coef, usm_lad_rma, & 564 land_category, pedestrant_category, wall_category, roof_category, & 565 write_svf_on_init, read_svf_on_init 718 PUBLIC usm_anthropogenic_heat, usm_material_model, ra_horiz_coef, & 719 usm_green_heat_model, usm_temperature_near_surface 720 566 721 567 722 568 723 CONTAINS 569 724 570 571 725 !------------------------------------------------------------------------------! 572 726 ! Description: … … 575 729 !> and plant canopy and it allocates the needed arrays for USM 576 730 !------------------------------------------------------------------------------! 577 SUBROUTINE usm_allocate_ urban_surface731 SUBROUTINE usm_allocate_surface 578 732 579 733 IMPLICIT NONE 580 734 581 INTEGER(iwp) :: i, j, k, d, l, ir, jr, ids, m 582 INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i) 583 INTEGER(iwp) :: k_topo2 !< vertical index indicating topography top for given (j,i) 584 INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb 585 INTEGER(iwp) :: procid 586 587 588 589 590 !-- auxiliary vars 591 ddxy2 = (/ddy2,ddy2,ddx2,ddx2/) !< 1/dx^2 or 1/dy^2 (in surface normal direction) 592 593 CALL location_message( '', .TRUE. ) 594 CALL location_message( ' allocation of needed arrays', .TRUE. ) 595 ! 596 !-- Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be 597 !-- removed later). The following contruct finds the lowest / largest index 598 !-- for any upward-facing wall (see bit 12). 599 nzubl = MINVAL( & 600 MAXLOC( & 601 MERGE( 1, 0, & 602 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 603 ), DIM = 1 & 604 ) - 1 & 605 ) 606 nzutl = MAXVAL( & 607 MAXLOC( & 608 MERGE( 1, 0, & 609 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 610 ), DIM = 1 & 611 ) - 1 & 612 ) 613 nzubl = max(nzubl,nzb) 614 615 616 IF ( plant_canopy ) THEN 617 !-- allocate needed arrays 618 ALLOCATE( pct(nys:nyn,nxl:nxr) ) 619 ALLOCATE( pch(nys:nyn,nxl:nxr) ) 620 621 !-- calculate plant canopy height 622 npcbl = 0 623 pct = 0.0_wp 624 pch = 0.0_wp 625 DO i = nxl, nxr 626 DO j = nys, nyn 627 ! 628 !-- Find topography top index 629 k_topo = get_topography_top_index( j, i, 's' ) 630 631 DO k = nzt+1, 0, -1 632 IF ( lad_s(k,j,i) /= 0.0_wp ) THEN 633 !-- we are at the top of the pcs 634 pct(j,i) = k + k_topo 635 pch(j,i) = k 636 npcbl = npcbl + pch(j,i) 637 EXIT 638 ENDIF 639 ENDDO 640 ENDDO 641 ENDDO 642 643 nzutl = max(nzutl, maxval(pct)) 644 !-- code of plant canopy model uses parameter pch_index 645 !-- we need to setup it here to right value 646 !-- (pch_index, lad_s and other arrays in PCM are defined flat) 647 pch_index = maxval(pch) 648 649 prototype_lad = maxval(lad_s) * .9_wp !< better be *1.0 if lad is either 0 or maxval(lad) everywhere 650 IF ( prototype_lad <= 0._wp ) prototype_lad = .3_wp 651 !WRITE(message_string, '(a,f6.3)') 'Precomputing effective box optical ' & 652 ! // 'depth using prototype leaf area density = ', prototype_lad 653 !CALL message('usm_init_urban_surface', 'PA0520', 0, 0, -1, 6, 0) 654 ENDIF 655 656 nzutl = min(nzutl+nzut_free, nzt) 657 658 #if defined( __parallel ) 659 CALL MPI_AllReduce(nzubl,nzub,1,MPI_INTEGER,MPI_MIN,comm2d,ierr); 660 CALL MPI_AllReduce(nzutl,nzut,1,MPI_INTEGER,MPI_MAX,comm2d,ierr); 661 #else 662 nzub = nzubl 663 nzut = nzutl 664 #endif 665 666 !-- global number of urban layers 667 nzu = nzut - nzub + 1 668 669 !-- allocate urban surfaces grid 670 !-- calc number of surfaces in local proc 671 CALL location_message( ' calculation of indices for surfaces', .TRUE. ) 672 nsurfl = 0 673 ! 674 !-- Number of land- and roof surfaces. Note, since horizontal surface elements 675 !-- are already counted in surface_mod, in case be simply reused here. 676 startland = 1 677 nsurfl = surf_usm_h%ns 678 endland = nsurfl 679 nlands = endland-startland+1 680 681 ! 682 !-- Number of vertical surfaces. As vertical surfaces are already 683 !-- counted in surface mod, it can be reused here. 684 startwall = nsurfl+1 685 nsurfl = nsurfl + surf_usm_v(0)%ns + surf_usm_v(1)%ns + & 686 surf_usm_v(2)%ns + surf_usm_v(3)%ns 687 endwall = nsurfl 688 nwalls = endwall-startwall+1 689 690 691 !-- range of energy balance surfaces ! will be treated separately by surf_usm_h and surf_usm_v 692 nenergy = 0 693 IF ( usm_energy_balance_land ) THEN 694 startenergy = startland 695 nenergy = nenergy + nlands 696 ELSE 697 startenergy = startwall 698 ENDIF 699 IF ( usm_energy_balance_wall ) THEN 700 endenergy = endwall 701 nenergy = nenergy + nwalls 702 ELSE 703 endenergy = endland 704 ENDIF 705 706 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 707 !-- block of virtual surfaces 708 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 709 !-- calculate sky surfaces ! not used so far! 710 startsky = nsurfl+1 711 nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1) 712 endsky = nsurfl 713 nskys = endsky-startsky+1 714 715 !-- border flags 716 #if defined( __parallel ) 717 isborder = (/ north_border_pe, south_border_pe, right_border_pe, left_border_pe /) 718 #else 719 isborder = (/.TRUE.,.TRUE.,.TRUE.,.TRUE./) 720 #endif 721 !-- fill array of the limits of the local domain borders 722 ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) ) 723 !-- calulation of the free borders of the domain 724 DO ids = 6,9 725 IF ( isborder(ids) ) THEN 726 !-- free border of the domain in direction ids 727 DO i = ijdb(1,ids), ijdb(2,ids) 728 DO j = ijdb(3,ids), ijdb(4,ids) 729 730 k_topo = get_topography_top_index( j, i, 's' ) 731 k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' ) 732 733 k = nzut - MAX( k_topo, k_topo2 ) 734 nsurfl = nsurfl + k 735 ENDDO 736 ENDDO 737 ENDIF 738 ENDDO 739 740 !-- fill gridpcbl and pcbl 741 IF ( plant_canopy ) THEN 742 ALLOCATE( pcbl(iz:ix, 1:npcbl) ) 743 ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) ) 744 gridpcbl(:,:,:) = 0 745 ipcgb = 0 746 DO i = nxl, nxr 747 DO j = nys, nyn 748 ! 749 !-- Find topography top index 750 k_topo = get_topography_top_index( j, i, 's' ) 751 752 DO k = k_topo + 1, pct(j,i) 753 ipcgb = ipcgb + 1 754 gridpcbl(k,j,i) = ipcgb 755 pcbl(:,ipcgb) = (/ k, j, i /) 756 ENDDO 757 ENDDO 758 ENDDO 759 760 ALLOCATE( pcbinsw( 1:npcbl ) ) 761 ALLOCATE( pcbinlw( 1:npcbl ) ) 762 ENDIF 763 764 !-- fill surfl 765 ALLOCATE(surfl(5,nsurfl)) 766 isurf = 0 767 768 !-- add land surfaces or roofs 769 DO i = nxl, nxr 770 DO j = nys, nyn 771 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 772 k = surf_usm_h%k(m) 773 774 isurf = isurf + 1 775 surfl(:,isurf) = (/iroof,k,j,i,m/) 776 ENDDO 777 ENDDO 778 ENDDO 779 780 !-- add walls 781 DO i = nxl, nxr 782 DO j = nys, nyn 783 l = 0 784 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 785 k = surf_usm_v(l)%k(m) 786 787 isurf = isurf + 1 788 surfl(:,isurf) = (/2,k,j,i,m/) 789 ENDDO 790 l = 1 791 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 792 k = surf_usm_v(l)%k(m) 793 794 isurf = isurf + 1 795 surfl(:,isurf) = (/1,k,j,i,m/) 796 ENDDO 797 l = 2 798 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 799 k = surf_usm_v(l)%k(m) 800 801 isurf = isurf + 1 802 surfl(:,isurf) = (/4,k,j,i,m/) 803 ENDDO 804 l = 3 805 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 806 k = surf_usm_v(l)%k(m) 807 808 isurf = isurf + 1 809 surfl(:,isurf) = (/3,k,j,i,m/) 810 ENDDO 811 ENDDO 812 ENDDO 813 814 !-- add sky 815 DO i = nxl, nxr 816 DO j = nys, nyn 817 isurf = isurf + 1 818 k = nzut 819 surfl(:,isurf) = (/isky,k,j,i,-1/) 820 ENDDO 821 ENDDO 822 823 !-- calulation of the free borders of the domain 824 DO ids = 6,9 825 IF ( isborder(ids) ) THEN 826 !-- free border of the domain in direction ids 827 DO i = ijdb(1,ids), ijdb(2,ids) 828 DO j = ijdb(3,ids), ijdb(4,ids) 829 k_topo = get_topography_top_index( j, i, 's' ) 830 k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' ) 831 832 DO k = MAX(k_topo,k_topo2)+1, nzut 833 isurf = isurf + 1 834 surfl(:,isurf) = (/ids,k,j,i,-1/) 835 ENDDO 836 ENDDO 837 ENDDO 838 ENDIF 839 ENDDO 840 841 !-- global array surf of indices of surfaces and displacement index array surfstart 842 ALLOCATE(nsurfs(0:numprocs-1)) 843 844 #if defined( __parallel ) 845 CALL MPI_Allgather(nsurfl,1,MPI_INTEGER,nsurfs,1,MPI_INTEGER,comm2d,ierr) 846 #else 847 nsurfs(0) = nsurfl 848 #endif 849 ALLOCATE(surfstart(0:numprocs)) 850 k = 0 851 DO i=0,numprocs-1 852 surfstart(i) = k 853 k = k+nsurfs(i) 854 ENDDO 855 surfstart(numprocs) = k 856 nsurf = k 857 ALLOCATE(surf(5,nsurf)) 858 859 #if defined( __parallel ) 860 CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, surfstart*5, MPI_INTEGER, comm2d, ierr) 861 #else 862 surf = surfl 863 #endif 864 865 !-- 866 !-- allocation of the arrays for direct and diffusion radiation 867 CALL location_message( ' allocation of radiation arrays', .TRUE. ) 868 !-- rad_sw_in, rad_lw_in are computed in radiation model, 869 !-- splitting of direct and diffusion part is done 870 !-- in usm_calc_diffusion_radiation for now 871 ALLOCATE( rad_sw_in_dir(nysg:nyng,nxlg:nxrg) ) 872 ALLOCATE( rad_sw_in_diff(nysg:nyng,nxlg:nxrg) ) 873 ALLOCATE( rad_lw_in_diff(nysg:nyng,nxlg:nxrg) ) 874 875 !-- allocate radiation arrays 876 ALLOCATE( surfins(nsurfl) ) 877 ALLOCATE( surfinl(nsurfl) ) 878 ALLOCATE( surfinsw(nsurfl) ) 879 ALLOCATE( surfinlw(nsurfl) ) 880 ALLOCATE( surfinswdir(nsurfl) ) 881 ALLOCATE( surfinswdif(nsurfl) ) 882 ALLOCATE( surfinlwdif(nsurfl) ) 883 ALLOCATE( surfoutsl(startenergy:endenergy) ) 884 ALLOCATE( surfoutll(startenergy:endenergy) ) 885 ALLOCATE( surfoutsw(startenergy:endenergy) ) 886 ALLOCATE( surfoutlw(startenergy:endenergy) ) 887 ALLOCATE( surfouts(nsurf) ) !TODO: global surfaces without virtual 888 ALLOCATE( surfoutl(nsurf) ) !TODO: global surfaces without virtual 889 890 735 INTEGER(iwp) :: l 891 736 892 737 ! … … 896 741 ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) ) 897 742 ! 898 !-- New899 ALLOCATE( surf_usm_h%rad_in_sw(1:surf_usm_h%ns) )900 ALLOCATE( surf_usm_h%rad_out_sw(1:surf_usm_h%ns) )901 ALLOCATE( surf_usm_h%rad_in_lw(1:surf_usm_h%ns) )902 ALLOCATE( surf_usm_h%rad_out_lw(1:surf_usm_h%ns) )903 !904 743 !-- For vertical surfaces 905 744 DO l = 0, 3 906 745 ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) ) 907 746 ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) ) 908 ALLOCATE( surf_usm_v(l)%rad_in_sw(1:surf_usm_v(l)%ns) )909 ALLOCATE( surf_usm_v(l)%rad_out_sw(1:surf_usm_v(l)%ns) )910 ALLOCATE( surf_usm_v(l)%rad_in_lw(1:surf_usm_v(l)%ns) )911 ALLOCATE( surf_usm_v(l)%rad_out_lw(1:surf_usm_v(l)%ns) )912 747 ENDDO 913 748 … … 920 755 ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) ) 921 756 ENDDO 922 923 !-- broadband albedo of the land, roof and wall surface 924 !-- for domain border and sky set artifically to 1.0 925 !-- what allows us to calculate heat flux leaving over 926 !-- side and top borders of the domain 927 ALLOCATE ( albedo_surf(nsurfl) ) 928 albedo_surf = 1.0_wp 929 ALLOCATE ( surf_usm_h%albedo_surf(1:surf_usm_h%ns) ) 757 ! 758 !-- Allocate albedo_type and albedo. Each surface element 759 !-- has 3 values, 0: wall fraction, 1: green fraction, 2: window fraction. 760 ALLOCATE( surf_usm_h%albedo_type(0:2,1:surf_usm_h%ns) ) 761 ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) ) 762 surf_usm_h%albedo_type = albedo_type 930 763 DO l = 0, 3 931 ALLOCATE( surf_usm_v(l)%albedo_surf(1:surf_usm_v(l)%ns) ) 932 ENDDO 764 ALLOCATE( surf_usm_v(l)%albedo_type(0:2,1:surf_usm_v(l)%ns) ) 765 ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) ) 766 surf_usm_v(l)%albedo_type = albedo_type 767 ENDDO 768 769 770 ! 771 !-- Allocate indoor target temperature for summer and winter 772 ALLOCATE( surf_usm_h%target_temp_summer(1:surf_usm_h%ns) ) 773 ALLOCATE( surf_usm_h%target_temp_winter(1:surf_usm_h%ns) ) 774 DO l = 0, 3 775 ALLOCATE( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) ) 776 ALLOCATE( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) ) 777 ENDDO 778 ! 779 !-- Allocate flag indicating ground floor level surface elements 780 ALLOCATE ( surf_usm_h%ground_level(1:surf_usm_h%ns) ) 781 DO l = 0, 3 782 ALLOCATE( surf_usm_v(l)%ground_level(1:surf_usm_v(l)%ns) ) 783 ENDDO 784 ! 785 !-- Allocate arrays for relative surface fraction. 786 !-- 0 - wall fraction, 1 - green fraction, 2 - window fraction 787 ALLOCATE( surf_usm_h%frac(0:2,1:surf_usm_h%ns) ) 788 surf_usm_h%frac = 0.0_wp 789 DO l = 0, 3 790 ALLOCATE( surf_usm_v(l)%frac(0:2,1:surf_usm_v(l)%ns) ) 791 surf_usm_v(l)%frac = 0.0_wp 792 ENDDO 933 793 934 794 !-- wall and roof surface parameters. First for horizontal surfaces 935 ALLOCATE ( emiss_surf(startenergy:endenergy) ) 936 937 ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) ) 938 ALLOCATE ( surf_usm_h%emiss_surf(1:surf_usm_h%ns) ) 939 ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) ) 940 ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) ) 941 ALLOCATE ( surf_usm_h%roughness_wall(1:surf_usm_h%ns) ) 795 ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) ) 796 ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) ) 797 ALLOCATE ( surf_usm_h%lambda_surf_window(1:surf_usm_h%ns) ) 798 ALLOCATE ( surf_usm_h%lambda_surf_green(1:surf_usm_h%ns) ) 799 ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) ) 800 ALLOCATE ( surf_usm_h%c_surface_window(1:surf_usm_h%ns) ) 801 ALLOCATE ( surf_usm_h%c_surface_green(1:surf_usm_h%ns) ) 802 ALLOCATE ( surf_usm_h%transmissivity(1:surf_usm_h%ns) ) 803 ALLOCATE ( surf_usm_h%lai(1:surf_usm_h%ns) ) 804 ALLOCATE ( surf_usm_h%emissivity(0:2,1:surf_usm_h%ns) ) 805 942 806 ! 943 807 !-- For vertical surfaces. 944 808 DO l = 0, 3 945 ALLOCATE ( surf_usm_v(l)%emiss_surf(1:surf_usm_v(l)%ns) ) 946 ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) ) 947 ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) ) 948 ALLOCATE ( surf_usm_v(l)%roughness_wall(1:surf_usm_v(l)%ns) ) 809 ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) ) 810 ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) ) 811 ALLOCATE ( surf_usm_v(l)%lambda_surf_window(1:surf_usm_v(l)%ns) ) 812 ALLOCATE ( surf_usm_v(l)%c_surface_window(1:surf_usm_v(l)%ns) ) 813 ALLOCATE ( surf_usm_v(l)%lambda_surf_green(1:surf_usm_v(l)%ns) ) 814 ALLOCATE ( surf_usm_v(l)%c_surface_green(1:surf_usm_v(l)%ns) ) 815 ALLOCATE ( surf_usm_v(l)%transmissivity(1:surf_usm_v(l)%ns) ) 816 ALLOCATE ( surf_usm_v(l)%lai(1:surf_usm_v(l)%ns) ) 817 818 ALLOCATE ( surf_usm_v(l)%emissivity(0:2,1:surf_usm_v(l)%ns) ) 949 819 ENDDO 950 820 821 ! 951 822 !-- allocate wall and roof material parameters. First for horizontal surfaces 952 823 ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) ) 824 ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns) ) 825 ALLOCATE ( surf_usm_h%thickness_green(1:surf_usm_h%ns) ) 953 826 ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 954 827 ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 828 ALLOCATE ( surf_usm_h%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 829 ALLOCATE ( surf_usm_h%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 830 ALLOCATE ( surf_usm_h%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 831 ALLOCATE ( surf_usm_h%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 832 955 833 ! 956 834 !-- For vertical surfaces. 957 835 DO l = 0, 3 958 836 ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) ) 837 ALLOCATE ( surf_usm_v(l)%thickness_window(1:surf_usm_v(l)%ns) ) 838 ALLOCATE ( surf_usm_v(l)%thickness_green(1:surf_usm_v(l)%ns) ) 959 839 ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 960 840 ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 841 ALLOCATE ( surf_usm_v(l)%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 842 ALLOCATE ( surf_usm_v(l)%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 843 ALLOCATE ( surf_usm_v(l)%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 844 ALLOCATE ( surf_usm_v(l)%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 961 845 ENDDO 962 846 … … 964 848 ALLOCATE ( zwn(nzb_wall:nzt_wall) ) 965 849 ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 850 ALLOCATE ( zwn_window(nzb_wall:nzt_wall) ) 851 ALLOCATE ( surf_usm_h%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 852 ALLOCATE ( zwn_green(nzb_wall:nzt_wall) ) 853 ALLOCATE ( surf_usm_h%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 966 854 ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 967 855 ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 968 856 ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 969 857 ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 858 ALLOCATE ( surf_usm_h%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 859 ALLOCATE ( surf_usm_h%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 860 ALLOCATE ( surf_usm_h%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 861 ALLOCATE ( surf_usm_h%zw_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 862 ALLOCATE ( surf_usm_h%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 863 ALLOCATE ( surf_usm_h%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 864 ALLOCATE ( surf_usm_h%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 865 ALLOCATE ( surf_usm_h%zw_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 970 866 ! 971 867 !-- For vertical surfaces. 972 868 DO l = 0, 3 973 869 ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 870 ALLOCATE ( surf_usm_v(l)%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 871 ALLOCATE ( surf_usm_v(l)%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 974 872 ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 975 873 ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 976 874 ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 977 875 ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 876 ALLOCATE ( surf_usm_v(l)%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 877 ALLOCATE ( surf_usm_v(l)%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 878 ALLOCATE ( surf_usm_v(l)%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 879 ALLOCATE ( surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 880 ALLOCATE ( surf_usm_v(l)%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 881 ALLOCATE ( surf_usm_v(l)%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 882 ALLOCATE ( surf_usm_v(l)%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 883 ALLOCATE ( surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 978 884 ENDDO 979 885 … … 988 894 IF ( .NOT. ALLOCATED( t_wall_h_p ) ) & 989 895 ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 896 IF ( .NOT. ALLOCATED( t_surf_window_h ) ) & 897 ALLOCATE ( t_surf_window_h(1:surf_usm_h%ns) ) 898 IF ( .NOT. ALLOCATED( t_surf_window_h_p ) ) & 899 ALLOCATE ( t_surf_window_h_p(1:surf_usm_h%ns) ) 900 IF ( .NOT. ALLOCATED( t_window_h ) ) & 901 ALLOCATE ( t_window_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 902 IF ( .NOT. ALLOCATED( t_window_h_p ) ) & 903 ALLOCATE ( t_window_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 904 IF ( .NOT. ALLOCATED( t_surf_green_h ) ) & 905 ALLOCATE ( t_surf_green_h(1:surf_usm_h%ns) ) 906 IF ( .NOT. ALLOCATED( t_surf_green_h_p ) ) & 907 ALLOCATE ( t_surf_green_h_p(1:surf_usm_h%ns) ) 908 IF ( .NOT. ALLOCATED( t_green_h ) ) & 909 ALLOCATE ( t_green_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 910 IF ( .NOT. ALLOCATED( t_green_h_p ) ) & 911 ALLOCATE ( t_green_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 912 IF ( .NOT. ALLOCATED( t_surf_whole_h ) ) & 913 ALLOCATE ( t_surf_whole_h(1:surf_usm_h%ns) ) 914 IF ( .NOT. ALLOCATED( t_surf_whole_h_p ) ) & 915 ALLOCATE ( t_surf_whole_h_p(1:surf_usm_h%ns) ) 916 IF ( .NOT. ALLOCATED( t_surf_10cm_h ) ) & 917 ALLOCATE ( t_surf_10cm_h(1:surf_usm_h%ns) ) 918 IF ( .NOT. ALLOCATED( t_surf_10cm_h_p ) ) & 919 ALLOCATE ( t_surf_10cm_h_p(1:surf_usm_h%ns) ) 990 920 #else 991 921 ! … … 1000 930 IF ( .NOT. ALLOCATED( t_wall_h_2 ) ) & 1001 931 ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 932 IF ( .NOT. ALLOCATED( t_surf_window_h_1 ) ) & 933 ALLOCATE ( t_surf_window_h_1(1:surf_usm_h%ns) ) 934 IF ( .NOT. ALLOCATED( t_surf_window_h_2 ) ) & 935 ALLOCATE ( t_surf_window_h_2(1:surf_usm_h%ns) ) 936 IF ( .NOT. ALLOCATED( t_window_h_1 ) ) & 937 ALLOCATE ( t_window_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 938 IF ( .NOT. ALLOCATED( t_window_h_2 ) ) & 939 ALLOCATE ( t_window_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 940 IF ( .NOT. ALLOCATED( t_surf_green_h_1 ) ) & 941 ALLOCATE ( t_surf_green_h_1(1:surf_usm_h%ns) ) 942 IF ( .NOT. ALLOCATED( t_surf_green_h_2 ) ) & 943 ALLOCATE ( t_surf_green_h_2(1:surf_usm_h%ns) ) 944 IF ( .NOT. ALLOCATED( t_green_h_1 ) ) & 945 ALLOCATE ( t_green_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 946 IF ( .NOT. ALLOCATED( t_green_h_2 ) ) & 947 ALLOCATE ( t_green_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 948 IF ( .NOT. ALLOCATED( t_surf_whole_h_1 ) ) & 949 ALLOCATE ( t_surf_whole_h_1(1:surf_usm_h%ns) ) 950 IF ( .NOT. ALLOCATED( t_surf_whole_h_2 ) ) & 951 ALLOCATE ( t_surf_whole_h_2(1:surf_usm_h%ns) ) 952 IF ( .NOT. ALLOCATED( t_surf_10cm_h_1 ) ) & 953 ALLOCATE ( t_surf_10cm_h_1(1:surf_usm_h%ns) ) 954 IF ( .NOT. ALLOCATED( t_surf_10cm_h_2 ) ) & 955 ALLOCATE ( t_surf_10cm_h_2(1:surf_usm_h%ns) ) 1002 956 ! 1003 957 !-- initial assignment of the pointers 1004 958 t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2 959 t_window_h => t_window_h_1; t_window_h_p => t_window_h_2 960 t_green_h => t_green_h_1; t_green_h_p => t_green_h_2 1005 961 t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2 962 t_surf_window_h => t_surf_window_h_1; t_surf_window_h_p => t_surf_window_h_2 963 t_surf_green_h => t_surf_green_h_1; t_surf_green_h_p => t_surf_green_h_2 964 t_surf_whole_h => t_surf_whole_h_1; t_surf_whole_h_p => t_surf_whole_h_2 965 t_surf_10cm_h => t_surf_10cm_h_1; t_surf_10cm_h_p => t_surf_10cm_h_2 966 1006 967 #endif 1007 968 … … 1017 978 IF ( .NOT. ALLOCATED( t_wall_v_p(l)%t ) ) & 1018 979 ALLOCATE ( t_wall_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 980 IF ( .NOT. ALLOCATED( t_surf_window_v(l)%t ) ) & 981 ALLOCATE ( t_surf_window_v(l)%t(1:surf_usm_v(l)%ns) ) 982 IF ( .NOT. ALLOCATED( t_surf_window_v_p(l)%t ) ) & 983 ALLOCATE ( t_surf_window_v_p(l)%t(1:surf_usm_v(l)%ns) ) 984 IF ( .NOT. ALLOCATED( t_window_v(l)%t ) ) & 985 ALLOCATE ( t_window_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 986 IF ( .NOT. ALLOCATED( t_window_v_p(l)%t ) ) & 987 ALLOCATE ( t_window_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 988 IF ( .NOT. ALLOCATED( t_green_v(l)%t ) ) & 989 ALLOCATE ( t_green_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 990 IF ( .NOT. ALLOCATED( t_green_v_p(l)%t ) ) & 991 ALLOCATE ( t_green_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 992 IF ( .NOT. ALLOCATED( t_surf_green_v(l)%t ) ) & 993 ALLOCATE ( t_surf_green_v(l)%t(1:surf_usm_v(l)%ns) ) 994 IF ( .NOT. ALLOCATED( t_surf_green_v_p(l)%t ) ) & 995 ALLOCATE ( t_surf_green_v_p(l)%t(1:surf_usm_v(l)%ns) ) 996 IF ( .NOT. ALLOCATED( t_surf_whole_v(l)%t ) ) & 997 ALLOCATE ( t_surf_whole_v(l)%t(1:surf_usm_v(l)%ns) ) 998 IF ( .NOT. ALLOCATED( t_surf_whole_v_p(l)%t ) ) & 999 ALLOCATE ( t_surf_whole_v_p(l)%t(1:surf_usm_v(l)%ns) ) 1000 IF ( .NOT. ALLOCATED( t_surf_10cm_v(l)%t ) ) & 1001 ALLOCATE ( t_surf_10cm_v(l)%t(1:surf_usm_v(l)%ns) ) 1002 IF ( .NOT. ALLOCATED( t_surf_10cm_v_p(l)%t ) ) & 1003 ALLOCATE ( t_surf_10cm_v_p(l)%t(1:surf_usm_v(l)%ns) ) 1019 1004 ENDDO 1020 1005 #else … … 1031 1016 IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) ) & 1032 1017 ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1018 IF ( .NOT. ALLOCATED( t_surf_window_v_1(l)%t ) ) & 1019 ALLOCATE ( t_surf_window_v_1(l)%t(1:surf_usm_v(l)%ns) ) 1020 IF ( .NOT. ALLOCATED( t_surf_window_v_2(l)%t ) ) & 1021 ALLOCATE ( t_surf_window_v_2(l)%t(1:surf_usm_v(l)%ns) ) 1022 IF ( .NOT. ALLOCATED( t_window_v_1(l)%t ) ) & 1023 ALLOCATE ( t_window_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1024 IF ( .NOT. ALLOCATED( t_window_v_2(l)%t ) ) & 1025 ALLOCATE ( t_window_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1026 IF ( .NOT. ALLOCATED( t_surf_green_v_1(l)%t ) ) & 1027 ALLOCATE ( t_surf_green_v_1(l)%t(1:surf_usm_v(l)%ns) ) 1028 IF ( .NOT. ALLOCATED( t_surf_green_v_2(l)%t ) ) & 1029 ALLOCATE ( t_surf_green_v_2(l)%t(1:surf_usm_v(l)%ns) ) 1030 IF ( .NOT. ALLOCATED( t_green_v_1(l)%t ) ) & 1031 ALLOCATE ( t_green_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1032 IF ( .NOT. ALLOCATED( t_green_v_2(l)%t ) ) & 1033 ALLOCATE ( t_green_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1034 IF ( .NOT. ALLOCATED( t_surf_whole_v_1(l)%t ) ) & 1035 ALLOCATE ( t_surf_whole_v_1(l)%t(1:surf_usm_v(l)%ns) ) 1036 IF ( .NOT. ALLOCATED( t_surf_whole_v_2(l)%t ) ) & 1037 ALLOCATE ( t_surf_whole_v_2(l)%t(1:surf_usm_v(l)%ns) ) 1038 IF ( .NOT. ALLOCATED( t_surf_10cm_v_1(l)%t ) ) & 1039 ALLOCATE ( t_surf_10cm_v_1(l)%t(1:surf_usm_v(l)%ns) ) 1040 IF ( .NOT. ALLOCATED( t_surf_10cm_v_2(l)%t ) ) & 1041 ALLOCATE ( t_surf_10cm_v_2(l)%t(1:surf_usm_v(l)%ns) ) 1033 1042 ENDDO 1034 1043 ! … … 1036 1045 t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2 1037 1046 t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2 1047 t_window_v => t_window_v_1; t_window_v_p => t_window_v_2 1048 t_green_v => t_green_v_1; t_green_v_p => t_green_v_2 1049 t_surf_window_v => t_surf_window_v_1; t_surf_window_v_p => t_surf_window_v_2 1050 t_surf_green_v => t_surf_green_v_1; t_surf_green_v_p => t_surf_green_v_2 1051 t_surf_whole_v => t_surf_whole_v_1; t_surf_whole_v_p => t_surf_whole_v_2 1052 t_surf_10cm_v => t_surf_10cm_v_1; t_surf_10cm_v_p => t_surf_10cm_v_2 1053 1038 1054 #endif 1039 1055 ! … … 1041 1057 ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns) ) 1042 1058 ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 1059 ALLOCATE ( surf_usm_h%tt_surface_window_m(1:surf_usm_h%ns) ) 1060 ALLOCATE ( surf_usm_h%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 1061 ALLOCATE ( surf_usm_h%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 1062 ALLOCATE ( surf_usm_h%tt_surface_green_m(1:surf_usm_h%ns) ) 1063 1043 1064 ! 1044 1065 !-- Set inital values for prognostic quantities 1045 1066 IF ( ALLOCATED( surf_usm_h%tt_surface_m ) ) surf_usm_h%tt_surface_m = 0.0_wp 1046 1067 IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp 1068 IF ( ALLOCATED( surf_usm_h%tt_surface_window_m ) ) surf_usm_h%tt_surface_window_m = 0.0_wp 1069 IF ( ALLOCATED( surf_usm_h%tt_window_m ) ) surf_usm_h%tt_window_m = 0.0_wp 1070 IF ( ALLOCATED( surf_usm_h%tt_green_m ) ) surf_usm_h%tt_green_m = 0.0_wp 1071 IF ( ALLOCATED( surf_usm_h%tt_surface_green_m ) ) surf_usm_h%tt_surface_green_m = 0.0_wp 1047 1072 ! 1048 1073 !-- Now, for vertical surfaces … … 1052 1077 IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) ) surf_usm_v(l)%tt_surface_m = 0.0_wp 1053 1078 IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp 1079 ALLOCATE ( surf_usm_v(l)%tt_surface_window_m(1:surf_usm_v(l)%ns) ) 1080 ALLOCATE ( surf_usm_v(l)%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1081 IF ( ALLOCATED( surf_usm_v(l)%tt_surface_window_m ) ) surf_usm_v(l)%tt_surface_window_m = 0.0_wp 1082 IF ( ALLOCATED( surf_usm_v(l)%tt_window_m ) ) surf_usm_v(l)%tt_window_m = 0.0_wp 1083 ALLOCATE ( surf_usm_v(l)%tt_surface_green_m(1:surf_usm_v(l)%ns) ) 1084 IF ( ALLOCATED( surf_usm_v(l)%tt_surface_green_m ) ) surf_usm_v(l)%tt_surface_green_m = 0.0_wp 1085 ALLOCATE ( surf_usm_v(l)%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1086 IF ( ALLOCATED( surf_usm_v(l)%tt_green_m ) ) surf_usm_v(l)%tt_green_m = 0.0_wp 1054 1087 ENDDO 1055 1088 … … 1058 1091 ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) ) 1059 1092 ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) ) 1093 ALLOCATE ( surf_usm_h%wghf_eb_window(1:surf_usm_h%ns) ) 1094 ALLOCATE ( surf_usm_h%wghf_eb_green(1:surf_usm_h%ns) ) 1095 ALLOCATE ( surf_usm_h%iwghf_eb(1:surf_usm_h%ns) ) 1096 ALLOCATE ( surf_usm_h%iwghf_eb_window(1:surf_usm_h%ns) ) 1060 1097 IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp 1061 1098 IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp 1062 1099 IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp 1100 IF ( ALLOCATED( surf_usm_h%wghf_eb_window ) ) surf_usm_h%wghf_eb_window = 0.0_wp 1101 IF ( ALLOCATED( surf_usm_h%wghf_eb_green ) ) surf_usm_h%wghf_eb_green = 0.0_wp 1102 IF ( ALLOCATED( surf_usm_h%iwghf_eb ) ) surf_usm_h%iwghf_eb = 0.0_wp 1103 IF ( ALLOCATED( surf_usm_h%iwghf_eb_window ) ) surf_usm_h%iwghf_eb_window = 0.0_wp 1063 1104 ! 1064 1105 !-- Now, for vertical surfaces … … 1067 1108 ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) ) 1068 1109 ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) ) 1110 ALLOCATE ( surf_usm_v(l)%wghf_eb_window(1:surf_usm_v(l)%ns) ) 1111 ALLOCATE ( surf_usm_v(l)%wghf_eb_green(1:surf_usm_v(l)%ns) ) 1112 ALLOCATE ( surf_usm_v(l)%iwghf_eb(1:surf_usm_v(l)%ns) ) 1113 ALLOCATE ( surf_usm_v(l)%iwghf_eb_window(1:surf_usm_v(l)%ns) ) 1069 1114 IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp 1070 1115 IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp 1071 1116 IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp 1117 IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_window ) ) surf_usm_v(l)%wghf_eb_window = 0.0_wp 1118 IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_green ) ) surf_usm_v(l)%wghf_eb_green = 0.0_wp 1119 IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb ) ) surf_usm_v(l)%iwghf_eb = 0.0_wp 1120 IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) ) surf_usm_v(l)%iwghf_eb_window = 0.0_wp 1072 1121 ENDDO 1073 1122 1074 END SUBROUTINE usm_allocate_urban_surface 1075 1123 END SUBROUTINE usm_allocate_surface 1076 1124 1077 1125 … … 1115 1163 ELSE 1116 1164 !-- wrong wall layer index 1165 RETURN 1166 ENDIF 1167 ENDIF 1168 IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN 1169 !-- wall layers 1170 READ(var(14:14), '(I1)', iostat=istat ) iwl 1171 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN 1172 var = var(1:12) 1173 ELSE 1174 !-- wrong window layer index 1175 RETURN 1176 ENDIF 1177 ENDIF 1178 IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN 1179 !-- wall layers 1180 READ(var(13:13), '(I1)', iostat=istat ) iwl 1181 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN 1182 var = var(1:11) 1183 ELSE 1184 !-- wrong green layer index 1117 1185 RETURN 1118 1186 ENDIF … … 1264 1332 ENDDO 1265 1333 1334 CASE ( 'usm_wghf_window' ) 1335 !-- array of heat flux from window ground (wall, roof, land) 1336 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN 1337 ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) ) 1338 surf_usm_h%wghf_eb_window_av = 0.0_wp 1339 ENDIF 1340 DO l = 0, 3 1341 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) ) THEN 1342 ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) ) 1343 surf_usm_v(l)%wghf_eb_window_av = 0.0_wp 1344 ENDIF 1345 ENDDO 1346 1347 CASE ( 'usm_wghf_green' ) 1348 !-- array of heat flux from green ground (wall, roof, land) 1349 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN 1350 ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) ) 1351 surf_usm_h%wghf_eb_green_av = 0.0_wp 1352 ENDIF 1353 DO l = 0, 3 1354 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) ) THEN 1355 ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) ) 1356 surf_usm_v(l)%wghf_eb_green_av = 0.0_wp 1357 ENDIF 1358 ENDDO 1359 1360 CASE ( 'usm_iwghf' ) 1361 !-- array of heat flux from indoor ground (wall, roof, land) 1362 IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN 1363 ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) ) 1364 surf_usm_h%iwghf_eb_av = 0.0_wp 1365 ENDIF 1366 DO l = 0, 3 1367 IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_av) ) THEN 1368 ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) ) 1369 surf_usm_v(l)%iwghf_eb_av = 0.0_wp 1370 ENDIF 1371 ENDDO 1372 1373 CASE ( 'usm_iwghf_window' ) 1374 !-- array of heat flux from indoor window ground (wall, roof, land) 1375 IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN 1376 ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) ) 1377 surf_usm_h%iwghf_eb_window_av = 0.0_wp 1378 ENDIF 1379 DO l = 0, 3 1380 IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) ) THEN 1381 ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) ) 1382 surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp 1383 ENDIF 1384 ENDDO 1385 1266 1386 CASE ( 'usm_t_surf' ) 1267 1387 !-- surface temperature for surfaces … … 1277 1397 ENDDO 1278 1398 1399 CASE ( 'usm_t_surf_window' ) 1400 !-- surface temperature for window surfaces 1401 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN 1402 ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) ) 1403 surf_usm_h%t_surf_window_av = 0.0_wp 1404 ENDIF 1405 DO l = 0, 3 1406 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_window_av) ) THEN 1407 ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) ) 1408 surf_usm_v(l)%t_surf_window_av = 0.0_wp 1409 ENDIF 1410 ENDDO 1411 1412 CASE ( 'usm_t_surf_green' ) 1413 !-- surface temperature for green surfaces 1414 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN 1415 ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) ) 1416 surf_usm_h%t_surf_green_av = 0.0_wp 1417 ENDIF 1418 DO l = 0, 3 1419 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_green_av) ) THEN 1420 ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) ) 1421 surf_usm_v(l)%t_surf_green_av = 0.0_wp 1422 ENDIF 1423 ENDDO 1424 1425 CASE ( 'usm_t_surf_whole' ) 1426 !-- surface temperature for whole surfaces 1427 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_whole_av) ) THEN 1428 ALLOCATE( surf_usm_h%t_surf_whole_av(1:surf_usm_h%ns) ) 1429 surf_usm_h%t_surf_whole_av = 0.0_wp 1430 ENDIF 1431 DO l = 0, 3 1432 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_whole_av) ) THEN 1433 ALLOCATE( surf_usm_v(l)%t_surf_whole_av(1:surf_usm_v(l)%ns) ) 1434 surf_usm_v(l)%t_surf_whole_av = 0.0_wp 1435 ENDIF 1436 ENDDO 1437 1438 CASE ( 'usm_t_surf_10cm' ) 1439 !-- near surface temperature for whole surfaces 1440 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_10cm_av) ) THEN 1441 ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) ) 1442 surf_usm_h%t_surf_10cm_av = 0.0_wp 1443 ENDIF 1444 DO l = 0, 3 1445 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) ) THEN 1446 ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) ) 1447 surf_usm_v(l)%t_surf_10cm_av = 0.0_wp 1448 ENDIF 1449 ENDDO 1450 1279 1451 CASE ( 'usm_t_wall' ) 1280 1452 !-- wall temperature for iwl layer of walls and land … … 1290 1462 ENDDO 1291 1463 1464 CASE ( 'usm_t_window' ) 1465 !-- window temperature for iwl layer of walls and land 1466 IF ( .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN 1467 ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 1468 surf_usm_h%t_window_av = 0.0_wp 1469 ENDIF 1470 DO l = 0, 3 1471 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_window_av) ) THEN 1472 ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 1473 surf_usm_v(l)%t_window_av = 0.0_wp 1474 ENDIF 1475 ENDDO 1476 1477 CASE ( 'usm_t_green' ) 1478 !-- green temperature for iwl layer of walls and land 1479 IF ( .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN 1480 ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 1481 surf_usm_h%t_green_av = 0.0_wp 1482 ENDIF 1483 DO l = 0, 3 1484 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_green_av) ) THEN 1485 ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 1486 surf_usm_v(l)%t_green_av = 0.0_wp 1487 ENDIF 1488 ENDDO 1489 1292 1490 CASE DEFAULT 1293 1491 CONTINUE … … 1451 1649 ENDDO 1452 1650 1651 CASE ( 'usm_wghf_window' ) 1652 !-- array of heat flux from window ground (wall, roof, land) 1653 DO m = 1, surf_usm_h%ns 1654 surf_usm_h%wghf_eb_window_av(m) = & 1655 surf_usm_h%wghf_eb_window_av(m) + & 1656 surf_usm_h%wghf_eb_window(m) 1657 ENDDO 1658 DO l = 0, 3 1659 DO m = 1, surf_usm_v(l)%ns 1660 surf_usm_v(l)%wghf_eb_window_av(m) = & 1661 surf_usm_v(l)%wghf_eb_window_av(m) + & 1662 surf_usm_v(l)%wghf_eb_window(m) 1663 ENDDO 1664 ENDDO 1665 1666 CASE ( 'usm_wghf_green' ) 1667 !-- array of heat flux from green ground (wall, roof, land) 1668 DO m = 1, surf_usm_h%ns 1669 surf_usm_h%wghf_eb_green_av(m) = & 1670 surf_usm_h%wghf_eb_green_av(m) + & 1671 surf_usm_h%wghf_eb_green(m) 1672 ENDDO 1673 DO l = 0, 3 1674 DO m = 1, surf_usm_v(l)%ns 1675 surf_usm_v(l)%wghf_eb_green_av(m) = & 1676 surf_usm_v(l)%wghf_eb_green_av(m) + & 1677 surf_usm_v(l)%wghf_eb_green(m) 1678 ENDDO 1679 ENDDO 1680 1681 CASE ( 'usm_iwghf' ) 1682 !-- array of heat flux from indoor ground (wall, roof, land) 1683 DO m = 1, surf_usm_h%ns 1684 surf_usm_h%iwghf_eb_av(m) = & 1685 surf_usm_h%iwghf_eb_av(m) + & 1686 surf_usm_h%iwghf_eb(m) 1687 ENDDO 1688 DO l = 0, 3 1689 DO m = 1, surf_usm_v(l)%ns 1690 surf_usm_v(l)%iwghf_eb_av(m) = & 1691 surf_usm_v(l)%iwghf_eb_av(m) + & 1692 surf_usm_v(l)%iwghf_eb(m) 1693 ENDDO 1694 ENDDO 1695 1696 CASE ( 'usm_iwghf_window' ) 1697 !-- array of heat flux from indoor window ground (wall, roof, land) 1698 DO m = 1, surf_usm_h%ns 1699 surf_usm_h%iwghf_eb_window_av(m) = & 1700 surf_usm_h%iwghf_eb_window_av(m) + & 1701 surf_usm_h%iwghf_eb_window(m) 1702 ENDDO 1703 DO l = 0, 3 1704 DO m = 1, surf_usm_v(l)%ns 1705 surf_usm_v(l)%iwghf_eb_window_av(m) = & 1706 surf_usm_v(l)%iwghf_eb_window_av(m) + & 1707 surf_usm_v(l)%iwghf_eb_window(m) 1708 ENDDO 1709 ENDDO 1710 1453 1711 CASE ( 'usm_t_surf' ) 1454 1712 !-- surface temperature for surfaces … … 1465 1723 ENDDO 1466 1724 ENDDO 1725 1726 CASE ( 'usm_t_surf_window' ) 1727 !-- surface temperature for window surfaces 1728 DO m = 1, surf_usm_h%ns 1729 surf_usm_h%t_surf_window_av(m) = & 1730 surf_usm_h%t_surf_window_av(m) + & 1731 t_surf_window_h(m) 1732 ENDDO 1733 DO l = 0, 3 1734 DO m = 1, surf_usm_v(l)%ns 1735 surf_usm_v(l)%t_surf_window_av(m) = & 1736 surf_usm_v(l)%t_surf_window_av(m) + & 1737 t_surf_window_v(l)%t(m) 1738 ENDDO 1739 ENDDO 1740 1741 CASE ( 'usm_t_surf_green' ) 1742 !-- surface temperature for green surfaces 1743 DO m = 1, surf_usm_h%ns 1744 surf_usm_h%t_surf_green_av(m) = & 1745 surf_usm_h%t_surf_green_av(m) + & 1746 t_surf_green_h(m) 1747 ENDDO 1748 DO l = 0, 3 1749 DO m = 1, surf_usm_v(l)%ns 1750 surf_usm_v(l)%t_surf_green_av(m) = & 1751 surf_usm_v(l)%t_surf_green_av(m) + & 1752 t_surf_green_v(l)%t(m) 1753 ENDDO 1754 ENDDO 1755 1756 CASE ( 'usm_t_surf_whole' ) 1757 !-- surface temperature for whole surfaces 1758 DO m = 1, surf_usm_h%ns 1759 surf_usm_h%t_surf_whole_av(m) = & 1760 surf_usm_h%t_surf_whole_av(m) + & 1761 t_surf_whole_h(m) 1762 ENDDO 1763 DO l = 0, 3 1764 DO m = 1, surf_usm_v(l)%ns 1765 surf_usm_v(l)%t_surf_whole_av(m) = & 1766 surf_usm_v(l)%t_surf_whole_av(m) + & 1767 t_surf_whole_v(l)%t(m) 1768 ENDDO 1769 ENDDO 1770 1771 CASE ( 'usm_t_surf_10cm' ) 1772 !-- near surface temperature for whole surfaces 1773 DO m = 1, surf_usm_h%ns 1774 surf_usm_h%t_surf_10cm_av(m) = & 1775 surf_usm_h%t_surf_10cm_av(m) + & 1776 t_surf_10cm_h(m) 1777 ENDDO 1778 DO l = 0, 3 1779 DO m = 1, surf_usm_v(l)%ns 1780 surf_usm_v(l)%t_surf_10cm_av(m) = & 1781 surf_usm_v(l)%t_surf_10cm_av(m) + & 1782 t_surf_10cm_v(l)%t(m) 1783 ENDDO 1784 ENDDO 1785 1467 1786 1468 1787 CASE ( 'usm_t_wall' ) … … 1481 1800 ENDDO 1482 1801 1802 CASE ( 'usm_t_window' ) 1803 !-- window temperature for iwl layer of walls and land 1804 DO m = 1, surf_usm_h%ns 1805 surf_usm_h%t_window_av(iwl,m) = & 1806 surf_usm_h%t_window_av(iwl,m) + & 1807 t_window_h(iwl,m) 1808 ENDDO 1809 DO l = 0, 3 1810 DO m = 1, surf_usm_v(l)%ns 1811 surf_usm_v(l)%t_window_av(iwl,m) = & 1812 surf_usm_v(l)%t_window_av(iwl,m) + & 1813 t_window_v(l)%t(iwl,m) 1814 ENDDO 1815 ENDDO 1816 1817 CASE ( 'usm_t_green' ) 1818 !-- green temperature for iwl layer of walls and land 1819 DO m = 1, surf_usm_h%ns 1820 surf_usm_h%t_green_av(iwl,m) = & 1821 surf_usm_h%t_green_av(iwl,m) + & 1822 t_green_h(iwl,m) 1823 ENDDO 1824 DO l = 0, 3 1825 DO m = 1, surf_usm_v(l)%ns 1826 surf_usm_v(l)%t_green_av(iwl,m) = & 1827 surf_usm_v(l)%t_green_av(iwl,m) + & 1828 t_green_v(l)%t(iwl,m) 1829 ENDDO 1830 ENDDO 1831 1483 1832 CASE DEFAULT 1484 1833 CONTINUE … … 1638 1987 ENDDO 1639 1988 1989 CASE ( 'usm_wghf_window' ) 1990 !-- array of heat flux from window ground (wall, roof, land) 1991 DO m = 1, surf_usm_h%ns 1992 surf_usm_h%wghf_eb_window_av(m) = & 1993 surf_usm_h%wghf_eb_window_av(m) / & 1994 REAL( average_count_3d, kind=wp ) 1995 ENDDO 1996 DO l = 0, 3 1997 DO m = 1, surf_usm_v(l)%ns 1998 surf_usm_v(l)%wghf_eb_window_av(m) = & 1999 surf_usm_v(l)%wghf_eb_window_av(m) / & 2000 REAL( average_count_3d, kind=wp ) 2001 ENDDO 2002 ENDDO 2003 2004 CASE ( 'usm_wghf_green' ) 2005 !-- array of heat flux from green ground (wall, roof, land) 2006 DO m = 1, surf_usm_h%ns 2007 surf_usm_h%wghf_eb_green_av(m) = & 2008 surf_usm_h%wghf_eb_green_av(m) / & 2009 REAL( average_count_3d, kind=wp ) 2010 ENDDO 2011 DO l = 0, 3 2012 DO m = 1, surf_usm_v(l)%ns 2013 surf_usm_v(l)%wghf_eb_green_av(m) = & 2014 surf_usm_v(l)%wghf_eb_green_av(m) / & 2015 REAL( average_count_3d, kind=wp ) 2016 ENDDO 2017 ENDDO 2018 2019 CASE ( 'usm_iwghf' ) 2020 !-- array of heat flux from indoor ground (wall, roof, land) 2021 DO m = 1, surf_usm_h%ns 2022 surf_usm_h%iwghf_eb_av(m) = & 2023 surf_usm_h%iwghf_eb_av(m) / & 2024 REAL( average_count_3d, kind=wp ) 2025 ENDDO 2026 DO l = 0, 3 2027 DO m = 1, surf_usm_v(l)%ns 2028 surf_usm_v(l)%iwghf_eb_av(m) = & 2029 surf_usm_v(l)%iwghf_eb_av(m) / & 2030 REAL( average_count_3d, kind=wp ) 2031 ENDDO 2032 ENDDO 2033 2034 CASE ( 'usm_iwghf_window' ) 2035 !-- array of heat flux from indoor window ground (wall, roof, land) 2036 DO m = 1, surf_usm_h%ns 2037 surf_usm_h%iwghf_eb_window_av(m) = & 2038 surf_usm_h%iwghf_eb_window_av(m) / & 2039 REAL( average_count_3d, kind=wp ) 2040 ENDDO 2041 DO l = 0, 3 2042 DO m = 1, surf_usm_v(l)%ns 2043 surf_usm_v(l)%iwghf_eb_window_av(m) = & 2044 surf_usm_v(l)%iwghf_eb_window_av(m) / & 2045 REAL( average_count_3d, kind=wp ) 2046 ENDDO 2047 ENDDO 2048 1640 2049 CASE ( 'usm_t_surf' ) 1641 2050 !-- surface temperature for surfaces … … 1649 2058 surf_usm_v(l)%t_surf_av(m) = & 1650 2059 surf_usm_v(l)%t_surf_av(m) / & 2060 REAL( average_count_3d, kind=wp ) 2061 ENDDO 2062 ENDDO 2063 2064 CASE ( 'usm_t_surf_window' ) 2065 !-- surface temperature for window surfaces 2066 DO m = 1, surf_usm_h%ns 2067 surf_usm_h%t_surf_window_av(m) = & 2068 surf_usm_h%t_surf_window_av(m) / & 2069 REAL( average_count_3d, kind=wp ) 2070 ENDDO 2071 DO l = 0, 3 2072 DO m = 1, surf_usm_v(l)%ns 2073 surf_usm_v(l)%t_surf_window_av(m) = & 2074 surf_usm_v(l)%t_surf_window_av(m) / & 2075 REAL( average_count_3d, kind=wp ) 2076 ENDDO 2077 ENDDO 2078 2079 CASE ( 'usm_t_surf_green' ) 2080 !-- surface temperature for green surfaces 2081 DO m = 1, surf_usm_h%ns 2082 surf_usm_h%t_surf_green_av(m) = & 2083 surf_usm_h%t_surf_green_av(m) / & 2084 REAL( average_count_3d, kind=wp ) 2085 ENDDO 2086 DO l = 0, 3 2087 DO m = 1, surf_usm_v(l)%ns 2088 surf_usm_v(l)%t_surf_green_av(m) = & 2089 surf_usm_v(l)%t_surf_green_av(m) / & 2090 REAL( average_count_3d, kind=wp ) 2091 ENDDO 2092 ENDDO 2093 2094 CASE ( 'usm_t_surf_whole' ) 2095 !-- surface temperature for whole surfaces 2096 DO m = 1, surf_usm_h%ns 2097 surf_usm_h%t_surf_whole_av(m) = & 2098 surf_usm_h%t_surf_whole_av(m) / & 2099 REAL( average_count_3d, kind=wp ) 2100 ENDDO 2101 DO l = 0, 3 2102 DO m = 1, surf_usm_v(l)%ns 2103 surf_usm_v(l)%t_surf_whole_av(m) = & 2104 surf_usm_v(l)%t_surf_whole_av(m) / & 2105 REAL( average_count_3d, kind=wp ) 2106 ENDDO 2107 ENDDO 2108 2109 CASE ( 'usm_t_surf_10cm' ) 2110 !-- near surface temperature for whole surfaces 2111 DO m = 1, surf_usm_h%ns 2112 surf_usm_h%t_surf_10cm_av(m) = & 2113 surf_usm_h%t_surf_10cm_av(m) / & 2114 REAL( average_count_3d, kind=wp ) 2115 ENDDO 2116 DO l = 0, 3 2117 DO m = 1, surf_usm_v(l)%ns 2118 surf_usm_v(l)%t_surf_10cm_av(m) = & 2119 surf_usm_v(l)%t_surf_10cm_av(m) / & 1651 2120 REAL( average_count_3d, kind=wp ) 1652 2121 ENDDO … … 1668 2137 ENDDO 1669 2138 2139 CASE ( 'usm_t_window' ) 2140 !-- window temperature for iwl layer of walls and land 2141 DO m = 1, surf_usm_h%ns 2142 surf_usm_h%t_window_av(iwl,m) = & 2143 surf_usm_h%t_window_av(iwl,m) / & 2144 REAL( average_count_3d, kind=wp ) 2145 ENDDO 2146 DO l = 0, 3 2147 DO m = 1, surf_usm_v(l)%ns 2148 surf_usm_v(l)%t_window_av(iwl,m) = & 2149 surf_usm_v(l)%t_window_av(iwl,m) / & 2150 REAL( average_count_3d, kind=wp ) 2151 ENDDO 2152 ENDDO 2153 2154 CASE ( 'usm_t_green' ) 2155 !-- green temperature for iwl layer of walls and land 2156 DO m = 1, surf_usm_h%ns 2157 surf_usm_h%t_green_av(iwl,m) = & 2158 surf_usm_h%t_green_av(iwl,m) / & 2159 REAL( average_count_3d, kind=wp ) 2160 ENDDO 2161 DO l = 0, 3 2162 DO m = 1, surf_usm_v(l)%ns 2163 surf_usm_v(l)%t_green_av(iwl,m) = & 2164 surf_usm_v(l)%t_green_av(iwl,m) / & 2165 REAL( average_count_3d, kind=wp ) 2166 ENDDO 2167 ENDDO 2168 2169 1670 2170 END SELECT 1671 2171 … … 1675 2175 1676 2176 1677 !------------------------------------------------------------------------------! 1678 !> Calculates radiation absorbed by box with given size and LAD. 1679 !> 1680 !> Simulates resol**2 rays (by equally spacing a bounding horizontal square 1681 !> conatining all possible rays that would cross the box) and calculates 1682 !> average transparency per ray. Returns fraction of absorbed radiation flux 1683 !> and area for which this fraction is effective. 1684 !------------------------------------------------------------------------------! 1685 PURE SUBROUTINE usm_box_absorb(boxsize, resol, dens, uvec, area, absorb) 1686 IMPLICIT NONE 1687 1688 REAL(wp), DIMENSION(3), INTENT(in) :: & 1689 boxsize, & !< z, y, x size of box in m 1690 uvec !< z, y, x unit vector of incoming flux 1691 INTEGER(iwp), INTENT(in) :: & 1692 resol !< No. of rays in x and y dimensions 1693 REAL(wp), INTENT(in) :: & 1694 dens !< box density (e.g. Leaf Area Density) 1695 REAL(wp), INTENT(out) :: & 1696 area, & !< horizontal area for flux absorbtion 1697 absorb !< fraction of absorbed flux 1698 REAL(wp) :: & 1699 xshift, yshift, & 1700 xmin, xmax, ymin, ymax, & 1701 xorig, yorig, & 1702 dx1, dy1, dz1, dx2, dy2, dz2, & 1703 crdist, & 1704 transp 1705 INTEGER(iwp) :: & 1706 i, j 1707 1708 xshift = uvec(3) / uvec(1) * boxsize(1) 1709 xmin = min(0._wp, -xshift) 1710 xmax = boxsize(3) + max(0._wp, -xshift) 1711 yshift = uvec(2) / uvec(1) * boxsize(1) 1712 ymin = min(0._wp, -yshift) 1713 ymax = boxsize(2) + max(0._wp, -yshift) 1714 1715 transp = 0._wp 1716 DO i = 1, resol 1717 xorig = xmin + (xmax-xmin) * (i-.5_wp) / resol 1718 DO j = 1, resol 1719 yorig = ymin + (ymax-ymin) * (j-.5_wp) / resol 1720 1721 dz1 = 0._wp 1722 dz2 = boxsize(1)/uvec(1) 1723 1724 IF ( uvec(2) > 0._wp ) THEN 1725 dy1 = -yorig / uvec(2) !< crossing with y=0 1726 dy2 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2) 1727 ELSE IF ( uvec(2) < 0._wp ) THEN 1728 dy1 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2) 1729 dy2 = -yorig / uvec(2) !< crossing with y=0 1730 ELSE !uvec(2)==0 1731 dy1 = -huge(1._wp) 1732 dy2 = huge(1._wp) 1733 ENDIF 1734 1735 IF ( uvec(3) > 0._wp ) THEN 1736 dx1 = -xorig / uvec(3) !< crossing with x=0 1737 dx2 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3) 1738 ELSE IF ( uvec(3) < 0._wp ) THEN 1739 dx1 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3) 1740 dx2 = -xorig / uvec(3) !< crossing with x=0 1741 ELSE !uvec(1)==0 1742 dx1 = -huge(1._wp) 1743 dx2 = huge(1._wp) 1744 ENDIF 1745 1746 crdist = max(0._wp, (min(dz2, dy2, dx2) - max(dz1, dy1, dx1))) 1747 transp = transp + exp(-ext_coef * dens * crdist) 1748 ENDDO 1749 ENDDO 1750 transp = transp / resol**2 1751 area = (boxsize(3)+xshift)*(boxsize(2)+yshift) 1752 absorb = 1._wp - transp 1753 1754 END SUBROUTINE usm_box_absorb 1755 1756 2177 1757 2178 !------------------------------------------------------------------------------! 1758 2179 ! Description: 1759 2180 ! ------------ 1760 !> This subroutine splits direct and diffusion dw radiation 1761 !> It sould not be called in case the radiation model already does it 1762 !> It follows <CITATION> 2181 !> Set internal Neumann boundary condition at outer soil grid points 2182 !> for temperature and humidity. 1763 2183 !------------------------------------------------------------------------------! 1764 SUBROUTINE usm_calc_diffusion_radiation 1765 1766 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 1767 INTEGER(iwp) :: i, j 1768 REAL(wp) :: year_angle !< angle 1769 REAL(wp) :: etr !< extraterestrial radiation 1770 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 1771 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 1772 REAL(wp) :: clearnessIndex !< clearness index 1773 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 1774 1775 1776 !-- Calculate current day and time based on the initial values and simulation time 1777 year_angle = ( (day_of_year_init * 86400) + time_utc_init & 1778 + time_since_reference_point ) * d_seconds_year & 1779 * 2.0_wp * pi 1780 1781 etr = solar_constant * (1.00011_wp + & 1782 0.034221_wp * cos(year_angle) + & 1783 0.001280_wp * sin(year_angle) + & 1784 0.000719_wp * cos(2.0_wp * year_angle) + & 1785 0.000077_wp * sin(2.0_wp * year_angle)) 1786 1787 !-- 1788 !-- Under a very low angle, we keep extraterestrial radiation at 1789 !-- the last small value, therefore the clearness index will be pushed 1790 !-- towards 0 while keeping full continuity. 1791 !-- 1792 IF ( zenith(0) <= lowest_solarUp ) THEN 1793 corrected_solarUp = lowest_solarUp 1794 ELSE 1795 corrected_solarUp = zenith(0) 1796 ENDIF 1797 1798 horizontalETR = etr * corrected_solarUp 1799 1800 DO i = nxlg, nxrg 1801 DO j = nysg, nyng 1802 clearnessIndex = rad_sw_in(0,j,i) / horizontalETR 1803 diff_frac = 1.0_wp / (1.0_wp + exp(-5.0033_wp + 8.6025_wp * clearnessIndex)) 1804 rad_sw_in_diff(j,i) = rad_sw_in(0,j,i) * diff_frac 1805 rad_sw_in_dir(j,i) = rad_sw_in(0,j,i) * (1.0_wp - diff_frac) 1806 rad_lw_in_diff(j,i) = rad_lw_in(0,j,i) 1807 ENDDO 1808 ENDDO 1809 1810 END SUBROUTINE usm_calc_diffusion_radiation 1811 1812 1813 !------------------------------------------------------------------------------! 1814 ! Description: 1815 ! ------------ 1816 !> Calculates shape view factors SVF and plant sink canopy factors PSCF 1817 !> !!!!!DESCRIPTION!!!!!!!!!! 1818 !------------------------------------------------------------------------------! 1819 SUBROUTINE usm_calc_svf 1820 1821 IMPLICIT NONE 1822 1823 INTEGER(iwp) :: i, j, k, l, d, ip, jp 1824 INTEGER(iwp) :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrtt, imrtf 1825 INTEGER(iwp) :: sd, td, ioln, iproc 1826 REAL(wp), DIMENSION(0:9) :: facearea 1827 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzterrl, planthl 1828 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csflt, pcsflt 1829 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: kcsflt,kpcsflt 1830 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: icsflt,dcsflt,ipcsflt,dpcsflt 1831 REAL(wp), DIMENSION(3) :: uv 1832 LOGICAL :: visible 1833 REAL(wp), DIMENSION(3) :: sa, ta !< real coordinates z,y,x of source and target 1834 REAL(wp) :: transparency, rirrf, sqdist, svfsum 1835 INTEGER(iwp) :: isurflt, isurfs, isurflt_prev 1836 INTEGER(iwp) :: itx, ity, itz 1837 CHARACTER(len=7) :: pid_char = '' 1838 INTEGER(iwp) :: win_lad, minfo 1839 REAL(wp), DIMENSION(:,:,:), POINTER :: lad_s_rma !< fortran pointer, but lower bounds are 1 1840 TYPE(c_ptr) :: lad_s_rma_p !< allocated c pointer 1841 #if defined( __parallel ) 1842 INTEGER(kind=MPI_ADDRESS_KIND) :: size_lad_rma 1843 #endif 1844 ! 1845 !-- calculation of the SVF 1846 CALL location_message( ' calculation of SVF and CSF', .TRUE. ) 1847 ! 1848 !-- precalculate face areas for different face directions using normal vector 1849 DO d = 0, 9 1850 facearea(d) = 1._wp 1851 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 1852 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 1853 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 1854 ENDDO 1855 1856 !-- initialize variables and temporary arrays for calculation of svf and csf 1857 nsvfl = 0 1858 ncsfl = 0 1859 nsvfla = gasize 1860 msvf = 1 1861 ALLOCATE( asvf1(nsvfla) ) 1862 asvf => asvf1 1863 IF ( plant_canopy ) THEN 1864 ncsfla = gasize 1865 mcsf = 1 1866 ALLOCATE( acsf1(ncsfla) ) 1867 acsf => acsf1 1868 ENDIF 1869 1870 !-- initialize temporary terrain and plant canopy height arrays (global 2D array!) 1871 ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) ) 1872 #if defined( __parallel ) 1873 ALLOCATE( nzterrl(nys:nyn,nxl:nxr) ) 1874 nzterrl = MAXLOC( & 1875 MERGE( 1, 0, & 1876 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 1877 ), DIM = 1 & 1878 ) - 1 ! = nzb_s_inner(nys:nyn,nxl:nxr) 1879 CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, & 1880 nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr ) 1881 DEALLOCATE(nzterrl) 1882 #else 1883 nzterr = RESHAPE( MAXLOC( & 1884 MERGE( 1, 0, & 1885 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 1886 ), DIM = 1 & 1887 ) - 1, & 1888 (/(nx+1)*(ny+1)/) & 1889 ) 1890 #endif 1891 IF ( plant_canopy ) THEN 1892 ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) ) 1893 maxboxesg = nx + ny + nzu + 1 1894 !-- temporary arrays storing values for csf calculation during raytracing 1895 ALLOCATE( boxes(3, maxboxesg) ) 1896 ALLOCATE( crlens(maxboxesg) ) 1897 1898 #if defined( __parallel ) 1899 ALLOCATE( planthl(nys:nyn,nxl:nxr) ) 1900 planthl = pch(nys:nyn,nxl:nxr) 1901 1902 CALL MPI_AllGather( planthl, nnx*nny, MPI_INTEGER, & 1903 plantt, nnx*nny, MPI_INTEGER, comm2d, ierr ) 1904 DEALLOCATE( planthl ) 1905 1906 !-- temporary arrays storing values for csf calculation during raytracing 1907 ALLOCATE( lad_ip(maxboxesg) ) 1908 ALLOCATE( lad_disp(maxboxesg) ) 1909 1910 IF ( usm_lad_rma ) THEN 1911 ALLOCATE( lad_s_ray(maxboxesg) ) 1912 1913 ! set conditions for RMA communication 1914 CALL MPI_Info_create(minfo, ierr) 1915 CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr) 1916 CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr) 1917 CALL MPI_Info_set(minfo, 'same_size', 'true', ierr) 1918 CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr) 1919 1920 !-- Allocate and initialize the MPI RMA window 1921 !-- must be in accordance with allocation of lad_s in plant_canopy_model 1922 !-- optimization of memory should be done 1923 !-- Argument X of function c_sizeof(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now 1924 size_lad_rma = c_sizeof(1.0_wp)*nnx*nny*nzu 1925 CALL MPI_Win_allocate(size_lad_rma, c_sizeof(1.0_wp), minfo, comm2d, & 1926 lad_s_rma_p, win_lad, ierr) 1927 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzu, nny, nnx /)) 1928 usm_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:) 1929 ELSE 1930 ALLOCATE(usm_lad(nzub:nzut, nys:nyn, nxl:nxr)) 1931 ENDIF 1932 #else 1933 plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) ) 1934 ALLOCATE(usm_lad(nzub:nzut, nys:nyn, nxl:nxr)) 1935 #endif 1936 usm_lad(:,:,:) = 0._wp 1937 DO i = nxl, nxr 1938 DO j = nys, nyn 1939 k = get_topography_top_index( j, i, 's' ) 1940 1941 usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i) 1942 ENDDO 1943 ENDDO 1944 1945 #if defined( __parallel ) 1946 IF ( usm_lad_rma ) THEN 1947 CALL MPI_Info_free(minfo, ierr) 1948 CALL MPI_Win_lock_all(0, win_lad, ierr) 1949 ELSE 1950 ALLOCATE( usm_lad_g(0:(nx+1)*(ny+1)*nzu-1) ) 1951 CALL MPI_AllGather( usm_lad, nnx*nny*nzu, MPI_REAL, & 1952 usm_lad_g, nnx*nny*nzu, MPI_REAL, comm2d, ierr ) 1953 ENDIF 1954 #endif 1955 ENDIF 1956 1957 IF ( mrt_factors ) THEN 1958 OPEN(153, file='MRT_TARGETS', access='SEQUENTIAL', & 1959 action='READ', status='OLD', form='FORMATTED', err=524) 1960 OPEN(154, file='MRT_FACTORS'//myid_char, access='DIRECT', recl=(5*4+2*8), & 1961 action='WRITE', status='REPLACE', form='UNFORMATTED', err=525) 1962 imrtf = 1 1963 DO 1964 READ(153, *, end=526, err=524) imrtt, i, j, k 1965 IF ( i < nxl .OR. i > nxr & 1966 .OR. j < nys .OR. j > nyn ) CYCLE 1967 ta = (/ REAL(k), REAL(j), REAL(i) /) 1968 1969 DO isurfs = 1, nsurf 1970 IF ( .NOT. usm_facing(i, j, k, -1, & 1971 surf(ix, isurfs), surf(iy, isurfs), & 1972 surf(iz, isurfs), surf(id, isurfs)) ) THEN 1973 CYCLE 1974 ENDIF 1975 1976 sd = surf(id, isurfs) 1977 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 1978 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 1979 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 1980 1981 !-- unit vector source -> target 1982 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 1983 sqdist = SUM(uv(:)**2) 1984 uv = uv / SQRT(sqdist) 1985 1986 !-- irradiance factor - see svf. Here we consider that target face is always normal, 1987 !-- i.e. the second dot product equals 1 1988 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & 1989 / (pi * sqdist) * facearea(sd) 1990 1991 !-- raytrace while not creating any canopy sink factors 1992 CALL usm_raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., & 1993 visible, transparency, win_lad) 1994 IF ( .NOT. visible ) CYCLE 1995 1996 !rsvf = rirrf * transparency 1997 WRITE(154, rec=imrtf, err=525) INT(imrtt, kind=4), & 1998 INT(surf(id, isurfs), kind=4), & 1999 INT(surf(iz, isurfs), kind=4), & 2000 INT(surf(iy, isurfs), kind=4), & 2001 INT(surf(ix, isurfs), kind=4), & 2002 REAL(rirrf, kind=8), REAL(transparency, kind=8) 2003 imrtf = imrtf + 1 2004 2005 ENDDO !< isurfs 2006 ENDDO !< MRT_TARGETS record 2007 2008 524 message_string = 'error reading file MRT_TARGETS' 2009 CALL message( 'usm_calc_svf', 'PA0524', 1, 2, 0, 6, 0 ) 2010 2011 525 message_string = 'error writing file MRT_FACTORS'//myid_char 2012 CALL message( 'usm_calc_svf', 'PA0525', 1, 2, 0, 6, 0 ) 2013 2014 526 CLOSE(153) 2015 CLOSE(154) 2016 ENDIF !< mrt_factors 2017 2018 2019 DO isurflt = 1, nsurfl 2020 !-- determine face centers 2021 td = surfl(id, isurflt) 2022 IF ( td >= isky .AND. .NOT. plant_canopy ) CYCLE 2023 ta = (/ REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td), & 2024 REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td), & 2025 REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td) /) 2026 DO isurfs = 1, nsurf 2027 IF ( .NOT. usm_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 2028 surfl(iz, isurflt), surfl(id, isurflt), & 2029 surf(ix, isurfs), surf(iy, isurfs), & 2030 surf(iz, isurfs), surf(id, isurfs)) ) THEN 2031 CYCLE 2032 ENDIF 2033 2034 sd = surf(id, isurfs) 2035 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 2036 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 2037 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 2038 2039 !-- unit vector source -> target 2040 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 2041 sqdist = SUM(uv(:)**2) 2042 uv = uv / SQRT(sqdist) 2043 2044 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 2045 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 2046 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 2047 / (pi * sqdist) & ! square of distance between centers 2048 * facearea(sd) 2049 2050 !-- raytrace + process plant canopy sinks within 2051 CALL usm_raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., & 2052 visible, transparency, win_lad) 2053 2054 IF ( .NOT. visible ) CYCLE 2055 IF ( td >= isky ) CYCLE !< we calculated these only for raytracing 2056 !< to find plant canopy sinks, we don't need svf for them 2057 ! rsvf = rirrf * transparency 2058 2059 !-- write to the svf array 2060 nsvfl = nsvfl + 1 2061 !-- check dimmension of asvf array and enlarge it if needed 2062 IF ( nsvfla < nsvfl ) THEN 2063 k = nsvfla * 2 2064 IF ( msvf == 0 ) THEN 2065 msvf = 1 2066 ALLOCATE( asvf1(k) ) 2067 asvf => asvf1 2068 asvf1(1:nsvfla) = asvf2 2069 DEALLOCATE( asvf2 ) 2070 ELSE 2071 msvf = 0 2072 ALLOCATE( asvf2(k) ) 2073 asvf => asvf2 2074 asvf2(1:nsvfla) = asvf1 2075 DEALLOCATE( asvf1 ) 2076 ENDIF 2077 nsvfla = k 2078 ENDIF 2079 !-- write svf values into the array 2080 asvf(nsvfl)%isurflt = isurflt 2081 asvf(nsvfl)%isurfs = isurfs 2082 asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency 2083 asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor 2084 ENDDO 2085 ENDDO 2086 2087 CALL location_message( ' waiting for completion of SVF and CSF calculation in all processes', .TRUE. ) 2088 !-- deallocate temporary global arrays 2089 DEALLOCATE(nzterr) 2090 2091 IF ( plant_canopy ) THEN 2092 !-- finalize mpi_rma communication and deallocate temporary arrays 2093 #if defined( __parallel ) 2094 IF ( usm_lad_rma ) THEN 2095 CALL MPI_Win_flush_all(win_lad, ierr) 2096 !-- unlock MPI window 2097 CALL MPI_Win_unlock_all(win_lad, ierr) 2098 !-- free MPI window 2099 CALL MPI_Win_free(win_lad, ierr) 2100 2101 !-- deallocate temporary arrays storing values for csf calculation during raytracing 2102 DEALLOCATE( lad_s_ray ) 2103 !-- usm_lad is the pointer to lad_s_rma in case of usm_lad_rma 2104 !-- and must not be deallocated here 2105 ELSE 2106 DEALLOCATE(usm_lad) 2107 DEALLOCATE(usm_lad_g) 2108 ENDIF 2109 #else 2110 DEALLOCATE(usm_lad) 2111 #endif 2112 DEALLOCATE( boxes ) 2113 DEALLOCATE( crlens ) 2114 DEALLOCATE( plantt ) 2115 ENDIF 2116 2117 CALL location_message( ' calculation of the complete SVF array', .TRUE. ) 2118 2119 !-- sort svf ( a version of quicksort ) 2120 CALL quicksort_svf(asvf,1,nsvfl) 2121 2122 ALLOCATE( svf(ndsvf,nsvfl) ) 2123 ALLOCATE( svfsurf(idsvf,nsvfl) ) 2124 2125 !< load svf from the structure array to plain arrays 2126 isurflt_prev = -1 2127 ksvf = 1 2128 svfsum = 0._wp 2129 DO isvf = 1, nsvfl 2130 !-- normalize svf per target face 2131 IF ( asvf(ksvf)%isurflt /= isurflt_prev ) THEN 2132 IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN 2133 !-- TODO detect and log when normalization differs too much from 1 2134 svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum 2135 ENDIF 2136 isurflt_prev = asvf(ksvf)%isurflt 2137 isvf_surflt = isvf 2138 svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp 2139 ELSE 2140 svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp 2141 ENDIF 2142 2143 svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /) 2144 svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /) 2145 2146 !-- next element 2147 ksvf = ksvf + 1 2148 ENDDO 2149 2150 IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN 2151 !-- TODO detect and log when normalization differs too much from 1 2152 svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum 2153 ENDIF 2154 2155 !-- deallocate temporary asvf array 2156 !-- DEALLOCATE(asvf) - ifort has a problem with deallocation of allocatable target 2157 !-- via pointing pointer - we need to test original targets 2158 IF ( ALLOCATED(asvf1) ) THEN 2159 DEALLOCATE(asvf1) 2160 ENDIF 2161 IF ( ALLOCATED(asvf2) ) THEN 2162 DEALLOCATE(asvf2) 2163 ENDIF 2164 2165 npcsfl = 0 2166 IF ( plant_canopy ) THEN 2167 2168 CALL location_message( ' calculation of the complete CSF array', .TRUE. ) 2169 2170 !-- sort and merge csf for the last time, keeping the array size to minimum 2171 CALL usm_merge_and_grow_csf(-1) 2172 2173 !-- aggregate csb among processors 2174 !-- allocate necessary arrays 2175 ALLOCATE( csflt(ndcsf,max(ncsfl,ndcsf)) ) 2176 ALLOCATE( kcsflt(kdcsf,max(ncsfl,kdcsf)) ) 2177 ALLOCATE( icsflt(0:numprocs-1) ) 2178 ALLOCATE( dcsflt(0:numprocs-1) ) 2179 ALLOCATE( ipcsflt(0:numprocs-1) ) 2180 ALLOCATE( dpcsflt(0:numprocs-1) ) 2181 2182 !-- fill out arrays of csf values and 2183 !-- arrays of number of elements and displacements 2184 !-- for particular precessors 2185 icsflt = 0 2186 dcsflt = 0 2187 ip = -1 2188 j = -1 2189 d = 0 2190 DO kcsf = 1, ncsfl 2191 j = j+1 2192 IF ( acsf(kcsf)%ip /= ip ) THEN 2193 !-- new block of the processor 2194 !-- number of elements of previous block 2195 IF ( ip>=0) icsflt(ip) = j 2196 d = d+j 2197 !-- blank blocks 2198 DO jp = ip+1, acsf(kcsf)%ip-1 2199 !-- number of elements is zero, displacement is equal to previous 2200 icsflt(jp) = 0 2201 dcsflt(jp) = d 2202 ENDDO 2203 !-- the actual block 2204 ip = acsf(kcsf)%ip 2205 dcsflt(ip) = d 2206 j = 0 2207 ENDIF 2208 !-- fill out real values of rsvf, rtransp 2209 csflt(1,kcsf) = acsf(kcsf)%rsvf 2210 csflt(2,kcsf) = acsf(kcsf)%rtransp 2211 !-- fill out integer values of itz,ity,itx,isurfs 2212 kcsflt(1,kcsf) = acsf(kcsf)%itz 2213 kcsflt(2,kcsf) = acsf(kcsf)%ity 2214 kcsflt(3,kcsf) = acsf(kcsf)%itx 2215 kcsflt(4,kcsf) = acsf(kcsf)%isurfs 2216 ENDDO 2217 !-- last blank blocks at the end of array 2218 j = j+1 2219 IF ( ip>=0 ) icsflt(ip) = j 2220 d = d+j 2221 DO jp = ip+1, numprocs-1 2222 !-- number of elements is zero, displacement is equal to previous 2223 icsflt(jp) = 0 2224 dcsflt(jp) = d 2225 ENDDO 2226 2227 !-- deallocate temporary acsf array 2228 !-- DEALLOCATE(acsf) - ifort has a problem with deallocation of allocatable target 2229 !-- via pointing pointer - we need to test original targets 2230 IF ( ALLOCATED(acsf1) ) THEN 2231 DEALLOCATE(acsf1) 2232 ENDIF 2233 IF ( ALLOCATED(acsf2) ) THEN 2234 DEALLOCATE(acsf2) 2235 ENDIF 2236 2237 #if defined( __parallel ) 2238 !-- scatter and gather the number of elements to and from all processor 2239 !-- and calculate displacements 2240 CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr) 2241 2242 npcsfl = SUM(ipcsflt) 2243 d = 0 2244 DO i = 0, numprocs-1 2245 dpcsflt(i) = d 2246 d = d + ipcsflt(i) 2247 ENDDO 2248 2249 !-- exchange csf fields between processors 2250 ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) ) 2251 ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) ) 2252 CALL MPI_AlltoAllv(csflt, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, & 2253 pcsflt, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr) 2254 CALL MPI_AlltoAllv(kcsflt, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, & 2255 kpcsflt, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr) 2256 2257 #else 2258 npcsfl = ncsfl 2259 ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) ) 2260 ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) ) 2261 pcsflt = csflt 2262 kpcsflt = kcsflt 2263 #endif 2264 2265 !-- deallocate temporary arrays 2266 DEALLOCATE( csflt ) 2267 DEALLOCATE( kcsflt ) 2268 DEALLOCATE( icsflt ) 2269 DEALLOCATE( dcsflt ) 2270 DEALLOCATE( ipcsflt ) 2271 DEALLOCATE( dpcsflt ) 2272 2273 !-- sort csf ( a version of quicksort ) 2274 CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl) 2275 2276 !-- aggregate canopy sink factor records with identical box & source 2277 !-- againg across all values from all processors 2278 IF ( npcsfl > 0 ) THEN 2279 icsf = 1 !< reading index 2280 kcsf = 1 !< writing index 2281 DO while (icsf < npcsfl) 2282 !-- here kpcsf(kcsf) already has values from kpcsf(icsf) 2283 IF ( kpcsflt(3,icsf) == kpcsflt(3,icsf+1) .AND. & 2284 kpcsflt(2,icsf) == kpcsflt(2,icsf+1) .AND. & 2285 kpcsflt(1,icsf) == kpcsflt(1,icsf+1) .AND. & 2286 kpcsflt(4,icsf) == kpcsflt(4,icsf+1) ) THEN 2287 !-- We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray 2288 !-- probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that 2289 !-- might mean that the traced beam passes longer through the canopy box. 2290 IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) ) THEN 2291 pcsflt(2,kcsf) = pcsflt(2,icsf+1) 2292 ENDIF 2293 pcsflt(1,kcsf) = pcsflt(1,kcsf) + pcsflt(1,icsf+1) 2294 2295 !-- advance reading index, keep writing index 2296 icsf = icsf + 1 2297 ELSE 2298 !-- not identical, just advance and copy 2299 icsf = icsf + 1 2300 kcsf = kcsf + 1 2301 kpcsflt(:,kcsf) = kpcsflt(:,icsf) 2302 pcsflt(:,kcsf) = pcsflt(:,icsf) 2303 ENDIF 2304 ENDDO 2305 !-- last written item is now also the last item in valid part of array 2306 npcsfl = kcsf 2307 ENDIF 2308 2309 ncsfl = npcsfl 2310 IF ( ncsfl > 0 ) THEN 2311 ALLOCATE( csf(ndcsf,ncsfl) ) 2312 ALLOCATE( csfsurf(idcsf,ncsfl) ) 2313 DO icsf = 1, ncsfl 2314 csf(:,icsf) = pcsflt(:,icsf) 2315 csfsurf(1,icsf) = gridpcbl(kpcsflt(1,icsf),kpcsflt(2,icsf),kpcsflt(3,icsf)) 2316 csfsurf(2,icsf) = kpcsflt(4,icsf) 2317 ENDDO 2318 ENDIF 2319 2320 !-- deallocation of temporary arrays 2321 DEALLOCATE( pcsflt ) 2322 DEALLOCATE( kpcsflt ) 2323 2324 ENDIF 2325 2326 RETURN 2327 2328 301 WRITE( message_string, * ) & 2329 'I/O error when processing shape view factors / ', & 2330 'plant canopy sink factors / direct irradiance factors.' 2331 CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 ) 2332 2333 END SUBROUTINE usm_calc_svf 2184 SUBROUTINE usm_boundary_condition 2185 2186 IMPLICIT NONE 2187 2188 INTEGER(iwp) :: i !< grid index x-direction 2189 INTEGER(iwp) :: ioff !< offset index x-direction indicating location of soil grid point 2190 INTEGER(iwp) :: j !< grid index y-direction 2191 INTEGER(iwp) :: joff !< offset index x-direction indicating location of soil grid point 2192 INTEGER(iwp) :: k !< grid index z-direction 2193 INTEGER(iwp) :: koff !< offset index x-direction indicating location of soil grid point 2194 INTEGER(iwp) :: l !< running index surface-orientation 2195 INTEGER(iwp) :: m !< running index surface elements 2196 2197 koff = surf_usm_h%koff 2198 DO m = 1, surf_usm_h%ns 2199 i = surf_usm_h%i(m) 2200 j = surf_usm_h%j(m) 2201 k = surf_usm_h%k(m) 2202 pt(k+koff,j,i) = pt(k,j,i) 2203 ENDDO 2204 2205 DO l = 0, 3 2206 ioff = surf_usm_v(l)%ioff 2207 joff = surf_usm_v(l)%joff 2208 DO m = 1, surf_usm_v(l)%ns 2209 i = surf_usm_v(l)%i(m) 2210 j = surf_usm_v(l)%j(m) 2211 k = surf_usm_v(l)%k(m) 2212 pt(k,j+joff,i+ioff) = pt(k,j,i) 2213 ENDDO 2214 ENDDO 2215 2216 END SUBROUTINE usm_boundary_condition 2334 2217 2335 2218 … … 2339 2222 ! ------------ 2340 2223 !> Subroutine checks variables and assigns units. 2341 !> It is ca aled out from subroutine check_parameters.2224 !> It is called out from subroutine check_parameters. 2342 2225 !------------------------------------------------------------------------------! 2343 2226 SUBROUTINE usm_check_data_output( variable, unit ) … … 2358 2241 var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. & 2359 2242 var(1:11) == 'usm_rad_hf_' .OR. & 2360 var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' ) THEN 2243 var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. & 2244 var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. & 2245 var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' ) THEN 2361 2246 unit = 'W/m2' 2362 ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' ) THEN 2247 ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. & 2248 var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR. & 2249 var(1:16) == 'usm_t_surf_green' .OR. & 2250 var(1:16) == 'usm_t_surf_whole' .OR. var(1:11) == 'usm_t_green' .OR. & 2251 var(1:15) == 'usm_t_surf_10cm') THEN 2363 2252 unit = 'K' 2364 2253 ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. & … … 2414 2303 ENDIF 2415 2304 2416 2417 2305 END SUBROUTINE usm_check_parameters 2418 2306 … … 2443 2331 INTEGER(iwp), PARAMETER :: nd = 5 2444 2332 CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 2445 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ i roof, isouth, inorth, iwest, ieast/)2333 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /) 2446 2334 INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart 2447 2335 INTEGER(iwp), DIMENSION(0:nd-1) :: dirend … … 2477 2365 ENDIF 2478 2366 ENDIF 2367 IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN 2368 !-- window layers 2369 READ(var(14:14), '(I1)', iostat=istat ) iwl 2370 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN 2371 var = var(1:12) 2372 ENDIF 2373 ENDIF 2374 IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN 2375 !-- green layers 2376 READ(var(13:13), '(I1)', iostat=istat ) iwl 2377 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN 2378 var = var(1:11) 2379 ENDIF 2380 ENDIF 2479 2381 IF ( (var(1:8) == 'usm_svf_' .OR. var(1:8) == 'usm_dif_') .AND. len(TRIM(var)) >= 13 ) THEN 2480 2382 !-- svf values to particular surface … … 2531 2433 2532 2434 CASE ( 'usm_surfalb' ) 2533 !-- surface albedo 2435 !-- surface albedo, weighted average 2534 2436 DO m = 1, surf_usm_h%ns 2535 2437 i = surf_usm_h%i(m) 2536 2438 j = surf_usm_h%j(m) 2537 2439 k = surf_usm_h%k(m) 2538 temp_pf(k,j,i) = surf_usm_h%albedo_surf(m) 2440 temp_pf(k,j,i) = surf_usm_h%frac(0,m) * & 2441 surf_usm_h%albedo(0,m) + & 2442 surf_usm_h%frac(1,m) * & 2443 surf_usm_h%albedo(1,m) + & 2444 surf_usm_h%frac(2,m) * & 2445 surf_usm_h%albedo(2,m) 2539 2446 ENDDO 2540 2447 DO l = 0, 3 … … 2543 2450 j = surf_usm_v(l)%j(m) 2544 2451 k = surf_usm_v(l)%k(m) 2545 temp_pf(k,j,i) = surf_usm_v(l)%albedo_surf(m) 2452 temp_pf(k,j,i) = surf_usm_v(l)%frac(0,m) * & 2453 surf_usm_v(l)%albedo(0,m) + & 2454 surf_usm_v(l)%frac(1,m) * & 2455 surf_usm_v(l)%albedo(1,m) + & 2456 surf_usm_v(l)%frac(2,m) * & 2457 surf_usm_v(l)%albedo(2,m) 2546 2458 ENDDO 2547 2459 ENDDO 2548 2460 2549 2461 CASE ( 'usm_surfemis' ) 2550 !-- surface albedo2462 !-- surface emissivity, weighted average 2551 2463 DO m = 1, surf_usm_h%ns 2552 2464 i = surf_usm_h%i(m) 2553 2465 j = surf_usm_h%j(m) 2554 2466 k = surf_usm_h%k(m) 2555 temp_pf(k,j,i) = surf_usm_h%emiss_surf(m) 2467 temp_pf(k,j,i) = surf_usm_h%frac(0,m) * & 2468 surf_usm_h%emissivity(0,m) + & 2469 surf_usm_h%frac(1,m) * & 2470 surf_usm_h%emissivity(1,m) + & 2471 surf_usm_h%frac(2,m) * & 2472 surf_usm_h%emissivity(2,m) 2556 2473 ENDDO 2557 2474 DO l = 0, 3 … … 2560 2477 j = surf_usm_v(l)%j(m) 2561 2478 k = surf_usm_v(l)%k(m) 2562 temp_pf(k,j,i) = surf_usm_v(l)%emiss_surf(m) 2479 temp_pf(k,j,i) = surf_usm_v(l)%frac(0,m) * & 2480 surf_usm_v(l)%emissivity(0,m) + & 2481 surf_usm_v(l)%frac(1,m) * & 2482 surf_usm_v(l)%emissivity(1,m) + & 2483 surf_usm_v(l)%frac(2,m) * & 2484 surf_usm_v(l)%emissivity(2,m) 2563 2485 ENDDO 2564 2486 ENDDO 2487 2488 CASE ( 'usm_surfwintrans' ) 2489 !-- transmissivity window tiles 2490 DO m = 1, surf_usm_h%ns 2491 i = surf_usm_h%i(m) 2492 j = surf_usm_h%j(m) 2493 k = surf_usm_h%k(m) 2494 temp_pf(k,j,i) = surf_usm_h%transmissivity(m) 2495 ENDDO 2496 DO l = 0, 3 2497 DO m = 1, surf_usm_v(l)%ns 2498 i = surf_usm_v(l)%i(m) 2499 j = surf_usm_v(l)%j(m) 2500 k = surf_usm_v(l)%k(m) 2501 temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m) 2502 ENDDO 2503 2504 ENDDO 2505 2565 2506 ! 2566 2507 !-- Not adjusted so far … … 2841 2782 ENDIF 2842 2783 2784 CASE ( 'usm_wghf_window' ) 2785 !-- array of heat flux from window ground (land, wall, roof) 2786 2787 IF ( av == 0 ) THEN 2788 DO m = 1, surf_usm_h%ns 2789 i = surf_usm_h%i(m) 2790 j = surf_usm_h%j(m) 2791 k = surf_usm_h%k(m) 2792 temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m) 2793 ENDDO 2794 DO l = 0, 3 2795 DO m = 1, surf_usm_v(l)%ns 2796 i = surf_usm_v(l)%i(m) 2797 j = surf_usm_v(l)%j(m) 2798 k = surf_usm_v(l)%k(m) 2799 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m) 2800 ENDDO 2801 ENDDO 2802 ELSE 2803 DO m = 1, surf_usm_h%ns 2804 i = surf_usm_h%i(m) 2805 j = surf_usm_h%j(m) 2806 k = surf_usm_h%k(m) 2807 temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m) 2808 ENDDO 2809 DO l = 0, 3 2810 DO m = 1, surf_usm_v(l)%ns 2811 i = surf_usm_v(l)%i(m) 2812 j = surf_usm_v(l)%j(m) 2813 k = surf_usm_v(l)%k(m) 2814 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m) 2815 ENDDO 2816 ENDDO 2817 ENDIF 2818 2819 CASE ( 'usm_wghf_green' ) 2820 !-- array of heat flux from green ground (land, wall, roof) 2821 2822 IF ( av == 0 ) THEN 2823 DO m = 1, surf_usm_h%ns 2824 i = surf_usm_h%i(m) 2825 j = surf_usm_h%j(m) 2826 k = surf_usm_h%k(m) 2827 temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m) 2828 ENDDO 2829 DO l = 0, 3 2830 DO m = 1, surf_usm_v(l)%ns 2831 i = surf_usm_v(l)%i(m) 2832 j = surf_usm_v(l)%j(m) 2833 k = surf_usm_v(l)%k(m) 2834 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m) 2835 ENDDO 2836 ENDDO 2837 ELSE 2838 DO m = 1, surf_usm_h%ns 2839 i = surf_usm_h%i(m) 2840 j = surf_usm_h%j(m) 2841 k = surf_usm_h%k(m) 2842 temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m) 2843 ENDDO 2844 DO l = 0, 3 2845 DO m = 1, surf_usm_v(l)%ns 2846 i = surf_usm_v(l)%i(m) 2847 j = surf_usm_v(l)%j(m) 2848 k = surf_usm_v(l)%k(m) 2849 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m) 2850 ENDDO 2851 ENDDO 2852 ENDIF 2853 2854 CASE ( 'usm_iwghf' ) 2855 !-- array of heat flux from indoor ground (land, wall, roof) 2856 IF ( av == 0 ) THEN 2857 DO m = 1, surf_usm_h%ns 2858 i = surf_usm_h%i(m) 2859 j = surf_usm_h%j(m) 2860 k = surf_usm_h%k(m) 2861 temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m) 2862 ENDDO 2863 DO l = 0, 3 2864 DO m = 1, surf_usm_v(l)%ns 2865 i = surf_usm_v(l)%i(m) 2866 j = surf_usm_v(l)%j(m) 2867 k = surf_usm_v(l)%k(m) 2868 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m) 2869 ENDDO 2870 ENDDO 2871 ELSE 2872 DO m = 1, surf_usm_h%ns 2873 i = surf_usm_h%i(m) 2874 j = surf_usm_h%j(m) 2875 k = surf_usm_h%k(m) 2876 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m) 2877 ENDDO 2878 DO l = 0, 3 2879 DO m = 1, surf_usm_v(l)%ns 2880 i = surf_usm_v(l)%i(m) 2881 j = surf_usm_v(l)%j(m) 2882 k = surf_usm_v(l)%k(m) 2883 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m) 2884 ENDDO 2885 ENDDO 2886 ENDIF 2887 2888 CASE ( 'usm_iwghf_window' ) 2889 !-- array of heat flux from indoor window ground (land, wall, roof) 2890 2891 IF ( av == 0 ) THEN 2892 DO m = 1, surf_usm_h%ns 2893 i = surf_usm_h%i(m) 2894 j = surf_usm_h%j(m) 2895 k = surf_usm_h%k(m) 2896 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m) 2897 ENDDO 2898 DO l = 0, 3 2899 DO m = 1, surf_usm_v(l)%ns 2900 i = surf_usm_v(l)%i(m) 2901 j = surf_usm_v(l)%j(m) 2902 k = surf_usm_v(l)%k(m) 2903 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m) 2904 ENDDO 2905 ENDDO 2906 ELSE 2907 DO m = 1, surf_usm_h%ns 2908 i = surf_usm_h%i(m) 2909 j = surf_usm_h%j(m) 2910 k = surf_usm_h%k(m) 2911 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m) 2912 ENDDO 2913 DO l = 0, 3 2914 DO m = 1, surf_usm_v(l)%ns 2915 i = surf_usm_v(l)%i(m) 2916 j = surf_usm_v(l)%j(m) 2917 k = surf_usm_v(l)%k(m) 2918 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m) 2919 ENDDO 2920 ENDDO 2921 ENDIF 2922 2843 2923 CASE ( 'usm_t_surf' ) 2844 2924 !-- surface temperature for surfaces … … 2874 2954 ENDDO 2875 2955 ENDIF 2956 2957 CASE ( 'usm_t_surf_window' ) 2958 !-- surface temperature for window surfaces 2959 2960 IF ( av == 0 ) THEN 2961 DO m = 1, surf_usm_h%ns 2962 i = surf_usm_h%i(m) 2963 j = surf_usm_h%j(m) 2964 k = surf_usm_h%k(m) 2965 temp_pf(k,j,i) = t_surf_window_h(m) 2966 ENDDO 2967 DO l = 0, 3 2968 DO m = 1, surf_usm_v(l)%ns 2969 i = surf_usm_v(l)%i(m) 2970 j = surf_usm_v(l)%j(m) 2971 k = surf_usm_v(l)%k(m) 2972 temp_pf(k,j,i) = t_surf_window_v(l)%t(m) 2973 ENDDO 2974 ENDDO 2975 2976 ELSE 2977 DO m = 1, surf_usm_h%ns 2978 i = surf_usm_h%i(m) 2979 j = surf_usm_h%j(m) 2980 k = surf_usm_h%k(m) 2981 temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m) 2982 ENDDO 2983 DO l = 0, 3 2984 DO m = 1, surf_usm_v(l)%ns 2985 i = surf_usm_v(l)%i(m) 2986 j = surf_usm_v(l)%j(m) 2987 k = surf_usm_v(l)%k(m) 2988 temp_pf(k,j,i) = surf_usm_v(l)%t_surf_window_av(m) 2989 ENDDO 2990 2991 ENDDO 2992 2993 ENDIF 2994 2995 CASE ( 'usm_t_surf_green' ) 2996 !-- surface temperature for green surfaces 2997 2998 IF ( av == 0 ) THEN 2999 DO m = 1, surf_usm_h%ns 3000 i = surf_usm_h%i(m) 3001 j = surf_usm_h%j(m) 3002 k = surf_usm_h%k(m) 3003 temp_pf(k,j,i) = t_surf_green_h(m) 3004 ENDDO 3005 DO l = 0, 3 3006 DO m = 1, surf_usm_v(l)%ns 3007 i = surf_usm_v(l)%i(m) 3008 j = surf_usm_v(l)%j(m) 3009 k = surf_usm_v(l)%k(m) 3010 temp_pf(k,j,i) = t_surf_green_v(l)%t(m) 3011 ENDDO 3012 ENDDO 3013 3014 ELSE 3015 DO m = 1, surf_usm_h%ns 3016 i = surf_usm_h%i(m) 3017 j = surf_usm_h%j(m) 3018 k = surf_usm_h%k(m) 3019 temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m) 3020 ENDDO 3021 DO l = 0, 3 3022 DO m = 1, surf_usm_v(l)%ns 3023 i = surf_usm_v(l)%i(m) 3024 j = surf_usm_v(l)%j(m) 3025 k = surf_usm_v(l)%k(m) 3026 temp_pf(k,j,i) = surf_usm_v(l)%t_surf_green_av(m) 3027 ENDDO 3028 3029 ENDDO 3030 3031 ENDIF 3032 3033 CASE ( 'usm_t_surf_whole' ) 3034 !-- surface temperature for whole surfaces 3035 3036 IF ( av == 0 ) THEN 3037 DO m = 1, surf_usm_h%ns 3038 i = surf_usm_h%i(m) 3039 j = surf_usm_h%j(m) 3040 k = surf_usm_h%k(m) 3041 temp_pf(k,j,i) = t_surf_whole_h(m) 3042 ENDDO 3043 DO l = 0, 3 3044 DO m = 1, surf_usm_v(l)%ns 3045 i = surf_usm_v(l)%i(m) 3046 j = surf_usm_v(l)%j(m) 3047 k = surf_usm_v(l)%k(m) 3048 temp_pf(k,j,i) = t_surf_whole_v(l)%t(m) 3049 ENDDO 3050 ENDDO 3051 3052 ELSE 3053 DO m = 1, surf_usm_h%ns 3054 i = surf_usm_h%i(m) 3055 j = surf_usm_h%j(m) 3056 k = surf_usm_h%k(m) 3057 temp_pf(k,j,i) = surf_usm_h%t_surf_whole_av(m) 3058 ENDDO 3059 DO l = 0, 3 3060 DO m = 1, surf_usm_v(l)%ns 3061 i = surf_usm_v(l)%i(m) 3062 j = surf_usm_v(l)%j(m) 3063 k = surf_usm_v(l)%k(m) 3064 temp_pf(k,j,i) = surf_usm_v(l)%t_surf_whole_av(m) 3065 ENDDO 3066 3067 ENDDO 3068 3069 ENDIF 3070 3071 CASE ( 'usm_t_surf_10cm' ) 3072 !-- near surface temperature for whole surfaces 3073 3074 IF ( av == 0 ) THEN 3075 DO m = 1, surf_usm_h%ns 3076 i = surf_usm_h%i(m) 3077 j = surf_usm_h%j(m) 3078 k = surf_usm_h%k(m) 3079 temp_pf(k,j,i) = t_surf_10cm_h(m) 3080 ENDDO 3081 DO l = 0, 3 3082 DO m = 1, surf_usm_v(l)%ns 3083 i = surf_usm_v(l)%i(m) 3084 j = surf_usm_v(l)%j(m) 3085 k = surf_usm_v(l)%k(m) 3086 temp_pf(k,j,i) = t_surf_10cm_v(l)%t(m) 3087 ENDDO 3088 ENDDO 3089 3090 ELSE 3091 DO m = 1, surf_usm_h%ns 3092 i = surf_usm_h%i(m) 3093 j = surf_usm_h%j(m) 3094 k = surf_usm_h%k(m) 3095 temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m) 3096 ENDDO 3097 DO l = 0, 3 3098 DO m = 1, surf_usm_v(l)%ns 3099 i = surf_usm_v(l)%i(m) 3100 j = surf_usm_v(l)%j(m) 3101 k = surf_usm_v(l)%k(m) 3102 temp_pf(k,j,i) = surf_usm_v(l)%t_surf_10cm_av(m) 3103 ENDDO 3104 3105 ENDDO 3106 3107 ENDIF 3108 2876 3109 2877 3110 CASE ( 'usm_t_wall' ) … … 2909 3142 ENDIF 2910 3143 3144 CASE ( 'usm_t_window' ) 3145 !-- window temperature for iwl layer of walls and land 3146 IF ( av == 0 ) THEN 3147 DO m = 1, surf_usm_h%ns 3148 i = surf_usm_h%i(m) 3149 j = surf_usm_h%j(m) 3150 k = surf_usm_h%k(m) 3151 temp_pf(k,j,i) = t_window_h(iwl,m) 3152 ENDDO 3153 DO l = 0, 3 3154 DO m = 1, surf_usm_v(l)%ns 3155 i = surf_usm_v(l)%i(m) 3156 j = surf_usm_v(l)%j(m) 3157 k = surf_usm_v(l)%k(m) 3158 temp_pf(k,j,i) = t_window_v(l)%t(iwl,m) 3159 ENDDO 3160 ENDDO 3161 ELSE 3162 DO m = 1, surf_usm_h%ns 3163 i = surf_usm_h%i(m) 3164 j = surf_usm_h%j(m) 3165 k = surf_usm_h%k(m) 3166 temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m) 3167 ENDDO 3168 DO l = 0, 3 3169 DO m = 1, surf_usm_v(l)%ns 3170 i = surf_usm_v(l)%i(m) 3171 j = surf_usm_v(l)%j(m) 3172 k = surf_usm_v(l)%k(m) 3173 temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m) 3174 ENDDO 3175 ENDDO 3176 ENDIF 3177 3178 CASE ( 'usm_t_green' ) 3179 !-- green temperature for iwl layer of walls and land 3180 IF ( av == 0 ) THEN 3181 DO m = 1, surf_usm_h%ns 3182 i = surf_usm_h%i(m) 3183 j = surf_usm_h%j(m) 3184 k = surf_usm_h%k(m) 3185 temp_pf(k,j,i) = t_green_h(iwl,m) 3186 ENDDO 3187 DO l = 0, 3 3188 DO m = 1, surf_usm_v(l)%ns 3189 i = surf_usm_v(l)%i(m) 3190 j = surf_usm_v(l)%j(m) 3191 k = surf_usm_v(l)%k(m) 3192 temp_pf(k,j,i) = t_green_v(l)%t(iwl,m) 3193 ENDDO 3194 ENDDO 3195 ELSE 3196 DO m = 1, surf_usm_h%ns 3197 i = surf_usm_h%i(m) 3198 j = surf_usm_h%j(m) 3199 k = surf_usm_h%k(m) 3200 temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m) 3201 ENDDO 3202 DO l = 0, 3 3203 DO m = 1, surf_usm_v(l)%ns 3204 i = surf_usm_v(l)%i(m) 3205 j = surf_usm_v(l)%j(m) 3206 k = surf_usm_v(l)%k(m) 3207 temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m) 3208 ENDDO 3209 ENDDO 3210 ENDIF 3211 3212 2911 3213 CASE DEFAULT 2912 3214 found = .FALSE. 2913 3215 2914 3216 END SELECT 3217 2915 3218 ! 2916 3219 !-- Rearrange dimensions for NetCDF output … … 2954 3257 var(1:11) == 'usm_rad_hf_' .OR. & 2955 3258 var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. & 3259 var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. & 3260 var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. & 2956 3261 var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. & 3262 var(1:17) == 'usm_t_surf_window' .OR. var(1:12) == 'usm_t_window' .OR. & 3263 var(1:16) == 'usm_t_surf_green' .OR. var(1:16) == 'usm_t_surf_whole' .OR. & 3264 var(1:15) == 'usm_t_surf_10cm' .OR. & 2957 3265 var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. & 2958 3266 var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. & 2959 var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis' ) THEN 3267 var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis' .OR. & 3268 var(1:16) == 'usm_surfwintrans' ) THEN 2960 3269 2961 3270 found = .TRUE. … … 2972 3281 END SUBROUTINE usm_define_netcdf_grid 2973 3282 2974 2975 !------------------------------------------------------------------------------!2976 !> Finds first model boundary crossed by a ray2977 !------------------------------------------------------------------------------!2978 PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross)2979 2980 IMPLICIT NONE2981 2982 INTEGER(iwp) :: d !<2983 INTEGER(iwp) :: seldim !< found fist crossing index2984 2985 INTEGER(iwp), DIMENSION(3) :: bdyd !< boundary direction2986 INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross !< found boundary crossing (d, z, y, x)2987 2988 REAL(wp) :: bdydim !<2989 REAL(wp) :: dist !<2990 2991 REAL(wp), DIMENSION(3) :: crossdist !< crossing distance2992 REAL(wp), DIMENSION(3), INTENT(in) :: origin !< ray origin2993 REAL(wp), DIMENSION(3), INTENT(in) :: uvect !< ray unit vector2994 2995 2996 bdydim = nzut + .5_wp !< top boundary2997 bdyd(1) = isky2998 crossdist(1) = ( bdydim - origin(1) ) / uvect(1) !< subroutine called only when uvect(1)>02999 3000 IF ( uvect(2) == 0._wp ) THEN3001 crossdist(2) = huge(1._wp)3002 ELSE3003 IF ( uvect(2) >= 0._wp ) THEN3004 bdydim = ny + .5_wp !< north global boundary3005 bdyd(2) = inorthb3006 ELSE3007 bdydim = -.5_wp !< south global boundary3008 bdyd(2) = isouthb3009 ENDIF3010 crossdist(2) = ( bdydim - origin(2) ) / uvect(2)3011 ENDIF3012 3013 IF ( uvect(3) == 0._wp ) THEN3014 crossdist(3) = huge(1._wp)3015 ELSE3016 IF ( uvect(3) >= 0._wp ) THEN3017 bdydim = nx + .5_wp !< east global boundary3018 bdyd(3) = ieastb3019 ELSE3020 bdydim = -.5_wp !< west global boundary3021 bdyd(3) = iwestb3022 ENDIF3023 crossdist(3) = ( bdydim - origin(3) ) / uvect(3)3024 ENDIF3025 3026 seldim = minloc(crossdist, 1)3027 dist = crossdist(seldim)3028 d = bdyd(seldim)3029 3030 bdycross(1) = d3031 bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &3032 + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )3033 3034 END SUBROUTINE3035 3036 3037 !------------------------------------------------------------------------------!3038 !> Determines whether two faces are oriented towards each other3039 !------------------------------------------------------------------------------!3040 PURE LOGICAL FUNCTION usm_facing(x, y, z, d, x2, y2, z2, d2)3041 IMPLICIT NONE3042 INTEGER(iwp), INTENT(in) :: x, y, z, d, x2, y2, z2, d23043 3044 usm_facing = .FALSE.3045 IF ( d==iroof .AND. d2==iroof ) RETURN3046 IF ( d==isky .AND. d2==isky ) RETURN3047 IF ( (d==isouth .OR. d==inorthb) .AND. (d2==isouth.OR.d2==inorthb) ) RETURN3048 IF ( (d==inorth .OR. d==isouthb) .AND. (d2==inorth.OR.d2==isouthb) ) RETURN3049 IF ( (d==iwest .OR. d==ieastb) .AND. (d2==iwest.OR.d2==ieastb) ) RETURN3050 IF ( (d==ieast .OR. d==iwestb) .AND. (d2==ieast.OR.d2==iwestb) ) RETURN3051 3052 SELECT CASE (d)3053 CASE (iroof) !< ground, roof3054 IF ( z2 < z ) RETURN3055 CASE (isky) !< sky3056 IF ( z2 > z ) RETURN3057 CASE (isouth, inorthb) !< south facing3058 IF ( y2 > y ) RETURN3059 CASE (inorth, isouthb) !< north facing3060 IF ( y2 < y ) RETURN3061 CASE (iwest, ieastb) !< west facing3062 IF ( x2 > x ) RETURN3063 CASE (ieast, iwestb) !< east facing3064 IF ( x2 < x ) RETURN3065 END SELECT3066 3067 SELECT CASE (d2)3068 CASE (iroof) !< ground, roof3069 IF ( z < z2 ) RETURN3070 CASE (isky) !< sky3071 IF ( z > z2 ) RETURN3072 CASE (isouth, inorthb) !< south facing3073 IF ( y > y2 ) RETURN3074 CASE (inorth, isouthb) !< north facing3075 IF ( y < y2 ) RETURN3076 CASE (iwest, ieastb) !< west facing3077 IF ( x > x2 ) RETURN3078 CASE (ieast, iwestb) !< east facing3079 IF ( x < x2 ) RETURN3080 CASE (-1)3081 CONTINUE3082 END SELECT3083 3084 usm_facing = .TRUE.3085 3086 END FUNCTION usm_facing3087 3088 3283 3089 3284 !------------------------------------------------------------------------------! … … 3102 3297 !-- Calculate wall grid spacings. 3103 3298 !-- Temperature is defined at the center of the wall layers, 3104 !-- whereas gradients/fluxes are defined at the edges (_stag) 3105 DO k = nzb_wall, nzt_wall 3106 zwn(k) = zwn_default(k) 3107 ENDDO 3108 ! 3299 !-- whereas gradients/fluxes are defined at the edges (_stag) 3109 3300 !-- apply for all particular surface grids. First for horizontal surfaces 3110 3301 DO m = 1, surf_usm_h%ns 3111 surf_usm_h%zw(:,m) = zwn(:) * & 3112 surf_usm_h%thickness_wall(m) 3302 3113 3303 surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m) 3114 3304 DO k = nzb_wall+1, nzt_wall … … 3116 3306 surf_usm_h%zw(k-1,m) 3117 3307 ENDDO 3308 surf_usm_h%dz_window(nzb_wall,m) = surf_usm_h%zw_window(nzb_wall,m) 3309 DO k = nzb_wall+1, nzt_wall 3310 surf_usm_h%dz_window(k,m) = surf_usm_h%zw_window(k,m) - & 3311 surf_usm_h%zw_window(k-1,m) 3312 ENDDO 3313 surf_usm_h%dz_green(nzb_wall,m) = surf_usm_h%zw_green(nzb_wall,m) 3314 DO k = nzb_wall+1, nzt_wall 3315 surf_usm_h%dz_green(k,m) = surf_usm_h%zw_green(k,m) - & 3316 surf_usm_h%zw_green(k-1,m) 3317 ENDDO 3118 3318 3119 3319 surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m) … … 3124 3324 ENDDO 3125 3325 surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m) 3326 3327 surf_usm_h%dz_window(nzt_wall+1,m) = surf_usm_h%dz_window(nzt_wall,m) 3328 3329 DO k = nzb_wall, nzt_wall-1 3330 surf_usm_h%dz_window_stag(k,m) = 0.5 * ( & 3331 surf_usm_h%dz_window(k+1,m) + surf_usm_h%dz_window(k,m) ) 3332 ENDDO 3333 surf_usm_h%dz_window_stag(nzt_wall,m) = surf_usm_h%dz_window(nzt_wall,m) 3334 3335 surf_usm_h%dz_green(nzt_wall+1,m) = surf_usm_h%dz_green(nzt_wall,m) 3336 3337 DO k = nzb_wall, nzt_wall-1 3338 surf_usm_h%dz_green_stag(k,m) = 0.5 * ( & 3339 surf_usm_h%dz_green(k+1,m) + surf_usm_h%dz_green(k,m) ) 3340 ENDDO 3341 surf_usm_h%dz_green_stag(nzt_wall,m) = surf_usm_h%dz_green(nzt_wall,m) 3126 3342 ENDDO 3127 surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall 3128 surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag 3343 surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall 3344 surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag 3345 surf_usm_h%ddz_window = 1.0_wp / surf_usm_h%dz_window 3346 surf_usm_h%ddz_window_stag = 1.0_wp / surf_usm_h%dz_window_stag 3347 surf_usm_h%ddz_green = 1.0_wp / surf_usm_h%dz_green 3348 surf_usm_h%ddz_green_stag = 1.0_wp / surf_usm_h%dz_green_stag 3129 3349 ! 3130 3350 !-- For vertical surfaces 3131 3351 DO l = 0, 3 3132 3352 DO m = 1, surf_usm_v(l)%ns 3133 surf_usm_v(l)%zw(:,m) = zwn(:) * &3134 surf_usm_v(l)%thickness_wall(m)3135 3353 surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m) 3136 3354 DO k = nzb_wall+1, nzt_wall 3137 3355 surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - & 3138 3356 surf_usm_v(l)%zw(k-1,m) 3357 ENDDO 3358 surf_usm_v(l)%dz_window(nzb_wall,m) = surf_usm_v(l)%zw_window(nzb_wall,m) 3359 DO k = nzb_wall+1, nzt_wall 3360 surf_usm_v(l)%dz_window(k,m) = surf_usm_v(l)%zw_window(k,m) - & 3361 surf_usm_v(l)%zw_window(k-1,m) 3362 ENDDO 3363 surf_usm_v(l)%dz_green(nzb_wall,m) = surf_usm_v(l)%zw_green(nzb_wall,m) 3364 DO k = nzb_wall+1, nzt_wall 3365 surf_usm_v(l)%dz_green(k,m) = surf_usm_v(l)%zw_green(k,m) - & 3366 surf_usm_v(l)%zw_green(k-1,m) 3139 3367 ENDDO 3140 3368 … … 3149 3377 surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = & 3150 3378 surf_usm_v(l)%dz_wall(nzt_wall,m) 3379 surf_usm_v(l)%dz_window(nzt_wall+1,m) = & 3380 surf_usm_v(l)%dz_window(nzt_wall,m) 3381 3382 DO k = nzb_wall, nzt_wall-1 3383 surf_usm_v(l)%dz_window_stag(k,m) = 0.5 * ( & 3384 surf_usm_v(l)%dz_window(k+1,m) + & 3385 surf_usm_v(l)%dz_window(k,m) ) 3386 ENDDO 3387 surf_usm_v(l)%dz_window_stag(nzt_wall,m) = & 3388 surf_usm_v(l)%dz_window(nzt_wall,m) 3389 surf_usm_v(l)%dz_green(nzt_wall+1,m) = & 3390 surf_usm_v(l)%dz_green(nzt_wall,m) 3391 3392 DO k = nzb_wall, nzt_wall-1 3393 surf_usm_v(l)%dz_green_stag(k,m) = 0.5 * ( & 3394 surf_usm_v(l)%dz_green(k+1,m) + & 3395 surf_usm_v(l)%dz_green(k,m) ) 3396 ENDDO 3397 surf_usm_v(l)%dz_green_stag(nzt_wall,m) = & 3398 surf_usm_v(l)%dz_green(nzt_wall,m) 3151 3399 ENDDO 3152 surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall 3153 surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag 3400 surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall 3401 surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag 3402 surf_usm_v(l)%ddz_window = 1.0_wp / surf_usm_v(l)%dz_window 3403 surf_usm_v(l)%ddz_window_stag = 1.0_wp / surf_usm_v(l)%dz_window_stag 3404 surf_usm_v(l)%ddz_green = 1.0_wp / surf_usm_v(l)%dz_green 3405 surf_usm_v(l)%ddz_green_stag = 1.0_wp / surf_usm_v(l)%dz_green_stag 3154 3406 ENDDO 3155 3407 … … 3168 3420 !------------------------------------------------------------------------------! 3169 3421 SUBROUTINE usm_init_urban_surface 3422 3423 USE arrays_3d, & 3424 ONLY: zw 3425 3426 USE netcdf_data_input_mod, & 3427 ONLY: building_pars_f, building_type_f, terrain_height_f 3170 3428 3171 3429 IMPLICIT NONE 3172 3430 3173 INTEGER(iwp) :: i, j, k, l, m !< running indices 3174 REAL(wp) :: c, d, tin, exn 3431 INTEGER(iwp) :: i !< loop index x-dirction 3432 INTEGER(iwp) :: ind_emis_wall !< index in input list for wall emissivity 3433 INTEGER(iwp) :: ind_emis_green !< index in input list for green emissivity 3434 INTEGER(iwp) :: ind_emis_win !< index in input list for window emissivity 3435 INTEGER(iwp) :: ind_green_frac_w !< index in input list for green fraction on wall 3436 INTEGER(iwp) :: ind_green_frac_r !< index in input list for green fraction on roof 3437 INTEGER(iwp) :: ind_hc1 !< index in input list for heat capacity at first wall layer 3438 INTEGER(iwp) :: ind_hc2 !< index in input list for heat capacity at second wall layer 3439 INTEGER(iwp) :: ind_hc3 !< index in input list for heat capacity at third wall layer 3440 INTEGER(iwp) :: ind_lai_r !< index in input list for LAI on roof 3441 INTEGER(iwp) :: ind_lai_w !< index in input list for LAI on wall 3442 INTEGER(iwp) :: ind_tc1 !< index in input list for thermal conductivity at first wall layer 3443 INTEGER(iwp) :: ind_tc2 !< index in input list for thermal conductivity at second wall layer 3444 INTEGER(iwp) :: ind_tc3 !< index in input list for thermal conductivity at third wall layer 3445 INTEGER(iwp) :: ind_trans !< index in input list for window transmissivity 3446 INTEGER(iwp) :: ind_wall_frac !< index in input list for wall fraction 3447 INTEGER(iwp) :: ind_win_frac !< index in input list for window fraction 3448 INTEGER(iwp) :: ind_z0 !< index in input list for z0 3449 INTEGER(iwp) :: ind_z0qh !< index in input list for z0h / z0q 3450 INTEGER(iwp) :: j !< loop index y-dirction 3451 INTEGER(iwp) :: k !< loop index z-dirction 3452 INTEGER(iwp) :: l !< loop index surface orientation 3453 INTEGER(iwp) :: m !< loop index surface element 3454 INTEGER(iwp) :: st !< dummy 3455 3456 REAL(wp) :: c, d, tin, twin, exn 3457 REAL(wp) :: ground_floor_level_l !< local height of ground floor level 3458 REAL(wp) :: z_agl !< height above ground 3175 3459 3176 3460 ! … … 3187 3471 lsf_surf = .FALSE. 3188 3472 ENDIF 3189 3190 !-- init anthropogenic sources of heat 3191 CALL usm_allocate_urban_surface() 3192 3193 !-- read the surface_types array somewhere 3473 3474 ! 3475 !-- Flag surface elements belonging to the ground floor level. Therefore, 3476 !-- use terrain height array from file, if available. This flag is later used 3477 !-- to control initialization of surface attributes. 3478 surf_usm_h%ground_level = .FALSE. 3479 DO m = 1, surf_usm_h%ns 3480 i = surf_usm_h%i(m) 3481 j = surf_usm_h%j(m) 3482 k = surf_usm_h%k(m) 3483 ! 3484 !-- Get local ground level. If no ground level is given in input file, 3485 !-- use default value. 3486 ground_floor_level_l = ground_floor_level 3487 IF ( building_pars_f%from_file ) THEN 3488 IF ( building_pars_f%pars_xy(ind_gflh,j,i) /= & 3489 building_pars_f%fill ) & 3490 ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i) 3491 ENDIF 3492 ! 3493 !-- Determine height of surface element above ground level 3494 IF ( terrain_height_f%from_file ) THEN 3495 z_agl = zw(k) - terrain_height_f%var(j,i) 3496 ELSE 3497 z_agl = zw(k) 3498 ENDIF 3499 ! 3500 !-- Set flag for ground level 3501 IF ( z_agl <= ground_floor_level_l ) & 3502 surf_usm_h%ground_level(m) = .TRUE. 3503 ENDDO 3504 3505 DO l = 0, 3 3506 surf_usm_v(l)%ground_level = .FALSE. 3507 DO m = 1, surf_usm_v(l)%ns 3508 i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff 3509 j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff 3510 k = surf_usm_v(l)%k(m) 3511 ! 3512 !-- Get local ground level. If no ground level is given in input file, 3513 !-- use default value. 3514 ground_floor_level_l = ground_floor_level 3515 IF ( building_pars_f%from_file ) THEN 3516 IF ( building_pars_f%pars_xy(ind_gflh,j,i) /= & 3517 building_pars_f%fill ) & 3518 ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i) 3519 ENDIF 3520 ! 3521 !-- Determine height of surface element above ground level. Please 3522 !-- note, height of surface element is determined with respect to 3523 !-- its height of the adjoing atmospheric grid point. 3524 IF ( terrain_height_f%from_file ) THEN 3525 z_agl = zw(k) - terrain_height_f%var(j-surf_usm_v(l)%joff, & 3526 i-surf_usm_v(l)%ioff) 3527 ELSE 3528 z_agl = zw(k) 3529 ENDIF 3530 ! 3531 !-- Set flag for ground level 3532 IF ( z_agl <= ground_floor_level_l ) & 3533 surf_usm_v(l)%ground_level(m) = .TRUE. 3534 3535 ENDDO 3536 ENDDO 3537 ! 3538 !-- Initialize urban-type surface attribute. According to initialization in 3539 !-- land-surface model, follow a 3-level approach. 3540 !-- Level 1 - initialization via default attributes 3541 DO m = 1, surf_usm_h%ns 3542 ! 3543 !-- Now, all horizontal surfaces are roof surfaces (?) 3544 surf_usm_h%isroof_surf(m) = .TRUE. 3545 surf_usm_h%surface_types(m) = roof_category !< default category for root surface 3546 ! 3547 !-- In order to distinguish between ground floor level and 3548 !-- above-ground-floor level surfaces, set input indices. 3549 ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, & 3550 surf_usm_h%ground_level(m) ) 3551 ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, & 3552 surf_usm_h%ground_level(m) ) 3553 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, & 3554 surf_usm_h%ground_level(m) ) 3555 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, & 3556 surf_usm_h%ground_level(m) ) 3557 ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, & 3558 surf_usm_h%ground_level(m) ) 3559 ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, & 3560 surf_usm_h%ground_level(m) ) 3561 ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, & 3562 surf_usm_h%ground_level(m) ) 3563 ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, & 3564 surf_usm_h%ground_level(m) ) 3565 ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, & 3566 surf_usm_h%ground_level(m) ) 3567 ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, & 3568 surf_usm_h%ground_level(m) ) 3569 ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, & 3570 surf_usm_h%ground_level(m) ) 3571 ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, & 3572 surf_usm_h%ground_level(m) ) 3573 ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, & 3574 surf_usm_h%ground_level(m) ) 3575 ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, & 3576 surf_usm_h%ground_level(m) ) 3577 ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, & 3578 surf_usm_h%ground_level(m) ) 3579 ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, & 3580 surf_usm_h%ground_level(m) ) 3581 ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, & 3582 surf_usm_h%ground_level(m) ) 3583 ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, & 3584 surf_usm_h%ground_level(m) ) 3585 ! 3586 !-- Initialize relatvie wall- (0), green- (1) and window (2) fractions 3587 surf_usm_h%frac(0,m) = building_pars(ind_wall_frac,building_type) 3588 surf_usm_h%frac(1,m) = building_pars(ind_green_frac_r,building_type) 3589 surf_usm_h%frac(2,m) = building_pars(ind_win_frac,building_type) 3590 surf_usm_h%lai(m) = building_pars(ind_green_frac_r,building_type) 3591 3592 surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,building_type) 3593 surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type) 3594 surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type) 3595 surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type) 3596 surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1,building_type) 3597 surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 3598 surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type) 3599 surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type) 3600 surf_usm_h%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,building_type) 3601 surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,building_type) 3602 surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,building_type) 3603 surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,building_type) 3604 surf_usm_h%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,building_type) 3605 surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 3606 surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,building_type) 3607 surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,building_type) 3608 surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,building_type) 3609 surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,building_type) 3610 surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,building_type) 3611 surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,building_type) 3612 surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,building_type) 3613 surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 3614 surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,building_type) 3615 surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,building_type) 3616 3617 surf_usm_h%target_temp_summer(m) = building_pars(12,building_type) 3618 surf_usm_h%target_temp_winter(m) = building_pars(13,building_type) 3619 ! 3620 !-- emissivity of wall-, green- and window fraction 3621 surf_usm_h%emissivity(0,m) = building_pars(ind_emis_wall,building_type) 3622 surf_usm_h%emissivity(1,m) = building_pars(ind_emis_green,building_type) 3623 surf_usm_h%emissivity(2,m) = building_pars(ind_emis_win,building_type) 3624 3625 surf_usm_h%transmissivity(m) = building_pars(ind_trans,building_type) 3626 3627 surf_usm_h%z0(m) = building_pars(ind_z0,building_type) 3628 surf_usm_h%z0h(m) = building_pars(ind_z0qh,building_type) 3629 surf_usm_h%z0q(m) = building_pars(ind_z0qh,building_type) 3630 ! 3631 !-- albedo type for wall fraction, green fraction, window fraction 3632 surf_usm_h%albedo_type(0,m) = INT( building_pars(ind_alb_wall,building_type) ) 3633 surf_usm_h%albedo_type(1,m) = INT( building_pars(ind_alb_green,building_type) ) 3634 surf_usm_h%albedo_type(2,m) = INT( building_pars(ind_alb_win,building_type) ) 3635 3636 surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1,building_type) 3637 surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2,building_type) 3638 surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3,building_type) 3639 surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4,building_type) 3640 3641 surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1,building_type) 3642 surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,building_type) 3643 surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,building_type) 3644 surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,building_type) 3645 3646 surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1,building_type) 3647 surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,building_type) 3648 surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,building_type) 3649 surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,building_type) 3650 3651 surf_usm_h%c_surface(m) = building_pars(45,building_type) 3652 surf_usm_h%lambda_surf(m) = building_pars(46,building_type) 3653 surf_usm_h%c_surface_green(m) = building_pars(45,building_type) 3654 surf_usm_h%lambda_surf_green(m) = building_pars(46,building_type) 3655 surf_usm_h%c_surface_window(m) = building_pars(45,building_type) 3656 surf_usm_h%lambda_surf_window(m) = building_pars(46,building_type) 3657 3658 ENDDO 3659 3660 DO l = 0, 3 3661 DO m = 1, surf_usm_v(l)%ns 3662 3663 surf_usm_v(l)%surface_types(m) = wall_category !< default category for root surface 3664 ! 3665 !-- In order to distinguish between ground floor level and 3666 !-- above-ground-floor level surfaces, set input indices. 3667 ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, & 3668 surf_usm_v(l)%ground_level(m) ) 3669 ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, & 3670 surf_usm_v(l)%ground_level(m) ) 3671 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, & 3672 surf_usm_v(l)%ground_level(m) ) 3673 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, & 3674 surf_usm_v(l)%ground_level(m) ) 3675 ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, & 3676 surf_usm_v(l)%ground_level(m) ) 3677 ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, & 3678 surf_usm_v(l)%ground_level(m) ) 3679 ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, & 3680 surf_usm_v(l)%ground_level(m) ) 3681 ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, & 3682 surf_usm_v(l)%ground_level(m) ) 3683 ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, & 3684 surf_usm_v(l)%ground_level(m) ) 3685 ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, & 3686 surf_usm_v(l)%ground_level(m) ) 3687 ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, & 3688 surf_usm_v(l)%ground_level(m) ) 3689 ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, & 3690 surf_usm_v(l)%ground_level(m) ) 3691 ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, & 3692 surf_usm_v(l)%ground_level(m) ) 3693 ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, & 3694 surf_usm_v(l)%ground_level(m) ) 3695 ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, & 3696 surf_usm_v(l)%ground_level(m) ) 3697 ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, & 3698 surf_usm_v(l)%ground_level(m) ) 3699 ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, & 3700 surf_usm_v(l)%ground_level(m) ) 3701 ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, & 3702 surf_usm_v(l)%ground_level(m) ) 3703 3704 ! 3705 !-- Initialize relatvie wall- (0), green- (1) and window (2) fractions 3706 surf_usm_v(l)%frac(0,m) = building_pars(ind_wall_frac,building_type) 3707 surf_usm_v(l)%frac(1,m) = building_pars(ind_green_frac_w,building_type) 3708 surf_usm_v(l)%frac(2,m) = building_pars(ind_win_frac,building_type) 3709 surf_usm_v(l)%lai(m) = building_pars(ind_lai_w,building_type) 3710 3711 surf_usm_v(l)%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,building_type) 3712 surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type) 3713 surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type) 3714 surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type) 3715 3716 surf_usm_v(l)%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,building_type) 3717 surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,building_type) 3718 surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,building_type) 3719 surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,building_type) 3720 3721 surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,building_type) 3722 surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,building_type) 3723 surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,building_type) 3724 surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,building_type) 3725 3726 surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars(ind_tc1,building_type) 3727 surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 3728 surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type) 3729 surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type) 3730 3731 surf_usm_v(l)%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,building_type) 3732 surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 3733 surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,building_type) 3734 surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,building_type) 3735 3736 surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,building_type) 3737 surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 3738 surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,building_type) 3739 surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,building_type) 3740 3741 surf_usm_v(l)%target_temp_summer(m) = building_pars(12,building_type) 3742 surf_usm_v(l)%target_temp_winter(m) = building_pars(13,building_type) 3743 ! 3744 !-- emissivity of wall-, green- and window fraction 3745 surf_usm_v(l)%emissivity(0,m) = building_pars(ind_emis_wall,building_type) 3746 surf_usm_v(l)%emissivity(1,m) = building_pars(ind_emis_green,building_type) 3747 surf_usm_v(l)%emissivity(2,m) = building_pars(ind_emis_win,building_type) 3748 3749 surf_usm_v(l)%transmissivity(m) = building_pars(ind_trans,building_type) 3750 3751 surf_usm_v(l)%z0(m) = building_pars(ind_z0,building_type) 3752 surf_usm_v(l)%z0h(m) = building_pars(ind_z0qh,building_type) 3753 surf_usm_v(l)%z0q(m) = building_pars(ind_z0qh,building_type) 3754 3755 surf_usm_v(l)%albedo_type(0,m) = INT( building_pars(ind_alb_wall,building_type) ) 3756 surf_usm_v(l)%albedo_type(1,m) = INT( building_pars(ind_alb_green,building_type) ) 3757 surf_usm_v(l)%albedo_type(2,m) = INT( building_pars(ind_alb_win,building_type) ) 3758 3759 surf_usm_v(l)%zw(nzb_wall,m) = building_pars(ind_thick_1,building_type) 3760 surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars(ind_thick_2,building_type) 3761 surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars(ind_thick_3,building_type) 3762 surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars(ind_thick_4,building_type) 3763 3764 surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars(ind_thick_1,building_type) 3765 surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,building_type) 3766 surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,building_type) 3767 surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,building_type) 3768 3769 surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars(ind_thick_1,building_type) 3770 surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,building_type) 3771 surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,building_type) 3772 surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,building_type) 3773 3774 surf_usm_v(l)%c_surface(m) = building_pars(45,building_type) 3775 surf_usm_v(l)%lambda_surf(m) = building_pars(46,building_type) 3776 surf_usm_v(l)%c_surface_green(m) = building_pars(45,building_type) 3777 surf_usm_v(l)%lambda_surf_green(m) = building_pars(46,building_type) 3778 surf_usm_v(l)%c_surface_window(m) = building_pars(45,building_type) 3779 surf_usm_v(l)%lambda_surf_window(m) = building_pars(46,building_type) 3780 3781 ENDDO 3782 ENDDO 3783 ! 3784 !-- Level 2 - initialization via building type read from file 3785 IF ( building_type_f%from_file ) THEN 3786 DO m = 1, surf_usm_h%ns 3787 i = surf_usm_h%i(m) 3788 j = surf_usm_h%j(m) 3789 ! 3790 !-- For the moment, limit building type to 6 (to overcome errors in input file). 3791 st = building_type_f%var(j,i) 3792 IF ( st /= building_type_f%fill ) THEN 3793 3794 ! 3795 !-- In order to distinguish between ground floor level and 3796 !-- above-ground-floor level surfaces, set input indices. 3797 ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, & 3798 surf_usm_h%ground_level(m) ) 3799 ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, & 3800 surf_usm_h%ground_level(m) ) 3801 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, & 3802 surf_usm_h%ground_level(m) ) 3803 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, & 3804 surf_usm_h%ground_level(m) ) 3805 ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, & 3806 surf_usm_h%ground_level(m) ) 3807 ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, & 3808 surf_usm_h%ground_level(m) ) 3809 ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, & 3810 surf_usm_h%ground_level(m) ) 3811 ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, & 3812 surf_usm_h%ground_level(m) ) 3813 ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, & 3814 surf_usm_h%ground_level(m) ) 3815 ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, & 3816 surf_usm_h%ground_level(m) ) 3817 ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, & 3818 surf_usm_h%ground_level(m) ) 3819 ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, & 3820 surf_usm_h%ground_level(m) ) 3821 ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, & 3822 surf_usm_h%ground_level(m) ) 3823 ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, & 3824 surf_usm_h%ground_level(m) ) 3825 ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, & 3826 surf_usm_h%ground_level(m) ) 3827 ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, & 3828 surf_usm_h%ground_level(m) ) 3829 ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, & 3830 surf_usm_h%ground_level(m) ) 3831 ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, & 3832 surf_usm_h%ground_level(m) ) 3833 3834 ! 3835 !-- Initialize relatvie wall- (0), green- (1) and window (2) fractions 3836 surf_usm_h%frac(0,m) = building_pars(ind_wall_frac,st) 3837 surf_usm_h%frac(1,m) = building_pars(ind_green_frac_r,st) 3838 surf_usm_h%frac(2,m) = building_pars(ind_win_frac,st) 3839 surf_usm_h%lai(m) = building_pars(ind_green_frac_r,st) 3840 3841 surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,st) 3842 surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st) 3843 surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st) 3844 surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st) 3845 surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1,st) 3846 surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st) 3847 surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st) 3848 surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st) 3849 3850 surf_usm_h%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,st) 3851 surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,st) 3852 surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,st) 3853 surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,st) 3854 surf_usm_h%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,st) 3855 surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,st) 3856 surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,st) 3857 surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,st) 3858 3859 surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,st) 3860 surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,st) 3861 surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,st) 3862 surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,st) 3863 surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,st) 3864 surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,st) 3865 surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,st) 3866 surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,st) 3867 3868 surf_usm_h%target_temp_summer(m) = building_pars(12,st) 3869 surf_usm_h%target_temp_winter(m) = building_pars(13,st) 3870 ! 3871 !-- emissivity of wall-, green- and window fraction 3872 surf_usm_h%emissivity(0,m) = building_pars(ind_emis_wall,st) 3873 surf_usm_h%emissivity(1,m) = building_pars(ind_emis_green,st) 3874 surf_usm_h%emissivity(2,m) = building_pars(ind_emis_win,st) 3875 3876 surf_usm_h%transmissivity(m) = building_pars(ind_trans,st) 3877 3878 surf_usm_h%z0(m) = building_pars(ind_z0,st) 3879 surf_usm_h%z0h(m) = building_pars(ind_z0qh,st) 3880 surf_usm_h%z0q(m) = building_pars(ind_z0qh,st) 3881 ! 3882 !-- albedo type for wall fraction, green fraction, window fraction 3883 surf_usm_h%albedo_type(0,m) = INT( building_pars(ind_alb_wall,st) ) 3884 surf_usm_h%albedo_type(1,m) = INT( building_pars(ind_alb_green,st) ) 3885 surf_usm_h%albedo_type(2,m) = INT( building_pars(ind_alb_win,st) ) 3886 3887 surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1,st) 3888 surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2,st) 3889 surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3,st) 3890 surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4,st) 3891 3892 surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1,st) 3893 surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,st) 3894 surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,st) 3895 surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,st) 3896 3897 surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1,st) 3898 surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,st) 3899 surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,st) 3900 surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,st) 3901 3902 surf_usm_h%c_surface(m) = building_pars(45,st) 3903 surf_usm_h%lambda_surf(m) = building_pars(46,st) 3904 surf_usm_h%c_surface_green(m) = building_pars(45,st) 3905 surf_usm_h%lambda_surf_green(m) = building_pars(46,st) 3906 surf_usm_h%c_surface_window(m) = building_pars(45,st) 3907 surf_usm_h%lambda_surf_window(m) = building_pars(46,st) 3908 3909 ENDIF 3910 ENDDO 3911 3912 DO l = 0, 3 3913 DO m = 1, surf_usm_v(l)%ns 3914 i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff 3915 j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff 3916 ! 3917 !-- For the moment, limit building type to 6 (to overcome errors in input file). 3918 3919 st = building_type_f%var(j,i) 3920 IF ( st /= building_type_f%fill ) THEN 3921 3922 ! 3923 !-- In order to distinguish between ground floor level and 3924 !-- above-ground-floor level surfaces, set input indices. 3925 ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, & 3926 surf_usm_v(l)%ground_level(m) ) 3927 ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, & 3928 surf_usm_v(l)%ground_level(m) ) 3929 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, & 3930 surf_usm_v(l)%ground_level(m) ) 3931 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, & 3932 surf_usm_v(l)%ground_level(m) ) 3933 ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, & 3934 surf_usm_v(l)%ground_level(m) ) 3935 ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, & 3936 surf_usm_v(l)%ground_level(m) ) 3937 ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, & 3938 surf_usm_v(l)%ground_level(m) ) 3939 ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, & 3940 surf_usm_v(l)%ground_level(m) ) 3941 ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, & 3942 surf_usm_v(l)%ground_level(m) ) 3943 ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, & 3944 surf_usm_v(l)%ground_level(m) ) 3945 ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, & 3946 surf_usm_v(l)%ground_level(m) ) 3947 ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, & 3948 surf_usm_v(l)%ground_level(m) ) 3949 ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, & 3950 surf_usm_v(l)%ground_level(m) ) 3951 ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, & 3952 surf_usm_v(l)%ground_level(m) ) 3953 ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, & 3954 surf_usm_v(l)%ground_level(m) ) 3955 ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, & 3956 surf_usm_v(l)%ground_level(m) ) 3957 ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, & 3958 surf_usm_v(l)%ground_level(m) ) 3959 ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, & 3960 surf_usm_v(l)%ground_level(m) ) 3961 3962 ! 3963 !-- Initialize relatvie wall- (0), green- (1) and window (2) fractions 3964 surf_usm_v(l)%frac(0,m) = building_pars(ind_wall_frac,st) 3965 surf_usm_v(l)%frac(1,m) = building_pars(ind_green_frac_w,st) 3966 surf_usm_v(l)%frac(2,m) = building_pars(ind_win_frac,st) 3967 surf_usm_v(l)%lai(m) = building_pars(ind_lai_w,st) 3968 3969 surf_usm_v(l)%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,st) 3970 surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st) 3971 surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st) 3972 surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st) 3973 3974 surf_usm_v(l)%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,st) 3975 surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,st) 3976 surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,st) 3977 surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,st) 3978 3979 surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,st) 3980 surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,st) 3981 surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,st) 3982 surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,st) 3983 3984 surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars(ind_tc1,st) 3985 surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st) 3986 surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st) 3987 surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st) 3988 3989 surf_usm_v(l)%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,st) 3990 surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,st) 3991 surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,st) 3992 surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,st) 3993 3994 surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,st) 3995 surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,st) 3996 surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,st) 3997 surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,st) 3998 3999 surf_usm_v(l)%target_temp_summer(m) = building_pars(12,st) 4000 surf_usm_v(l)%target_temp_winter(m) = building_pars(13,st) 4001 ! 4002 !-- emissivity of wall-, green- and window fraction 4003 surf_usm_v(l)%emissivity(0,m) = building_pars(ind_emis_wall,st) 4004 surf_usm_v(l)%emissivity(1,m) = building_pars(ind_emis_green,st) 4005 surf_usm_v(l)%emissivity(2,m) = building_pars(ind_emis_win,st) 4006 4007 surf_usm_v(l)%transmissivity(m) = building_pars(ind_trans,st) 4008 4009 surf_usm_v(l)%z0(m) = building_pars(ind_z0,st) 4010 surf_usm_v(l)%z0h(m) = building_pars(ind_z0qh,st) 4011 surf_usm_v(l)%z0q(m) = building_pars(ind_z0qh,st) 4012 4013 surf_usm_v(l)%albedo_type(0,m) = INT( building_pars(ind_alb_wall,st) ) 4014 surf_usm_v(l)%albedo_type(1,m) = INT( building_pars(ind_alb_green,st) ) 4015 surf_usm_v(l)%albedo_type(2,m) = INT( building_pars(ind_alb_win,st) ) 4016 4017 surf_usm_v(l)%zw(nzb_wall,m) = building_pars(ind_thick_1,st) 4018 surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars(ind_thick_2,st) 4019 surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars(ind_thick_3,st) 4020 surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars(ind_thick_4,st) 4021 4022 surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars(ind_thick_1,st) 4023 surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,st) 4024 surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,st) 4025 surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,st) 4026 4027 surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars(ind_thick_1,st) 4028 surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,st) 4029 surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,st) 4030 surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,st) 4031 4032 surf_usm_v(l)%c_surface(m) = building_pars(45,st) 4033 surf_usm_v(l)%lambda_surf(m) = building_pars(46,st) 4034 surf_usm_v(l)%c_surface_green(m) = building_pars(45,st) 4035 surf_usm_v(l)%lambda_surf_green(m) = building_pars(46,st) 4036 surf_usm_v(l)%c_surface_window(m) = building_pars(45,st) 4037 surf_usm_v(l)%lambda_surf_window(m) = building_pars(46,st) 4038 4039 4040 ENDIF 4041 ENDDO 4042 ENDDO 4043 ENDIF 4044 4045 ! 4046 !-- Level 3 - initialization via building_pars read from file 4047 IF ( building_pars_f%from_file ) THEN 4048 DO m = 1, surf_usm_h%ns 4049 i = surf_usm_h%i(m) 4050 j = surf_usm_h%j(m) 4051 4052 ! 4053 !-- In order to distinguish between ground floor level and 4054 !-- above-ground-floor level surfaces, set input indices. 4055 ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, & 4056 surf_usm_h%ground_level(m) ) 4057 ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, & 4058 surf_usm_h%ground_level(m) ) 4059 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, & 4060 surf_usm_h%ground_level(m) ) 4061 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, & 4062 surf_usm_h%ground_level(m) ) 4063 ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, & 4064 surf_usm_h%ground_level(m) ) 4065 ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, & 4066 surf_usm_h%ground_level(m) ) 4067 ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, & 4068 surf_usm_h%ground_level(m) ) 4069 ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, & 4070 surf_usm_h%ground_level(m) ) 4071 ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, & 4072 surf_usm_h%ground_level(m) ) 4073 ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, & 4074 surf_usm_h%ground_level(m) ) 4075 ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, & 4076 surf_usm_h%ground_level(m) ) 4077 ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, & 4078 surf_usm_h%ground_level(m) ) 4079 ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, & 4080 surf_usm_h%ground_level(m) ) 4081 ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, & 4082 surf_usm_h%ground_level(m) ) 4083 ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, & 4084 surf_usm_h%ground_level(m) ) 4085 ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, & 4086 surf_usm_h%ground_level(m) ) 4087 ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, & 4088 surf_usm_h%ground_level(m) ) 4089 ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, & 4090 surf_usm_h%ground_level(m) ) 4091 4092 ! 4093 !-- Initialize relatvie wall- (0), green- (1) and window (2) fractions 4094 IF ( building_pars_f%pars_xy(ind_wall_frac,j,i) /= building_pars_f%fill ) & 4095 surf_usm_h%frac(0,m) = building_pars_f%pars_xy(ind_wall_frac,j,i) 4096 IF ( building_pars_f%pars_xy(ind_green_frac_r,j,i) /= building_pars_f%fill ) & 4097 surf_usm_h%frac(1,m) = building_pars_f%pars_xy(ind_green_frac_r,j,i) 4098 IF ( building_pars_f%pars_xy(ind_win_frac,j,i) /= building_pars_f%fill ) & 4099 surf_usm_h%frac(2,m) = building_pars_f%pars_xy(ind_win_frac,j,i) 4100 4101 4102 IF ( building_pars_f%pars_xy(ind_lai_r,j,i) /= building_pars_f%fill ) & 4103 surf_usm_h%lai(m) = building_pars_f%pars_xy(ind_lai_r,j,i) 4104 4105 IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN 4106 surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4107 surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4108 ENDIF 4109 IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) & 4110 surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i) 4111 IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) & 4112 surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i) 4113 IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN 4114 surf_usm_h%rho_c_green(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4115 surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4116 ENDIF 4117 IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) & 4118 surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i) 4119 IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) & 4120 surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i) 4121 IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN 4122 surf_usm_h%rho_c_window(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4123 surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4124 ENDIF 4125 IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) & 4126 surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i) 4127 IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) & 4128 surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i) 4129 4130 IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN 4131 surf_usm_h%lambda_h(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4132 surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4133 ENDIF 4134 IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) & 4135 surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i) 4136 IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) & 4137 surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i) 4138 IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN 4139 surf_usm_h%lambda_h_green(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4140 surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4141 ENDIF 4142 IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) & 4143 surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i) 4144 IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) & 4145 surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i) 4146 IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN 4147 surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4148 surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4149 ENDIF 4150 IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) & 4151 surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i) 4152 IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) & 4153 surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i) 4154 4155 IF ( building_pars_f%pars_xy(12,j,i) /= building_pars_f%fill ) & 4156 surf_usm_h%target_temp_summer(m) = building_pars_f%pars_xy(12,j,i) 4157 IF ( building_pars_f%pars_xy(13,j,i) /= building_pars_f%fill ) & 4158 surf_usm_h%target_temp_winter(m) = building_pars_f%pars_xy(13,j,i) 4159 4160 IF ( building_pars_f%pars_xy(ind_emis_wall,j,i) /= building_pars_f%fill ) & 4161 surf_usm_h%emissivity(0,m) = building_pars_f%pars_xy(ind_emis_wall,j,i) 4162 IF ( building_pars_f%pars_xy(ind_emis_green,j,i) /= building_pars_f%fill )& 4163 surf_usm_h%emissivity(1,m) = building_pars_f%pars_xy(ind_emis_green,j,i) 4164 IF ( building_pars_f%pars_xy(ind_emis_win,j,i) /= building_pars_f%fill ) & 4165 surf_usm_h%emissivity(2,m) = building_pars_f%pars_xy(ind_emis_win,j,i) 4166 4167 IF ( building_pars_f%pars_xy(ind_trans,j,i) /= building_pars_f%fill ) & 4168 surf_usm_h%transmissivity(m) = building_pars_f%pars_xy(ind_trans,j,i) 4169 4170 IF ( building_pars_f%pars_xy(ind_z0,j,i) /= building_pars_f%fill ) & 4171 surf_usm_h%z0(m) = building_pars_f%pars_xy(ind_z0,j,i) 4172 IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill ) & 4173 surf_usm_h%z0h(m) = building_pars_f%pars_xy(ind_z0qh,j,i) 4174 IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill ) & 4175 surf_usm_h%z0q(m) = building_pars_f%pars_xy(ind_z0qh,j,i) 4176 4177 IF ( building_pars_f%pars_xy(ind_alb_wall,j,i) /= building_pars_f%fill ) & 4178 surf_usm_h%albedo_type(0,m) = building_pars_f%pars_xy(ind_alb_wall,j,i) 4179 IF ( building_pars_f%pars_xy(ind_alb_green,j,i) /= building_pars_f%fill ) & 4180 surf_usm_h%albedo_type(1,m) = building_pars_f%pars_xy(ind_alb_green,j,i) 4181 IF ( building_pars_f%pars_xy(ind_alb_win,j,i) /= building_pars_f%fill ) & 4182 surf_usm_h%albedo_type(2,m) = building_pars_f%pars_xy(ind_alb_win,j,i) 4183 4184 IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) & 4185 surf_usm_h%zw(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i) 4186 IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) & 4187 surf_usm_h%zw(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i) 4188 IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) & 4189 surf_usm_h%zw(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i) 4190 IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) & 4191 surf_usm_h%zw(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i) 4192 IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) & 4193 surf_usm_h%zw_green(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i) 4194 IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) & 4195 surf_usm_h%zw_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i) 4196 IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) & 4197 surf_usm_h%zw_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i) 4198 IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) & 4199 surf_usm_h%zw_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i) 4200 IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) & 4201 surf_usm_h%zw_window(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i) 4202 IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) & 4203 surf_usm_h%zw_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i) 4204 IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) & 4205 surf_usm_h%zw_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i) 4206 IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) & 4207 surf_usm_h%zw_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i) 4208 4209 IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) & 4210 surf_usm_h%c_surface(m) = building_pars_f%pars_xy(45,j,i) 4211 IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) & 4212 surf_usm_h%lambda_surf(m) = building_pars_f%pars_xy(46,j,i) 4213 IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) & 4214 surf_usm_h%c_surface_green(m) = building_pars_f%pars_xy(45,j,i) 4215 IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) & 4216 surf_usm_h%lambda_surf_green(m) = building_pars_f%pars_xy(46,j,i) 4217 IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) & 4218 surf_usm_h%c_surface_window(m) = building_pars_f%pars_xy(45,j,i) 4219 IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) & 4220 surf_usm_h%lambda_surf_window(m) = building_pars_f%pars_xy(46,j,i) 4221 ENDDO 4222 4223 4224 4225 DO l = 0, 3 4226 DO m = 1, surf_usm_v(l)%ns 4227 i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff 4228 j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff 4229 4230 ! 4231 !-- In order to distinguish between ground floor level and 4232 !-- above-ground-floor level surfaces, set input indices. 4233 ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, & 4234 surf_usm_v(l)%ground_level(m) ) 4235 ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, & 4236 surf_usm_v(l)%ground_level(m) ) 4237 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, & 4238 surf_usm_v(l)%ground_level(m) ) 4239 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, & 4240 surf_usm_v(l)%ground_level(m) ) 4241 ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, & 4242 surf_usm_v(l)%ground_level(m) ) 4243 ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, & 4244 surf_usm_v(l)%ground_level(m) ) 4245 ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, & 4246 surf_usm_v(l)%ground_level(m) ) 4247 ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, & 4248 surf_usm_v(l)%ground_level(m) ) 4249 ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, & 4250 surf_usm_v(l)%ground_level(m) ) 4251 ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, & 4252 surf_usm_v(l)%ground_level(m) ) 4253 ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, & 4254 surf_usm_v(l)%ground_level(m) ) 4255 ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, & 4256 surf_usm_v(l)%ground_level(m) ) 4257 ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, & 4258 surf_usm_v(l)%ground_level(m) ) 4259 ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, & 4260 surf_usm_v(l)%ground_level(m) ) 4261 ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, & 4262 surf_usm_v(l)%ground_level(m) ) 4263 ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, & 4264 surf_usm_v(l)%ground_level(m) ) 4265 ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, & 4266 surf_usm_v(l)%ground_level(m) ) 4267 ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, & 4268 surf_usm_v(l)%ground_level(m) ) 4269 4270 ! 4271 !-- Initialize relatvie wall- (0), green- (1) and window (2) fractions 4272 IF ( building_pars_f%pars_xy(ind_wall_frac,j,i) /= & 4273 building_pars_f%fill ) & 4274 surf_usm_v(l)%frac(0,m) = building_pars_f%pars_xy(ind_wall_frac,j,i) 4275 IF ( building_pars_f%pars_xy(ind_green_frac_w,j,i) /= & 4276 building_pars_f%fill ) & 4277 surf_usm_v(l)%frac(1,m) = building_pars_f%pars_xy(ind_green_frac_w,j,i) 4278 IF ( building_pars_f%pars_xy(ind_win_frac,j,i) /= & 4279 building_pars_f%fill ) & 4280 surf_usm_v(l)%frac(2,m) = building_pars_f%pars_xy(ind_win_frac,j,i) 4281 4282 IF ( building_pars_f%pars_xy(ind_lai_w,j,i) /= building_pars_f%fill ) & 4283 surf_usm_v(l)%lai(m) = building_pars_f%pars_xy(ind_lai_w,j,i) 4284 4285 IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) & 4286 THEN 4287 surf_usm_v(l)%rho_c_wall(nzb_wall,m) = & 4288 building_pars_f%pars_xy(ind_hc1,j,i) 4289 surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = & 4290 building_pars_f%pars_xy(ind_hc1,j,i) 4291 ENDIF 4292 IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) & 4293 surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = & 4294 building_pars_f%pars_xy(ind_hc2,j,i) 4295 IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) & 4296 surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = & 4297 building_pars_f%pars_xy(ind_hc3,j,i) 4298 IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN 4299 surf_usm_v(l)%rho_c_green(nzb_wall,m) = & 4300 building_pars_f%pars_xy(ind_hc1,j,i) 4301 surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = & 4302 building_pars_f%pars_xy(ind_hc1,j,i) 4303 ENDIF 4304 IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) & 4305 surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i) 4306 IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) & 4307 surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i) 4308 IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN 4309 surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4310 surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i) 4311 ENDIF 4312 IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) & 4313 surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i) 4314 IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) & 4315 surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i) 4316 4317 IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN 4318 surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4319 surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4320 ENDIF 4321 IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) & 4322 surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i) 4323 IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) & 4324 surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i) 4325 IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN 4326 surf_usm_v(l)%lambda_h_green(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4327 surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4328 ENDIF 4329 IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) & 4330 surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i) 4331 IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) & 4332 surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i) 4333 IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN 4334 surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4335 surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i) 4336 ENDIF 4337 IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) & 4338 surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i) 4339 IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) & 4340 surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i) 4341 4342 IF ( building_pars_f%pars_xy(12,j,i) /= building_pars_f%fill ) & 4343 surf_usm_v(l)%target_temp_summer(m) = building_pars_f%pars_xy(12,j,i) 4344 IF ( building_pars_f%pars_xy(13,j,i) /= building_pars_f%fill ) & 4345 surf_usm_v(l)%target_temp_winter(m) = building_pars_f%pars_xy(13,j,i) 4346 4347 IF ( building_pars_f%pars_xy(ind_emis_wall,j,i) /= building_pars_f%fill ) & 4348 surf_usm_v(l)%emissivity(0,m) = building_pars_f%pars_xy(ind_emis_wall,j,i) 4349 IF ( building_pars_f%pars_xy(ind_emis_green,j,i) /= building_pars_f%fill )& 4350 surf_usm_v(l)%emissivity(1,m) = building_pars_f%pars_xy(ind_emis_green,j,i) 4351 IF ( building_pars_f%pars_xy(ind_emis_win,j,i) /= building_pars_f%fill ) & 4352 surf_usm_v(l)%emissivity(2,m) = building_pars_f%pars_xy(ind_emis_win,j,i) 4353 4354 IF ( building_pars_f%pars_xy(ind_trans,j,i) /= building_pars_f%fill ) & 4355 surf_usm_v(l)%transmissivity(m) = building_pars_f%pars_xy(ind_trans,j,i) 4356 4357 IF ( building_pars_f%pars_xy(ind_z0,j,i) /= building_pars_f%fill ) & 4358 surf_usm_v(l)%z0(m) = building_pars_f%pars_xy(ind_z0,j,i) 4359 IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill ) & 4360 surf_usm_v(l)%z0h(m) = building_pars_f%pars_xy(ind_z0qh,j,i) 4361 IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill ) & 4362 surf_usm_v(l)%z0q(m) = building_pars_f%pars_xy(ind_z0qh,j,i) 4363 4364 IF ( building_pars_f%pars_xy(ind_alb_wall,j,i) /= building_pars_f%fill ) & 4365 surf_usm_v(l)%albedo_type(0,m) = building_pars_f%pars_xy(ind_alb_wall,j,i) 4366 IF ( building_pars_f%pars_xy(ind_alb_green,j,i) /= building_pars_f%fill ) & 4367 surf_usm_v(l)%albedo_type(1,m) = building_pars_f%pars_xy(ind_alb_green,j,i) 4368 IF ( building_pars_f%pars_xy(ind_alb_win,j,i) /= building_pars_f%fill ) & 4369 surf_usm_v(l)%albedo_type(2,m) = building_pars_f%pars_xy(ind_alb_win,j,i) 4370 4371 IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) & 4372 surf_usm_v(l)%zw(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i) 4373 IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) & 4374 surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i) 4375 IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) & 4376 surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i) 4377 IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) & 4378 surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i) 4379 IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) & 4380 surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i) 4381 IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) & 4382 surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i) 4383 IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) & 4384 surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i) 4385 IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) & 4386 surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i) 4387 IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) & 4388 surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i) 4389 IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) & 4390 surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i) 4391 IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) & 4392 surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i) 4393 IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) & 4394 surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i) 4395 4396 IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) & 4397 surf_usm_v(l)%c_surface(m) = building_pars_f%pars_xy(45,j,i) 4398 IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) & 4399 surf_usm_v(l)%lambda_surf(m) = building_pars_f%pars_xy(46,j,i) 4400 IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) & 4401 surf_usm_v(l)%c_surface_green(m) = building_pars_f%pars_xy(45,j,i) 4402 IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) & 4403 surf_usm_v(l)%lambda_surf_green(m) = building_pars_f%pars_xy(46,j,i) 4404 IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) & 4405 surf_usm_v(l)%c_surface_window(m) = building_pars_f%pars_xy(45,j,i) 4406 IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) & 4407 surf_usm_v(l)%lambda_surf_window(m) = building_pars_f%pars_xy(46,j,i) 4408 4409 ENDDO 4410 ENDDO 4411 ENDIF 4412 ! 4413 !-- Read the surface_types array. 4414 !-- Please note, here also initialization of surface attributes is done as 4415 !-- long as _urbsurf and _surfpar files are available. Values from above 4416 !-- will be overwritten. This might be removed later, but is still in the 4417 !-- code to enable compatibility with older model version. 3194 4418 CALL usm_read_urban_surface_types() 3195 4419 … … 3197 4421 CALL usm_init_material_model() 3198 4422 4423 !-- init anthropogenic sources of heat 3199 4424 IF ( usm_anthropogenic_heat ) THEN 3200 4425 !-- init anthropogenic sources of heat (from transportation for now) 3201 4426 CALL usm_read_anthropogenic_heat() 3202 4427 ENDIF 3203 3204 IF ( read_svf_on_init ) THEN 3205 !-- read svf, csf, svfsurf and csfsurf data from file 3206 CALL location_message( ' Start reading SVF from file', .TRUE. ) 3207 CALL usm_read_svf_from_file() 3208 CALL location_message( ' Reading SVF from file has finished', .TRUE. ) 3209 ELSE 3210 !-- calculate SFV and CSF 3211 CALL location_message( ' Start calculation of SVF', .TRUE. ) 3212 CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'start' ) 3213 CALL usm_calc_svf() 3214 CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'stop' ) 3215 CALL location_message( ' Calculation of SVF has finished', .TRUE. ) 3216 ENDIF 3217 3218 IF ( write_svf_on_init ) THEN 3219 !-- write svf, csf svfsurf and csfsurf data to file 3220 CALL location_message( ' Store SVF and CSF to file', .TRUE. ) 3221 CALL usm_write_svf_to_file() 3222 ENDIF 3223 4428 3224 4429 IF ( plant_canopy ) THEN 3225 !-- gridpcbl was only necessary for initialization 3226 DEALLOCATE( gridpcbl ) 3227 IF ( .NOT. ALLOCATED(pc_heating_rate) ) THEN 4430 4431 IF ( .NOT. ALLOCATED( pc_heating_rate) ) THEN 3228 4432 !-- then pc_heating_rate is allocated in init_plant_canopy 3229 4433 !-- in case of cthf /= 0 => we need to allocate it for our use here … … 3251 4455 3252 4456 t_surf_h(m) = pt(k,j,i) * exn 4457 t_surf_window_h(m) = pt(k,j,i) * exn 4458 t_surf_green_h(m) = pt(k,j,i) * exn 4459 t_surf_whole_h(m) = pt(k,j,i) * exn 3253 4460 ENDDO 3254 4461 ! … … 3261 4468 3262 4469 t_surf_v(l)%t(m) = pt(k,j,i) * exn 4470 t_surf_window_v(l)%t(m) = pt(k,j,i) * exn 4471 t_surf_green_v(l)%t(m) = pt(k,j,i) * exn 4472 t_surf_whole_v(l)%t(m) = pt(k,j,i) * exn 3263 4473 ENDDO 3264 4474 ENDDO … … 3275 4485 IF ( surf_usm_h%isroof_surf(m) ) THEN 3276 4486 tin = roof_inner_temperature 4487 twin = window_inner_temperature 3277 4488 ! 3278 4489 !-- Normal land surface 3279 4490 ELSE 3280 4491 tin = soil_inner_temperature 4492 twin = window_inner_temperature 3281 4493 ENDIF 3282 4494 … … 3286 4498 3287 4499 t_wall_h(k,m) = ( 1.0_wp - c ) * t_surf_h(m) + c * tin 4500 t_window_h(k,m) = ( 1.0_wp - c ) * t_surf_window_h(m) + c * twin 4501 t_green_h(k,m) = t_surf_h(m) 3288 4502 ENDDO 3289 4503 ENDDO … … 3295 4509 !-- Inner wall 3296 4510 tin = wall_inner_temperature 4511 twin = window_inner_temperature 3297 4512 3298 4513 DO k = nzb_wall, nzt_wall+1 3299 4514 c = REAL( k - nzb_wall, wp ) / & 3300 4515 REAL( nzt_wall + 1 - nzb_wall , wp ) 3301 3302 t_w all_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) + &3303 c * tin4516 t_wall_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) + c * tin 4517 t_window_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_window_v(l)%t(m) + c * twin 4518 t_green_v(l)%t(k,m) = t_surf_v(l)%t(m) 3304 4519 ENDDO 3305 4520 ENDDO … … 3315 4530 t_surf_h_p = t_surf_h 3316 4531 t_surf_v_p = t_surf_v 4532 t_surf_window_h_p = t_surf_window_h 4533 t_surf_window_v_p = t_surf_window_v 4534 t_surf_whole_h_p = t_surf_whole_h 4535 t_surf_whole_v_p = t_surf_whole_v 4536 t_surf_green_h_p = t_surf_green_h 4537 t_surf_green_v_p = t_surf_green_v 4538 t_surf_10cm_h_p = t_surf_10cm_h 4539 t_surf_10cm_v_p = t_surf_10cm_v 3317 4540 3318 4541 t_wall_h_p = t_wall_h 3319 4542 t_wall_v_p = t_wall_v 3320 3321 !-- Adjust radiative fluxes for urban surface at model start 3322 CALL usm_radiation 4543 t_window_h_p = t_window_h 4544 t_window_v_p = t_window_v 4545 t_green_h_p = t_green_h 4546 t_green_v_p = t_green_v 3323 4547 3324 4548 CALL cpu_log( log_point_s(78), 'usm_init', 'stop' ) 3325 4549 3326 3327 4550 END SUBROUTINE usm_init_urban_surface 3328 4551 … … 3342 4565 INTEGER(iwp) :: i,j,k,l,kw, m !< running indices 3343 4566 3344 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend !< tendency4567 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend, wintend !< tendency 3345 4568 3346 4569 ! … … 3355 4578 !-- prognostic equation for ground/roof temperature t_wall_h 3356 4579 wtend(:) = 0.0_wp 3357 wtend(nzb_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m)) * & 3358 ( surf_usm_h%lambda_h(nzb_wall,m) * & 3359 ( t_wall_h(nzb_wall+1,m) & 3360 - t_wall_h(nzb_wall,m) ) * & 3361 surf_usm_h%ddz_wall(nzb_wall+1,m) & 3362 + surf_usm_h%wghf_eb(m) ) * & 3363 surf_usm_h%ddz_wall_stag(nzb_wall,m) 3364 3365 DO kw = nzb_wall+1, nzt_wall 3366 wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m)) & 3367 * ( surf_usm_h%lambda_h(kw,m) & 3368 * ( t_wall_h(kw+1,m) - t_wall_h(kw,m) ) & 3369 * surf_usm_h%ddz_wall(kw+1,m) & 3370 - surf_usm_h%lambda_h(kw-1,m) & 3371 * ( t_wall_h(kw,m) - t_wall_h(kw-1,m) ) & 3372 * surf_usm_h%ddz_wall(kw,m) & 3373 ) * surf_usm_h%ddz_wall_stag(kw,m) 3374 ENDDO 4580 wtend(nzb_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m)) * & 4581 ( surf_usm_h%lambda_h(nzb_wall,m) * & 4582 ( t_wall_h(nzb_wall+1,m) & 4583 - t_wall_h(nzb_wall,m) ) * & 4584 surf_usm_h%ddz_wall(nzb_wall+1,m) & 4585 + surf_usm_h%frac(0,m) & 4586 / (surf_usm_h%frac(0,m) & 4587 + surf_usm_h%frac(1,m) ) & 4588 * surf_usm_h%wghf_eb(m) & 4589 - surf_usm_h%frac(1,m) & 4590 / (surf_usm_h%frac(0,m) & 4591 + surf_usm_h%frac(1,m) ) & 4592 * ( surf_usm_h%lambda_h_green(nzt_wall,m) & 4593 * surf_usm_h%ddz_green(nzt_wall,m) & 4594 + surf_usm_h%lambda_h(nzb_wall,m) & 4595 * surf_usm_h%ddz_wall(nzb_wall,m) ) & 4596 / ( surf_usm_h%ddz_green(nzt_wall,m) & 4597 + surf_usm_h%ddz_wall(nzb_wall,m) ) & 4598 * ( t_wall_h(nzb_wall,m) & 4599 - t_green_h(nzt_wall,m) ) ) * & 4600 surf_usm_h%ddz_wall_stag(nzb_wall,m) 4601 4602 !dummy value for testing 4603 surf_usm_h%iwghf_eb(m) = 0. 4604 4605 IF ( indoor_model ) then 4606 DO kw = nzb_wall+1, nzt_wall-1 4607 wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m)) & 4608 * ( surf_usm_h%lambda_h(kw,m) & 4609 * ( t_wall_h(kw+1,m) - t_wall_h(kw,m) ) & 4610 * surf_usm_h%ddz_wall(kw+1,m) & 4611 - surf_usm_h%lambda_h(kw-1,m) & 4612 * ( t_wall_h(kw,m) - t_wall_h(kw-1,m) ) & 4613 * surf_usm_h%ddz_wall(kw,m) & 4614 ) * surf_usm_h%ddz_wall_stag(kw,m) 4615 ENDDO 4616 wtend(nzt_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzt_wall,m)) * & 4617 ( surf_usm_h%lambda_h(nzt_wall-1,m) * & 4618 ( t_wall_h(nzt_wall,m) & 4619 - t_wall_h(nzt_wall-1,m) ) * & 4620 surf_usm_h%ddz_wall(nzt_wall,m) & 4621 + surf_usm_h%iwghf_eb(m) ) * & 4622 surf_usm_h%ddz_wall_stag(nzt_wall,m) 4623 ELSE 4624 DO kw = nzb_wall+1, nzt_wall 4625 wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m)) &