Home

Context Navigation

urban_surface_mod.f90 @ 2544

Last change on this file since 2544 was 2544, checked in by maronga, 5 years ago
introduced new module date_and_time_mod
Property svn:keywords set to `Id`
File size: 264.8 KB

Line
1	!> @file urban_surface_mod.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2015-2017 Czech Technical University in Prague
18	! Copyright 1997-2017 Leibniz Universitaet Hannover
19	!--------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: urban_surface_mod.f90 2544 2017-10-13 18:09:32Z maronga $
28	! Date and time quantities are now read from date_and_time_mod. Solar constant is
29	! read from radiation_model_mod
30	!
31	! 2516 2017-10-04 11:03:04Z suehring
32	! Remove tabs
33	!
34	! 2514 2017-10-04 09:52:37Z suehring
35	! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
36	! no output of ghost layer data
37	!
38	! 2350 2017-08-15 11:48:26Z kanani
39	! Bugfix and error message for nopointer version.
40	! Additional "! defined(__nopointer)" as workaround to enable compilation of
41	! nopointer version.
42	!
43	! 2318 2017-07-20 17:27:44Z suehring
44	! Get topography top index via Function call
45	!
46	! 2317 2017-07-20 17:27:19Z suehring
47	! Bugfix: adjust output of shf. Added support for spinups
48	!
49	! 2287 2017-06-15 16:46:30Z suehring
50	! Bugfix in determination topography-top index
51	!
52	! 2269 2017-06-09 11:57:32Z suehring
53	! Enable restart runs with different number of PEs
54	! Bugfixes nopointer branch
55	!
56	! 2258 2017-06-08 07:55:13Z suehring
57	! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
58	!
59	! 2233 2017-05-30 18:08:54Z suehring
60	!
61	! 2232 2017-05-30 17:47:52Z suehring
62	! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
63	! insteat, heat fluxes are directly applied in diffusion_s.
64	!
65	! 2213 2017-04-24 15:10:35Z kanani
66	! Removal of output quantities usm_lad and usm_canopy_hr
67	!
68	! 2209 2017-04-19 09:34:46Z kanani
69	! cpp switch __mpi3 removed,
70	! minor formatting,
71	! small bugfix for division by zero (Krc)
72	!
73	! 2113 2017-01-12 13:40:46Z kanani
74	! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
75	!
76	! 2071 2016-11-17 11:22:14Z maronga
77	! Small bugfix (Resler)
78	!
79	! 2031 2016-10-21 15:11:58Z knoop
80	! renamed variable rho to rho_ocean
81	!
82	! 2024 2016-10-12 16:42:37Z kanani
83	! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
84	! optimization of MPI-RMA operations,
85	! declaration of pcbl as integer,
86	! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
87	! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
88	! use of new control parameter varnamelength,
89	! added output variables usm_rad_ressw, usm_rad_reslw,
90	! minor formatting changes,
91	! minor optimizations.
92	!
93	! 2011 2016-09-19 17:29:57Z kanani
94	! Major reformatting according to PALM coding standard (comments, blanks,
95	! alphabetical ordering, etc.),
96	! removed debug_prints,
97	! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
98	! defined in MODULE control_parameters (modules.f90) to avoid circular
99	! dependencies,
100	! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
101	!
102	! 2007 2016-08-24 15:47:17Z kanani
103	! Initial revision
104	!
105	!
106	! Description:
107	! ------------
108	! 2016/6/9 - Initial version of the USM (Urban Surface Model)
109	! authors: Jaroslav Resler, Pavel Krc
110	! (Czech Technical University in Prague and Institute of
111	! Computer Science of the Czech Academy of Sciences, Prague)
112	! with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
113	! partly inspired by PALM LSM (B. Maronga)
114	! parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
115	!> Module for Urban Surface Model (USM)
116	!> The module includes:
117	!> 1. radiation model with direct/diffuse radiation, shading, reflections
118	!> and integration with plant canopy
119	!> 2. wall and wall surface model
120	!> 3. surface layer energy balance
121	!> 4. anthropogenic heat (only from transportation so far)
122	!> 5. necessary auxiliary subroutines (reading inputs, writing outputs,
123	!> restart simulations, ...)
124	!> It also make use of standard radiation and integrates it into
125	!> urban surface model.
126	!>
127	!> Further work:
128	!> -------------
129	!> 1. Reduce number of shape view factors by merging factors for distant surfaces
130	!> under shallow angles. Idea: Iteratively select the smallest shape view
131	!> factor by value (among all sources and targets) which has a similarly
132	!> oriented source neighbor (or near enough) SVF and merge them by adding
133	!> value of the smaller SVF to the larger one and deleting the smaller one.
134	!> This will allow for better scaling at higher resolutions.
135	!>
136	!> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity
137	!> from surfaces that are visible from local surfaces (i.e. there is a SVF
138	!> where target is local). To do that, radiosity will be exchanged after each
139	!> reflection step using MPI_Alltoall instead of current MPI_Allgather.
140	!>
141	!> 3. Temporarily large values of surface heat flux can be observed, up to
142	!> 1.2 Km/s, which seem to be not realistic.
143	!>
144	!> @todo Revise flux conversion in energy-balance solver
145	!> @todo Bugfixing in nopointer branch
146	!> @todo Check optimizations for RMA operations
147	!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
148	!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
149	!> factor 3 between min and max time
150	!------------------------------------------------------------------------------!
151	MODULE urban_surface_mod
152
153	#if ! defined( __nopointer )
154	USE arrays_3d, &
155	ONLY: zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
156	#endif
157
158	USE cloud_parameters, &
159	ONLY: cp, r_d
160
161	USE constants, &
162	ONLY: pi
163
164	USE control_parameters, &
165	ONLY: coupling_start_time, dz, topography, dt_3d, &
166	intermediate_timestep_count, initializing_actions, &
167	intermediate_timestep_count_max, simulated_time, end_time, &
168	timestep_scheme, tsc, coupling_char, io_blocks, io_group, &
169	message_string, time_since_reference_point, surface_pressure, &
170	g, pt_surface, large_scale_forcing, lsf_surf, spinup, &
171	spinup_pt_mean, spinup_time, time_do3d, dt_do3d, &
172	average_count_3d, varnamelength, urban_surface
173
174	USE cpulog, &
175	ONLY: cpu_log, log_point, log_point_s
176
177	USE date_and_time_mod, &
178	ONLY: d_seconds_year, day_of_year_init, time_utc_init
179
180	USE grid_variables, &
181	ONLY: dx, dy, ddx, ddy, ddx2, ddy2
182
183	USE indices, &
184	ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, &
185	nysg, nzb, nzt, nbgp, wall_flags_0
186
187	USE, INTRINSIC :: iso_c_binding
188
189	USE kinds
190
191	USE pegrid
192
193	USE plant_canopy_model_mod, &
194	ONLY: plant_canopy, pch_index, &
195	pc_heating_rate, lad_s
196
197	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, &
200	sigma_sb, solar_constant, sun_direction, sun_dir_lat, &
201	sun_dir_lon, force_radiation_call
202
203	USE statistics, &
204	ONLY: hom, statistic_regions
205
206	USE surface_mod
207
208
209
210	IMPLICIT NONE
211
212	!-- configuration parameters (they can be setup in PALM config)
213	LOGICAL :: split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation
214	!< (.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 roofs
216	LOGICAL :: usm_energy_balance_wall = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs
217	LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used
218	LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used
219	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
227	INTEGER(iwp) :: land_category = 2 !< default category for land surface
228	INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone
229	INTEGER(iwp) :: pedestrant_category = 2 !< default category for wall surface in pedestrian zone
230	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
235	!< 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
315	!
316	!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
317	TYPE t_surf_vertical
318	REAL(wp), DIMENSION(:), ALLOCATABLE :: t
319	END TYPE t_surf_vertical
320	!
321	!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
322	TYPE t_wall_vertical
323	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t
324	END TYPE t_wall_vertical
325
326	!-- arrays for calculation of svf and csf
327	TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array
328	TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array
329	TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: asvf1, asvf2 !< realizations of svf array
330	TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET :: acsf1, acsf2 !< realizations of csf array
331	INTEGER(iwp) :: nsvfla !< dimmension of array allocated for storage of svf in local processor
332	INTEGER(iwp) :: ncsfla !< dimmension of array allocated for storage of csf in local processor
333	INTEGER(iwp) :: msvf, mcsf !< mod for swapping the growing array
334	INTEGER(iwp), PARAMETER :: gasize = 10000 !< initial size of growing arrays
335	!-- temporary arrays for calculation of csf in raytracing
336	INTEGER(iwp) :: maxboxesg !< max number of boxes ray can cross in the domain
337	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: boxes !< coordinates of gridboxes being crossed by ray
338	REAL(wp), DIMENSION(:), ALLOCATABLE :: crlens !< array of crossing lengths of ray for particular grid boxes
339	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored
340	#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_ip
343	#endif
344	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad_s_ray !< array of received lad_s for appropriate gridboxes crossed by ray
345
346	!-- arrays storing the values of USM
347	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 surfaces
349	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins !< array of sw radiation falling to local surface after i-th reflection
350	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl !< array of lw radiation for local surface after i-th reflection
351
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 reflections
354	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< array of lw radiation falling to local surface including radiation from reflections
355	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir !< array of direct sw radiation falling to local surface
356	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif !< array of diffuse sw radiation from sky and model boundary falling to local surface
357	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif !< array of diffuse lw radiation from sky and model boundary falling to local surface
358
359	!< Outward radiation is only valid for nonvirtual surfaces
360	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< array of reflected sw radiation for local surface in i-th reflection
361	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< array of reflected + emitted lw radiation for local surface in i-th reflection
362	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection
363	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
364	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
365	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
366	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf !< array of total radiation flux incoming to minus outgoing from local surface
367	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< local copy of rad_net (net radiation at surface)
368
369	!-- arrays for time averages
370	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l
371	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
372	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
373	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
374	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
375	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
376	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
377	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
378	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
379	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
380	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
381	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
382	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface
383
384	!-- block variables needed for calculation of the plant canopy model inside the urban surface model
385	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 surfaces
388	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 box
391	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw !< array of absorbed lw radiation for local plant canopy box
392	INTEGER(iwp) :: npcbl !< number of the plant canopy gridboxes in local processor
393	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy
394	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy
395	REAL(wp), DIMENSION(:,:,:), POINTER :: usm_lad !< subset of lad_s within urban surface, transformed to plain Z coordinate
396	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 depth
398	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing
399
400	!-- radiation related arrays (it should be better in interface of radiation module of PALM
401	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_dir !< direct sw radiation
402	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_diff !< diffusion sw radiation
403	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_diff !< diffusion lw radiation
404
405	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
406	!-- anthropogenic heat sources
407	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
408	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2)
409	REAL(wp), DIMENSION(:), ALLOCATABLE :: aheatprof !< diurnal profile of anthropogenic heat
410
411	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
412	!-- wall surface model
413	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
414	!-- wall surface model constants
415	INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched)
416	INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched)
417	INTEGER(iwp), PARAMETER :: nzw = 4 !< number of wall layers (fixed for now)
418
419	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
420	!< normalized soil, wall and roof layer depths (m/m)
421
422	REAL(wp) :: wall_inner_temperature = 296.0_wp !< temperature of the inner wall surface (~23 degrees C) (K)
423	REAL(wp) :: roof_inner_temperature = 296.0_wp !< temperature of the inner roof surface (~23 degrees C) (K)
424	REAL(wp) :: soil_inner_temperature = 283.0_wp !< temperature of the deep soil (~10 degrees C) (K)
425
426	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
427	!-- surface and material model variables for walls, ground, roofs
428	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
429	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m)
430
431	#if defined( __nopointer )
432	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h !< wall surface temperature (K) at horizontal walls
433	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_p !< progn. wall surface temperature (K) at horizontal walls
434
435	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v
436	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_p
437	#else
438	REAL(wp), DIMENSION(:), POINTER :: t_surf_h
439	REAL(wp), DIMENSION(:), POINTER :: t_surf_h_p
440
441	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_1
442	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_2
443
444	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v
445	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v_p
446
447	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_1
448	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_2
449	#endif
450	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_av !< average of wall surface temperature (K)
451
452	!-- Temporal tendencies for time stepping
453	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K)
454
455	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456	!-- Energy balance variables
457	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
458	!-- parameters of the land, roof and wall surfaces
459	REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface
460	!-- parameters of the wall surfaces
461	REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface
462
463	#if defined( __nopointer )
464	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h !< Wall temperature (K)
465	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av !< Average of t_wall
466	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_p !< Prog. wall temperature (K)
467
468	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v !< Wall temperature (K)
469	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av !< Average of t_wall
470	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_p !< Prog. wall temperature (K)
471	#else
472	REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p
473	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av, t_wall_h_1, t_wall_h_2
474
475	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p
476	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av, t_wall_v_1, t_wall_v_2
477	#endif
478
479	!-- Wall temporal tendencies for time stepping
480	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array
481
482	!-- Surface and material parameters classes (surface_type)
483	!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
484	INTEGER(iwp) :: n_surface_types !< number of the wall type categories
485	INTEGER(iwp), PARAMETER :: n_surface_params = 8 !< number of parameters for each type of the wall
486	INTEGER(iwp), PARAMETER :: ialbedo = 1 !< albedo of the surface
487	INTEGER(iwp), PARAMETER :: iemiss = 2 !< emissivity of the surface
488	INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 )
489	INTEGER(iwp), PARAMETER :: irough = 4 !< roughness relative to concrete
490	INTEGER(iwp), PARAMETER :: icsurf = 5 !< Surface skin layer heat capacity (J mâ2 Kâ1 )
491	INTEGER(iwp), PARAMETER :: ithick = 6 !< thickness of the surface (wall, roof, land) ( m )
492	INTEGER(iwp), PARAMETER :: irhoC = 7 !< volumetric heat capacity rho*C of the material ( J mâ3 Kâ1 )
493	INTEGER(iwp), PARAMETER :: ilambdah = 8 !< thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
494	CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports)
495	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types
496	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: surface_params !< parameters of wall types
497
498	CHARACTER(len=*), PARAMETER :: svf_file_name='usm_svf'
499
500	!-- interfaces of subroutines accessed from outside of this module
501	INTERFACE usm_check_data_output
502	MODULE PROCEDURE usm_check_data_output
503	END INTERFACE usm_check_data_output
504
505	INTERFACE usm_check_parameters
506	MODULE PROCEDURE usm_check_parameters
507	END INTERFACE usm_check_parameters
508
509	INTERFACE usm_data_output_3d
510	MODULE PROCEDURE usm_data_output_3d
511	END INTERFACE usm_data_output_3d
512
513	INTERFACE usm_define_netcdf_grid
514	MODULE PROCEDURE usm_define_netcdf_grid
515	END INTERFACE usm_define_netcdf_grid
516
517	INTERFACE usm_init_urban_surface
518	MODULE PROCEDURE usm_init_urban_surface
519	END INTERFACE usm_init_urban_surface
520
521	INTERFACE usm_material_heat_model
522	MODULE PROCEDURE usm_material_heat_model
523	END INTERFACE usm_material_heat_model
524
525	INTERFACE usm_parin
526	MODULE PROCEDURE usm_parin
527	END INTERFACE usm_parin
528
529	INTERFACE usm_radiation
530	MODULE PROCEDURE usm_radiation
531	END INTERFACE usm_radiation
532
533	INTERFACE usm_read_restart_data
534	MODULE PROCEDURE usm_read_restart_data
535	END INTERFACE usm_read_restart_data
536
537	INTERFACE usm_surface_energy_balance
538	MODULE PROCEDURE usm_surface_energy_balance
539	END INTERFACE usm_surface_energy_balance
540
541	INTERFACE usm_swap_timelevel
542	MODULE PROCEDURE usm_swap_timelevel
543	END INTERFACE usm_swap_timelevel
544
545	INTERFACE usm_write_restart_data
546	MODULE PROCEDURE usm_write_restart_data
547	END INTERFACE usm_write_restart_data
548
549	SAVE
550
551	PRIVATE
552
553	!-- 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
566
567
568	CONTAINS
569
570
571	!------------------------------------------------------------------------------!
572	! Description:
573	! ------------
574	!> This subroutine creates the necessary indices of the urban surfaces
575	!> and plant canopy and it allocates the needed arrays for USM
576	!------------------------------------------------------------------------------!
577	SUBROUTINE usm_allocate_urban_surface
578
579	IMPLICIT NONE
580
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, nsurfl5, MPI_INTEGER, surf, nsurfs5, 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
891
892	!
893	!-- Allocate radiation arrays which are part of the new data type.
894	!-- For horizontal surfaces.
895	ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns) )
896	ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
897	!
898	!-- New
899	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	!-- For vertical surfaces
905	DO l = 0, 3
906	ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) )
907	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	ENDDO
913
914	!-- Wall surface model
915	!-- allocate arrays for wall surface model and define pointers
916
917	!-- allocate array of wall types and wall parameters
918	ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
919	DO l = 0, 3
920	ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
921	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) )
930	DO l = 0, 3
931	ALLOCATE( surf_usm_v(l)%albedo_surf(1:surf_usm_v(l)%ns) )
932	ENDDO
933
934	!-- 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) )
942	!
943	!-- For vertical surfaces.
944	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) )
949	ENDDO
950
951	!-- allocate wall and roof material parameters. First for horizontal surfaces
952	ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) )
953	ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
954	ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
955	!
956	!-- For vertical surfaces.
957	DO l = 0, 3
958	ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) )
959	ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
960	ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
961	ENDDO
962
963	!-- allocate wall and roof layers sizes. For horizontal surfaces.
964	ALLOCATE ( zwn(nzb_wall:nzt_wall) )
965	ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
966	ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
967	ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
968	ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
969	ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
970	!
971	!-- For vertical surfaces.
972	DO l = 0, 3
973	ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
974	ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
975	ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
976	ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
977	ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
978	ENDDO
979
980	!-- allocate wall and roof temperature arrays, for horizontal walls
981	#if defined( __nopointer )
982	IF ( .NOT. ALLOCATED( t_surf_h ) ) &
983	ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
984	IF ( .NOT. ALLOCATED( t_surf_h_p ) ) &
985	ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
986	IF ( .NOT. ALLOCATED( t_wall_h ) ) &
987	ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
988	IF ( .NOT. ALLOCATED( t_wall_h_p ) ) &
989	ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
990	#else
991	!
992	!-- Allocate if required. Note, in case of restarts, some of these arrays
993	!-- might be already allocated.
994	IF ( .NOT. ALLOCATED( t_surf_h_1 ) ) &
995	ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) )
996	IF ( .NOT. ALLOCATED( t_surf_h_2 ) ) &
997	ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) )
998	IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) &
999	ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1000	IF ( .NOT. ALLOCATED( t_wall_h_2 ) ) &
1001	ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1002	!
1003	!-- initial assignment of the pointers
1004	t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2
1005	t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2
1006	#endif
1007
1008	!-- allocate wall and roof temperature arrays, for vertical walls if required
1009	#if defined( __nopointer )
1010	DO l = 0, 3
1011	IF ( .NOT. ALLOCATED( t_surf_v(l)%t ) ) &
1012	ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) )
1013	IF ( .NOT. ALLOCATED( t_surf_v_p(l)%t ) ) &
1014	ALLOCATE ( t_surf_v_p(l)%t(1:surf_usm_v(l)%ns) )
1015	IF ( .NOT. ALLOCATED( t_wall_v(l)%t ) ) &
1016	ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1017	IF ( .NOT. ALLOCATED( t_wall_v_p(l)%t ) ) &
1018	ALLOCATE ( t_wall_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1019	ENDDO
1020	#else
1021	!
1022	!-- Allocate if required. Note, in case of restarts, some of these arrays
1023	!-- might be already allocated.
1024	DO l = 0, 3
1025	IF ( .NOT. ALLOCATED( t_surf_v_1(l)%t ) ) &
1026	ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) )
1027	IF ( .NOT. ALLOCATED( t_surf_v_2(l)%t ) ) &
1028	ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) )
1029	IF ( .NOT. ALLOCATED( t_wall_v_1(l)%t ) ) &
1030	ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1031	IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) ) &
1032	ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1033	ENDDO
1034	!
1035	!-- initial assignment of the pointers
1036	t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2
1037	t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
1038	#endif
1039	!
1040	!-- Allocate intermediate timestep arrays. For horizontal surfaces.
1041	ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns) )
1042	ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1043	!
1044	!-- Set inital values for prognostic quantities
1045	IF ( ALLOCATED( surf_usm_h%tt_surface_m ) ) surf_usm_h%tt_surface_m = 0.0_wp
1046	IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp
1047	!
1048	!-- Now, for vertical surfaces
1049	DO l = 0, 3
1050	ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns) )
1051	ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1052	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) ) surf_usm_v(l)%tt_surface_m = 0.0_wp
1053	IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp
1054	ENDDO
1055
1056	!-- allocate wall heat flux output array and set initial values. For horizontal surfaces
1057	! ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns) ) !can be removed
1058	ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
1059	ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
1060	IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp
1061	IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp
1062	IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp
1063	!
1064	!-- Now, for vertical surfaces
1065	DO l = 0, 3
1066	! ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns) ) ! can be removed
1067	ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
1068	ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
1069	IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp
1070	IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp
1071	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp
1072	ENDDO
1073
1074	END SUBROUTINE usm_allocate_urban_surface
1075
1076
1077
1078	!------------------------------------------------------------------------------!
1079	! Description:
1080	! ------------
1081	!> Sum up and time-average urban surface output quantities as well as allocate
1082	!> the array necessary for storing the average.
1083	!------------------------------------------------------------------------------!
1084	SUBROUTINE usm_average_3d_data( mode, variable )
1085
1086	IMPLICIT NONE
1087
1088	CHARACTER (len=*), INTENT(IN) :: mode
1089	CHARACTER (len=*), INTENT(IN) :: variable
1090
1091	INTEGER(iwp) :: i, j, k, l, m, ids, iwl,istat
1092	CHARACTER (len=varnamelength) :: var, surfid
1093	INTEGER(iwp), PARAMETER :: nd = 5
1094	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
1095
1096	!-- find the real name of the variable
1097	var = TRIM(variable)
1098	DO i = 0, nd-1
1099	k = len(TRIM(var))
1100	j = len(TRIM(dirname(i)))
1101	IF ( var(k-j+1:k) == dirname(i) ) THEN
1102	ids = i
1103	var = var(:k-j)
1104	EXIT
1105	ENDIF
1106	ENDDO
1107	IF ( ids == -1 ) THEN
1108	var = TRIM(variable)
1109	ENDIF
1110	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
1111	!-- wall layers
1112	READ(var(12:12), '(I1)', iostat=istat ) iwl
1113	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1114	var = var(1:10)
1115	ELSE
1116	!-- wrong wall layer index
1117	RETURN
1118	ENDIF
1119	ENDIF
1120
1121	IF ( mode == 'allocate' ) THEN
1122
1123	SELECT CASE ( TRIM( var ) )
1124
1125	CASE ( 'usm_rad_net' )
1126	!-- array of complete radiation balance
1127	IF ( .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN
1128	ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
1129	surf_usm_h%rad_net_av = 0.0_wp
1130	ENDIF
1131	DO l = 0, 3
1132	IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN
1133	ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
1134	surf_usm_v(l)%rad_net_av = 0.0_wp
1135	ENDIF
1136	ENDDO
1137
1138	CASE ( 'usm_rad_insw' )
1139	!-- array of sw radiation falling to surface after i-th reflection
1140	IF ( .NOT. ALLOCATED(surf_usm_h%surfinsw_av) ) THEN
1141	ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
1142	surf_usm_h%surfinsw_av = 0.0_wp
1143	ENDIF
1144	DO l = 0, 3
1145	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinsw_av) ) THEN
1146	ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
1147	surf_usm_v(l)%surfinsw_av = 0.0_wp
1148	ENDIF
1149	ENDDO
1150
1151	CASE ( 'usm_rad_inlw' )
1152	!-- array of lw radiation falling to surface after i-th reflection
1153	IF ( .NOT. ALLOCATED(surf_usm_h%surfinlw_av) ) THEN
1154	ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
1155	surf_usm_h%surfinlw_av = 0.0_wp
1156	ENDIF
1157	DO l = 0, 3
1158	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinlw_av) ) THEN
1159	ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
1160	surf_usm_v(l)%surfinlw_av = 0.0_wp
1161	ENDIF
1162	ENDDO
1163
1164	CASE ( 'usm_rad_inswdir' )
1165	!-- array of direct sw radiation falling to surface from sun
1166	IF ( .NOT. ALLOCATED(surfinswdir_av) ) THEN
1167	ALLOCATE( surfinswdir_av(startenergy:endenergy) )
1168	surfinswdir_av = 0.0_wp
1169	ENDIF
1170
1171	CASE ( 'usm_rad_inswdif' )
1172	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1173	IF ( .NOT. ALLOCATED(surfinswdif_av) ) THEN
1174	ALLOCATE( surfinswdif_av(startenergy:endenergy) )
1175	surfinswdif_av = 0.0_wp
1176	ENDIF
1177
1178	CASE ( 'usm_rad_inswref' )
1179	!-- array of sw radiation falling to surface from reflections
1180	IF ( .NOT. ALLOCATED(surfinswref_av) ) THEN
1181	ALLOCATE( surfinswref_av(startenergy:endenergy) )
1182	surfinswref_av = 0.0_wp
1183	ENDIF
1184
1185	CASE ( 'usm_rad_inlwdif' )
1186	!-- array of sw radiation falling to surface after i-th reflection
1187	IF ( .NOT. ALLOCATED(surfinlwdif_av) ) THEN
1188	ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
1189	surfinlwdif_av = 0.0_wp
1190	ENDIF
1191
1192	CASE ( 'usm_rad_inlwref' )
1193	!-- array of lw radiation falling to surface from reflections
1194	IF ( .NOT. ALLOCATED(surfinlwref_av) ) THEN
1195	ALLOCATE( surfinlwref_av(startenergy:endenergy) )
1196	surfinlwref_av = 0.0_wp
1197	ENDIF
1198
1199	CASE ( 'usm_rad_outsw' )
1200	!-- array of sw radiation emitted from surface after i-th reflection
1201	IF ( .NOT. ALLOCATED(surfoutsw_av) ) THEN
1202	ALLOCATE( surfoutsw_av(startenergy:endenergy) )
1203	surfoutsw_av = 0.0_wp
1204	ENDIF
1205
1206	CASE ( 'usm_rad_outlw' )
1207	!-- array of lw radiation emitted from surface after i-th reflection
1208	IF ( .NOT. ALLOCATED(surfoutlw_av) ) THEN
1209	ALLOCATE( surfoutlw_av(startenergy:endenergy) )
1210	surfoutlw_av = 0.0_wp
1211	ENDIF
1212	CASE ( 'usm_rad_ressw' )
1213	!-- array of residua of sw radiation absorbed in surface after last reflection
1214	IF ( .NOT. ALLOCATED(surfins_av) ) THEN
1215	ALLOCATE( surfins_av(startenergy:endenergy) )
1216	surfins_av = 0.0_wp
1217	ENDIF
1218
1219	CASE ( 'usm_rad_reslw' )
1220	!-- array of residua of lw radiation absorbed in surface after last reflection
1221	IF ( .NOT. ALLOCATED(surfinl_av) ) THEN
1222	ALLOCATE( surfinl_av(startenergy:endenergy) )
1223	surfinl_av = 0.0_wp
1224	ENDIF
1225
1226	CASE ( 'usm_rad_hf' )
1227	!-- array of heat flux from radiation for surfaces after i-th reflection
1228	IF ( .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN
1229	ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
1230	surf_usm_h%surfhf_av = 0.0_wp
1231	ENDIF
1232	DO l = 0, 3
1233	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfhf_av) ) THEN
1234	ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
1235	surf_usm_v(l)%surfhf_av = 0.0_wp
1236	ENDIF
1237	ENDDO
1238
1239	CASE ( 'usm_wshf' )
1240	!-- array of sensible heat flux from surfaces
1241	!-- land surfaces
1242	IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN
1243	ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1244	surf_usm_h%wshf_eb_av = 0.0_wp
1245	ENDIF
1246	DO l = 0, 3
1247	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN
1248	ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1249	surf_usm_v(l)%wshf_eb_av = 0.0_wp
1250	ENDIF
1251	ENDDO
1252
1253	CASE ( 'usm_wghf' )
1254	!-- array of heat flux from ground (wall, roof, land)
1255	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN
1256	ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1257	surf_usm_h%wghf_eb_av = 0.0_wp
1258	ENDIF
1259	DO l = 0, 3
1260	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN
1261	ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1262	surf_usm_v(l)%wghf_eb_av = 0.0_wp
1263	ENDIF
1264	ENDDO
1265
1266	CASE ( 'usm_t_surf' )
1267	!-- surface temperature for surfaces
1268	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_av) ) THEN
1269	ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
1270	surf_usm_h%t_surf_av = 0.0_wp
1271	ENDIF
1272	DO l = 0, 3
1273	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_av) ) THEN
1274	ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
1275	surf_usm_v(l)%t_surf_av = 0.0_wp
1276	ENDIF
1277	ENDDO
1278
1279	CASE ( 'usm_t_wall' )
1280	!-- wall temperature for iwl layer of walls and land
1281	IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN
1282	ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1283	surf_usm_h%t_wall_av = 0.0_wp
1284	ENDIF
1285	DO l = 0, 3
1286	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN
1287	ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1288	surf_usm_v(l)%t_wall_av = 0.0_wp
1289	ENDIF
1290	ENDDO
1291
1292	CASE DEFAULT
1293	CONTINUE
1294
1295	END SELECT
1296
1297	ELSEIF ( mode == 'sum' ) THEN
1298
1299	SELECT CASE ( TRIM( var ) )
1300
1301	CASE ( 'usm_rad_net' )
1302	!-- array of complete radiation balance
1303	DO m = 1, surf_usm_h%ns
1304	surf_usm_h%rad_net_av(m) = &
1305	surf_usm_h%rad_net_av(m) + &
1306	surf_usm_h%rad_net_l(m)
1307	ENDDO
1308	DO l = 0, 3
1309	DO m = 1, surf_usm_v(l)%ns
1310	surf_usm_v(l)%rad_net_av(m) = &
1311	surf_usm_v(l)%rad_net_av(m) + &
1312	surf_usm_v(l)%rad_net_l(m)
1313	ENDDO
1314	ENDDO
1315
1316	CASE ( 'usm_rad_insw' )
1317	!-- array of sw radiation falling to surface after i-th reflection
1318	DO l = startenergy, endenergy
1319	IF ( surfl(id,l) == ids ) THEN
1320	surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
1321	ENDIF
1322	ENDDO
1323
1324	CASE ( 'usm_rad_inlw' )
1325	!-- array of lw radiation falling to surface after i-th reflection
1326	DO l = startenergy, endenergy
1327	IF ( surfl(id,l) == ids ) THEN
1328	surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
1329	ENDIF
1330	ENDDO
1331
1332	CASE ( 'usm_rad_inswdir' )
1333	!-- array of direct sw radiation falling to surface from sun
1334	DO l = startenergy, endenergy
1335	IF ( surfl(id,l) == ids ) THEN
1336	surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
1337	ENDIF
1338	ENDDO
1339
1340	CASE ( 'usm_rad_inswdif' )
1341	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1342	DO l = startenergy, endenergy
1343	IF ( surfl(id,l) == ids ) THEN
1344	surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
1345	ENDIF
1346	ENDDO
1347
1348	CASE ( 'usm_rad_inswref' )
1349	!-- array of sw radiation falling to surface from reflections
1350	DO l = startenergy, endenergy
1351	IF ( surfl(id,l) == ids ) THEN
1352	surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1353	surfinswdir(l) - surfinswdif(l)
1354	ENDIF
1355	ENDDO
1356
1357
1358	CASE ( 'usm_rad_inlwdif' )
1359	!-- array of sw radiation falling to surface after i-th reflection
1360	DO l = startenergy, endenergy
1361	IF ( surfl(id,l) == ids ) THEN
1362	surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1363	surfinswdir(l) - surfinswdif(l)
1364	ENDIF
1365	ENDDO
1366	!
1367	CASE ( 'usm_rad_inlwref' )
1368	!-- array of lw radiation falling to surface from reflections
1369	DO l = startenergy, endenergy
1370	IF ( surfl(id,l) == ids ) THEN
1371	surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
1372	ENDIF
1373	ENDDO
1374
1375	CASE ( 'usm_rad_outsw' )
1376	!-- array of sw radiation emitted from surface after i-th reflection
1377	DO l = startenergy, endenergy
1378	IF ( surfl(id,l) == ids ) THEN
1379	surfinlwref_av(l) = surfinlwref_av(l) + &
1380	surfinlw(l) - surfinlwdif(l)
1381	ENDIF
1382	ENDDO
1383
1384	CASE ( 'usm_rad_outlw' )
1385	!-- array of lw radiation emitted from surface after i-th reflection
1386	DO l = startenergy, endenergy
1387	IF ( surfl(id,l) == ids ) THEN
1388	surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
1389	ENDIF
1390	ENDDO
1391
1392	CASE ( 'usm_rad_ressw' )
1393	!-- array of residua of sw radiation absorbed in surface after last reflection
1394	DO l = startenergy, endenergy
1395	IF ( surfl(id,l) == ids ) THEN
1396	surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
1397	ENDIF
1398	ENDDO
1399
1400	CASE ( 'usm_rad_reslw' )
1401	!-- array of residua of lw radiation absorbed in surface after last reflection
1402	DO l = startenergy, endenergy
1403	IF ( surfl(id,l) == ids ) THEN
1404	surfins_av(l) = surfins_av(l) + surfins(l)
1405	ENDIF
1406	ENDDO
1407
1408	CASE ( 'usm_rad_hf' )
1409	!-- array of heat flux from radiation for surfaces after i-th reflection
1410	DO m = 1, surf_usm_h%ns
1411	surf_usm_h%surfhf_av(m) = &
1412	surf_usm_h%surfhf_av(m) + &
1413	surf_usm_h%surfhf(m)
1414	ENDDO
1415	DO l = 0, 3
1416	DO m = 1, surf_usm_v(l)%ns
1417	surf_usm_v(l)%surfhf_av(m) = &
1418	surf_usm_v(l)%surfhf_av(m) + &
1419	surf_usm_v(l)%surfhf(m)
1420	ENDDO
1421	ENDDO
1422
1423	CASE ( 'usm_wshf' )
1424	!-- array of sensible heat flux from surfaces (land, roof, wall)
1425	DO m = 1, surf_usm_h%ns
1426	surf_usm_h%wshf_eb_av(m) = &
1427	surf_usm_h%wshf_eb_av(m) + &
1428	surf_usm_h%wshf_eb(m)
1429	ENDDO
1430	DO l = 0, 3
1431	DO m = 1, surf_usm_v(l)%ns
1432	surf_usm_v(l)%wshf_eb_av(m) = &
1433	surf_usm_v(l)%wshf_eb_av(m) + &
1434	surf_usm_v(l)%wshf_eb(m)
1435	ENDDO
1436	ENDDO
1437
1438	CASE ( 'usm_wghf' )
1439	!-- array of heat flux from ground (wall, roof, land)
1440	DO m = 1, surf_usm_h%ns
1441	surf_usm_h%wghf_eb_av(m) = &
1442	surf_usm_h%wghf_eb_av(m) + &
1443	surf_usm_h%wghf_eb(m)
1444	ENDDO
1445	DO l = 0, 3
1446	DO m = 1, surf_usm_v(l)%ns
1447	surf_usm_v(l)%wghf_eb_av(m) = &
1448	surf_usm_v(l)%wghf_eb_av(m) + &
1449	surf_usm_v(l)%wghf_eb(m)
1450	ENDDO
1451	ENDDO
1452
1453	CASE ( 'usm_t_surf' )
1454	!-- surface temperature for surfaces
1455	DO m = 1, surf_usm_h%ns
1456	surf_usm_h%t_surf_av(m) = &
1457	surf_usm_h%t_surf_av(m) + &
1458	t_surf_h(m)
1459	ENDDO
1460	DO l = 0, 3
1461	DO m = 1, surf_usm_v(l)%ns
1462	surf_usm_v(l)%t_surf_av(m) = &
1463	surf_usm_v(l)%t_surf_av(m) + &
1464	t_surf_v(l)%t(m)
1465	ENDDO
1466	ENDDO
1467
1468	CASE ( 'usm_t_wall' )
1469	!-- wall temperature for iwl layer of walls and land
1470	DO m = 1, surf_usm_h%ns
1471	surf_usm_h%t_wall_av(iwl,m) = &
1472	surf_usm_h%t_wall_av(iwl,m) + &
1473	t_wall_h(iwl,m)
1474	ENDDO
1475	DO l = 0, 3
1476	DO m = 1, surf_usm_v(l)%ns
1477	surf_usm_v(l)%t_wall_av(iwl,m) = &
1478	surf_usm_v(l)%t_wall_av(iwl,m) + &
1479	t_wall_v(l)%t(iwl,m)
1480	ENDDO
1481	ENDDO
1482
1483	CASE DEFAULT
1484	CONTINUE
1485
1486	END SELECT
1487
1488	ELSEIF ( mode == 'average' ) THEN
1489
1490	SELECT CASE ( TRIM( var ) )
1491
1492	CASE ( 'usm_rad_net' )
1493	!-- array of complete radiation balance
1494	DO m = 1, surf_usm_h%ns
1495	surf_usm_h%rad_net_av(m) = &
1496	surf_usm_h%rad_net_av(m) / &
1497	REAL( average_count_3d, kind=wp )
1498	ENDDO
1499	DO l = 0, 3
1500	DO m = 1, surf_usm_v(l)%ns
1501	surf_usm_v(l)%rad_net_av(m) = &
1502	surf_usm_v(l)%rad_net_av(m) / &
1503	REAL( average_count_3d, kind=wp )
1504	ENDDO
1505	ENDDO
1506
1507	CASE ( 'usm_rad_insw' )
1508	!-- array of sw radiation falling to surface after i-th reflection
1509	DO l = startenergy, endenergy
1510	IF ( surfl(id,l) == ids ) THEN
1511	surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
1512	ENDIF
1513	ENDDO
1514
1515	CASE ( 'usm_rad_inlw' )
1516	!-- array of lw radiation falling to surface after i-th reflection
1517	DO l = startenergy, endenergy
1518	IF ( surfl(id,l) == ids ) THEN
1519	surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
1520	ENDIF
1521	ENDDO
1522
1523	CASE ( 'usm_rad_inswdir' )
1524	!-- array of direct sw radiation falling to surface from sun
1525	DO l = startenergy, endenergy
1526	IF ( surfl(id,l) == ids ) THEN
1527	surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
1528	ENDIF
1529	ENDDO
1530
1531	CASE ( 'usm_rad_inswdif' )
1532	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1533	DO l = startenergy, endenergy
1534	IF ( surfl(id,l) == ids ) THEN
1535	surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
1536	ENDIF
1537	ENDDO
1538
1539	CASE ( 'usm_rad_inswref' )
1540	!-- array of sw radiation falling to surface from reflections
1541	DO l = startenergy, endenergy
1542	IF ( surfl(id,l) == ids ) THEN
1543	surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
1544	ENDIF
1545	ENDDO
1546
1547	CASE ( 'usm_rad_inlwdif' )
1548	!-- array of sw radiation falling to surface after i-th reflection
1549	DO l = startenergy, endenergy
1550	IF ( surfl(id,l) == ids ) THEN
1551	surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
1552	ENDIF
1553	ENDDO
1554
1555	CASE ( 'usm_rad_inlwref' )
1556	!-- array of lw radiation falling to surface from reflections
1557	DO l = startenergy, endenergy
1558	IF ( surfl(id,l) == ids ) THEN
1559	surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
1560	ENDIF
1561	ENDDO
1562
1563	CASE ( 'usm_rad_outsw' )
1564	!-- array of sw radiation emitted from surface after i-th reflection
1565	DO l = startenergy, endenergy
1566	IF ( surfl(id,l) == ids ) THEN
1567	surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
1568	ENDIF
1569	ENDDO
1570
1571	CASE ( 'usm_rad_outlw' )
1572	!-- array of lw radiation emitted from surface after i-th reflection
1573	DO l = startenergy, endenergy
1574	IF ( surfl(id,l) == ids ) THEN
1575	surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
1576	ENDIF
1577	ENDDO
1578
1579	CASE ( 'usm_rad_ressw' )
1580	!-- array of residua of sw radiation absorbed in surface after last reflection
1581	DO l = startenergy, endenergy
1582	IF ( surfl(id,l) == ids ) THEN
1583	surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
1584	ENDIF
1585	ENDDO
1586
1587	CASE ( 'usm_rad_reslw' )
1588	!-- array of residua of lw radiation absorbed in surface after last reflection
1589	DO l = startenergy, endenergy
1590	IF ( surfl(id,l) == ids ) THEN
1591	surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
1592	ENDIF
1593	ENDDO
1594
1595	CASE ( 'usm_rad_hf' )
1596	!-- array of heat flux from radiation for surfaces after i-th reflection
1597	DO m = 1, surf_usm_h%ns
1598	surf_usm_h%surfhf_av(m) = &
1599	surf_usm_h%surfhf_av(m) / &
1600	REAL( average_count_3d, kind=wp )
1601	ENDDO
1602	DO l = 0, 3
1603	DO m = 1, surf_usm_v(l)%ns
1604	surf_usm_v(l)%surfhf_av(m) = &
1605	surf_usm_v(l)%surfhf_av(m) / &
1606	REAL( average_count_3d, kind=wp )
1607	ENDDO
1608	ENDDO
1609
1610	CASE ( 'usm_wshf' )
1611	!-- array of sensible heat flux from surfaces (land, roof, wall)
1612	DO m = 1, surf_usm_h%ns
1613	surf_usm_h%wshf_eb_av(m) = &
1614	surf_usm_h%wshf_eb_av(m) / &
1615	REAL( average_count_3d, kind=wp )
1616	ENDDO
1617	DO l = 0, 3
1618	DO m = 1, surf_usm_v(l)%ns
1619	surf_usm_v(l)%wshf_eb_av(m) = &
1620	surf_usm_v(l)%wshf_eb_av(m) / &
1621	REAL( average_count_3d, kind=wp )
1622	ENDDO
1623	ENDDO
1624
1625	CASE ( 'usm_wghf' )
1626	!-- array of heat flux from ground (wall, roof, land)
1627	DO m = 1, surf_usm_h%ns
1628	surf_usm_h%wghf_eb_av(m) = &
1629	surf_usm_h%wghf_eb_av(m) / &
1630	REAL( average_count_3d, kind=wp )
1631	ENDDO
1632	DO l = 0, 3
1633	DO m = 1, surf_usm_v(l)%ns
1634	surf_usm_v(l)%wghf_eb_av(m) = &
1635	surf_usm_v(l)%wghf_eb_av(m) / &
1636	REAL( average_count_3d, kind=wp )
1637	ENDDO
1638	ENDDO
1639
1640	CASE ( 'usm_t_surf' )
1641	!-- surface temperature for surfaces
1642	DO m = 1, surf_usm_h%ns
1643	surf_usm_h%t_surf_av(m) = &
1644	surf_usm_h%t_surf_av(m) / &
1645	REAL( average_count_3d, kind=wp )
1646	ENDDO
1647	DO l = 0, 3
1648	DO m = 1, surf_usm_v(l)%ns
1649	surf_usm_v(l)%t_surf_av(m) = &
1650	surf_usm_v(l)%t_surf_av(m) / &
1651	REAL( average_count_3d, kind=wp )
1652	ENDDO
1653	ENDDO
1654
1655	CASE ( 'usm_t_wall' )
1656	!-- wall temperature for iwl layer of walls and land
1657	DO m = 1, surf_usm_h%ns
1658	surf_usm_h%t_wall_av(iwl,m) = &
1659	surf_usm_h%t_wall_av(iwl,m) / &
1660	REAL( average_count_3d, kind=wp )
1661	ENDDO
1662	DO l = 0, 3
1663	DO m = 1, surf_usm_v(l)%ns
1664	surf_usm_v(l)%t_wall_av(iwl,m) = &
1665	surf_usm_v(l)%t_wall_av(iwl,m) / &
1666	REAL( average_count_3d, kind=wp )
1667	ENDDO
1668	ENDDO
1669
1670	END SELECT
1671
1672	ENDIF
1673
1674	END SUBROUTINE usm_average_3d_data
1675
1676
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
1757	!------------------------------------------------------------------------------!
1758	! Description:
1759	! ------------
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>
1763	!------------------------------------------------------------------------------!
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)nnxnny*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, nnxnnynzu, MPI_REAL, &
1952	usm_lad_g, nnxnnynzu, 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=(54+28), &
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, ndcsficsflt, ndcsfdcsflt, MPI_REAL, &
2253	pcsflt, ndcsfipcsflt, ndcsfdpcsflt, MPI_REAL, comm2d, ierr)
2254	CALL MPI_AlltoAllv(kcsflt, kdcsficsflt, kdcsfdcsflt, MPI_INTEGER, &
2255	kpcsflt, kdcsfipcsflt, kdcsfdpcsflt, 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
2334
2335
2336	!------------------------------------------------------------------------------!
2337	!
2338	! Description:
2339	! ------------
2340	!> Subroutine checks variables and assigns units.
2341	!> It is caaled out from subroutine check_parameters.
2342	!------------------------------------------------------------------------------!
2343	SUBROUTINE usm_check_data_output( variable, unit )
2344
2345	IMPLICIT NONE
2346
2347	CHARACTER (len=*),INTENT(IN) :: variable !:
2348	CHARACTER (len=*),INTENT(OUT) :: unit !:
2349
2350	CHARACTER (len=varnamelength) :: var
2351
2352	var = TRIM(variable)
2353	IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. &
2354	var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. &
2355	var(1:16) == 'usm_rad_inswdif_' .OR. var(1:16) == 'usm_rad_inswref_' .OR. &
2356	var(1:16) == 'usm_rad_inlwdif_' .OR. var(1:16) == 'usm_rad_inlwref_' .OR. &
2357	var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. &
2358	var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. &
2359	var(1:11) == 'usm_rad_hf_' .OR. &
2360	var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' ) THEN
2361	unit = 'W/m2'
2362	ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' ) THEN
2363	unit = 'K'
2364	ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. &
2365	var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. &
2366	var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis') THEN
2367	unit = '1'
2368	ELSE
2369	unit = 'illegal'
2370	ENDIF
2371
2372	END SUBROUTINE usm_check_data_output
2373
2374
2375	!------------------------------------------------------------------------------!
2376	! Description:
2377	! ------------
2378	!> Check parameters routine for urban surface model
2379	!------------------------------------------------------------------------------!
2380	SUBROUTINE usm_check_parameters
2381
2382	USE control_parameters, &
2383	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing, &
2384	lsf_surf, topography
2385
2386	!
2387	!-- Dirichlet boundary conditions are required as the surface fluxes are
2388	!-- calculated from the temperature/humidity gradients in the urban surface
2389	!-- model
2390	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
2391	message_string = 'urban surface model requires setting of '// &
2392	'bc_pt_b = "dirichlet" and '// &
2393	'bc_q_b = "dirichlet"'
2394	CALL message( 'check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
2395	ENDIF
2396
2397	IF ( .NOT. constant_flux_layer ) THEN
2398	message_string = 'urban surface model requires '// &
2399	'constant_flux_layer = .T.'
2400	CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
2401	ENDIF
2402	!
2403	!-- Surface forcing has to be disabled for LSF in case of enabled
2404	!-- urban surface module
2405	IF ( large_scale_forcing ) THEN
2406	lsf_surf = .FALSE.
2407	ENDIF
2408	!
2409	!-- Topography
2410	IF ( topography == 'flat' ) THEN
2411	message_string = 'topography /= "flat" is required '// &
2412	'when using the urban surface model'
2413	CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2414	ENDIF
2415
2416
2417	END SUBROUTINE usm_check_parameters
2418
2419
2420	!------------------------------------------------------------------------------!
2421	!
2422	! Description:
2423	! ------------
2424	!> Output of the 3D-arrays in netCDF and/or AVS format
2425	!> for variables of urban_surface model.
2426	!> It resorts the urban surface module output quantities from surf style
2427	!> indexing into temporary 3D array with indices (i,j,k).
2428	!> It is called from subroutine data_output_3d.
2429	!------------------------------------------------------------------------------!
2430	SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2431
2432	IMPLICIT NONE
2433
2434	INTEGER(iwp), INTENT(IN) :: av !<
2435	CHARACTER (len=*), INTENT(IN) :: variable !<
2436	INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0)
2437	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
2438	LOGICAL, INTENT(OUT) :: found !<
2439	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< sp - it has to correspond to module data_output_3d
2440	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: temp_pf !< temp array for urban surface output procedure
2441
2442	CHARACTER (len=varnamelength) :: var, surfid
2443	INTEGER(iwp), PARAMETER :: nd = 5
2444	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
2445	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iroof, isouth, inorth, iwest, ieast /)
2446	INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart
2447	INTEGER(iwp), DIMENSION(0:nd-1) :: dirend
2448	INTEGER(iwp) :: ids,isurf,isvf,isurfs,isurflt
2449	INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat, l, m
2450	INTEGER(iwp) :: k_topo !< topography top index
2451
2452	dirstart = (/ startland, startwall, startwall, startwall, startwall /)
2453	dirend = (/ endland, endwall, endwall, endwall, endwall /)
2454
2455	found = .TRUE.
2456	temp_pf = -1._wp
2457
2458	ids = -1
2459	var = TRIM(variable)
2460	DO i = 0, nd-1
2461	k = len(TRIM(var))
2462	j = len(TRIM(dirname(i)))
2463	IF ( var(k-j+1:k) == dirname(i) ) THEN
2464	ids = i
2465	var = var(:k-j)
2466	EXIT
2467	ENDIF
2468	ENDDO
2469	IF ( ids == -1 ) THEN
2470	var = TRIM(variable)
2471	ENDIF
2472	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
2473	!-- wall layers
2474	READ(var(12:12), '(I1)', iostat=istat ) iwl
2475	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2476	var = var(1:10)
2477	ENDIF
2478	ENDIF
2479	IF ( (var(1:8) == 'usm_svf_' .OR. var(1:8) == 'usm_dif_') .AND. len(TRIM(var)) >= 13 ) THEN
2480	!-- svf values to particular surface
2481	surfid = var(9:)
2482	i = index(surfid,'_')
2483	j = index(surfid(i+1:),'_')
2484	READ(surfid(1:i-1),*, iostat=istat ) is
2485	IF ( istat == 0 ) THEN
2486	READ(surfid(i+1:i+j-1),*, iostat=istat ) js
2487	ENDIF
2488	IF ( istat == 0 ) THEN
2489	READ(surfid(i+j+1:),*, iostat=istat ) ks
2490	ENDIF
2491	IF ( istat == 0 ) THEN
2492	var = var(1:7)
2493	ENDIF
2494	ENDIF
2495
2496	SELECT CASE ( TRIM(var) )
2497
2498	CASE ( 'usm_surfz' )
2499	!-- array of lw radiation falling to local surface after i-th reflection
2500	DO m = 1, surf_usm_h%ns
2501	i = surf_usm_h%i(m)
2502	j = surf_usm_h%j(m)
2503	k = surf_usm_h%k(m)
2504	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
2505	ENDDO
2506	DO l = 0, 3
2507	DO m = 1, surf_usm_v(l)%ns
2508	i = surf_usm_v(l)%i(m)
2509	j = surf_usm_v(l)%j(m)
2510	k = surf_usm_v(l)%k(m)
2511	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
2512	ENDDO
2513	ENDDO
2514
2515	CASE ( 'usm_surfcat' )
2516	!-- surface category
2517	DO m = 1, surf_usm_h%ns
2518	i = surf_usm_h%i(m)
2519	j = surf_usm_h%j(m)
2520	k = surf_usm_h%k(m)
2521	temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2522	ENDDO
2523	DO l = 0, 3
2524	DO m = 1, surf_usm_v(l)%ns
2525	i = surf_usm_v(l)%i(m)
2526	j = surf_usm_v(l)%j(m)
2527	k = surf_usm_v(l)%k(m)
2528	temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
2529	ENDDO
2530	ENDDO
2531
2532	CASE ( 'usm_surfalb' )
2533	!-- surface albedo
2534	DO m = 1, surf_usm_h%ns
2535	i = surf_usm_h%i(m)
2536	j = surf_usm_h%j(m)
2537	k = surf_usm_h%k(m)
2538	temp_pf(k,j,i) = surf_usm_h%albedo_surf(m)
2539	ENDDO
2540	DO l = 0, 3
2541	DO m = 1, surf_usm_v(l)%ns
2542	i = surf_usm_v(l)%i(m)
2543	j = surf_usm_v(l)%j(m)
2544	k = surf_usm_v(l)%k(m)
2545	temp_pf(k,j,i) = surf_usm_v(l)%albedo_surf(m)
2546	ENDDO
2547	ENDDO
2548
2549	CASE ( 'usm_surfemis' )
2550	!-- surface albedo
2551	DO m = 1, surf_usm_h%ns
2552	i = surf_usm_h%i(m)
2553	j = surf_usm_h%j(m)
2554	k = surf_usm_h%k(m)
2555	temp_pf(k,j,i) = surf_usm_h%emiss_surf(m)
2556	ENDDO
2557	DO l = 0, 3
2558	DO m = 1, surf_usm_v(l)%ns
2559	i = surf_usm_v(l)%i(m)
2560	j = surf_usm_v(l)%j(m)
2561	k = surf_usm_v(l)%k(m)
2562	temp_pf(k,j,i) = surf_usm_v(l)%emiss_surf(m)
2563	ENDDO
2564	ENDDO
2565	!
2566	!-- Not adjusted so far
2567	CASE ( 'usm_svf', 'usm_dif' )
2568	!-- shape view factors or iradiance factors to selected surface
2569	IF ( TRIM(var)=='usm_svf' ) THEN
2570	k = 1
2571	ELSE
2572	k = 2
2573	ENDIF
2574	DO isvf = 1, nsvfl
2575	isurflt = svfsurf(1, isvf)
2576	isurfs = svfsurf(2, isvf)
2577
2578	IF ( surf(ix,isurfs) == is .AND. surf(iy,isurfs) == js .AND. &
2579	surf(iz,isurfs) == ks .AND. surf(id,isurfs) == ids ) THEN
2580	!-- correct source surface
2581	temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
2582	ENDIF
2583	ENDDO
2584
2585	CASE ( 'usm_rad_net' )
2586	!-- array of complete radiation balance
2587	IF ( av == 0 ) THEN
2588	DO m = 1, surf_usm_h%ns
2589	i = surf_usm_h%i(m)
2590	j = surf_usm_h%j(m)
2591	k = surf_usm_h%k(m)
2592	temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
2593	ENDDO
2594	DO l = 0, 3
2595	DO m = 1, surf_usm_v(l)%ns
2596	i = surf_usm_v(l)%i(m)
2597	j = surf_usm_v(l)%j(m)
2598	k = surf_usm_v(l)%k(m)
2599	temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
2600	ENDDO
2601	ENDDO
2602	ELSE
2603	DO m = 1, surf_usm_h%ns
2604	i = surf_usm_h%i(m)
2605	j = surf_usm_h%j(m)
2606	k = surf_usm_h%k(m)
2607	temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
2608	ENDDO
2609	DO l = 0, 3
2610	DO m = 1, surf_usm_v(l)%ns
2611	i = surf_usm_v(l)%i(m)
2612	j = surf_usm_v(l)%j(m)
2613	k = surf_usm_v(l)%k(m)
2614	temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
2615	ENDDO
2616	ENDDO
2617	ENDIF
2618
2619	CASE ( 'usm_rad_insw' )
2620	!-- array of sw radiation falling to surface after i-th reflection
2621	DO isurf = dirstart(ids), dirend(ids)
2622	IF ( surfl(id,isurf) == ids ) THEN
2623	IF ( av == 0 ) THEN
2624	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
2625	ELSE
2626	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
2627	ENDIF
2628	ENDIF
2629	ENDDO
2630
2631	CASE ( 'usm_rad_inlw' )
2632	!-- array of lw radiation falling to surface after i-th reflection
2633	DO isurf = dirstart(ids), dirend(ids)
2634	IF ( surfl(id,isurf) == ids ) THEN
2635	IF ( av == 0 ) THEN
2636	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
2637	ELSE
2638	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
2639	ENDIF
2640	ENDIF
2641	ENDDO
2642
2643	CASE ( 'usm_rad_inswdir' )
2644	!-- array of direct sw radiation falling to surface from sun
2645	DO isurf = dirstart(ids), dirend(ids)
2646	IF ( surfl(id,isurf) == ids ) THEN
2647	IF ( av == 0 ) THEN
2648	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
2649	ELSE
2650	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
2651	ENDIF
2652	ENDIF
2653	ENDDO
2654
2655	CASE ( 'usm_rad_inswdif' )
2656	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
2657	DO isurf = dirstart(ids), dirend(ids)
2658	IF ( surfl(id,isurf) == ids ) THEN
2659	IF ( av == 0 ) THEN
2660	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
2661	ELSE
2662	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
2663	ENDIF
2664	ENDIF
2665	ENDDO
2666
2667	CASE ( 'usm_rad_inswref' )
2668	!-- array of sw radiation falling to surface from reflections
2669	DO isurf = dirstart(ids), dirend(ids)
2670	IF ( surfl(id,isurf) == ids ) THEN
2671	IF ( av == 0 ) THEN
2672	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
2673	surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
2674	ELSE
2675	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
2676	ENDIF
2677	ENDIF
2678	ENDDO
2679
2680	CASE ( 'usm_rad_inlwref' )
2681	!-- array of lw radiation falling to surface from reflections
2682	DO isurf = dirstart(ids), dirend(ids)
2683	IF ( surfl(id,isurf) == ids ) THEN
2684	IF ( av == 0 ) THEN
2685	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
2686	ELSE
2687	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
2688	ENDIF
2689	ENDIF
2690	ENDDO
2691
2692	CASE ( 'usm_rad_outsw' )
2693	!-- array of sw radiation emitted from surface after i-th reflection
2694	DO isurf = dirstart(ids), dirend(ids)
2695	IF ( surfl(id,isurf) == ids ) THEN
2696	IF ( av == 0 ) THEN
2697	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
2698	ELSE
2699	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
2700	ENDIF
2701	ENDIF
2702	ENDDO
2703
2704	CASE ( 'usm_rad_outlw' )
2705	!-- array of lw radiation emitted from surface after i-th reflection
2706	DO isurf = dirstart(ids), dirend(ids)
2707	IF ( surfl(id,isurf) == ids ) THEN
2708	IF ( av == 0 ) THEN
2709	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
2710	ELSE
2711	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
2712	ENDIF
2713	ENDIF
2714	ENDDO
2715
2716	CASE ( 'usm_rad_ressw' )
2717	!-- average of array of residua of sw radiation absorbed in surface after last reflection
2718	DO isurf = dirstart(ids), dirend(ids)
2719	IF ( surfl(id,isurf) == ids ) THEN
2720	IF ( av == 0 ) THEN
2721	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
2722	ELSE
2723	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
2724	ENDIF
2725	ENDIF
2726	ENDDO
2727
2728	CASE ( 'usm_rad_reslw' )
2729	!-- average of array of residua of lw radiation absorbed in surface after last reflection
2730	DO isurf = dirstart(ids), dirend(ids)
2731	IF ( surfl(id,isurf) == ids ) THEN
2732	IF ( av == 0 ) THEN
2733	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
2734	ELSE
2735	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
2736	ENDIF
2737	ENDIF
2738	ENDDO
2739
2740	CASE ( 'usm_rad_hf' )
2741	!-- array of heat flux from radiation for surfaces after all reflections
2742	IF ( av == 0 ) THEN
2743	DO m = 1, surf_usm_h%ns
2744	i = surf_usm_h%i(m)
2745	j = surf_usm_h%j(m)
2746	k = surf_usm_h%k(m)
2747	temp_pf(k,j,i) = surf_usm_h%surfhf(m)
2748	ENDDO
2749	DO l = 0, 3
2750	DO m = 1, surf_usm_v(l)%ns
2751	i = surf_usm_v(l)%i(m)
2752	j = surf_usm_v(l)%j(m)
2753	k = surf_usm_v(l)%k(m)
2754	temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
2755	ENDDO
2756	ENDDO
2757	ELSE
2758	DO m = 1, surf_usm_h%ns
2759	i = surf_usm_h%i(m)
2760	j = surf_usm_h%j(m)
2761	k = surf_usm_h%k(m)
2762	temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
2763	ENDDO
2764	DO l = 0, 3
2765	DO m = 1, surf_usm_v(l)%ns
2766	i = surf_usm_v(l)%i(m)
2767	j = surf_usm_v(l)%j(m)
2768	k = surf_usm_v(l)%k(m)
2769	temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
2770	ENDDO
2771	ENDDO
2772	ENDIF
2773
2774	CASE ( 'usm_wshf' )
2775	!-- array of sensible heat flux from surfaces
2776	IF ( av == 0 ) THEN
2777	DO m = 1, surf_usm_h%ns
2778	i = surf_usm_h%i(m)
2779	j = surf_usm_h%j(m)
2780	k = surf_usm_h%k(m)
2781	temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2782	ENDDO
2783	DO l = 0, 3
2784	DO m = 1, surf_usm_v(l)%ns
2785	i = surf_usm_v(l)%i(m)
2786	j = surf_usm_v(l)%j(m)
2787	k = surf_usm_v(l)%k(m)
2788	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2789	ENDDO
2790	ENDDO
2791	ELSE
2792	DO m = 1, surf_usm_h%ns
2793	i = surf_usm_h%i(m)
2794	j = surf_usm_h%j(m)
2795	k = surf_usm_h%k(m)
2796	temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2797	ENDDO
2798	DO l = 0, 3
2799	DO m = 1, surf_usm_v(l)%ns
2800	i = surf_usm_v(l)%i(m)
2801	j = surf_usm_v(l)%j(m)
2802	k = surf_usm_v(l)%k(m)
2803	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2804	ENDDO
2805	ENDDO
2806	ENDIF
2807
2808
2809	CASE ( 'usm_wghf' )
2810	!-- array of heat flux from ground (land, wall, roof)
2811	IF ( av == 0 ) THEN
2812	DO m = 1, surf_usm_h%ns
2813	i = surf_usm_h%i(m)
2814	j = surf_usm_h%j(m)
2815	k = surf_usm_h%k(m)
2816	temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2817	ENDDO
2818	DO l = 0, 3
2819	DO m = 1, surf_usm_v(l)%ns
2820	i = surf_usm_v(l)%i(m)
2821	j = surf_usm_v(l)%j(m)
2822	k = surf_usm_v(l)%k(m)
2823	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2824	ENDDO
2825	ENDDO
2826	ELSE
2827	DO m = 1, surf_usm_h%ns
2828	i = surf_usm_h%i(m)
2829	j = surf_usm_h%j(m)
2830	k = surf_usm_h%k(m)
2831	temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2832	ENDDO
2833	DO l = 0, 3
2834	DO m = 1, surf_usm_v(l)%ns
2835	i = surf_usm_v(l)%i(m)
2836	j = surf_usm_v(l)%j(m)
2837	k = surf_usm_v(l)%k(m)
2838	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2839	ENDDO
2840	ENDDO
2841	ENDIF
2842
2843	CASE ( 'usm_t_surf' )
2844	!-- surface temperature for surfaces
2845	IF ( av == 0 ) THEN
2846	DO m = 1, surf_usm_h%ns
2847	i = surf_usm_h%i(m)
2848	j = surf_usm_h%j(m)
2849	k = surf_usm_h%k(m)
2850	temp_pf(k,j,i) = t_surf_h(m)
2851	ENDDO
2852	DO l = 0, 3
2853	DO m = 1, surf_usm_v(l)%ns
2854	i = surf_usm_v(l)%i(m)
2855	j = surf_usm_v(l)%j(m)
2856	k = surf_usm_v(l)%k(m)
2857	temp_pf(k,j,i) = t_surf_v(l)%t(m)
2858	ENDDO
2859	ENDDO
2860	ELSE
2861	DO m = 1, surf_usm_h%ns
2862	i = surf_usm_h%i(m)
2863	j = surf_usm_h%j(m)
2864	k = surf_usm_h%k(m)
2865	temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
2866	ENDDO
2867	DO l = 0, 3
2868	DO m = 1, surf_usm_v(l)%ns
2869	i = surf_usm_v(l)%i(m)
2870	j = surf_usm_v(l)%j(m)
2871	k = surf_usm_v(l)%k(m)
2872	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
2873	ENDDO
2874	ENDDO
2875	ENDIF
2876
2877	CASE ( 'usm_t_wall' )
2878	!-- wall temperature for iwl layer of walls and land
2879	IF ( av == 0 ) THEN
2880	DO m = 1, surf_usm_h%ns
2881	i = surf_usm_h%i(m)
2882	j = surf_usm_h%j(m)
2883	k = surf_usm_h%k(m)
2884	temp_pf(k,j,i) = t_wall_h(iwl,m)
2885	ENDDO
2886	DO l = 0, 3
2887	DO m = 1, surf_usm_v(l)%ns
2888	i = surf_usm_v(l)%i(m)
2889	j = surf_usm_v(l)%j(m)
2890	k = surf_usm_v(l)%k(m)
2891	temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
2892	ENDDO
2893	ENDDO
2894	ELSE
2895	DO m = 1, surf_usm_h%ns
2896	i = surf_usm_h%i(m)
2897	j = surf_usm_h%j(m)
2898	k = surf_usm_h%k(m)
2899	temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
2900	ENDDO
2901	DO l = 0, 3
2902	DO m = 1, surf_usm_v(l)%ns
2903	i = surf_usm_v(l)%i(m)
2904	j = surf_usm_v(l)%j(m)
2905	k = surf_usm_v(l)%k(m)
2906	temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
2907	ENDDO
2908	ENDDO
2909	ENDIF
2910
2911	CASE DEFAULT
2912	found = .FALSE.
2913
2914	END SELECT
2915	!
2916	!-- Rearrange dimensions for NetCDF output
2917	DO j = nys, nyn
2918	DO i = nxl, nxr
2919	DO k = nzb_do, nzt_do
2920	local_pf(i,j,k) = temp_pf(k,j,i)
2921	ENDDO
2922	ENDDO
2923	ENDDO
2924
2925	END SUBROUTINE usm_data_output_3d
2926
2927
2928	!------------------------------------------------------------------------------!
2929	!
2930	! Description:
2931	! ------------
2932	!> Soubroutine defines appropriate grid for netcdf variables.
2933	!> It is called out from subroutine netcdf.
2934	!------------------------------------------------------------------------------!
2935	SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
2936
2937	IMPLICIT NONE
2938
2939	CHARACTER (len=*), INTENT(IN) :: variable !<
2940	LOGICAL, INTENT(OUT) :: found !<
2941	CHARACTER (len=*), INTENT(OUT) :: grid_x !<
2942	CHARACTER (len=*), INTENT(OUT) :: grid_y !<
2943	CHARACTER (len=*), INTENT(OUT) :: grid_z !<
2944
2945	CHARACTER (len=varnamelength) :: var
2946
2947	var = TRIM(variable)
2948	IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. &
2949	var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. &
2950	var(1:16) == 'usm_rad_inswdif_' .OR. var(1:16) == 'usm_rad_inswref_' .OR. &
2951	var(1:16) == 'usm_rad_inlwdif_' .OR. var(1:16) == 'usm_rad_inlwref_' .OR. &
2952	var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. &
2953	var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. &
2954	var(1:11) == 'usm_rad_hf_' .OR. &
2955	var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
2956	var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. &
2957	var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. &
2958	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
2960
2961	found = .TRUE.
2962	grid_x = 'x'
2963	grid_y = 'y'
2964	grid_z = 'zu'
2965	ELSE
2966	found = .FALSE.
2967	grid_x = 'none'
2968	grid_y = 'none'
2969	grid_z = 'none'
2970	ENDIF
2971
2972	END SUBROUTINE usm_define_netcdf_grid
2973
2974
2975	!------------------------------------------------------------------------------!
2976	!> Finds first model boundary crossed by a ray
2977	!------------------------------------------------------------------------------!
2978	PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross)
2979
2980	IMPLICIT NONE
2981
2982	INTEGER(iwp) :: d !<
2983	INTEGER(iwp) :: seldim !< found fist crossing index
2984
2985	INTEGER(iwp), DIMENSION(3) :: bdyd !< boundary direction
2986	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 distance
2992	REAL(wp), DIMENSION(3), INTENT(in) :: origin !< ray origin
2993	REAL(wp), DIMENSION(3), INTENT(in) :: uvect !< ray unit vector
2994
2995
2996	bdydim = nzut + .5_wp !< top boundary
2997	bdyd(1) = isky
2998	crossdist(1) = ( bdydim - origin(1) ) / uvect(1) !< subroutine called only when uvect(1)>0
2999
3000	IF ( uvect(2) == 0._wp ) THEN
3001	crossdist(2) = huge(1._wp)
3002	ELSE
3003	IF ( uvect(2) >= 0._wp ) THEN
3004	bdydim = ny + .5_wp !< north global boundary
3005	bdyd(2) = inorthb
3006	ELSE
3007	bdydim = -.5_wp !< south global boundary
3008	bdyd(2) = isouthb
3009	ENDIF
3010	crossdist(2) = ( bdydim - origin(2) ) / uvect(2)
3011	ENDIF
3012
3013	IF ( uvect(3) == 0._wp ) THEN
3014	crossdist(3) = huge(1._wp)
3015	ELSE
3016	IF ( uvect(3) >= 0._wp ) THEN
3017	bdydim = nx + .5_wp !< east global boundary
3018	bdyd(3) = ieastb
3019	ELSE
3020	bdydim = -.5_wp !< west global boundary
3021	bdyd(3) = iwestb
3022	ENDIF
3023	crossdist(3) = ( bdydim - origin(3) ) / uvect(3)
3024	ENDIF
3025
3026	seldim = minloc(crossdist, 1)
3027	dist = crossdist(seldim)
3028	d = bdyd(seldim)
3029
3030	bdycross(1) = d
3031	bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &
3032	+ .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
3033
3034	END SUBROUTINE
3035
3036
3037	!------------------------------------------------------------------------------!
3038	!> Determines whether two faces are oriented towards each other
3039	!------------------------------------------------------------------------------!
3040	PURE LOGICAL FUNCTION usm_facing(x, y, z, d, x2, y2, z2, d2)
3041	IMPLICIT NONE
3042	INTEGER(iwp), INTENT(in) :: x, y, z, d, x2, y2, z2, d2
3043
3044	usm_facing = .FALSE.
3045	IF ( d==iroof .AND. d2==iroof ) RETURN
3046	IF ( d==isky .AND. d2==isky ) RETURN
3047	IF ( (d==isouth .OR. d==inorthb) .AND. (d2==isouth.OR.d2==inorthb) ) RETURN
3048	IF ( (d==inorth .OR. d==isouthb) .AND. (d2==inorth.OR.d2==isouthb) ) RETURN
3049	IF ( (d==iwest .OR. d==ieastb) .AND. (d2==iwest.OR.d2==ieastb) ) RETURN
3050	IF ( (d==ieast .OR. d==iwestb) .AND. (d2==ieast.OR.d2==iwestb) ) RETURN
3051
3052	SELECT CASE (d)
3053	CASE (iroof) !< ground, roof
3054	IF ( z2 < z ) RETURN
3055	CASE (isky) !< sky
3056	IF ( z2 > z ) RETURN
3057	CASE (isouth, inorthb) !< south facing
3058	IF ( y2 > y ) RETURN
3059	CASE (inorth, isouthb) !< north facing
3060	IF ( y2 < y ) RETURN
3061	CASE (iwest, ieastb) !< west facing
3062	IF ( x2 > x ) RETURN
3063	CASE (ieast, iwestb) !< east facing
3064	IF ( x2 < x ) RETURN
3065	END SELECT
3066
3067	SELECT CASE (d2)
3068	CASE (iroof) !< ground, roof
3069	IF ( z < z2 ) RETURN
3070	CASE (isky) !< sky
3071	IF ( z > z2 ) RETURN
3072	CASE (isouth, inorthb) !< south facing
3073	IF ( y > y2 ) RETURN
3074	CASE (inorth, isouthb) !< north facing
3075	IF ( y < y2 ) RETURN
3076	CASE (iwest, ieastb) !< west facing
3077	IF ( x > x2 ) RETURN
3078	CASE (ieast, iwestb) !< east facing
3079	IF ( x < x2 ) RETURN
3080	CASE (-1)
3081	CONTINUE
3082	END SELECT
3083
3084	usm_facing = .TRUE.
3085
3086	END FUNCTION usm_facing
3087
3088
3089	!------------------------------------------------------------------------------!
3090	! Description:
3091	! ------------
3092	!> Initialization of the wall surface model
3093	!------------------------------------------------------------------------------!
3094	SUBROUTINE usm_init_material_model
3095
3096	IMPLICIT NONE
3097
3098	INTEGER(iwp) :: k, l, m !< running indices
3099
3100	CALL location_message( ' initialization of wall surface model', .TRUE. )
3101
3102	!-- Calculate wall grid spacings.
3103	!-- 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	!
3109	!-- apply for all particular surface grids. First for horizontal surfaces
3110	DO m = 1, surf_usm_h%ns
3111	surf_usm_h%zw(:,m) = zwn(:) * &
3112	surf_usm_h%thickness_wall(m)
3113	surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
3114	DO k = nzb_wall+1, nzt_wall
3115	surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) - &
3116	surf_usm_h%zw(k-1,m)
3117	ENDDO
3118
3119	surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
3120
3121	DO k = nzb_wall, nzt_wall-1
3122	surf_usm_h%dz_wall_stag(k,m) = 0.5 * ( &
3123	surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
3124	ENDDO
3125	surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
3126	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
3129	!
3130	!-- For vertical surfaces
3131	DO l = 0, 3
3132	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	surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
3136	DO k = nzb_wall+1, nzt_wall
3137	surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - &
3138	surf_usm_v(l)%zw(k-1,m)
3139	ENDDO
3140
3141	surf_usm_v(l)%dz_wall(nzt_wall+1,m) = &
3142	surf_usm_v(l)%dz_wall(nzt_wall,m)
3143
3144	DO k = nzb_wall, nzt_wall-1
3145	surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * ( &
3146	surf_usm_v(l)%dz_wall(k+1,m) + &
3147	surf_usm_v(l)%dz_wall(k,m) )
3148	ENDDO
3149	surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = &
3150	surf_usm_v(l)%dz_wall(nzt_wall,m)
3151	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
3154	ENDDO
3155
3156
3157	CALL location_message( ' wall structures filed out', .TRUE. )
3158
3159	CALL location_message( ' initialization of wall surface model finished', .TRUE. )
3160
3161	END SUBROUTINE usm_init_material_model
3162
3163
3164	!------------------------------------------------------------------------------!
3165	! Description:
3166	! ------------
3167	!> Initialization of the urban surface model
3168	!------------------------------------------------------------------------------!
3169	SUBROUTINE usm_init_urban_surface
3170
3171	IMPLICIT NONE
3172
3173	INTEGER(iwp) :: i, j, k, l, m !< running indices
3174	REAL(wp) :: c, d, tin, exn
3175
3176	!
3177	!-- NOPOINTER version not implemented yet
3178	#if defined( __nopointer )
3179	message_string = 'The urban surface module only runs with POINTER version'
3180	CALL message( 'urban_surface_mod', 'PA0452', 1, 2, 0, 6, 0 )
3181	#endif
3182
3183	CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
3184	!-- surface forcing have to be disabled for LSF
3185	!-- in case of enabled urban surface module
3186	IF ( large_scale_forcing ) THEN
3187	lsf_surf = .FALSE.
3188	ENDIF
3189
3190	!-- init anthropogenic sources of heat
3191	CALL usm_allocate_urban_surface()
3192
3193	!-- read the surface_types array somewhere
3194	CALL usm_read_urban_surface_types()
3195
3196	!-- init material heat model
3197	CALL usm_init_material_model()
3198
3199	IF ( usm_anthropogenic_heat ) THEN
3200	!-- init anthropogenic sources of heat (from transportation for now)
3201	CALL usm_read_anthropogenic_heat()
3202	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
3224	IF ( plant_canopy ) THEN
3225	!-- gridpcbl was only necessary for initialization
3226	DEALLOCATE( gridpcbl )
3227	IF ( .NOT. ALLOCATED(pc_heating_rate) ) THEN
3228	!-- then pc_heating_rate is allocated in init_plant_canopy
3229	!-- in case of cthf /= 0 => we need to allocate it for our use here
3230	ALLOCATE( pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3231	ENDIF
3232	ENDIF
3233
3234	!-- Intitialization of the surface and wall/ground/roof temperature
3235
3236	!-- Initialization for restart runs
3237	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
3238	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
3239
3240	!-- Calculate initial surface temperature from pt of adjacent gridbox
3241	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
3242
3243	!
3244	!-- At horizontal surfaces. Please note, t_surf_h is defined on a
3245	!-- different data type, but with the same dimension.
3246	#if ! defined( __nopointer )
3247	DO m = 1, surf_usm_h%ns
3248	i = surf_usm_h%i(m)
3249	j = surf_usm_h%j(m)
3250	k = surf_usm_h%k(m)
3251
3252	t_surf_h(m) = pt(k,j,i) * exn
3253	ENDDO
3254	!
3255	!-- At vertical surfaces.
3256	DO l = 0, 3
3257	DO m = 1, surf_usm_v(l)%ns
3258	i = surf_usm_v(l)%i(m)
3259	j = surf_usm_v(l)%j(m)
3260	k = surf_usm_v(l)%k(m)
3261
3262	t_surf_v(l)%t(m) = pt(k,j,i) * exn
3263	ENDDO
3264	ENDDO
3265	#endif
3266
3267	!-- initial values for t_wall
3268	!-- outer value is set to surface temperature
3269	!-- inner value is set to wall_inner_temperature
3270	!-- and profile is logaritmic (linear in nz).
3271	!-- Horizontal surfaces
3272	DO m = 1, surf_usm_h%ns
3273	!
3274	!-- Roof
3275	IF ( surf_usm_h%isroof_surf(m) ) THEN
3276	tin = roof_inner_temperature
3277	!
3278	!-- Normal land surface
3279	ELSE
3280	tin = soil_inner_temperature
3281	ENDIF
3282
3283	DO k = nzb_wall, nzt_wall+1
3284	c = REAL( k - nzb_wall, wp ) / &
3285	REAL( nzt_wall + 1 - nzb_wall , wp )
3286
3287	t_wall_h(k,m) = ( 1.0_wp - c ) * t_surf_h(m) + c * tin
3288	ENDDO
3289	ENDDO
3290	!
3291	!-- Vertical surfaces
3292	DO l = 0, 3
3293	DO m = 1, surf_usm_v(l)%ns
3294	!
3295	!-- Inner wall
3296	tin = wall_inner_temperature
3297
3298	DO k = nzb_wall, nzt_wall+1
3299	c = REAL( k - nzb_wall, wp ) / &
3300	REAL( nzt_wall + 1 - nzb_wall , wp )
3301
3302	t_wall_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) + &
3303	c * tin
3304	ENDDO
3305	ENDDO
3306	ENDDO
3307
3308	ENDIF
3309
3310	!--
3311	!-- Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
3312	CALL user_init_urban_surface
3313
3314	!-- initialize prognostic values for the first timestep
3315	t_surf_h_p = t_surf_h
3316	t_surf_v_p = t_surf_v
3317
3318	t_wall_h_p = t_wall_h
3319	t_wall_v_p = t_wall_v
3320
3321	!-- Adjust radiative fluxes for urban surface at model start
3322	CALL usm_radiation
3323
3324	CALL cpu_log( log_point_s(78), 'usm_init', 'stop' )
3325
3326
3327	END SUBROUTINE usm_init_urban_surface
3328
3329
3330	!------------------------------------------------------------------------------!
3331	! Description:
3332	! ------------
3333	!
3334	!> Wall model as part of the urban surface model. The model predicts wall
3335	!> temperature.
3336	!------------------------------------------------------------------------------!
3337	SUBROUTINE usm_material_heat_model
3338
3339
3340	IMPLICIT NONE
3341
3342	INTEGER(iwp) :: i,j,k,l,kw, m !< running indices
3343
3344	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend !< tendency
3345
3346	!
3347	!-- For horizontal surfaces
3348	DO m = 1, surf_usm_h%ns
3349	!
3350	!-- Obtain indices
3351	i = surf_usm_h%i(m)
3352	j = surf_usm_h%j(m)
3353	k = surf_usm_h%k(m)
3354	!
3355	!-- prognostic equation for ground/roof temperature t_wall_h
3356	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
3375
3376	t_wall_h_p(nzb_wall:nzt_wall,m) = t_wall_h(nzb_wall:nzt_wall,m) &
3377	+ dt_3d * ( tsc(2) &
3378	* wtend(nzb_wall:nzt_wall) + tsc(3) &
3379	* surf_usm_h%tt_wall_m(nzb_wall:nzt_wall,m) )
3380
3381	!
3382	!-- calculate t_wall tendencies for the next Runge-Kutta step
3383	IF ( timestep_scheme(1:5) == 'runge' ) THEN
3384	IF ( intermediate_timestep_count == 1 ) THEN
3385	DO kw = nzb_wall, nzt_wall
3386	surf_usm_h%tt_wall_m(kw,m) = wtend(kw)
3387	ENDDO
3388	ELSEIF ( intermediate_timestep_count < &
3389	intermediate_timestep_count_max ) THEN
3390	DO kw = nzb_wall, nzt_wall
3391	surf_usm_h%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) + &
3392	5.3125_wp * surf_usm_h%tt_wall_m(kw,m)
3393	ENDDO
3394	ENDIF
3395	ENDIF
3396	ENDDO
3397	!
3398	!-- For vertical surfaces
3399	DO l = 0, 3
3400	DO m = 1, surf_usm_v(l)%ns
3401	!
3402	!-- Obtain indices
3403	i = surf_usm_v(l)%i(m)
3404	j = surf_usm_v(l)%j(m)
3405	k = surf_usm_v(l)%k(m)
3406	!
3407	!-- prognostic equation for wall temperature t_wall_v
3408	wtend(:) = 0.0_wp
3409	wtend(nzb_wall) = (1.0_wp / surf_usm_v(l)%rho_c_wall(nzb_wall,m)) * &
3410	( surf_usm_v(l)%lambda_h(nzb_wall,m) * &
3411	( t_wall_v(l)%t(nzb_wall+1,m) &
3412	- t_wall_v(l)%t(nzb_wall,m) ) * &
3413	surf_usm_v(l)%ddz_wall(nzb_wall+1,m) &
3414	+ surf_usm_v(l)%wghf_eb(m) ) * &
3415	surf_usm_v(l)%ddz_wall_stag(nzb_wall,m)
3416
3417	DO kw = nzb_wall+1, nzt_wall
3418	wtend(kw) = (1.0_wp / surf_usm_v(l)%rho_c_wall(kw,m)) &
3419	* ( surf_usm_v(l)%lambda_h(kw,m) &
3420	* ( t_wall_v(l)%t(kw+1,m) - t_wall_v(l)%t(kw,m) )&
3421	* surf_usm_v(l)%ddz_wall(kw+1,m) &
3422	- surf_usm_v(l)%lambda_h(kw-1,m) &
3423	* ( t_wall_v(l)%t(kw,m) - t_wall_v(l)%t(kw-1,m) )&
3424	* surf_usm_v(l)%ddz_wall(kw,m) &
3425	) * surf_usm_v(l)%ddz_wall_stag(kw,m)
3426	ENDDO
3427
3428	t_wall_v_p(l)%t(nzb_wall:nzt_wall,m) = &
3429	t_wall_v(l)%t(nzb_wall:nzt_wall,m) &
3430	+ dt_3d * ( tsc(2) &
3431	* wtend(nzb_wall:nzt_wall) + tsc(3) &
3432	* surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall,m) )
3433
3434	!
3435	!-- calculate t_wall tendencies for the next Runge-Kutta step
3436	IF ( timestep_scheme(1:5) == 'runge' ) THEN
3437	IF ( intermediate_timestep_count == 1 ) THEN
3438	DO kw = nzb_wall, nzt_wall
3439	surf_usm_v(l)%tt_wall_m(kw,m) = wtend(kw)
3440	ENDDO
3441	ELSEIF ( intermediate_timestep_count < &
3442	intermediate_timestep_count_max ) THEN
3443	DO kw = nzb_wall, nzt_wall
3444	surf_usm_v(l)%tt_wall_m(kw,m) = &
3445	- 9.5625_wp * wtend(kw) + &
3446	5.3125_wp * surf_usm_v(l)%tt_wall_m(kw,m)
3447	ENDDO
3448	ENDIF
3449	ENDIF
3450	ENDDO
3451	ENDDO
3452
3453	END SUBROUTINE usm_material_heat_model
3454
3455
3456	!------------------------------------------------------------------------------!
3457	! Description:
3458	! ------------
3459	!> Parin for &usm_par for urban surface model
3460	!------------------------------------------------------------------------------!
3461	SUBROUTINE usm_parin
3462
3463	IMPLICIT NONE
3464
3465	CHARACTER (LEN=80) :: line !< string containing current line of file PARIN
3466
3467	NAMELIST /urban_surface_par/ &
3468	land_category, &
3469	mrt_factors, &
3470	nrefsteps, &
3471	pedestrant_category, &
3472	ra_horiz_coef, &
3473	read_svf_on_init, &
3474	roof_category, &
3475	split_diffusion_radiation, &
3476	urban_surface, &
3477	usm_anthropogenic_heat, &
3478	usm_energy_balance_land, &
3479	usm_energy_balance_wall, &
3480	usm_material_model, &
3481	usm_lad_rma, &
3482	wall_category, &
3483	write_svf_on_init
3484
3485	line = ' '
3486
3487	!
3488	!-- Try to find urban surface model package
3489	REWIND ( 11 )
3490	line = ' '
3491	DO WHILE ( INDEX( line, '&urban_surface_par' ) == 0 )
3492	READ ( 11, '(A)', END=10 ) line
3493	ENDDO
3494	BACKSPACE ( 11 )
3495
3496	!
3497	!-- Read user-defined namelist
3498	READ ( 11, urban_surface_par )
3499	!
3500	!-- Set flag that indicates that the land surface model is switched on
3501	urban_surface = .TRUE.
3502
3503	!
3504	!-- Activate spinup
3505	IF ( spinup_time > 0.0_wp ) THEN
3506	coupling_start_time = spinup_time
3507	end_time = end_time + spinup_time
3508	IF ( spinup_pt_mean == 9999999.9_wp ) THEN
3509	spinup_pt_mean = pt_surface
3510	ENDIF
3511	spinup = .TRUE.
3512	ENDIF
3513
3514	10 CONTINUE
3515
3516	END SUBROUTINE usm_parin
3517
3518
3519	!------------------------------------------------------------------------------!
3520	! Description:
3521	! ------------
3522	!> This subroutine calculates interaction of the solar radiation
3523	!> with urban surface and updates surface, roofs and walls heatfluxes.
3524	!> It also updates rad_sw_out and rad_lw_out.
3525	!------------------------------------------------------------------------------!
3526	SUBROUTINE usm_radiation
3527
3528	IMPLICIT NONE
3529
3530	INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
3531	INTEGER(iwp) :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
3532	INTEGER(iwp), DIMENSION(4) :: bdycross
3533	REAL(wp), DIMENSION(3,3) :: mrot !< grid rotation matrix (xyz)
3534	REAL(wp), DIMENSION(3,0:9) :: vnorm !< face direction normal vectors (xyz)
3535	REAL(wp), DIMENSION(3) :: sunorig !< grid rotated solar direction unit vector (xyz)
3536	REAL(wp), DIMENSION(3) :: sunorig_grid !< grid squashed solar direction unit vector (zyx)
3537	REAL(wp), DIMENSION(0:9) :: costheta !< direct irradiance factor of solar angle
3538	REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temp tendency
3539	REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove)
3540	REAL(wp) :: rx, ry, rz
3541	REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff
3542	INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy
3543	INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /)
3544
3545	#if ! defined( __nopointer )
3546	IF ( plant_canopy ) THEN
3547	pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
3548	/ (cp * hyp(nzub:nzut) * dxdydz) !< equals to 1 / (rho * c_p * Vbox * T)
3549	ENDIF
3550	#endif
3551	sun_direction = .TRUE.
3552	CALL calc_zenith !< required also for diffusion radiation
3553
3554	!-- prepare rotated normal vectors and irradiance factor
3555	vnorm(1,:) = idir(:)
3556	vnorm(2,:) = jdir(:)
3557	vnorm(3,:) = kdir(:)
3558	mrot(1, :) = (/ cos(alpha), -sin(alpha), 0._wp /)
3559	mrot(2, :) = (/ sin(alpha), cos(alpha), 0._wp /)
3560	mrot(3, :) = (/ 0._wp, 0._wp, 1._wp /)
3561	sunorig = (/ sun_dir_lon, sun_dir_lat, zenith(0) /)
3562	sunorig = matmul(mrot, sunorig)
3563	DO d = 0, 9
3564	costheta(d) = dot_product(sunorig, vnorm(:,d))
3565	ENDDO
3566
3567	IF ( zenith(0) > 0 ) THEN
3568	!-- now we will "squash" the sunorig vector by grid box size in
3569	!-- each dimension, so that this new direction vector will allow us
3570	!-- to traverse the ray path within grid coordinates directly
3571	sunorig_grid = (/ sunorig(3)/dz, sunorig(2)/dy, sunorig(1)/dx /)
3572	!-- sunorig_grid = sunorig_grid / norm2(sunorig_grid)
3573	sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
3574
3575	IF ( plant_canopy ) THEN
3576	!-- precompute effective box depth with prototype Leaf Area Density
3577	pc_box_dimshift = maxloc(sunorig, 1) - 1
3578	CALL usm_box_absorb(cshift((/dx,dy,dz/), pc_box_dimshift), &
3579	60, prototype_lad, &
3580	cshift(sunorig, pc_box_dimshift), &
3581	pc_box_area, pc_abs_frac)
3582	pc_box_area = pc_box_area * sunorig(pc_box_dimshift+1) / sunorig(3)
3583	pc_abs_eff = log(1._wp - pc_abs_frac) / prototype_lad
3584	ENDIF
3585	ENDIF
3586
3587	!-- split diffusion and direct part of the solar downward radiation
3588	!-- comming from radiation model and store it in 2D arrays
3589	!-- rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff
3590	IF ( split_diffusion_radiation ) THEN
3591	CALL usm_calc_diffusion_radiation
3592	ELSE
3593	rad_sw_in_diff = 0.0_wp
3594	rad_sw_in_dir(:,:) = rad_sw_in(0,:,:)
3595	rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
3596	ENDIF
3597
3598	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3599	!-- First pass: direct + diffuse irradiance
3600	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3601	surfinswdir = 0._wp !nsurfl
3602	surfinswdif = 0._wp !nsurfl
3603	surfinlwdif = 0._wp !nsurfl
3604	surfins = 0._wp !nsurfl
3605	surfinl = 0._wp !nsurfl
3606	surfoutsl(:) = 0.0_wp !start-end
3607	surfoutll(:) = 0.0_wp !start-end
3608
3609	!-- Set up thermal radiation from surfaces
3610	!-- emiss_surf is defined only for surfaces for which energy balance is calculated
3611	!-- Workaround: reorder surface data type back on 1D array including all surfaces,
3612	!-- which implies to reorder horizontal and vertical surfaces
3613	!
3614	!-- Horizontal walls
3615	mm = 1
3616	DO i = nxl, nxr
3617	DO j = nys, nyn
3618
3619	DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
3620	surfoutll(mm) = surf_usm_h%emiss_surf(m) * sigma_sb &
3621	* t_surf_h(m)**4
3622	albedo_surf(mm) = surf_usm_h%albedo_surf(m)
3623	emiss_surf(mm) = surf_usm_h%emiss_surf(m)
3624	mm = mm + 1
3625	ENDDO
3626	ENDDO
3627	ENDDO
3628	!
3629	!-- Vertical walls
3630	DO i = nxl, nxr
3631	DO j = nys, nyn
3632	DO ll = 0, 3
3633	l = reorder(ll)
3634	DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
3635	surfoutll(mm) = surf_usm_v(l)%emiss_surf(m) * sigma_sb &
3636	* t_surf_v(l)%t(m)**4
3637	albedo_surf(mm) = surf_usm_v(l)%albedo_surf(m)
3638	emiss_surf(mm) = surf_usm_v(l)%emiss_surf(m)
3639	mm = mm + 1
3640	ENDDO
3641	ENDDO
3642	ENDDO
3643	ENDDO
3644
3645	#if defined( __parallel )
3646	!-- might be optimized and gather only values relevant for current processor
3647
3648	CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, &
3649	surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
3650	#else
3651	surfoutl(:) = surfoutll(:) !nsurf global
3652	#endif
3653
3654	isurf1 = -1 !< previous processed surface
3655	DO isvf = 1, nsvfl
3656	isurf = svfsurf(1, isvf)
3657	k = surfl(iz, isurf)
3658	j = surfl(iy, isurf)
3659	i = surfl(ix, isurf)
3660	isurfsrc = svfsurf(2, isvf)
3661	IF ( zenith(0) > 0 .AND. isurf /= isurf1 ) THEN
3662	!-- locate the virtual surface where the direct solar ray crosses domain boundary
3663	!-- (once per target surface)
3664	d = surfl(id, isurf)
3665	rz = REAL(k, wp) - 0.5_wp * kdir(d)
3666	ry = REAL(j, wp) - 0.5_wp * jdir(d)
3667	rx = REAL(i, wp) - 0.5_wp * idir(d)
3668
3669	CALL usm_find_boundary_face( (/ rz, ry, rx /), sunorig_grid, bdycross)
3670
3671	isurf1 = isurf
3672	ENDIF
3673
3674	IF ( surf(id, isurfsrc) >= isky ) THEN
3675	!-- diffuse rad from boundary surfaces. Since it is a simply
3676	!-- calculated value, it is not assigned to surfref(s/l),
3677	!-- instead it is used directly here
3678	!-- we consider the radiation from the radiation model falling on surface
3679	!-- as the radiation falling on the top of urban layer into the place of the source surface
3680	!-- we consider it as a very reasonable simplification which allow as avoid
3681	!-- necessity of other global range arrays and some all to all mpi communication
3682	surfinswdif(isurf) = surfinswdif(isurf) + rad_sw_in_diff(j,i) * svf(1,isvf) * svf(2,isvf)
3683	!< canopy shading is applied only to shortwave
3684	surfinlwdif(isurf) = surfinlwdif(isurf) + rad_lw_in_diff(j,i) * svf(1,isvf)
3685	ELSE
3686	!-- for surface-to-surface factors we calculate thermal radiation in 1st pass
3687	surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
3688	ENDIF
3689
3690	IF ( zenith(0) > 0 .AND. all( surf(1:4,isurfsrc) == bdycross ) ) THEN
3691	!-- found svf between model boundary and the face => face isn't shaded
3692	surfinswdir(isurf) = rad_sw_in_dir(j,i) &
3693	* costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0)
3694
3695	ENDIF
3696	ENDDO
3697
3698	IF ( plant_canopy ) THEN
3699
3700	pcbinsw(:) = 0._wp
3701	pcbinlw(:) = 0._wp !< will stay always 0 since we don't absorb lw anymore
3702	!
3703	!-- pcsf first pass
3704	isurf1 = -1 !< previous processed pcgb
3705	DO icsf = 1, ncsfl
3706	ipcgb = csfsurf(1, icsf)
3707	i = pcbl(ix,ipcgb)
3708	j = pcbl(iy,ipcgb)
3709	k = pcbl(iz,ipcgb)
3710	isurfsrc = csfsurf(2, icsf)
3711
3712	IF ( zenith(0) > 0 .AND. ipcgb /= isurf1 ) THEN
3713	!-- locate the virtual surface where the direct solar ray crosses domain boundary
3714	!-- (once per target PC gridbox)
3715	rz = REAL(k, wp)
3716	ry = REAL(j, wp)
3717	rx = REAL(i, wp)
3718	CALL usm_find_boundary_face( (/ rz, ry, rx /), &
3719	sunorig_grid, bdycross)
3720
3721	isurf1 = ipcgb
3722	ENDIF
3723
3724	IF ( surf(id, isurfsrc) >= isky ) THEN
3725	!-- Diffuse rad from boundary surfaces. See comments for svf above.
3726	pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * rad_sw_in_diff(j,i)
3727	!-- canopy shading is applied only to shortwave, therefore no absorbtion for lw
3728	!-- pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * rad_lw_in_diff(j,i)
3729	!ELSE
3730	!-- Thermal radiation in 1st pass
3731	!-- pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc)
3732	ENDIF
3733
3734	IF ( zenith(0) > 0 .AND. all( surf(1:4,isurfsrc) == bdycross ) ) THEN
3735	!-- found svf between model boundary and the pcgb => pcgb isn't shaded
3736	pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
3737	pcbinsw(ipcgb) = pcbinsw(ipcgb) &
3738	+ rad_sw_in_dir(j, i) * pc_box_area * csf(2,icsf) * pc_abs_frac
3739	ENDIF
3740	ENDDO
3741	ENDIF
3742
3743	surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy)
3744	surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy)
3745	surfinsw(:) = surfins(:)
3746	surfinlw(:) = surfinl(:)
3747	surfoutsw(:) = 0.0_wp
3748	surfoutlw(:) = surfoutll(:)
3749	! surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
3750	! - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
3751
3752	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3753	!-- Next passes - reflections
3754	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3755	DO refstep = 1, nrefsteps
3756
3757	surfoutsl(startenergy:endenergy) = albedo_surf(startenergy:endenergy) * surfins(startenergy:endenergy)
3758	!-- for non-transparent surfaces, longwave albedo is 1 - emissivity
3759	surfoutll(startenergy:endenergy) = (1._wp - emiss_surf(startenergy:endenergy)) * surfinl(startenergy:endenergy)
3760
3761	#if defined( __parallel )
3762	CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
3763	surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
3764	CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
3765	surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
3766	#else
3767	surfouts(:) = surfoutsl(:)
3768	surfoutl(:) = surfoutll(:)
3769	#endif
3770
3771	!-- reset for next pass input
3772	surfins(:) = 0._wp
3773	surfinl(:) = 0._wp
3774
3775	!-- reflected radiation
3776	DO isvf = 1, nsvfl
3777	isurf = svfsurf(1, isvf)
3778	isurfsrc = svfsurf(2, isvf)
3779
3780	!-- TODO: to remove if, use start+end for isvf
3781	IF ( surf(id, isurfsrc) < isky ) THEN
3782	surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
3783	surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
3784	ENDIF
3785	ENDDO
3786
3787	!-- radiation absorbed by plant canopy
3788	DO icsf = 1, ncsfl
3789	ipcgb = csfsurf(1, icsf)
3790	isurfsrc = csfsurf(2, icsf)
3791
3792	IF ( surf(id, isurfsrc) < isky ) THEN
3793	pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
3794	!-- pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc)
3795	ENDIF
3796	ENDDO
3797
3798	surfinsw(:) = surfinsw(:) + surfins(:)
3799	surfinlw(:) = surfinlw(:) + surfinl(:)
3800	surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy)
3801	surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy)
3802	! surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
3803	! - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
3804
3805	ENDDO
3806
3807	!-- push heat flux absorbed by plant canopy to respective 3D arrays
3808	#if ! defined( __nopointer )
3809	IF ( plant_canopy ) THEN
3810	pc_heating_rate(:,:,:) = 0._wp
3811	DO ipcgb = 1, npcbl
3812	j = pcbl(iy, ipcgb)
3813	i = pcbl(ix, ipcgb)
3814	k = pcbl(iz, ipcgb)
3815	!
3816	!-- Following expression equals former kk = k - nzb_s_inner(j,i)
3817	kk = k - ( get_topography_top_index( j, i, 's' ) ) !- lad arrays are defined flat
3818	pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
3819	* pchf_prep(k) * pt(k, j, i) !-- = dT/dt
3820	ENDDO
3821	ENDIF
3822	#endif
3823	!
3824	!-- Transfer radiation arrays required for energy balance to the respective data types
3825	DO i = startenergy, endenergy
3826	m = surfl(5,i)
3827	!
3828	!-- upward-facing
3829	IF ( surfl(1,i) == 0 ) THEN
3830	surf_usm_h%rad_in_sw(m) = surfinsw(i)
3831	surf_usm_h%rad_out_sw(m) = surfoutsw(i)
3832	surf_usm_h%rad_in_lw(m) = surfinlw(i)
3833	surf_usm_h%rad_out_lw(m) = surfoutlw(i)
3834	!
3835	!-- southward-facding
3836	ELSEIF ( surfl(1,i) == 1 ) THEN
3837	surf_usm_v(1)%rad_in_sw(m) = surfinsw(i)
3838	surf_usm_v(1)%rad_out_sw(m) = surfoutsw(i)
3839	surf_usm_v(1)%rad_in_lw(m) = surfinlw(i)
3840	surf_usm_v(1)%rad_out_lw(m) = surfoutlw(i)
3841	!
3842	!-- northward-facding
3843	ELSEIF ( surfl(1,i) == 2 ) THEN
3844	surf_usm_v(0)%rad_in_sw(m) = surfinsw(i)
3845	surf_usm_v(0)%rad_out_sw(m) = surfoutsw(i)
3846	surf_usm_v(0)%rad_in_lw(m) = surfinlw(i)
3847	surf_usm_v(0)%rad_out_lw(m) = surfoutlw(i)
3848	!
3849	!-- westward-facding
3850	ELSEIF ( surfl(1,i) == 3 ) THEN
3851	surf_usm_v(3)%rad_in_sw(m) = surfinsw(i)
3852	surf_usm_v(3)%rad_out_sw(m) = surfoutsw(i)
3853	surf_usm_v(3)%rad_in_lw(m) = surfinlw(i)
3854	surf_usm_v(3)%rad_out_lw(m) = surfoutlw(i)
3855	!
3856	!-- eastward-facing
3857	ELSEIF ( surfl(1,i) == 4 ) THEN
3858	surf_usm_v(2)%rad_in_sw(m) = surfinsw(i)
3859	surf_usm_v(2)%rad_out_sw(m) = surfoutsw(i)
3860	surf_usm_v(2)%rad_in_lw(m) = surfinlw(i)
3861	surf_usm_v(2)%rad_out_lw(m) = surfoutlw(i)
3862	ENDIF
3863
3864	ENDDO
3865
3866
3867	DO m = 1, surf_usm_h%ns
3868	surf_usm_h%surfhf(m) = surf_usm_h%rad_in_sw(m) + &
3869	surf_usm_h%rad_in_lw(m) - &
3870	surf_usm_h%rad_out_sw(m) - &
3871	surf_usm_h%rad_out_lw(m)
3872	ENDDO
3873
3874	DO l = 0, 3
3875	DO m = 1, surf_usm_v(l)%ns
3876	surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_in_sw(m) + &
3877	surf_usm_v(l)%rad_in_lw(m) - &
3878	surf_usm_v(l)%rad_out_sw(m) - &
3879	surf_usm_v(l)%rad_out_lw(m)
3880	ENDDO
3881	ENDDO
3882
3883	!-- return surface radiation to horizontal surfaces
3884	!-- to rad_sw_in, rad_lw_in and rad_net for outputs
3885	!!!!!!!!!!
3886	!-- we need the original radiation on urban top layer
3887	!-- for calculation of MRT so we can't do adjustment here for now
3888	!!!!!!!!!!
3889	!!!DO isurf = 1, nsurfl
3890	!!! i = surfl(ix,isurf)
3891	!!! j = surfl(iy,isurf)
3892	!!! k = surfl(iz,isurf)
3893	!!! d = surfl(id,isurf)
3894	!!! IF ( d==iroof ) THEN
3895	!!! rad_sw_in(:,j,i) = surfinsw(isurf)
3896	!!! rad_lw_in(:,j,i) = surfinlw(isurf)
3897	!!! rad_net(j,i) = rad_sw_in(k,j,i) - rad_sw_out(k,j,i) + rad_lw_in(k,j,i) - rad_lw_out(k,j,i)
3898	!!! ENDIF
3899	!!!ENDDO
3900
3901	END SUBROUTINE usm_radiation
3902
3903
3904	!------------------------------------------------------------------------------!
3905	! Description:
3906	! ------------
3907	!> Raytracing for detecting obstacles and calculating compound canopy sink
3908	!> factors. (A simple obstacle detection would only need to process faces in
3909	!> 3 dimensions without any ordering.)
3910	!> Assumtions:
3911	!> -----------
3912	!> 1. The ray always originates from a face midpoint (only one coordinate equals
3913	!> *.5, i.e. wall) and doesn't travel parallel to the surface (that would mean
3914	!> shape factor=0). Therefore, the ray may never travel exactly along a face
3915	!> or an edge.
3916	!> 2. From grid bottom to urban surface top the grid has to be equidistant
3917	!> within each of the dimensions, including vertical (but the resolution
3918	!> doesn't need to be the same in all three dimensions).
3919	!------------------------------------------------------------------------------!
3920	SUBROUTINE usm_raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency, win_lad)
3921	IMPLICIT NONE
3922
3923	REAL(wp), DIMENSION(3), INTENT(in) :: src, targ !< real coordinates z,y,x
3924	INTEGER(iwp), INTENT(in) :: isrc !< index of source face for csf
3925	REAL(wp), INTENT(in) :: rirrf !< irradiance factor for csf
3926	REAL(wp), INTENT(in) :: atarg !< target surface area for csf
3927	LOGICAL, INTENT(in) :: create_csf !< whether to generate new CSFs during raytracing
3928	LOGICAL, INTENT(out) :: visible
3929	REAL(wp), INTENT(out) :: transparency !< along whole path
3930	INTEGER(iwp), INTENT(in) :: win_lad
3931	INTEGER(iwp) :: i, j, k, d
3932	INTEGER(iwp) :: seldim !< dimension to be incremented
3933	INTEGER(iwp) :: ncsb !< no of written plant canopy sinkboxes
3934	INTEGER(iwp) :: maxboxes !< max no of gridboxes visited
3935	REAL(wp) :: distance !< euclidean along path
3936	REAL(wp) :: crlen !< length of gridbox crossing
3937	REAL(wp) :: lastdist !< beginning of current crossing
3938	REAL(wp) :: nextdist !< end of current crossing
3939	REAL(wp) :: realdist !< distance in meters per unit distance
3940	REAL(wp) :: crmid !< midpoint of crossing
3941	REAL(wp) :: cursink !< sink factor for current canopy box
3942	REAL(wp), DIMENSION(3) :: delta !< path vector
3943	REAL(wp), DIMENSION(3) :: uvect !< unit vector
3944	REAL(wp), DIMENSION(3) :: dimnextdist !< distance for each dimension increments
3945	INTEGER(iwp), DIMENSION(3) :: box !< gridbox being crossed
3946	INTEGER(iwp), DIMENSION(3) :: dimnext !< next dimension increments along path
3947	INTEGER(iwp), DIMENSION(3) :: dimdelta !< dimension direction = +- 1
3948	INTEGER(iwp) :: px, py !< number of processors in x and y dir before
3949	!< the processor in the question
3950	INTEGER(iwp) :: ip !< number of processor where gridbox reside
3951	INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array
3952	REAL(wp) :: lad_s_target !< recieved lad_s of particular grid box
3953	REAL(wp), PARAMETER :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
3954
3955	!
3956	!-- Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
3957	!-- the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
3958	maxboxes = SUM(ABS(NINT(targ) - NINT(src))) + 1
3959	IF ( plant_canopy .AND. ncsfl + maxboxes > ncsfla ) THEN
3960	!-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
3961	!-- k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
3962	!-- / log(grow_factor)), kind=wp))
3963	!-- or use this code to simply always keep some extra space after growing
3964	k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
3965
3966	CALL usm_merge_and_grow_csf(k)
3967	ENDIF
3968
3969	transparency = 1._wp
3970	ncsb = 0
3971
3972	delta(:) = targ(:) - src(:)
3973	distance = SQRT(SUM(delta(:)**2))
3974	IF ( distance == 0._wp ) THEN
3975	visible = .TRUE.
3976	RETURN
3977	ENDIF
3978	uvect(:) = delta(:) / distance
3979	realdist = SQRT(SUM( (uvect(:)(/dz,dy,dx/))*2 ))
3980
3981	lastdist = 0._wp
3982
3983	!-- Since all face coordinates have values *.5 and we'd like to use
3984	!-- integers, all these have .5 added
3985	DO d = 1, 3
3986	IF ( uvect(d) == 0._wp ) THEN
3987	dimnext(d) = 999999999
3988	dimdelta(d) = 999999999
3989	dimnextdist(d) = 1.0E20_wp
3990	ELSE IF ( uvect(d) > 0._wp ) THEN
3991	dimnext(d) = CEILING(src(d) + .5_wp)
3992	dimdelta(d) = 1
3993	dimnextdist(d) = (dimnext(d) - .5_wp - src(d)) / uvect(d)
3994	ELSE
3995	dimnext(d) = FLOOR(src(d) + .5_wp)
3996	dimdelta(d) = -1
3997	dimnextdist(d) = (dimnext(d) - .5_wp - src(d)) / uvect(d)
3998	ENDIF
3999	ENDDO
4000
4001	DO
4002	!-- along what dimension will the next wall crossing be?
4003	seldim = minloc(dimnextdist, 1)
4004	nextdist = dimnextdist(seldim)
4005	IF ( nextdist > distance ) nextdist = distance
4006
4007	crlen = nextdist - lastdist
4008	IF ( crlen > .001_wp ) THEN
4009	crmid = (lastdist + nextdist) * .5_wp
4010	box = NINT(src(:) + uvect(:) * crmid)
4011
4012	!-- calculate index of the grid with global indices (box(2),box(3))
4013	!-- in the array nzterr and plantt and id of the coresponding processor
4014	px = box(3)/nnx
4015	py = box(2)/nny
4016	ip = px*pdims(2)+py
4017	ig = ipnnxnny + (box(3)-pxnnx)nny + box(2)-py*nny
4018	IF ( box(1) <= nzterr(ig) ) THEN
4019	visible = .FALSE.
4020	RETURN
4021	ENDIF
4022
4023	IF ( plant_canopy ) THEN
4024	IF ( box(1) <= plantt(ig) ) THEN
4025	ncsb = ncsb + 1
4026	boxes(:,ncsb) = box
4027	crlens(ncsb) = crlen
4028	#if defined( __parallel )
4029	lad_ip(ncsb) = ip
4030	lad_disp(ncsb) = (box(3)-pxnnx)(nnynzu) + (box(2)-pynny)*nzu + box(1)-nzub
4031	#endif
4032	ENDIF
4033	ENDIF
4034	ENDIF
4035
4036	IF ( nextdist >= distance ) EXIT
4037	lastdist = nextdist
4038	dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
4039	dimnextdist(seldim) = (dimnext(seldim) - .5_wp - src(seldim)) / uvect(seldim)
4040	ENDDO
4041
4042	IF ( plant_canopy ) THEN
4043	#if defined( __parallel )
4044	IF ( usm_lad_rma ) THEN
4045	!-- send requests for lad_s to appropriate processor
4046	CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
4047	DO i = 1, ncsb
4048	CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
4049	1, MPI_REAL, win_lad, ierr)
4050	IF ( ierr /= 0 ) THEN
4051	WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
4052	CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
4053	ENDIF
4054	ENDDO
4055
4056	!-- wait for all pending local requests complete
4057	CALL MPI_Win_flush_local_all(win_lad, ierr)
4058	IF ( ierr /= 0 ) THEN
4059	WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
4060	CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
4061	ENDIF
4062	CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
4063
4064	ENDIF
4065	#endif
4066
4067	!-- calculate csf and transparency
4068	DO i = 1, ncsb
4069	#if defined( __parallel )
4070	IF ( usm_lad_rma ) THEN
4071	lad_s_target = lad_s_ray(i)
4072	ELSE
4073	lad_s_target = usm_lad_g(lad_ip(i)nnxnny*nzu + lad_disp(i))
4074	ENDIF
4075	#else
4076	lad_s_target = usm_lad(boxes(1,i),boxes(2,i),boxes(3,i))
4077	#endif
4078	cursink = 1._wp - exp(-ext_coef * lad_s_target * crlens(i)*realdist)
4079
4080	IF ( create_csf ) THEN
4081	!-- write svf values into the array
4082	ncsfl = ncsfl + 1
4083	acsf(ncsfl)%ip = lad_ip(i)
4084	acsf(ncsfl)%itx = boxes(3,i)
4085	acsf(ncsfl)%ity = boxes(2,i)
4086	acsf(ncsfl)%itz = boxes(1,i)
4087	acsf(ncsfl)%isurfs = isrc
4088	acsf(ncsfl)%rsvf = REAL(cursinkrirrfatarg, wp) !-- we postpone multiplication by transparency
4089	acsf(ncsfl)%rtransp = REAL(transparency, wp)
4090	ENDIF !< create_csf
4091
4092	transparency = transparency * (1._wp - cursink)
4093
4094	ENDDO
4095	ENDIF
4096
4097	visible = .TRUE.
4098
4099	END SUBROUTINE usm_raytrace
4100
4101
4102	!------------------------------------------------------------------------------!
4103	! Description:
4104	! ------------
4105	!
4106	!> This subroutine is part of the urban surface model.
4107	!> It reads daily heat produced by anthropogenic sources
4108	!> and the diurnal cycle of the heat.
4109	!------------------------------------------------------------------------------!
4110	SUBROUTINE usm_read_anthropogenic_heat
4111
4112	INTEGER(iwp) :: i,j,ii
4113	REAL(wp) :: heat
4114
4115	!-- allocation of array of sources of anthropogenic heat and their diural profile
4116	ALLOCATE( aheat(nys:nyn,nxl:nxr) )
4117	ALLOCATE( aheatprof(0:24) )
4118
4119	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4120	!-- read daily amount of heat and its daily cycle
4121	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4122	aheat = 0.0_wp
4123	DO ii = 0, io_blocks-1
4124	IF ( ii == io_group ) THEN
4125
4126	!-- open anthropogenic heat file
4127	OPEN( 151, file='ANTHROPOGENIC_HEAT'//TRIM(coupling_char), action='read', &
4128	status='old', form='formatted', err=11 )
4129	i = 0
4130	j = 0
4131	DO
4132	READ( 151, *, err=12, end=13 ) i, j, heat
4133	IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn ) THEN
4134	!-- write heat into the array
4135	aheat(j,i) = heat
4136	ENDIF
4137	CYCLE
4138	12 WRITE(message_string,'(a,2i4)') 'error in file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' after line ',i,j
4139	CALL message( 'usm_read_anthropogenic_heat', 'PA0515', 0, 1, 0, 6, 0 )
4140	ENDDO
4141	13 CLOSE(151)
4142	CYCLE
4143	11 message_string = 'file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' does not exist'
4144	CALL message( 'usm_read_anthropogenic_heat', 'PA0516', 1, 2, 0, 6, 0 )
4145	ENDIF
4146
4147	#if defined( __parallel ) && ! defined ( __check )
4148	CALL MPI_BARRIER( comm2d, ierr )
4149	#endif
4150	ENDDO
4151
4152	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4153	!-- read diurnal profiles of heat sources
4154	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4155	aheatprof = 0.0_wp
4156	DO ii = 0, io_blocks-1
4157	IF ( ii == io_group ) THEN
4158
4159	!-- open anthropogenic heat profile file
4160	OPEN( 151, file='ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char), action='read', &
4161	status='old', form='formatted', err=21 )
4162	i = 0
4163	DO
4164	READ( 151, *, err=22, end=23 ) i, heat
4165	IF ( i >= 0 .AND. i <= 24 ) THEN
4166	!-- write heat into the array
4167	aheatprof(i) = heat
4168	ENDIF
4169	CYCLE
4170	22 WRITE(message_string,'(a,i4)') 'error in file ANTHROPOGENIC_HEAT_PROFILE'// &
4171	TRIM(coupling_char)//' after line ',i
4172	CALL message( 'usm_read_anthropogenic_heat', 'PA0517', 0, 1, 0, 6, 0 )
4173	ENDDO
4174	aheatprof(24) = aheatprof(0)
4175	23 CLOSE(151)
4176	CYCLE
4177	21 message_string = 'file ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char)//' does not exist'
4178	CALL message( 'usm_read_anthropogenic_heat', 'PA0518', 1, 2, 0, 6, 0 )
4179	ENDIF
4180
4181	#if defined( __parallel ) && ! defined ( __check )
4182	CALL MPI_BARRIER( comm2d, ierr )
4183	#endif
4184	ENDDO
4185
4186	END SUBROUTINE usm_read_anthropogenic_heat
4187
4188
4189	!------------------------------------------------------------------------------!
4190	!
4191	! Description:
4192	! ------------
4193	!> Soubroutine reads t_surf and t_wall data from restart files
4194	!kanani: Renamed this routine according to corresponging routines in PALM
4195	!kanani: Modified the routine to match read_var_list, from where usm_read_restart_data
4196	! shall be called in the future. This part has not been tested yet. (see virtual_flight_mod)
4197	! Also, I had some trouble with the allocation of t_surf, since this is a pointer.
4198	! So, I added some directives here.
4199	!------------------------------------------------------------------------------!
4200	SUBROUTINE usm_read_restart_data( ii, &
4201	nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
4202	nynfa, nyn_on_file, nysfa, nys_on_file, &
4203	offset_xa, offset_ya, overlap_count )
4204
4205
4206	USE pegrid, &
4207	ONLY: numprocs_previous_run
4208
4209	IMPLICIT NONE
4210
4211	CHARACTER (LEN=1) :: dum !< dummy to create correct string for reading input variable
4212	CHARACTER (LEN=30) :: field_chr !< input variable
4213
4214	INTEGER(iwp) :: l !< index variable for surface type
4215	INTEGER(iwp) :: ii !< running index over input files
4216	INTEGER(iwp) :: kk !< running index over previous input files covering current local domain
4217	INTEGER(iwp) :: ns_h_on_file_usm !< number of horizontal surface elements (urban type) on file
4218	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
4219	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
4220	INTEGER(iwp) :: nxl_on_file !< index of left boundary on former local domain
4221	INTEGER(iwp) :: nxrc !< index of right boundary on current subdomain
4222	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
4223	INTEGER(iwp) :: nxr_on_file !< index of right boundary on former local domain
4224	INTEGER(iwp) :: nync !< index of north boundary on current subdomain
4225	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
4226	INTEGER(iwp) :: nyn_on_file !< index of norht boundary on former local domain
4227	INTEGER(iwp) :: nysc !< index of south boundary on current subdomain
4228	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
4229	INTEGER(iwp) :: nys_on_file !< index of south boundary on former local domain
4230	INTEGER(iwp) :: overlap_count !< number of overlaps
4231
4232	INTEGER(iwp) :: ns_v_on_file_usm(0:3) !< number of vertical surface elements (urban type) on file
4233
4234	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !<
4235	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !<
4236	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !<
4237	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !<
4238	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !<
4239	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !<
4240
4241	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: start_index_on_file
4242	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: end_index_on_file
4243
4244	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp_surf_h
4245	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tmp_wall_h
4246
4247	TYPE( t_surf_vertical ), DIMENSION(0:3) :: tmp_surf_v
4248	TYPE( t_wall_vertical ), DIMENSION(0:3) :: tmp_wall_v
4249
4250
4251	IF ( initializing_actions == 'read_restart_data' .OR. &
4252	initializing_actions == 'cyclic_fill' ) THEN
4253
4254	!
4255	!-- Read number of respective surface elements on file
4256	READ ( 13 ) field_chr
4257	IF ( TRIM( field_chr ) /= 'ns_h_on_file_usm' ) THEN
4258	!
4259	!-- Add a proper error message
4260	ENDIF
4261	READ ( 13 ) ns_h_on_file_usm
4262
4263	READ ( 13 ) field_chr
4264	IF ( TRIM( field_chr ) /= 'ns_v_on_file_usm' ) THEN
4265	!
4266	!-- Add a proper error message
4267	ENDIF
4268	READ ( 13 ) ns_v_on_file_usm
4269	!
4270	!-- Allocate temporary arrays for reading data on file. Note, the
4271	!-- size of allocated surface elements do not necessarily need to match
4272	!-- the size of present surface elements on current processor, as the
4273	!-- number of processors between restarts can change.
4274	ALLOCATE( tmp_surf_h(1:ns_h_on_file_usm) )
4275	ALLOCATE( tmp_wall_h(nzb_wall:nzt_wall+1,1:ns_h_on_file_usm) )
4276
4277	DO l = 0, 3
4278	ALLOCATE( tmp_surf_v(l)%t(1:ns_v_on_file_usm(l)) )
4279	ALLOCATE( tmp_wall_v(l)%t(nzb_wall:nzt_wall+1,1:ns_v_on_file_usm(l) ) )
4280	ENDDO
4281
4282
4283	READ ( 13 ) field_chr
4284
4285	DO WHILE ( TRIM( field_chr ) /= '* end usm *' )
4286	!
4287	!-- Map data on file as often as needed (data are read only for k=1)
4288	DO kk = 1, overlap_count
4289	!
4290	!-- Get the index range of the subdomain on file which overlap with the
4291	!-- current subdomain
4292	nxlf = nxlfa(ii,kk)
4293	nxlc = nxlfa(ii,kk) + offset_xa(ii,kk)
4294	nxrf = nxrfa(ii,kk)
4295	nxrc = nxrfa(ii,kk) + offset_xa(ii,kk)
4296	nysf = nysfa(ii,kk)
4297	nysc = nysfa(ii,kk) + offset_ya(ii,kk)
4298	nynf = nynfa(ii,kk)
4299	nync = nynfa(ii,kk) + offset_ya(ii,kk)
4300
4301	SELECT CASE ( TRIM( field_chr ) )
4302
4303	CASE ( 'usm_start_index_h', 'usm_start_index_v' )
4304	IF ( kk == 1 ) &
4305	READ ( 13 ) start_index_on_file
4306
4307	CASE ( 'usm_end_index_h', 'usm_end_index_v' )
4308	IF ( kk == 1 ) &
4309	READ ( 13 ) end_index_on_file
4310
4311	CASE ( 't_surf_h' )
4312	#if defined( __nopointer )
4313	IF ( kk == 1 ) THEN
4314	IF ( .NOT. ALLOCATED( t_surf_h ) ) &
4315	ALLOCATE( t_surf_h(1:surf_usm_h%ns) )
4316	READ ( 13 ) tmp_surf_h
4317	ENDIF
4318	CALL restore_surface_elements_usm_1d( &
4319	t_surf_h, tmp_surf_h, &
4320	surf_usm_h%start_index )
4321	#else
4322	IF ( kk == 1 ) THEN
4323	IF ( .NOT. ALLOCATED( t_surf_h_1 ) ) &
4324	ALLOCATE( t_surf_h_1(1:surf_usm_h%ns) )
4325	READ ( 13 ) tmp_surf_h
4326	ENDIF
4327	CALL restore_surface_elements_usm_1d( &
4328	t_surf_h_1, tmp_surf_h, &
4329	surf_usm_h%start_index )
4330	#endif
4331
4332	CASE ( 't_surf_v(0)' )
4333	#if defined( __nopointer )
4334	IF ( kk == 1 ) THEN
4335	IF ( .NOT. ALLOCATED( t_surf_v(0)%t ) ) &
4336	ALLOCATE( t_surf_v(0)%t(1:surf_usm_v(0)%ns) )
4337	READ ( 13 ) tmp_surf_v(0)%t
4338	ENDIF
4339	CALL restore_surface_elements_usm_1d( &
4340	t_surf_v(0)%t, tmp_surf_v(0)%t, &
4341	surf_usm_v(0)%start_index )
4342	#else
4343	IF ( kk == 1 ) THEN
4344	IF ( .NOT. ALLOCATED( t_surf_v_1(0)%t ) ) &
4345	ALLOCATE( t_surf_v_1(0)%t(1:surf_usm_v(0)%ns) )
4346	READ ( 13 ) tmp_surf_v(0)%t
4347	ENDIF
4348	CALL restore_surface_elements_usm_1d( &
4349	t_surf_v_1(0)%t, tmp_surf_v(0)%t,&
4350	surf_usm_v(0)%start_index )
4351	#endif
4352
4353	CASE ( 't_surf_v(1)' )
4354	#if defined( __nopointer )
4355	IF ( kk == 1 ) THEN
4356	IF ( .NOT. ALLOCATED( t_surf_v(1)%t ) ) &
4357	ALLOCATE( t_surf_v(1)%t(1:surf_usm_v(1)%ns) )
4358	READ ( 13 ) tmp_surf_v(1)%t
4359	ENDIF
4360	CALL restore_surface_elements_usm_1d( &
4361	t_surf_v(1)%t, tmp_surf_v(1)%t, &
4362	surf_usm_v(1)%start_index )
4363	#else
4364	IF ( kk == 1 ) THEN
4365	IF ( .NOT. ALLOCATED( t_surf_v_1(1)%t ) ) &
4366	ALLOCATE( t_surf_v_1(1)%t(1:surf_usm_v(1)%ns) )
4367	READ ( 13 ) tmp_surf_v(1)%t
4368	ENDIF
4369	CALL restore_surface_elements_usm_1d( &
4370	t_surf_v_1(1)%t, tmp_surf_v(1)%t,&
4371	surf_usm_v(1)%start_index )
4372	#endif
4373
4374	CASE ( 't_surf_v(2)' )
4375	#if defined( __nopointer )
4376	IF ( kk == 1 ) THEN
4377	IF ( .NOT. ALLOCATED( t_surf_v(2)%t ) ) &
4378	ALLOCATE( t_surf_v(2)%t(1:surf_usm_v(2)%ns) )
4379	READ ( 13 ) tmp_surf_v(2)%t
4380	ENDIF
4381	CALL restore_surface_elements_usm_1d( &
4382	t_surf_v(2)%t, tmp_surf_v(2)%t, &
4383	surf_usm_v(2)%start_index )
4384	#else
4385	IF ( kk == 1 ) THEN
4386	IF ( .NOT. ALLOCATED( t_surf_v_1(2)%t ) ) &
4387	ALLOCATE( t_surf_v_1(2)%t(1:surf_usm_v(2)%ns) )
4388	READ ( 13 ) tmp_surf_v(2)%t
4389	ENDIF
4390	CALL restore_surface_elements_usm_1d( &
4391	t_surf_v_1(2)%t, tmp_surf_v(2)%t,&
4392	surf_usm_v(2)%start_index )
4393	#endif
4394
4395	CASE ( 't_surf_v(3)' )
4396	#if defined( __nopointer )
4397	IF ( kk == 1 ) THEN
4398	IF ( .NOT. ALLOCATED( t_surf_v(3)%t ) ) &
4399	ALLOCATE( t_surf_v(3)%t(1:surf_usm_v(3)%ns) )
4400	READ ( 13 ) tmp_surf_v(3)%t
4401	ENDIF
4402	CALL restore_surface_elements_usm_1d( &
4403	t_surf_v(3)%t, tmp_surf_v(3)%t, &
4404	surf_usm_v(3)%start_index )
4405	#else
4406	IF ( kk == 1 ) THEN
4407	IF ( .NOT. ALLOCATED( t_surf_v_1(3)%t ) ) &
4408	ALLOCATE( t_surf_v_1(3)%t(1:surf_usm_v(3)%ns) )
4409	READ ( 13 ) tmp_surf_v(3)%t
4410	ENDIF
4411	CALL restore_surface_elements_usm_1d( &
4412	t_surf_v_1(3)%t, tmp_surf_v(3)%t,&
4413	surf_usm_v(3)%start_index )
4414	#endif
4415	CASE ( 't_wall_h' )
4416	#if defined( __nopointer )
4417	IF ( kk == 1 ) THEN
4418	IF ( .NOT. ALLOCATED( t_wall_h ) ) &
4419	ALLOCATE( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
4420	READ ( 13 ) tmp_wall_h
4421	ENDIF
4422	CALL restore_surface_elements_usm_2d( &
4423	t_wall_h, tmp_wall_h, &
4424	surf_usm_h%start_index )
4425	#else
4426	IF ( kk == 1 ) THEN
4427	IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) &
4428	ALLOCATE( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
4429	READ ( 13 ) tmp_wall_h
4430	ENDIF
4431	CALL restore_surface_elements_usm_2d( &
4432	t_wall_h_1, tmp_wall_h, &
4433	surf_usm_h%start_index )
4434	#endif
4435	CASE ( 't_wall_v(0)' )
4436	#if defined( __nopointer )
4437	IF ( kk == 1 ) THEN
4438	IF ( .NOT. ALLOCATED( t_wall_v(0)%t ) ) &
4439	ALLOCATE( t_wall_v(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
4440	READ ( 13 ) tmp_wall_v(0)%t
4441	ENDIF
4442	CALL restore_surface_elements_usm_2d( &
4443	t_wall_v(0)%t, tmp_wall_v(0)%t, &
4444	surf_usm_v(0)%start_index )
4445	#else
4446	IF ( kk == 1 ) THEN
4447	IF ( .NOT. ALLOCATED( t_wall_v_1(0)%t ) ) &
4448	ALLOCATE( t_wall_v_1(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
4449	READ ( 13 ) tmp_wall_v(0)%t
4450	ENDIF
4451	CALL restore_surface_elements_usm_2d( &
4452	t_wall_v_1(0)%t, tmp_wall_v(0)%t,&
4453	surf_usm_v(0)%start_index )
4454	#endif
4455	CASE ( 't_wall_v(1)' )
4456	#if defined( __nopointer )
4457	IF ( kk == 1 ) THEN
4458	IF ( .NOT. ALLOCATED( t_wall_v(1)%t ) ) &
4459	ALLOCATE( t_wall_v(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
4460	READ ( 13 ) tmp_wall_v(1)%t
4461	ENDIF
4462	CALL restore_surface_elements_usm_2d( &
4463	t_wall_v(1)%t, tmp_wall_v(1)%t, &
4464	surf_usm_v(1)%start_index )
4465	#else
4466	IF ( kk == 1 ) THEN
4467	IF ( .NOT. ALLOCATED( t_wall_v_1(1)%t ) ) &
4468	ALLOCATE( t_wall_v_1(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
4469	READ ( 13 ) tmp_wall_v(1)%t
4470	ENDIF
4471	CALL restore_surface_elements_usm_2d( &
4472	t_wall_v_1(1)%t, tmp_wall_v(1)%t,&
4473	surf_usm_v(1)%start_index )
4474	#endif
4475	CASE ( 't_wall_v(2)' )
4476	#if defined( __nopointer )
4477	IF ( kk == 1 ) THEN
4478	IF ( .NOT. ALLOCATED( t_wall_v(2)%t ) ) &
4479	ALLOCATE( t_wall_v(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
4480	READ ( 13 ) tmp_wall_v(2)%t
4481	ENDIF
4482	CALL restore_surface_elements_usm_2d( &
4483	t_wall_v(2)%t, tmp_wall_v(2)%t, &
4484	surf_usm_v(2)%start_index )
4485	#else
4486	IF ( kk == 1 ) THEN
4487	IF ( .NOT. ALLOCATED( t_wall_v_1(2)%t ) ) &
4488	ALLOCATE( t_wall_v_1(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
4489	READ ( 13 ) tmp_wall_v(2)%t
4490	ENDIF
4491	CALL restore_surface_elements_usm_2d( &
4492	t_wall_v_1(2)%t, tmp_wall_v(2)%t,&
4493	surf_usm_v(2)%start_index )
4494	#endif
4495	CASE ( 't_wall_v(3)' )
4496	#if defined( __nopointer )
4497	IF ( kk == 1 ) THEN
4498	IF ( .NOT. ALLOCATED( t_wall_v(3)%t ) ) &
4499	ALLOCATE( t_wall_v(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
4500	READ ( 13 ) tmp_wall_v(3)%t
4501	ENDIF
4502	CALL restore_surface_elements_usm_2d( &
4503	t_wall_v(3)%t, tmp_wall_v(3)%t, &
4504	surf_usm_v(3)%start_index )
4505	#else
4506	IF ( kk == 1 ) THEN
4507	IF ( .NOT. ALLOCATED( t_wall_v_1(3)%t ) ) &
4508	ALLOCATE( t_wall_v_1(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
4509	READ ( 13 ) tmp_wall_v(3)%t
4510	ENDIF
4511	CALL restore_surface_elements_usm_2d( &
4512	t_wall_v_1(3)%t, tmp_wall_v(3)%t,&
4513	surf_usm_v(3)%start_index )
4514	#endif
4515
4516	CASE DEFAULT
4517	WRITE ( message_string, * ) 'unknown variable named "', &
4518	TRIM( field_chr ), '" found in', &
4519	'&data from prior run on PE ', myid
4520	CALL message( 'user_read_restart_data', 'UI0012', 1, 2, 0, 6, 0 )
4521
4522	END SELECT
4523
4524	ENDDO
4525
4526	READ ( 13 ) field_chr
4527
4528	ENDDO
4529
4530	ENDIF
4531
4532	CONTAINS
4533
4534	SUBROUTINE restore_surface_elements_usm_1d( surf_target, surf_file, start_index_c )
4535
4536	IMPLICIT NONE
4537
4538	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
4539	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
4540	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
4541	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
4542	INTEGER(iwp) :: m !< surface-element index on file
4543	INTEGER(iwp) :: mm !< surface-element index on current subdomain
4544
4545	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: start_index_c
4546
4547	REAL(wp), DIMENSION(:) :: surf_target !< target surface type
4548	REAL(wp), DIMENSION(:) :: surf_file !< surface type on file
4549
4550	ic = nxlc
4551	DO i = nxlf, nxrf
4552	jc = nysc
4553	DO j = nysf, nynf
4554
4555	mm = start_index_c(jc,ic)
4556	DO m = start_index_on_file(j,i), end_index_on_file(j,i)
4557	surf_target(mm) = surf_file(m)
4558	mm = mm + 1
4559	ENDDO
4560
4561	jc = jc + 1
4562	ENDDO
4563	ic = ic + 1
4564	ENDDO
4565
4566
4567	END SUBROUTINE restore_surface_elements_usm_1d
4568
4569	SUBROUTINE restore_surface_elements_usm_2d( surf_target, surf_file, start_index_c )
4570
4571	IMPLICIT NONE
4572
4573	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
4574	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
4575	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
4576	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
4577	INTEGER(iwp) :: m !< surface-element index on file
4578	INTEGER(iwp) :: mm !< surface-element index on current subdomain
4579
4580	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: start_index_c
4581
4582	REAL(wp), DIMENSION(:,:) :: surf_target !< target surface type
4583	REAL(wp), DIMENSION(:,:) :: surf_file !< surface type on file
4584
4585	ic = nxlc
4586	DO i = nxlf, nxrf
4587	jc = nysc
4588	DO j = nysf, nynf
4589
4590	mm = start_index_c(jc,ic)
4591	DO m = start_index_on_file(j,i), end_index_on_file(j,i)
4592	surf_target(:,mm) = surf_file(:,m)
4593	mm = mm + 1
4594	ENDDO
4595
4596	jc = jc + 1
4597	ENDDO
4598	ic = ic + 1
4599	ENDDO
4600
4601	END SUBROUTINE restore_surface_elements_usm_2d
4602
4603	END SUBROUTINE usm_read_restart_data
4604
4605
4606
4607	!------------------------------------------------------------------------------!
4608	!
4609	! Description:
4610	! ------------
4611	!> Soubroutine reads svf and svfsurf data from saved file
4612	!------------------------------------------------------------------------------!
4613	SUBROUTINE usm_read_svf_from_file
4614
4615	IMPLICIT NONE
4616	INTEGER(iwp) :: fsvf = 89
4617	INTEGER(iwp) :: i
4618	CHARACTER(usm_version_len) :: usm_version_field
4619	CHARACTER(svf_code_len) :: svf_code_field
4620
4621	DO i = 0, io_blocks-1
4622	IF ( i == io_group ) THEN
4623	OPEN ( fsvf, file=TRIM(svf_file_name)//TRIM(coupling_char)//myid_char, &
4624	form='unformatted', status='old' )
4625
4626	!-- read and check version
4627	READ ( fsvf ) usm_version_field
4628	IF ( TRIM(usm_version_field) /= TRIM(usm_version) ) THEN
4629	WRITE( message_string, * ) 'Version of binary SVF file "', &
4630	TRIM(usm_version_field), '" does not match ', &
4631	'the version of model "', TRIM(usm_version), '"'
4632	CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
4633	ENDIF
4634
4635	!-- read nsvfl, ncsfl
4636	READ ( fsvf ) nsvfl, ncsfl
4637	IF ( nsvfl <= 0 .OR. ncsfl < 0 ) THEN
4638	WRITE( message_string, * ) 'Wrong number of SVF or CSF'
4639	CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
4640	ELSE
4641	WRITE(message_string,*) ' Number of SVF and CSF to read', nsvfl, ncsfl
4642	CALL location_message( message_string, .TRUE. )
4643	ENDIF
4644
4645	ALLOCATE(svf(ndsvf,nsvfl))
4646	ALLOCATE(svfsurf(idsvf,nsvfl))
4647	READ(fsvf) svf
4648	READ(fsvf) svfsurf
4649	IF ( plant_canopy ) THEN
4650	ALLOCATE(csf(ndcsf,ncsfl))
4651	ALLOCATE(csfsurf(idcsf,ncsfl))
4652	READ(fsvf) csf
4653	READ(fsvf) csfsurf
4654	ENDIF
4655	READ ( fsvf ) svf_code_field
4656
4657	IF ( TRIM(svf_code_field) /= TRIM(svf_code) ) THEN
4658	WRITE( message_string, * ) 'Wrong structure of binary svf file'
4659	CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
4660	ENDIF
4661
4662	CLOSE (fsvf)
4663
4664	ENDIF
4665	#if defined( __parallel )
4666	CALL MPI_BARRIER( comm2d, ierr )
4667	#endif
4668	ENDDO
4669
4670	END SUBROUTINE usm_read_svf_from_file
4671
4672
4673	!------------------------------------------------------------------------------!
4674	! Description:
4675	! ------------
4676	!
4677	!> This subroutine reads walls, roofs and land categories and it parameters
4678	!> from input files.
4679	!------------------------------------------------------------------------------!
4680	SUBROUTINE usm_read_urban_surface_types
4681
4682	CHARACTER(12) :: wtn
4683	INTEGER(iwp) :: wtc
4684	REAL(wp), DIMENSION(n_surface_params) :: wtp
4685
4686	INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg) :: usm_par
4687	REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg) :: usm_val
4688	INTEGER(iwp) :: k, l, d, iw, jw, kw, it, ip, ii, ij, m
4689	INTEGER(iwp) :: i, j
4690	INTEGER(iwp) :: nz, roof, dirwe, dirsn
4691	INTEGER(iwp) :: category
4692	INTEGER(iwp) :: weheight1, wecat1, snheight1, sncat1
4693	INTEGER(iwp) :: weheight2, wecat2, snheight2, sncat2
4694	INTEGER(iwp) :: weheight3, wecat3, snheight3, sncat3
4695	REAL(wp) :: height, albedo, thick
4696	REAL(wp) :: wealbedo1, wethick1, snalbedo1, snthick1
4697	REAL(wp) :: wealbedo2, wethick2, snalbedo2, snthick2
4698	REAL(wp) :: wealbedo3, wethick3, snalbedo3, snthick3
4699
4700	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4701	!-- read categories of walls and their parameters
4702	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4703	DO ii = 0, io_blocks-1
4704	IF ( ii == io_group ) THEN
4705
4706	!-- open urban surface file
4707	OPEN( 151, file='SURFACE_PARAMETERS'//coupling_char, action='read', &
4708	status='old', form='formatted', err=15 )
4709	!-- first test and get n_surface_types
4710	k = 0
4711	l = 0
4712	DO
4713	l = l+1
4714	READ( 151, *, err=11, end=12 ) wtc, wtp, wtn
4715	k = k+1
4716	CYCLE
4717	11 CONTINUE
4718	ENDDO
4719	12 n_surface_types = k
4720	ALLOCATE( surface_type_names(n_surface_types) )
4721	ALLOCATE( surface_type_codes(n_surface_types) )
4722	ALLOCATE( surface_params(n_surface_params, n_surface_types) )
4723	!-- real reading
4724	rewind( 151 )
4725	k = 0
4726	DO
4727	READ( 151, *, err=13, end=14 ) wtc, wtp, wtn
4728	k = k+1
4729	surface_type_codes(k) = wtc
4730	surface_params(:,k) = wtp
4731	surface_type_names(k) = wtn
4732	CYCLE
4733	13 WRITE(6,'(i3,a,2i5)') myid, 'readparams2 error k=', k
4734	FLUSH(6)
4735	CONTINUE
4736	ENDDO
4737	14 CLOSE(151)
4738	CYCLE
4739	15 message_string = 'file SURFACE_PARAMETERS'//TRIM(coupling_char)//' does not exist'
4740	CALL message( 'usm_read_urban_surface_types', 'PA0513', 1, 2, 0, 6, 0 )
4741	ENDIF
4742	ENDDO
4743
4744	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4745	!-- read types of surfaces
4746	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4747	usm_par = 0
4748	DO ii = 0, io_blocks-1
4749	IF ( ii == io_group ) THEN
4750
4751	!
4752	!-- open csv urban surface file
4753	OPEN(