Home

Context Navigation

source: palm/trunk/SOURCE/urban_surface_mod.f90 @ 2233

Last change on this file since 2233 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 249.5 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 2233 2017-05-30 18:08:54Z suehring $
28	!
29	! 2232 2017-05-30 17:47:52Z suehring
30	! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
31	! insteat, heat fluxes are directly applied in diffusion_s.
32	!
33	! 2213 2017-04-24 15:10:35Z kanani
34	! Removal of output quantities usm_lad and usm_canopy_hr
35	!
36	! 2209 2017-04-19 09:34:46Z kanani
37	! cpp switch __mpi3 removed,
38	! minor formatting,
39	! small bugfix for division by zero (Krc)
40	!
41	! 2113 2017-01-12 13:40:46Z kanani
42	! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
43	!
44	! 2071 2016-11-17 11:22:14Z maronga
45	! Small bugfix (Resler)
46	!
47	! 2031 2016-10-21 15:11:58Z knoop
48	! renamed variable rho to rho_ocean
49	!
50	! 2024 2016-10-12 16:42:37Z kanani
51	! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
52	! optimization of MPI-RMA operations,
53	! declaration of pcbl as integer,
54	! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
55	! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
56	! use of new control parameter varnamelength,
57	! added output variables usm_rad_ressw, usm_rad_reslw,
58	! minor formatting changes,
59	! minor optimizations.
60	!
61	! 2011 2016-09-19 17:29:57Z kanani
62	! Major reformatting according to PALM coding standard (comments, blanks,
63	! alphabetical ordering, etc.),
64	! removed debug_prints,
65	! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
66	! defined in MODULE control_parameters (modules.f90) to avoid circular
67	! dependencies,
68	! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
69	!
70	! 2007 2016-08-24 15:47:17Z kanani
71	! Initial revision
72	!
73	!
74	! Description:
75	! ------------
76	! 2016/6/9 - Initial version of the USM (Urban Surface Model)
77	! authors: Jaroslav Resler, Pavel Krc
78	! (Czech Technical University in Prague and Institute of
79	! Computer Science of the Czech Academy of Sciences, Prague)
80	! with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
81	! partly inspired by PALM LSM (B. Maronga)
82	! parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
83	!> Module for Urban Surface Model (USM)
84	!> The module includes:
85	!> 1. radiation model with direct/diffuse radiation, shading, reflections
86	!> and integration with plant canopy
87	!> 2. wall and wall surface model
88	!> 3. surface layer energy balance
89	!> 4. anthropogenic heat (only from transportation so far)
90	!> 5. necessary auxiliary subroutines (reading inputs, writing outputs,
91	!> restart simulations, ...)
92	!> It also make use of standard radiation and integrates it into
93	!> urban surface model.
94	!>
95	!> Further work:
96	!> -------------
97	!> 1. Reduce number of shape view factors by merging factors for distant surfaces
98	!> under shallow angles. Idea: Iteratively select the smallest shape view
99	!> factor by value (among all sources and targets) which has a similarly
100	!> oriented source neighbor (or near enough) SVF and merge them by adding
101	!> value of the smaller SVF to the larger one and deleting the smaller one.
102	!> This will allow for better scaling at higher resolutions.
103	!>
104	!> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity
105	!> from surfaces that are visible from local surfaces (i.e. there is a SVF
106	!> where target is local). To do that, radiosity will be exchanged after each
107	!> reflection step using MPI_Alltoall instead of current MPI_Allgather.
108	!>
109	!> 3. Temporarily large values of surface heat flux can be observed, up to
110	!> 1.2 Km/s, which seem to be not realistic.
111	!>
112	!> @todo Check optimizations for RMA operations
113	!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
114	!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
115	!> factor 3 between min and max time
116	!------------------------------------------------------------------------------!
117	MODULE urban_surface_mod
118
119	USE arrays_3d, &
120	ONLY: zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
121
122	USE cloud_parameters, &
123	ONLY: cp, r_d
124
125	USE constants, &
126	ONLY: pi
127
128	USE control_parameters, &
129	ONLY: dz, topography, dt_3d, intermediate_timestep_count, &
130	initializing_actions, intermediate_timestep_count_max, &
131	simulated_time, end_time, timestep_scheme, tsc, &
132	coupling_char, io_blocks, io_group, message_string, &
133	time_since_reference_point, surface_pressure, &
134	g, pt_surface, large_scale_forcing, lsf_surf, &
135	time_do3d, dt_do3d, average_count_3d, varnamelength, &
136	urban_surface
137
138	USE cpulog, &
139	ONLY: cpu_log, log_point, log_point_s
140
141	USE grid_variables, &
142	ONLY: dx, dy, ddx, ddy, ddx2, ddy2
143
144	USE indices, &
145	ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, &
146	nysg, nzb, nzt, nbgp, wall_flags_0
147
148	USE, INTRINSIC :: iso_c_binding
149
150	USE kinds
151
152	USE pegrid
153
154	USE plant_canopy_model_mod, &
155	ONLY: plant_canopy, pch_index, &
156	pc_heating_rate, lad_s
157
158	USE radiation_model_mod, &
159	ONLY: radiation, calc_zenith, zenith, day_init, time_utc_init, &
160	rad_net, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, &
161	sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon, &
162	force_radiation_call
163
164	USE statistics, &
165	ONLY: hom, statistic_regions
166
167	USE surface_mod
168
169
170
171	IMPLICIT NONE
172
173	!-- configuration parameters (they can be setup in PALM config)
174	LOGICAL :: split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation
175	!< (.F. in case the radiation model already does it)
176	LOGICAL :: usm_energy_balance_land = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs
177	LOGICAL :: usm_energy_balance_wall = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs
178	LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used
179	LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used
180	LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls
181	LOGICAL :: mrt_factors = .FALSE. !< whether to generate MRT factor files during init
182	LOGICAL :: write_svf_on_init = .FALSE.
183	LOGICAL :: read_svf_on_init = .FALSE.
184	LOGICAL :: usm_lad_rma = .TRUE. !< use MPI RMA to access LAD for raytracing (instead of global array)
185
186	INTEGER(iwp) :: nrefsteps = 0 !< number of reflection steps to perform
187
188	INTEGER(iwp) :: land_category = 2 !< default category for land surface
189	INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone
190	INTEGER(iwp) :: pedestrant_category = 2 !< default category for wall surface in pedestrian zone
191	INTEGER(iwp) :: roof_category = 2 !< default category for root surface
192	REAL(wp) :: roof_height_limit = 4._wp !< height for distinguish between land surfaces and roofs
193
194	REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha)
195	REAL(wp) :: ra_horiz_coef = 5.0_wp !< mysterious coefficient for correction of overestimation
196	!< of r_a for horizontal surfaces -> TODO
197
198	!-- parameters of urban surface model
199	INTEGER(iwp), PARAMETER :: usm_version_len = 10 !< length of identification string of usm version
200	CHARACTER(usm_version_len), PARAMETER :: usm_version = 'USM v. 1.0' !< identification of version of binary svf and restart files
201	INTEGER(iwp), PARAMETER :: svf_code_len = 15 !< length of code for verification of the end of svf file
202	CHARACTER(svf_code_len), PARAMETER :: svf_code = '* end svf *' !< code for verification of the end of svf file
203	INTEGER(iwp) :: nzu !< number of layers of urban surface (will be calculated)
204	INTEGER(iwp) :: nzub,nzut !< bottom and top layer of urban surface (will be calculated)
205	INTEGER(iwp), PARAMETER :: nzut_free = 3 !< number of free layers in urban surface layer above top of buildings
206	INTEGER(iwp), PARAMETER :: ndsvf = 2 !< number of dimensions of real values in SVF
207	INTEGER(iwp), PARAMETER :: idsvf = 2 !< number of dimensions of integer values in SVF
208	INTEGER(iwp), PARAMETER :: ndcsf = 2 !< number of dimensions of real values in CSF
209	INTEGER(iwp), PARAMETER :: idcsf = 2 !< number of dimensions of integer values in CSF
210	INTEGER(iwp), PARAMETER :: kdcsf = 4 !< number of dimensions of integer values in CSF calculation array
211	INTEGER(iwp), PARAMETER :: id = 1 !< position of d-index in surfl and surf
212	INTEGER(iwp), PARAMETER :: iz = 2 !< position of k-index in surfl and surf
213	INTEGER(iwp), PARAMETER :: iy = 3 !< position of j-index in surfl and surf
214	INTEGER(iwp), PARAMETER :: ix = 4 !< position of i-index in surfl and surf
215	INTEGER(iwp), PARAMETER :: iroof = 0 !< 0 - index of ground or roof
216	INTEGER(iwp), PARAMETER :: isouth = 1 !< 1 - index of south facing wall
217	INTEGER(iwp), PARAMETER :: inorth = 2 !< 2 - index of north facing wall
218	INTEGER(iwp), PARAMETER :: iwest = 3 !< 3 - index of west facing wall
219	INTEGER(iwp), PARAMETER :: ieast = 4 !< 4 - index of east facing wall
220	INTEGER(iwp), PARAMETER :: isky = 5 !< 5 - index of top border of the urban surface layer ("urban sky")
221	INTEGER(iwp), PARAMETER :: inorthb = 6 !< 6 - index of free north border of the domain (south facing)
222	INTEGER(iwp), PARAMETER :: isouthb = 7 !< 7 - index of north south border of the domain (north facing)
223	INTEGER(iwp), PARAMETER :: ieastb = 8 !< 8 - index of east border of the domain (west facing)
224	INTEGER(iwp), PARAMETER :: iwestb = 9 !< 9 - index of wast border of the domain (east facing)
225	INTEGER(iwp), DIMENSION(0:9), PARAMETER :: idir = (/0,0,0,-1,1,0,0,0,-1,1/) !< surface normal direction x indices
226	INTEGER(iwp), DIMENSION(0:9), PARAMETER :: jdir = (/0,-1,1,0,0,0,-1,1,0,0/) !< surface normal direction y indices
227	INTEGER(iwp), DIMENSION(0:9), PARAMETER :: kdir = (/1,0,0,0,0,-1,0,0,0,0/) !< surface normal direction z indices
228	REAL(wp), DIMENSION(1:4) :: ddxy2 !< 1/dx^2 or 1/dy^2 (in surface normal direction)
229	INTEGER(iwp), DIMENSION(1:4,6:9) :: ijdb !< start and end of the local domain border coordinates (set in code)
230	LOGICAL, DIMENSION(6:9) :: isborder !< is PE on the border of the domain in four corresponding directions
231	!< parameter but set in the code
232
233	!-- indices and sizes of urban surface model
234	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surfl !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
235	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surf !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
236	INTEGER(iwp) :: nsurfl !< number of all surfaces in local processor
237	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nsurfs !< array of number of all surfaces in individual processors
238	INTEGER(iwp) :: startsky !< start index of block of sky
239	INTEGER(iwp) :: endsky !< end index of block of sky
240	INTEGER(iwp) :: nskys !< number of sky surfaces in local processor
241	INTEGER(iwp) :: startland !< start index of block of land and roof surfaces
242	INTEGER(iwp) :: endland !< end index of block of land and roof surfaces
243	INTEGER(iwp) :: nlands !< number of land and roof surfaces in local processor
244	INTEGER(iwp) :: startwall !< start index of block of wall surfaces
245	INTEGER(iwp) :: endwall !< end index of block of wall surfaces
246	INTEGER(iwp) :: nwalls !< number of wall surfaces in local processor
247	INTEGER(iwp) :: startenergy !< start index of block of real surfaces (land, walls and roofs)
248	INTEGER(iwp) :: endenergy !< end index of block of real surfaces (land, walls and roofs)
249	INTEGER(iwp) :: nenergy !< number of real surfaces in local processor
250	INTEGER(iwp) :: nsurf !< global number of surfaces in index array of surfaces (nsurf = Î£proc nsurfs)
251	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surfstart !< starts of blocks of surfaces for individual processors in array surf
252	!< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1]
253	INTEGER(iwp) :: nsvfl !< number of svf for local processor
254	INTEGER(iwp) :: ncsfl !< no. of csf in local processor
255	!< needed only during calc_svf but must be here because it is
256	!< shared between subroutines usm_calc_svf and usm_raytrace
257
258	!-- type for calculation of svf
259	TYPE t_svf
260	INTEGER(iwp) :: isurflt !<
261	INTEGER(iwp) :: isurfs !<
262	REAL(wp) :: rsvf !<
263	REAL(wp) :: rtransp !<
264	END TYPE
265
266	!-- type for calculation of csf
267	TYPE t_csf
268	INTEGER(iwp) :: ip !<
269	INTEGER(iwp) :: itx !<
270	INTEGER(iwp) :: ity !<
271	INTEGER(iwp) :: itz !<
272	INTEGER(iwp) :: isurfs !<
273	REAL(wp) :: rsvf !<
274	REAL(wp) :: rtransp !<
275	END TYPE
276	!
277	!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
278	TYPE t_surf_vertical
279	REAL(wp), DIMENSION(:), ALLOCATABLE :: t
280	END TYPE t_surf_vertical
281	!
282	!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
283	TYPE t_wall_vertical
284	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t
285	END TYPE t_wall_vertical
286
287	!-- arrays for calculation of svf and csf
288	TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array
289	TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array
290	TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: asvf1, asvf2 !< realizations of svf array
291	TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET :: acsf1, acsf2 !< realizations of csf array
292	INTEGER(iwp) :: nsvfla !< dimmension of array allocated for storage of svf in local processor
293	INTEGER(iwp) :: ncsfla !< dimmension of array allocated for storage of csf in local processor
294	INTEGER(iwp) :: msvf, mcsf !< mod for swapping the growing array
295	INTEGER(iwp), PARAMETER :: gasize = 10000 !< initial size of growing arrays
296	!-- temporary arrays for calculation of csf in raytracing
297	INTEGER(iwp) :: maxboxesg !< max number of boxes ray can cross in the domain
298	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: boxes !< coordinates of gridboxes being crossed by ray
299	REAL(wp), DIMENSION(:), ALLOCATABLE :: crlens !< array of crossing lengths of ray for particular grid boxes
300	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored
301	INTEGER(kind=MPI_ADDRESS_KIND), &
302	DIMENSION(:), ALLOCATABLE :: lad_disp !< array of displaycements of lad in local array of proc lad_ip
303	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad_s_ray !< array of received lad_s for appropriate gridboxes crossed by ray
304
305	!-- arrays storing the values of USM
306	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of source and target surface for svf[isvf]
307	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: svf !< array of shape view factors+direct irradiation factors for local surfaces
308	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins !< array of sw radiation falling to local surface after i-th reflection
309	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl !< array of lw radiation for local surface after i-th reflection
310
311	!< Inward radiation is also valid for virtual surfaces (radiation leaving domain)
312	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< array of sw radiation falling to local surface including radiation from reflections
313	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< array of lw radiation falling to local surface including radiation from reflections
314	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir !< array of direct sw radiation falling to local surface
315	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif !< array of diffuse sw radiation from sky and model boundary falling to local surface
316	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif !< array of diffuse lw radiation from sky and model boundary falling to local surface
317
318	!< Outward radiation is only valid for nonvirtual surfaces
319	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< array of reflected sw radiation for local surface in i-th reflection
320	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< array of reflected + emitted lw radiation for local surface in i-th reflection
321	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection
322	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
323	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
324	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
325	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf !< array of total radiation flux incoming to minus outgoing from local surface
326	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< local copy of rad_net (net radiation at surface)
327
328	!-- arrays for time averages
329	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l
330	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
331	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
332	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
333	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
334	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
335	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
336	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
337	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
338	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
339	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
340	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
341	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface
342
343	!-- block variables needed for calculation of the plant canopy model inside the urban surface model
344	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: csfsurf !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
345	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csf !< array of plant canopy sink fators + direct irradiation factors (transparency)
346	!< for local surfaces
347	INTEGER(wp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
348	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< index of local pcb[k,j,i]
349	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw !< array of absorbed sw radiation for local plant canopy box
350	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw !< array of absorbed lw radiation for local plant canopy box
351	INTEGER(iwp) :: npcbl !< number of the plant canopy gridboxes in local processor
352	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy
353	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy
354	REAL(wp), DIMENSION(:,:,:), POINTER :: usm_lad !< subset of lad_s within urban surface, transformed to plain Z coordinate
355	REAL(wp), DIMENSION(:), POINTER :: usm_lad_g !< usm_lad globalized (used to avoid MPI RMA calls in raytracing)
356	REAL(wp) :: prototype_lad !< prototype leaf area density for computing effective optical depth
357	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing
358
359	!-- radiation related arrays (it should be better in interface of radiation module of PALM
360	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_dir !< direct sw radiation
361	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_diff !< diffusion sw radiation
362	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_diff !< diffusion lw radiation
363
364	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
365	!-- anthropogenic heat sources
366	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
367	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2)
368	REAL(wp), DIMENSION(:), ALLOCATABLE :: aheatprof !< diurnal profile of anthropogenic heat
369
370	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
371	!-- wall surface model
372	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
373	!-- wall surface model constants
374	INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched)
375	INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched)
376	INTEGER(iwp), PARAMETER :: nzw = 4 !< number of wall layers (fixed for now)
377
378	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
379	!< normalized soil, wall and roof layer depths (m/m)
380
381	REAL(wp) :: wall_inner_temperature = 296.0_wp !< temperature of the inner wall surface (~23 degrees C) (K)
382	REAL(wp) :: roof_inner_temperature = 296.0_wp !< temperature of the inner roof surface (~23 degrees C) (K)
383	REAL(wp) :: soil_inner_temperature = 283.0_wp !< temperature of the deep soil (~10 degrees C) (K)
384
385	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
386	!-- surface and material model variables for walls, ground, roofs
387	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
388	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m)
389
390	#if defined( __nopointer )
391	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h !< wall surface temperature (K) at horizontal walls
392	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_p !< progn. wall surface temperature (K) at horizontal walls
393
394	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v
395	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_p
396	#else
397	REAL(wp), DIMENSION(:), POINTER :: t_surf_h
398	REAL(wp), DIMENSION(:), POINTER :: t_surf_h_p
399
400	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_1
401	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_2
402
403	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v
404	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v_p
405
406	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_1
407	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_2
408	#endif
409	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_av !< average of wall surface temperature (K)
410
411	!-- Temporal tendencies for time stepping
412	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K)
413
414	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415	!-- Energy balance variables
416	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417	!-- parameters of the land, roof and wall surfaces
418	REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface
419	!-- parameters of the wall surfaces
420	REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface
421
422	#if defined( __nopointer )
423	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_wall_h !< Wall temperature (K)
424	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_wall_h_av !< Average of t_wall
425	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_wall_h_p !< Prog. wall temperature (K)
426
427	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v !< Wall temperature (K)
428	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av !< Average of t_wall
429	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_p !< Prog. wall temperature (K)
430	#else
431	REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p
432	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av, t_wall_h_1, t_wall_h_2
433
434	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p
435	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av, t_wall_v_1, t_wall_v_2
436	#endif
437
438	!-- Wall temporal tendencies for time stepping
439	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array
440
441	!-- Surface and material parameters classes (surface_type)
442	!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
443	INTEGER(iwp) :: n_surface_types !< number of the wall type categories
444	INTEGER(iwp), PARAMETER :: n_surface_params = 8 !< number of parameters for each type of the wall
445	INTEGER(iwp), PARAMETER :: ialbedo = 1 !< albedo of the surface
446	INTEGER(iwp), PARAMETER :: iemiss = 2 !< emissivity of the surface
447	INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 )
448	INTEGER(iwp), PARAMETER :: irough = 4 !< roughness relative to concrete
449	INTEGER(iwp), PARAMETER :: icsurf = 5 !< Surface skin layer heat capacity (J mâ2 Kâ1 )
450	INTEGER(iwp), PARAMETER :: ithick = 6 !< thickness of the surface (wall, roof, land) ( m )
451	INTEGER(iwp), PARAMETER :: irhoC = 7 !< volumetric heat capacity rho*C of the material ( J mâ3 Kâ1 )
452	INTEGER(iwp), PARAMETER :: ilambdah = 8 !< thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
453	CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports)
454	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types
455	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: surface_params !< parameters of wall types
456
457	CHARACTER(len=*), PARAMETER :: svf_file_name='usm_svf'
458
459	!-- interfaces of subroutines accessed from outside of this module
460	INTERFACE usm_check_data_output
461	MODULE PROCEDURE usm_check_data_output
462	END INTERFACE usm_check_data_output
463
464	INTERFACE usm_check_parameters
465	MODULE PROCEDURE usm_check_parameters
466	END INTERFACE usm_check_parameters
467
468	INTERFACE usm_data_output_3d
469	MODULE PROCEDURE usm_data_output_3d
470	END INTERFACE usm_data_output_3d
471
472	INTERFACE usm_define_netcdf_grid
473	MODULE PROCEDURE usm_define_netcdf_grid
474	END INTERFACE usm_define_netcdf_grid
475
476	INTERFACE usm_init_urban_surface
477	MODULE PROCEDURE usm_init_urban_surface
478	END INTERFACE usm_init_urban_surface
479
480	INTERFACE usm_material_heat_model
481	MODULE PROCEDURE usm_material_heat_model
482	END INTERFACE usm_material_heat_model
483
484	INTERFACE usm_parin
485	MODULE PROCEDURE usm_parin
486	END INTERFACE usm_parin
487
488	INTERFACE usm_radiation
489	MODULE PROCEDURE usm_radiation
490	END INTERFACE usm_radiation
491
492	INTERFACE usm_read_restart_data
493	MODULE PROCEDURE usm_read_restart_data
494	END INTERFACE usm_read_restart_data
495
496	INTERFACE usm_surface_energy_balance
497	MODULE PROCEDURE usm_surface_energy_balance
498	END INTERFACE usm_surface_energy_balance
499
500	INTERFACE usm_swap_timelevel
501	MODULE PROCEDURE usm_swap_timelevel
502	END INTERFACE usm_swap_timelevel
503
504	INTERFACE usm_write_restart_data
505	MODULE PROCEDURE usm_write_restart_data
506	END INTERFACE usm_write_restart_data
507
508	SAVE
509
510	PRIVATE
511
512	!-- Public parameters, constants and initial values
513	PUBLIC split_diffusion_radiation, &
514	usm_anthropogenic_heat, usm_material_model, mrt_factors, &
515	usm_check_parameters, &
516	usm_energy_balance_land, usm_energy_balance_wall, nrefsteps, &
517	usm_init_urban_surface, usm_radiation, usm_read_restart_data, &
518	usm_surface_energy_balance, usm_material_heat_model, &
519	usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, &
520	usm_data_output_3d, usm_define_netcdf_grid, usm_parin, &
521	usm_write_restart_data, &
522	nzub, nzut, ra_horiz_coef, usm_lad_rma, &
523	land_category, pedestrant_category, wall_category, roof_category, &
524	write_svf_on_init, read_svf_on_init
525
526
527	CONTAINS
528
529
530	!------------------------------------------------------------------------------!
531	! Description:
532	! ------------
533	!> This subroutine creates the necessary indices of the urban surfaces
534	!> and plant canopy and it allocates the needed arrays for USM
535	!------------------------------------------------------------------------------!
536	SUBROUTINE usm_allocate_urban_surface
537
538	IMPLICIT NONE
539
540	INTEGER(iwp) :: i, j, k, d, l, ir, jr, ids, m
541	INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i)
542	INTEGER(iwp) :: k_topo2 !< vertical index indicating topography top for given (j,i)
543	INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb
544	INTEGER(iwp) :: procid
545
546
547
548
549	!-- auxiliary vars
550	ddxy2 = (/ddy2,ddy2,ddx2,ddx2/) !< 1/dx^2 or 1/dy^2 (in surface normal direction)
551
552	CALL location_message( '', .TRUE. )
553	CALL location_message( ' allocation of needed arrays', .TRUE. )
554	!
555	!-- Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
556	!-- removed later). The following contruct finds the lowest / largest index
557	!-- for any upward-facing wall (see bit 12).
558	nzubl = MINVAL( &
559	MAXLOC( &
560	MERGE( 1, 0, &
561	BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) &
562	), DIM = 1 &
563	) - 1 &
564	)
565	nzutl = MAXVAL( &
566	MAXLOC( &
567	MERGE( 1, 0, &
568	BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) &
569	), DIM = 1 &
570	) - 1 &
571	)
572	nzubl = max(nzubl,nzb)
573
574
575	IF ( plant_canopy ) THEN
576	!-- allocate needed arrays
577	ALLOCATE( pct(nys:nyn,nxl:nxr) )
578	ALLOCATE( pch(nys:nyn,nxl:nxr) )
579
580	!-- calculate plant canopy height
581	npcbl = 0
582	pct = 0.0_wp
583	pch = 0.0_wp
584	DO i = nxl, nxr
585	DO j = nys, nyn
586	!
587	!-- Find topography top index
588	k_topo = MAXLOC( MERGE( 1, 0, &
589	BTEST( wall_flags_0(:,j,i), 12 ) &
590	), DIM = 1 &
591	) - 1
592	DO k = nzt+1, 0, -1
593	IF ( lad_s(k,j,i) /= 0.0_wp ) THEN
594	!-- we are at the top of the pcs
595	pct(j,i) = k + k_topo
596	pch(j,i) = k
597	npcbl = npcbl + pch(j,i)
598	EXIT
599	ENDIF
600	ENDDO
601	ENDDO
602	ENDDO
603
604	nzutl = max(nzutl, maxval(pct))
605	!-- code of plant canopy model uses parameter pch_index
606	!-- we need to setup it here to right value
607	!-- (pch_index, lad_s and other arrays in PCM are defined flat)
608	pch_index = maxval(pch)
609
610	prototype_lad = maxval(lad_s) * .9_wp !< better be *1.0 if lad is either 0 or maxval(lad) everywhere
611	IF ( prototype_lad <= 0._wp ) prototype_lad = .3_wp
612	!WRITE(message_string, '(a,f6.3)') 'Precomputing effective box optical ' &
613	! // 'depth using prototype leaf area density = ', prototype_lad
614	!CALL message('usm_init_urban_surface', 'PA0520', 0, 0, -1, 6, 0)
615	ENDIF
616
617	nzutl = min(nzutl+nzut_free, nzt)
618
619	#if defined( __parallel )
620	CALL MPI_AllReduce(nzubl,nzub,1,MPI_INTEGER,MPI_MIN,comm2d,ierr);
621	CALL MPI_AllReduce(nzutl,nzut,1,MPI_INTEGER,MPI_MAX,comm2d,ierr);
622	#else
623	nzub = nzubl
624	nzut = nzutl
625	#endif
626
627	!-- global number of urban layers
628	nzu = nzut - nzub + 1
629
630	!-- allocate urban surfaces grid
631	!-- calc number of surfaces in local proc
632	CALL location_message( ' calculation of indices for surfaces', .TRUE. )
633	nsurfl = 0
634	!
635	!-- Number of land- and roof surfaces. Note, since horizontal surface elements
636	!-- are already counted in surface_mod, in case be simply reused here.
637	startland = 1
638	nsurfl = surf_usm_h%ns
639	endland = nsurfl
640	nlands = endland-startland+1
641
642	!
643	!-- Number of vertical surfaces. As vertical surfaces are already
644	!-- counted in surface mod, it can be reused here.
645	startwall = nsurfl+1
646	nsurfl = nsurfl + surf_usm_v(0)%ns + surf_usm_v(1)%ns + &
647	surf_usm_v(2)%ns + surf_usm_v(3)%ns
648	endwall = nsurfl
649	nwalls = endwall-startwall+1
650
651
652	!-- range of energy balance surfaces ! will be treated separately by surf_usm_h and surf_usm_v
653	nenergy = 0
654	IF ( usm_energy_balance_land ) THEN
655	startenergy = startland
656	nenergy = nenergy + nlands
657	ELSE
658	startenergy = startwall
659	ENDIF
660	IF ( usm_energy_balance_wall ) THEN
661	endenergy = endwall
662	nenergy = nenergy + nwalls
663	ELSE
664	endenergy = endland
665	ENDIF
666
667	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
668	!-- block of virtual surfaces
669	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
670	!-- calculate sky surfaces ! not used so far!
671	startsky = nsurfl+1
672	nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1)
673	endsky = nsurfl
674	nskys = endsky-startsky+1
675
676	!-- border flags
677	#if defined( __parallel )
678	isborder = (/ north_border_pe, south_border_pe, right_border_pe, left_border_pe /)
679	#else
680	isborder = (/.TRUE.,.TRUE.,.TRUE.,.TRUE./)
681	#endif
682	!-- fill array of the limits of the local domain borders
683	ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) )
684	!-- calulation of the free borders of the domain
685	DO ids = 6,9
686	IF ( isborder(ids) ) THEN
687	!-- free border of the domain in direction ids
688	DO i = ijdb(1,ids), ijdb(2,ids)
689	DO j = ijdb(3,ids), ijdb(4,ids)
690
691	k_topo = MAXLOC( MERGE( 1, 0, &
692	BTEST( wall_flags_0(:,j,i), 12 ) &
693	), DIM = 1 &
694	) - 1
695	k_topo2 = MAXLOC( MERGE( 1, 0, &
696	BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 ) &
697	), DIM = 1 &
698	) - 1
699
700	k = nzut - MAX( k_topo, k_topo2 )
701	nsurfl = nsurfl + k
702	ENDDO
703	ENDDO
704	ENDIF
705	ENDDO
706
707	!-- fill gridpcbl and pcbl
708	IF ( plant_canopy ) THEN
709	ALLOCATE( pcbl(iz:ix, 1:npcbl) )
710	ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) )
711	gridpcbl(:,:,:) = 0
712	ipcgb = 0
713	DO i = nxl, nxr
714	DO j = nys, nyn
715	!
716	!-- Find topography top index
717	k_topo = MAXLOC( MERGE( 1, 0, &
718	BTEST( wall_flags_0(:,j,i), 12 ) &
719	), DIM = 1 &
720	) - 1
721	DO k = k_topo + 1, pct(j,i)
722	ipcgb = ipcgb + 1
723	gridpcbl(k,j,i) = ipcgb
724	pcbl(:,ipcgb) = (/ k, j, i /)
725	ENDDO
726	ENDDO
727	ENDDO
728
729	ALLOCATE( pcbinsw( 1:npcbl ) )
730	ALLOCATE( pcbinlw( 1:npcbl ) )
731	ENDIF
732
733	!-- fill surfl
734	ALLOCATE(surfl(5,nsurfl))
735	isurf = 0
736
737	!-- add land surfaces or roofs
738	DO i = nxl, nxr
739	DO j = nys, nyn
740	DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
741	k = surf_usm_h%k(m)
742
743	isurf = isurf + 1
744	surfl(:,isurf) = (/iroof,k,j,i,m/)
745	ENDDO
746	ENDDO
747	ENDDO
748
749	!-- add walls
750	DO i = nxl, nxr
751	DO j = nys, nyn
752	l = 0
753	DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
754	k = surf_usm_v(l)%k(m)
755
756	isurf = isurf + 1
757	surfl(:,isurf) = (/2,k,j,i,m/)
758	ENDDO
759	l = 1
760	DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
761	k = surf_usm_v(l)%k(m)
762
763	isurf = isurf + 1
764	surfl(:,isurf) = (/1,k,j,i,m/)
765	ENDDO
766	l = 2
767	DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
768	k = surf_usm_v(l)%k(m)
769
770	isurf = isurf + 1
771	surfl(:,isurf) = (/4,k,j,i,m/)
772	ENDDO
773	l = 3
774	DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
775	k = surf_usm_v(l)%k(m)
776
777	isurf = isurf + 1
778	surfl(:,isurf) = (/3,k,j,i,m/)
779	ENDDO
780	ENDDO
781	ENDDO
782
783	!-- add sky
784	DO i = nxl, nxr
785	DO j = nys, nyn
786	isurf = isurf + 1
787	k = nzut
788	surfl(:,isurf) = (/isky,k,j,i,-1/)
789	ENDDO
790	ENDDO
791
792	!-- calulation of the free borders of the domain
793	DO ids = 6,9
794	IF ( isborder(ids) ) THEN
795	!-- free border of the domain in direction ids
796	DO i = ijdb(1,ids), ijdb(2,ids)
797	DO j = ijdb(3,ids), ijdb(4,ids)
798	k_topo = MAXLOC( MERGE( 1, 0, &
799	BTEST( wall_flags_0(:,j,i), 12 ) &
800	), DIM = 1 &
801	) - 1
802	k_topo2 = MAXLOC( MERGE( 1, 0, &
803	BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 ) &
804	), DIM = 1 &
805	) - 1
806	DO k = MAX(k_topo,k_topo2)+1, nzut
807	isurf = isurf + 1
808	surfl(:,isurf) = (/ids,k,j,i,-1/)
809	ENDDO
810	ENDDO
811	ENDDO
812	ENDIF
813	ENDDO
814
815	!-- global array surf of indices of surfaces and displacement index array surfstart
816	ALLOCATE(nsurfs(0:numprocs-1))
817
818	#if defined( __parallel )
819	CALL MPI_Allgather(nsurfl,1,MPI_INTEGER,nsurfs,1,MPI_INTEGER,comm2d,ierr)
820	#else
821	nsurfs(0) = nsurfl
822	#endif
823	ALLOCATE(surfstart(0:numprocs))
824	k = 0
825	DO i=0,numprocs-1
826	surfstart(i) = k
827	k = k+nsurfs(i)
828	ENDDO
829	surfstart(numprocs) = k
830	nsurf = k
831	ALLOCATE(surf(5,nsurf))
832
833	#if defined( __parallel )
834	CALL MPI_AllGatherv(surfl, nsurfl5, MPI_INTEGER, surf, nsurfs5, surfstart*5, MPI_INTEGER, comm2d, ierr)
835	#else
836	surf = surfl
837	#endif
838
839	!--
840	!-- allocation of the arrays for direct and diffusion radiation
841	CALL location_message( ' allocation of radiation arrays', .TRUE. )
842	!-- rad_sw_in, rad_lw_in are computed in radiation model,
843	!-- splitting of direct and diffusion part is done
844	!-- in usm_calc_diffusion_radiation for now
845	ALLOCATE( rad_sw_in_dir(nysg:nyng,nxlg:nxrg) )
846	ALLOCATE( rad_sw_in_diff(nysg:nyng,nxlg:nxrg) )
847	ALLOCATE( rad_lw_in_diff(nysg:nyng,nxlg:nxrg) )
848
849	!-- allocate radiation arrays
850	ALLOCATE( surfins(nsurfl) )
851	ALLOCATE( surfinl(nsurfl) )
852	ALLOCATE( surfinsw(nsurfl) )
853	ALLOCATE( surfinlw(nsurfl) )
854	ALLOCATE( surfinswdir(nsurfl) )
855	ALLOCATE( surfinswdif(nsurfl) )
856	ALLOCATE( surfinlwdif(nsurfl) )
857	ALLOCATE( surfoutsl(startenergy:endenergy) )
858	ALLOCATE( surfoutll(startenergy:endenergy) )
859	ALLOCATE( surfoutsw(startenergy:endenergy) )
860	ALLOCATE( surfoutlw(startenergy:endenergy) )
861	ALLOCATE( surfouts(nsurf) ) !TODO: global surfaces without virtual
862	ALLOCATE( surfoutl(nsurf) ) !TODO: global surfaces without virtual
863
864
865
866	!
867	!-- Allocate radiation arrays which are part of the new data type.
868	!-- For horizontal surfaces.
869	ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns) )
870	ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
871	!
872	!-- New
873	ALLOCATE( surf_usm_h%rad_in_sw(1:surf_usm_h%ns) )
874	ALLOCATE( surf_usm_h%rad_out_sw(1:surf_usm_h%ns) )
875	ALLOCATE( surf_usm_h%rad_in_lw(1:surf_usm_h%ns) )
876	ALLOCATE( surf_usm_h%rad_out_lw(1:surf_usm_h%ns) )
877	!
878	!-- For vertical surfaces
879	DO l = 0, 3
880	ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) )
881	ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
882	ALLOCATE( surf_usm_v(l)%rad_in_sw(1:surf_usm_v(l)%ns) )
883	ALLOCATE( surf_usm_v(l)%rad_out_sw(1:surf_usm_v(l)%ns) )
884	ALLOCATE( surf_usm_v(l)%rad_in_lw(1:surf_usm_v(l)%ns) )
885	ALLOCATE( surf_usm_v(l)%rad_out_lw(1:surf_usm_v(l)%ns) )
886	ENDDO
887
888	!-- Wall surface model
889	!-- allocate arrays for wall surface model and define pointers
890
891	!-- allocate array of wall types and wall parameters
892	ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
893	DO l = 0, 3
894	ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
895	ENDDO
896
897	!-- broadband albedo of the land, roof and wall surface
898	!-- for domain border and sky set artifically to 1.0
899	!-- what allows us to calculate heat flux leaving over
900	!-- side and top borders of the domain
901	ALLOCATE ( albedo_surf(nsurfl) )
902	albedo_surf = 1.0_wp
903	ALLOCATE ( surf_usm_h%albedo_surf(1:surf_usm_h%ns) )
904	DO l = 0, 3
905	ALLOCATE( surf_usm_v(l)%albedo_surf(1:surf_usm_v(l)%ns) )
906	ENDDO
907
908	!-- wall and roof surface parameters. First for horizontal surfaces
909	ALLOCATE ( emiss_surf(startenergy:endenergy) )
910
911	ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) )
912	ALLOCATE ( surf_usm_h%emiss_surf(1:surf_usm_h%ns) )
913	ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) )
914	ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) )
915	ALLOCATE ( surf_usm_h%roughness_wall(1:surf_usm_h%ns) )
916	!
917	!-- For vertical surfaces.
918	DO l = 0, 3
919	ALLOCATE ( surf_usm_v(l)%emiss_surf(1:surf_usm_v(l)%ns) )
920	ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) )
921	ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) )
922	ALLOCATE ( surf_usm_v(l)%roughness_wall(1:surf_usm_v(l)%ns) )
923	ENDDO
924
925	!-- allocate wall and roof material parameters. First for horizontal surfaces
926	ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) )
927	ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
928	ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
929	!
930	!-- For vertical surfaces.
931	DO l = 0, 3
932	ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) )
933	ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
934	ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
935	ENDDO
936
937	!-- allocate wall and roof layers sizes. For horizontal surfaces.
938	ALLOCATE ( zwn(nzb_wall:nzt_wall) )
939	ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
940	ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
941	ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
942	ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
943	ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
944	!
945	!-- For vertical surfaces.
946	DO l = 0, 3
947	ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
948	ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
949	ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
950	ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
951	ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
952	ENDDO
953
954	!-- allocate wall and roof temperature arrays, for horizontal walls
955	#if defined( __nopointer )
956	ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
957	ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
958	ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
959	ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
960
961	ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
962	ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
963	ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
964	ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
965	#else
966	ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) )
967	ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) )
968	ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
969	ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
970
971	!-- initial assignment of the pointers
972	t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2
973	t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2
974	#endif
975
976	!-- allocate wall and roof temperature arrays, for vertical walls
977	#if defined( __nopointer )
978	DO l = 0, 3
979	ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) )
980	ALLOCATE ( t_surf_v(l)%t_p(1:surf_usm_v(l)%ns) )
981	ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
982	ALLOCATE ( t_wall_v(l)%t_p(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
983	ENDDO
984	#else
985	DO l = 0, 3
986	ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) )
987	ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) )
988	ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
989	ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
990	ENDDO
991	!
992	!-- initial assignment of the pointers
993	t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2
994	t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
995	#endif
996	!
997	!-- Allocate intermediate timestep arrays. For horizontal surfaces.
998	ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns) )
999	ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1000	!
1001	!-- Set inital values for prognostic quantities
1002	IF ( ALLOCATED( surf_usm_h%tt_surface_m ) ) surf_usm_h%tt_surface_m = 0.0_wp
1003	IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp
1004	!
1005	!-- Now, for vertical surfaces
1006	DO l = 0, 3
1007	ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns) )
1008	ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1009	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) ) surf_usm_v(l)%tt_surface_m = 0.0_wp
1010	IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp
1011	ENDDO
1012
1013	!-- allocate wall heat flux output array and set initial values. For horizontal surfaces
1014	! ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns) ) !can be removed
1015	ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
1016	ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
1017	IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp
1018	IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp
1019	IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp
1020	!
1021	!-- Now, for vertical surfaces
1022	DO l = 0, 3
1023	! ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns) ) ! can be removed
1024	ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
1025	ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
1026	IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp
1027	IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp
1028	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp
1029	ENDDO
1030
1031	END SUBROUTINE usm_allocate_urban_surface
1032
1033
1034
1035	!------------------------------------------------------------------------------!
1036	! Description:
1037	! ------------
1038	!> Sum up and time-average urban surface output quantities as well as allocate
1039	!> the array necessary for storing the average.
1040	!------------------------------------------------------------------------------!
1041	SUBROUTINE usm_average_3d_data( mode, variable )
1042
1043	IMPLICIT NONE
1044
1045	CHARACTER (len=*), INTENT(IN) :: mode
1046	CHARACTER (len=*), INTENT(IN) :: variable
1047
1048	INTEGER(iwp) :: i, j, k, l, m, ids, iwl,istat
1049	CHARACTER (len=varnamelength) :: var, surfid
1050	INTEGER(iwp), PARAMETER :: nd = 5
1051	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
1052
1053	!-- find the real name of the variable
1054	var = TRIM(variable)
1055	DO i = 0, nd-1
1056	k = len(TRIM(var))
1057	j = len(TRIM(dirname(i)))
1058	IF ( var(k-j+1:k) == dirname(i) ) THEN
1059	ids = i
1060	var = var(:k-j)
1061	EXIT
1062	ENDIF
1063	ENDDO
1064	IF ( ids == -1 ) THEN
1065	var = TRIM(variable)
1066	ENDIF
1067	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
1068	!-- wall layers
1069	READ(var(12:12), '(I1)', iostat=istat ) iwl
1070	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1071	var = var(1:10)
1072	ELSE
1073	!-- wrong wall layer index
1074	RETURN
1075	ENDIF
1076	ENDIF
1077
1078	IF ( mode == 'allocate' ) THEN
1079
1080	SELECT CASE ( TRIM( var ) )
1081
1082	CASE ( 'usm_rad_net' )
1083	!-- array of complete radiation balance
1084	IF ( .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN
1085	ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
1086	surf_usm_h%rad_net_av = 0.0_wp
1087	ENDIF
1088	DO l = 0, 3
1089	IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN
1090	ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
1091	surf_usm_v(l)%rad_net_av = 0.0_wp
1092	ENDIF
1093	ENDDO
1094
1095	CASE ( 'usm_rad_insw' )
1096	!-- array of sw radiation falling to surface after i-th reflection
1097	IF ( .NOT. ALLOCATED(surf_usm_h%surfinsw_av) ) THEN
1098	ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
1099	surf_usm_h%surfinsw_av = 0.0_wp
1100	ENDIF
1101	DO l = 0, 3
1102	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinsw_av) ) THEN
1103	ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
1104	surf_usm_v(l)%surfinsw_av = 0.0_wp
1105	ENDIF
1106	ENDDO
1107
1108	CASE ( 'usm_rad_inlw' )
1109	!-- array of lw radiation falling to surface after i-th reflection
1110	IF ( .NOT. ALLOCATED(surf_usm_h%surfinlw_av) ) THEN
1111	ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
1112	surf_usm_h%surfinlw_av = 0.0_wp
1113	ENDIF
1114	DO l = 0, 3
1115	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinlw_av) ) THEN
1116	ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
1117	surf_usm_v(l)%surfinlw_av = 0.0_wp
1118	ENDIF
1119	ENDDO
1120
1121	CASE ( 'usm_rad_inswdir' )
1122	!-- array of direct sw radiation falling to surface from sun
1123	IF ( .NOT. ALLOCATED(surfinswdir_av) ) THEN
1124	ALLOCATE( surfinswdir_av(startenergy:endenergy) )
1125	surfinswdir_av = 0.0_wp
1126	ENDIF
1127
1128	CASE ( 'usm_rad_inswdif' )
1129	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1130	IF ( .NOT. ALLOCATED(surfinswdif_av) ) THEN
1131	ALLOCATE( surfinswdif_av(startenergy:endenergy) )
1132	surfinswdif_av = 0.0_wp
1133	ENDIF
1134
1135	CASE ( 'usm_rad_inswref' )
1136	!-- array of sw radiation falling to surface from reflections
1137	IF ( .NOT. ALLOCATED(surfinswref_av) ) THEN
1138	ALLOCATE( surfinswref_av(startenergy:endenergy) )
1139	surfinswref_av = 0.0_wp
1140	ENDIF
1141
1142	CASE ( 'usm_rad_inlwdif' )
1143	!-- array of sw radiation falling to surface after i-th reflection
1144	IF ( .NOT. ALLOCATED(surfinlwdif_av) ) THEN
1145	ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
1146	surfinlwdif_av = 0.0_wp
1147	ENDIF
1148
1149	CASE ( 'usm_rad_inlwref' )
1150	!-- array of lw radiation falling to surface from reflections
1151	IF ( .NOT. ALLOCATED(surfinlwref_av) ) THEN
1152	ALLOCATE( surfinlwref_av(startenergy:endenergy) )
1153	surfinlwref_av = 0.0_wp
1154	ENDIF
1155
1156	CASE ( 'usm_rad_outsw' )
1157	!-- array of sw radiation emitted from surface after i-th reflection
1158	IF ( .NOT. ALLOCATED(surfoutsw_av) ) THEN
1159	ALLOCATE( surfoutsw_av(startenergy:endenergy) )
1160	surfoutsw_av = 0.0_wp
1161	ENDIF
1162
1163	CASE ( 'usm_rad_outlw' )
1164	!-- array of lw radiation emitted from surface after i-th reflection
1165	IF ( .NOT. ALLOCATED(surfoutlw_av) ) THEN
1166	ALLOCATE( surfoutlw_av(startenergy:endenergy) )
1167	surfoutlw_av = 0.0_wp
1168	ENDIF
1169	CASE ( 'usm_rad_ressw' )
1170	!-- array of residua of sw radiation absorbed in surface after last reflection
1171	IF ( .NOT. ALLOCATED(surfins_av) ) THEN
1172	ALLOCATE( surfins_av(startenergy:endenergy) )
1173	surfins_av = 0.0_wp
1174	ENDIF
1175
1176	CASE ( 'usm_rad_reslw' )
1177	!-- array of residua of lw radiation absorbed in surface after last reflection
1178	IF ( .NOT. ALLOCATED(surfinl_av) ) THEN
1179	ALLOCATE( surfinl_av(startenergy:endenergy) )
1180	surfinl_av = 0.0_wp
1181	ENDIF
1182
1183	CASE ( 'usm_rad_hf' )
1184	!-- array of heat flux from radiation for surfaces after i-th reflection
1185	IF ( .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN
1186	ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
1187	surf_usm_h%surfhf_av = 0.0_wp
1188	ENDIF
1189	DO l = 0, 3
1190	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfhf_av) ) THEN
1191	ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
1192	surf_usm_v(l)%surfhf_av = 0.0_wp
1193	ENDIF
1194	ENDDO
1195
1196	CASE ( 'usm_wshf' )
1197	!-- array of sensible heat flux from surfaces
1198	!-- land surfaces
1199	IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN
1200	ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1201	surf_usm_h%wshf_eb_av = 0.0_wp
1202	ENDIF
1203	DO l = 0, 3
1204	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN
1205	ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1206	surf_usm_v(l)%wshf_eb_av = 0.0_wp
1207	ENDIF
1208	ENDDO
1209
1210	CASE ( 'usm_wghf' )
1211	!-- array of heat flux from ground (wall, roof, land)
1212	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN
1213	ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1214	surf_usm_h%wghf_eb_av = 0.0_wp
1215	ENDIF
1216	DO l = 0, 3
1217	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN
1218	ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1219	surf_usm_v(l)%wghf_eb_av = 0.0_wp
1220	ENDIF
1221	ENDDO
1222
1223	CASE ( 'usm_t_surf' )
1224	!-- surface temperature for surfaces
1225	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_av) ) THEN
1226	ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
1227	surf_usm_h%t_surf_av = 0.0_wp
1228	ENDIF
1229	DO l = 0, 3
1230	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_av) ) THEN
1231	ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
1232	surf_usm_v(l)%t_surf_av = 0.0_wp
1233	ENDIF
1234	ENDDO
1235
1236	CASE ( 'usm_t_wall' )
1237	!-- wall temperature for iwl layer of walls and land
1238	IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN
1239	ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1240	surf_usm_h%t_wall_av = 0.0_wp
1241	ENDIF
1242	DO l = 0, 3
1243	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN
1244	ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1245	surf_usm_v(l)%t_wall_av = 0.0_wp
1246	ENDIF
1247	ENDDO
1248
1249	CASE DEFAULT
1250	CONTINUE
1251
1252	END SELECT
1253
1254	ELSEIF ( mode == 'sum' ) THEN
1255
1256	SELECT CASE ( TRIM( var ) )
1257
1258	CASE ( 'usm_rad_net' )
1259	!-- array of complete radiation balance
1260	DO m = 1, surf_usm_h%ns
1261	surf_usm_h%rad_net_av(m) = &
1262	surf_usm_h%rad_net_av(m) + &
1263	surf_usm_h%rad_net_l(m)
1264	ENDDO
1265	DO l = 0, 3
1266	DO m = 1, surf_usm_v(l)%ns
1267	surf_usm_v(l)%rad_net_av(m) = &
1268	surf_usm_v(l)%rad_net_av(m) + &
1269	surf_usm_v(l)%rad_net_l(m)
1270	ENDDO
1271	ENDDO
1272
1273	CASE ( 'usm_rad_insw' )
1274	!-- array of sw radiation falling to surface after i-th reflection
1275	DO l = startenergy, endenergy
1276	IF ( surfl(id,l) == ids ) THEN
1277	surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
1278	ENDIF
1279	ENDDO
1280
1281	CASE ( 'usm_rad_inlw' )
1282	!-- array of lw radiation falling to surface after i-th reflection
1283	DO l = startenergy, endenergy
1284	IF ( surfl(id,l) == ids ) THEN
1285	surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
1286	ENDIF
1287	ENDDO
1288
1289	CASE ( 'usm_rad_inswdir' )
1290	!-- array of direct sw radiation falling to surface from sun
1291	DO l = startenergy, endenergy
1292	IF ( surfl(id,l) == ids ) THEN
1293	surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
1294	ENDIF
1295	ENDDO
1296
1297	CASE ( 'usm_rad_inswdif' )
1298	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1299	DO l = startenergy, endenergy
1300	IF ( surfl(id,l) == ids ) THEN
1301	surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
1302	ENDIF
1303	ENDDO
1304
1305	CASE ( 'usm_rad_inswref' )
1306	!-- array of sw radiation falling to surface from reflections
1307	DO l = startenergy, endenergy
1308	IF ( surfl(id,l) == ids ) THEN
1309	surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1310	surfinswdir(l) - surfinswdif(l)
1311	ENDIF
1312	ENDDO
1313
1314
1315	CASE ( 'usm_rad_inlwdif' )
1316	!-- array of sw radiation falling to surface after i-th reflection
1317	DO l = startenergy, endenergy
1318	IF ( surfl(id,l) == ids ) THEN
1319	surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1320	surfinswdir(l) - surfinswdif(l)
1321	ENDIF
1322	ENDDO
1323	!
1324	CASE ( 'usm_rad_inlwref' )
1325	!-- array of lw radiation falling to surface from reflections
1326	DO l = startenergy, endenergy
1327	IF ( surfl(id,l) == ids ) THEN
1328	surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
1329	ENDIF
1330	ENDDO
1331
1332	CASE ( 'usm_rad_outsw' )
1333	!-- array of sw radiation emitted from surface after i-th reflection
1334	DO l = startenergy, endenergy
1335	IF ( surfl(id,l) == ids ) THEN
1336	surfinlwref_av(l) = surfinlwref_av(l) + &
1337	surfinlw(l) - surfinlwdif(l)
1338	ENDIF
1339	ENDDO
1340
1341	CASE ( 'usm_rad_outlw' )
1342	!-- array of lw radiation emitted from surface after i-th reflection
1343	DO l = startenergy, endenergy
1344	IF ( surfl(id,l) == ids ) THEN
1345	surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
1346	ENDIF
1347	ENDDO
1348
1349	CASE ( 'usm_rad_ressw' )
1350	!-- array of residua of sw radiation absorbed in surface after last reflection
1351	DO l = startenergy, endenergy
1352	IF ( surfl(id,l) == ids ) THEN
1353	surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
1354	ENDIF
1355	ENDDO
1356
1357	CASE ( 'usm_rad_reslw' )
1358	!-- array of residua of lw radiation absorbed in surface after last reflection
1359	DO l = startenergy, endenergy
1360	IF ( surfl(id,l) == ids ) THEN
1361	surfins_av(l) = surfins_av(l) + surfins(l)
1362	ENDIF
1363	ENDDO
1364
1365	CASE ( 'usm_rad_hf' )
1366	!-- array of heat flux from radiation for surfaces after i-th reflection
1367	DO m = 1, surf_usm_h%ns
1368	surf_usm_h%surfhf_av(m) = &
1369	surf_usm_h%surfhf_av(m) + &
1370	surf_usm_h%surfhf(m)
1371	ENDDO
1372	DO l = 0, 3
1373	DO m = 1, surf_usm_v(l)%ns
1374	surf_usm_v(l)%surfhf_av(m) = &
1375	surf_usm_v(l)%surfhf_av(m) + &
1376	surf_usm_v(l)%surfhf(m)
1377	ENDDO
1378	ENDDO
1379
1380	CASE ( 'usm_wshf' )
1381	!-- array of sensible heat flux from surfaces (land, roof, wall)
1382	DO m = 1, surf_usm_h%ns
1383	surf_usm_h%wshf_eb_av(m) = &
1384	surf_usm_h%wshf_eb_av(m) + &
1385	surf_usm_h%wshf_eb(m)
1386	ENDDO
1387	DO l = 0, 3
1388	DO m = 1, surf_usm_v(l)%ns
1389	surf_usm_v(l)%wshf_eb_av(m) = &
1390	surf_usm_v(l)%wshf_eb_av(m) + &
1391	surf_usm_v(l)%wshf_eb(m)
1392	ENDDO
1393	ENDDO
1394
1395	CASE ( 'usm_wghf' )
1396	!-- array of heat flux from ground (wall, roof, land)
1397	DO m = 1, surf_usm_h%ns
1398	surf_usm_h%wghf_eb_av(m) = &
1399	surf_usm_h%wghf_eb_av(m) + &
1400	surf_usm_h%wghf_eb(m)
1401	ENDDO
1402	DO l = 0, 3
1403	DO m = 1, surf_usm_v(l)%ns
1404	surf_usm_v(l)%wghf_eb_av(m) = &
1405	surf_usm_v(l)%wghf_eb_av(m) + &
1406	surf_usm_v(l)%wghf_eb(m)
1407	ENDDO
1408	ENDDO
1409
1410	CASE ( 'usm_t_surf' )
1411	!-- surface temperature for surfaces
1412	DO m = 1, surf_usm_h%ns
1413	surf_usm_h%t_surf_av(m) = &
1414	surf_usm_h%t_surf_av(m) + &
1415	t_surf_h(m)
1416	ENDDO
1417	DO l = 0, 3
1418	DO m = 1, surf_usm_v(l)%ns
1419	surf_usm_v(l)%t_surf_av(m) = &
1420	surf_usm_v(l)%t_surf_av(m) + &
1421	t_surf_v(l)%t(m)
1422	ENDDO
1423	ENDDO
1424
1425	CASE ( 'usm_t_wall' )
1426	!-- wall temperature for iwl layer of walls and land
1427	DO m = 1, surf_usm_h%ns
1428	surf_usm_h%t_wall_av(iwl,m) = &
1429	surf_usm_h%t_wall_av(iwl,m) + &
1430	t_wall_h(iwl,m)
1431	ENDDO
1432	DO l = 0, 3
1433	DO m = 1, surf_usm_v(l)%ns
1434	surf_usm_v(l)%t_wall_av(iwl,m) = &
1435	surf_usm_v(l)%t_wall_av(iwl,m) + &
1436	t_wall_v(l)%t(iwl,m)
1437	ENDDO
1438	ENDDO
1439
1440	CASE DEFAULT
1441	CONTINUE
1442
1443	END SELECT
1444
1445	ELSEIF ( mode == 'average' ) THEN
1446
1447	SELECT CASE ( TRIM( var ) )
1448
1449	CASE ( 'usm_rad_net' )
1450	!-- array of complete radiation balance
1451	DO m = 1, surf_usm_h%ns
1452	surf_usm_h%rad_net_av(m) = &
1453	surf_usm_h%rad_net_av(m) / &
1454	REAL( average_count_3d, kind=wp )
1455	ENDDO
1456	DO l = 0, 3
1457	DO m = 1, surf_usm_v(l)%ns
1458	surf_usm_v(l)%rad_net_av(m) = &
1459	surf_usm_v(l)%rad_net_av(m) / &
1460	REAL( average_count_3d, kind=wp )
1461	ENDDO
1462	ENDDO
1463
1464	CASE ( 'usm_rad_insw' )
1465	!-- array of sw radiation falling to surface after i-th reflection
1466	DO l = startenergy, endenergy
1467	IF ( surfl(id,l) == ids ) THEN
1468	surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
1469	ENDIF
1470	ENDDO
1471
1472	CASE ( 'usm_rad_inlw' )
1473	!-- array of lw radiation falling to surface after i-th reflection
1474	DO l = startenergy, endenergy
1475	IF ( surfl(id,l) == ids ) THEN
1476	surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
1477	ENDIF
1478	ENDDO
1479
1480	CASE ( 'usm_rad_inswdir' )
1481	!-- array of direct sw radiation falling to surface from sun
1482	DO l = startenergy, endenergy
1483	IF ( surfl(id,l) == ids ) THEN
1484	surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
1485	ENDIF
1486	ENDDO
1487
1488	CASE ( 'usm_rad_inswdif' )
1489	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1490	DO l = startenergy, endenergy
1491	IF ( surfl(id,l) == ids ) THEN
1492	surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
1493	ENDIF
1494	ENDDO
1495
1496	CASE ( 'usm_rad_inswref' )
1497	!-- array of sw radiation falling to surface from reflections
1498	DO l = startenergy, endenergy
1499	IF ( surfl(id,l) == ids ) THEN
1500	surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
1501	ENDIF
1502	ENDDO
1503
1504	CASE ( 'usm_rad_inlwdif' )
1505	!-- array of sw radiation falling to surface after i-th reflection
1506	DO l = startenergy, endenergy
1507	IF ( surfl(id,l) == ids ) THEN
1508	surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
1509	ENDIF
1510	ENDDO
1511
1512	CASE ( 'usm_rad_inlwref' )
1513	!-- array of lw radiation falling to surface from reflections
1514	DO l = startenergy, endenergy
1515	IF ( surfl(id,l) == ids ) THEN
1516	surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
1517	ENDIF
1518	ENDDO
1519
1520	CASE ( 'usm_rad_outsw' )
1521	!-- array of sw radiation emitted from surface after i-th reflection
1522	DO l = startenergy, endenergy
1523	IF ( surfl(id,l) == ids ) THEN
1524	surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
1525	ENDIF
1526	ENDDO
1527
1528	CASE ( 'usm_rad_outlw' )
1529	!-- array of lw radiation emitted from surface after i-th reflection
1530	DO l = startenergy, endenergy
1531	IF ( surfl(id,l) == ids ) THEN
1532	surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
1533	ENDIF
1534	ENDDO
1535
1536	CASE ( 'usm_rad_ressw' )
1537	!-- array of residua of sw radiation absorbed in surface after last reflection
1538	DO l = startenergy, endenergy
1539	IF ( surfl(id,l) == ids ) THEN
1540	surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
1541	ENDIF
1542	ENDDO
1543
1544	CASE ( 'usm_rad_reslw' )
1545	!-- array of residua of lw radiation absorbed in surface after last reflection
1546	DO l = startenergy, endenergy
1547	IF ( surfl(id,l) == ids ) THEN
1548	surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
1549	ENDIF
1550	ENDDO
1551
1552	CASE ( 'usm_rad_hf' )
1553	!-- array of heat flux from radiation for surfaces after i-th reflection
1554	DO m = 1, surf_usm_h%ns
1555	surf_usm_h%surfhf_av(m) = &
1556	surf_usm_h%surfhf_av(m) / &
1557	REAL( average_count_3d, kind=wp )
1558	ENDDO
1559	DO l = 0, 3
1560	DO m = 1, surf_usm_v(l)%ns
1561	surf_usm_v(l)%surfhf_av(m) = &
1562	surf_usm_v(l)%surfhf_av(m) / &
1563	REAL( average_count_3d, kind=wp )
1564	ENDDO
1565	ENDDO
1566
1567	CASE ( 'usm_wshf' )
1568	!-- array of sensible heat flux from surfaces (land, roof, wall)
1569	DO m = 1, surf_usm_h%ns
1570	surf_usm_h%wshf_eb_av(m) = &
1571	surf_usm_h%wshf_eb_av(m) / &
1572	REAL( average_count_3d, kind=wp )
1573	ENDDO
1574	DO l = 0, 3
1575	DO m = 1, surf_usm_v(l)%ns
1576	surf_usm_v(l)%wshf_eb_av(m) = &
1577	surf_usm_v(l)%wshf_eb_av(m) / &
1578	REAL( average_count_3d, kind=wp )
1579	ENDDO
1580	ENDDO
1581
1582	CASE ( 'usm_wghf' )
1583	!-- array of heat flux from ground (wall, roof, land)
1584	DO m = 1, surf_usm_h%ns
1585	surf_usm_h%wghf_eb_av(m) = &
1586	surf_usm_h%wghf_eb_av(m) / &
1587	REAL( average_count_3d, kind=wp )
1588	ENDDO
1589	DO l = 0, 3
1590	DO m = 1, surf_usm_v(l)%ns
1591	surf_usm_v(l)%wghf_eb_av(m) = &
1592	surf_usm_v(l)%wghf_eb_av(m) / &
1593	REAL( average_count_3d, kind=wp )
1594	ENDDO
1595	ENDDO
1596
1597	CASE ( 'usm_t_surf' )
1598	!-- surface temperature for surfaces
1599	DO m = 1, surf_usm_h%ns
1600	surf_usm_h%t_surf_av(m) = &
1601	surf_usm_h%t_surf_av(m) / &
1602	REAL( average_count_3d, kind=wp )
1603	ENDDO
1604	DO l = 0, 3
1605	DO m = 1, surf_usm_v(l)%ns
1606	surf_usm_v(l)%t_surf_av(m) = &
1607	surf_usm_v(l)%t_surf_av(m) / &
1608	REAL( average_count_3d, kind=wp )
1609	ENDDO
1610	ENDDO
1611
1612	CASE ( 'usm_t_wall' )
1613	!-- wall temperature for iwl layer of walls and land
1614	DO m = 1, surf_usm_h%ns
1615	surf_usm_h%t_wall_av(iwl,m) = &
1616	surf_usm_h%t_wall_av(iwl,m) / &
1617	REAL( average_count_3d, kind=wp )
1618	ENDDO
1619	DO l = 0, 3
1620	DO m = 1, surf_usm_v(l)%ns
1621	surf_usm_v(l)%t_wall_av(iwl,m) = &
1622	surf_usm_v(l)%t_wall_av(iwl,m) / &
1623	REAL( average_count_3d, kind=wp )
1624	ENDDO
1625	ENDDO
1626
1627	END SELECT
1628
1629	ENDIF
1630
1631	END SUBROUTINE usm_average_3d_data
1632
1633
1634	!------------------------------------------------------------------------------!
1635	!> Calculates radiation absorbed by box with given size and LAD.
1636	!>
1637	!> Simulates resol**2 rays (by equally spacing a bounding horizontal square
1638	!> conatining all possible rays that would cross the box) and calculates
1639	!> average transparency per ray. Returns fraction of absorbed radiation flux
1640	!> and area for which this fraction is effective.
1641	!------------------------------------------------------------------------------!
1642	PURE SUBROUTINE usm_box_absorb(boxsize, resol, dens, uvec, area, absorb)
1643	IMPLICIT NONE
1644
1645	REAL(wp), DIMENSION(3), INTENT(in) :: &
1646	boxsize, & !< z, y, x size of box in m
1647	uvec !< z, y, x unit vector of incoming flux
1648	INTEGER(iwp), INTENT(in) :: &
1649	resol !< No. of rays in x and y dimensions
1650	REAL(wp), INTENT(in) :: &
1651	dens !< box density (e.g. Leaf Area Density)
1652	REAL(wp), INTENT(out) :: &
1653	area, & !< horizontal area for flux absorbtion
1654	absorb !< fraction of absorbed flux
1655	REAL(wp) :: &
1656	xshift, yshift, &
1657	xmin, xmax, ymin, ymax, &
1658	xorig, yorig, &
1659	dx1, dy1, dz1, dx2, dy2, dz2, &
1660	crdist, &
1661	transp
1662	INTEGER(iwp) :: &
1663	i, j
1664
1665	xshift = uvec(3) / uvec(1) * boxsize(1)
1666	xmin = min(0._wp, -xshift)
1667	xmax = boxsize(3) + max(0._wp, -xshift)
1668	yshift = uvec(2) / uvec(1) * boxsize(1)
1669	ymin = min(0._wp, -yshift)
1670	ymax = boxsize(2) + max(0._wp, -yshift)
1671
1672	transp = 0._wp
1673	DO i = 1, resol
1674	xorig = xmin + (xmax-xmin) * (i-.5_wp) / resol
1675	DO j = 1, resol
1676	yorig = ymin + (ymax-ymin) * (j-.5_wp) / resol
1677
1678	dz1 = 0._wp
1679	dz2 = boxsize(1)/uvec(1)
1680
1681	IF ( uvec(2) > 0._wp ) THEN
1682	dy1 = -yorig / uvec(2) !< crossing with y=0
1683	dy2 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2)
1684	ELSE IF ( uvec(2) < 0._wp ) THEN
1685	dy1 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2)
1686	dy2 = -yorig / uvec(2) !< crossing with y=0
1687	ELSE !uvec(2)==0
1688	dy1 = -huge(1._wp)
1689	dy2 = huge(1._wp)
1690	ENDIF
1691
1692	IF ( uvec(3) > 0._wp ) THEN
1693	dx1 = -xorig / uvec(3) !< crossing with x=0
1694	dx2 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3)
1695	ELSE IF ( uvec(3) < 0._wp ) THEN
1696	dx1 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3)
1697	dx2 = -xorig / uvec(3) !< crossing with x=0
1698	ELSE !uvec(1)==0
1699	dx1 = -huge(1._wp)
1700	dx2 = huge(1._wp)
1701	ENDIF
1702
1703	crdist = max(0._wp, (min(dz2, dy2, dx2) - max(dz1, dy1, dx1)))
1704	transp = transp + exp(-ext_coef * dens * crdist)
1705	ENDDO
1706	ENDDO
1707	transp = transp / resol**2
1708	area = (boxsize(3)+xshift)*(boxsize(2)+yshift)
1709	absorb = 1._wp - transp
1710
1711	END SUBROUTINE usm_box_absorb
1712
1713
1714	!------------------------------------------------------------------------------!
1715	! Description:
1716	! ------------
1717	!> This subroutine splits direct and diffusion dw radiation
1718	!> It sould not be called in case the radiation model already does it
1719	!> It follows <CITATION>
1720	!------------------------------------------------------------------------------!
1721	SUBROUTINE usm_calc_diffusion_radiation
1722
1723	REAL(wp), PARAMETER :: sol_const = 1367.0_wp !< solar conbstant
1724	REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation
1725	INTEGER(iwp) :: i, j
1726	REAL(wp), PARAMETER :: year_seconds = 86400._wp * 365._wp
1727	REAL(wp) :: year_angle !< angle
1728	REAL(wp) :: etr !< extraterestrial radiation
1729	REAL(wp) :: corrected_solarUp !< corrected solar up radiation
1730	REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation
1731	REAL(wp) :: clearnessIndex !< clearness index
1732	REAL(wp) :: diff_frac !< diffusion fraction of the radiation
1733
1734
1735	!-- Calculate current day and time based on the initial values and simulation time
1736	year_angle = ((day_init*86400) + time_utc_init+time_since_reference_point) &
1737	/ year_seconds * 2.0_wp * pi
1738
1739	etr = sol_const * (1.00011_wp + &
1740	0.034221_wp * cos(year_angle) + &
1741	0.001280_wp * sin(year_angle) + &
1742	0.000719_wp * cos(2.0_wp * year_angle) + &
1743	0.000077_wp * sin(2.0_wp * year_angle))
1744
1745	!--
1746	!-- Under a very low angle, we keep extraterestrial radiation at
1747	!-- the last small value, therefore the clearness index will be pushed
1748	!-- towards 0 while keeping full continuity.
1749	!--
1750	IF ( zenith(0) <= lowest_solarUp ) THEN
1751	corrected_solarUp = lowest_solarUp
1752	ELSE
1753	corrected_solarUp = zenith(0)
1754	ENDIF
1755
1756	horizontalETR = etr * corrected_solarUp
1757
1758	DO i = nxlg, nxrg
1759	DO j = nysg, nyng
1760	clearnessIndex = rad_sw_in(0,j,i) / horizontalETR
1761	diff_frac = 1.0_wp / (1.0_wp + exp(-5.0033_wp + 8.6025_wp * clearnessIndex))
1762	rad_sw_in_diff(j,i) = rad_sw_in(0,j,i) * diff_frac
1763	rad_sw_in_dir(j,i) = rad_sw_in(0,j,i) * (1.0_wp - diff_frac)
1764	rad_lw_in_diff(j,i) = rad_lw_in(0,j,i)
1765	ENDDO
1766	ENDDO
1767
1768	END SUBROUTINE usm_calc_diffusion_radiation
1769
1770
1771	!------------------------------------------------------------------------------!
1772	! Description:
1773	! ------------
1774	!> Calculates shape view factors SVF and plant sink canopy factors PSCF
1775	!> !!!!!DESCRIPTION!!!!!!!!!!
1776	!------------------------------------------------------------------------------!
1777	SUBROUTINE usm_calc_svf
1778
1779	IMPLICIT NONE
1780
1781	INTEGER(iwp) :: i, j, k, l, d, ip, jp
1782	INTEGER(iwp) :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrtt, imrtf
1783	INTEGER(iwp) :: sd, td, ioln, iproc
1784	REAL(wp), DIMENSION(0:9) :: facearea
1785	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzterrl, planthl
1786	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csflt, pcsflt
1787	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: kcsflt,kpcsflt
1788	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: icsflt,dcsflt,ipcsflt,dpcsflt
1789	REAL(wp), DIMENSION(3) :: uv
1790	LOGICAL :: visible
1791	REAL(wp), DIMENSION(3) :: sa, ta !< real coordinates z,y,x of source and target
1792	REAL(wp) :: transparency, rirrf, sqdist, svfsum
1793	INTEGER(iwp) :: isurflt, isurfs, isurflt_prev
1794	INTEGER(iwp) :: itx, ity, itz
1795	CHARACTER(len=7) :: pid_char = ''
1796	INTEGER(iwp) :: win_lad, minfo
1797	REAL(wp), DIMENSION(:,:,:), POINTER :: lad_s_rma !< fortran pointer, but lower bounds are 1
1798	TYPE(c_ptr) :: lad_s_rma_p !< allocated c pointer
1799	INTEGER(kind=MPI_ADDRESS_KIND) :: size_lad_rma
1800	!
1801	!-- calculation of the SVF
1802	CALL location_message( ' calculation of SVF and CSF', .TRUE. )
1803	!
1804	!-- precalculate face areas for different face directions using normal vector
1805	DO d = 0, 9
1806	facearea(d) = 1._wp
1807	IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
1808	IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
1809	IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
1810	ENDDO
1811
1812	!-- initialize variables and temporary arrays for calculation of svf and csf
1813	nsvfl = 0
1814	ncsfl = 0
1815	nsvfla = gasize
1816	msvf = 1
1817	ALLOCATE( asvf1(nsvfla) )
1818	asvf => asvf1
1819	IF ( plant_canopy ) THEN
1820	ncsfla = gasize
1821	mcsf = 1
1822	ALLOCATE( acsf1(ncsfla) )
1823	acsf => acsf1
1824	ENDIF
1825
1826	!-- initialize temporary terrain and plant canopy height arrays (global 2D array!)
1827	ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) )
1828	#if defined( __parallel )
1829	ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
1830	nzterrl = MAXLOC( &
1831	MERGE( 1, 0, &
1832	BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) &
1833	), DIM = 1 &
1834	) - 1 ! = nzb_s_inner(nys:nyn,nxl:nxr)
1835	CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, &
1836	nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
1837	DEALLOCATE(nzterrl)
1838	#else
1839	nzterr = RESHAPE( MAXLOC( &
1840	MERGE( 1, 0, &
1841	BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) &
1842	), DIM = 1 &
1843	) - 1, &
1844	(/(nx+1)*(ny+1)/) &
1845	)
1846	#endif
1847	IF ( plant_canopy ) THEN
1848	ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) )
1849	maxboxesg = nx + ny + nzu + 1
1850	!-- temporary arrays storing values for csf calculation during raytracing
1851	ALLOCATE( boxes(3, maxboxesg) )
1852	ALLOCATE( crlens(maxboxesg) )
1853
1854	#if defined( __parallel )
1855	ALLOCATE( planthl(nys:nyn,nxl:nxr) )
1856	planthl = pch(nys:nyn,nxl:nxr)
1857
1858	CALL MPI_AllGather( planthl, nnx*nny, MPI_INTEGER, &
1859	plantt, nnx*nny, MPI_INTEGER, comm2d, ierr )
1860	DEALLOCATE( planthl )
1861
1862	!-- temporary arrays storing values for csf calculation during raytracing
1863	ALLOCATE( lad_ip(maxboxesg) )
1864	ALLOCATE( lad_disp(maxboxesg) )
1865
1866	IF ( usm_lad_rma ) THEN
1867	ALLOCATE( lad_s_ray(maxboxesg) )
1868
1869	! set conditions for RMA communication
1870	CALL MPI_Info_create(minfo, ierr)
1871	CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
1872	CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr)
1873	CALL MPI_Info_set(minfo, 'same_size', 'true', ierr)
1874	CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr)
1875
1876	!-- Allocate and initialize the MPI RMA window
1877	!-- must be in accordance with allocation of lad_s in plant_canopy_model
1878	!-- optimization of memory should be done
1879	!-- Argument X of function c_sizeof(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now
1880	size_lad_rma = c_sizeof(1.0_wp)nnxnny*nzu
1881	CALL MPI_Win_allocate(size_lad_rma, c_sizeof(1.0_wp), minfo, comm2d, &
1882	lad_s_rma_p, win_lad, ierr)
1883	CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzu, nny, nnx /))
1884	usm_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:)
1885	ELSE
1886	ALLOCATE(usm_lad(nzub:nzut, nys:nyn, nxl:nxr))
1887	ENDIF
1888	#else
1889	plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) )
1890	ALLOCATE(usm_lad(nzub:nzut, nys:nyn, nxl:nxr))
1891	#endif
1892	usm_lad(:,:,:) = 0._wp
1893	DO i = nxl, nxr
1894	DO j = nys, nyn
1895	k = MAXLOC( &
1896	MERGE( 1, 0, &
1897	BTEST( wall_flags_0(:,j,i), 12 ) &
1898	), DIM = 1 &
1899	) - 1
1900
1901	usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
1902	ENDDO
1903	ENDDO
1904
1905	#if defined( __parallel )
1906	IF ( usm_lad_rma ) THEN
1907	CALL MPI_Info_free(minfo, ierr)
1908	CALL MPI_Win_lock_all(0, win_lad, ierr)
1909	ELSE
1910	ALLOCATE( usm_lad_g(0:(nx+1)(ny+1)nzu-1) )
1911	CALL MPI_AllGather( usm_lad, nnxnnynzu, MPI_REAL, &
1912	usm_lad_g, nnxnnynzu, MPI_REAL, comm2d, ierr )
1913	ENDIF
1914	#endif
1915	ENDIF
1916
1917	IF ( mrt_factors ) THEN
1918	OPEN(153, file='MRT_TARGETS', access='SEQUENTIAL', &
1919	action='READ', status='OLD', form='FORMATTED', err=524)
1920	OPEN(154, file='MRT_FACTORS'//myid_char, access='DIRECT', recl=(54+28), &
1921	action='WRITE', status='REPLACE', form='UNFORMATTED', err=525)
1922	imrtf = 1
1923	DO
1924	READ(153, *, end=526, err=524) imrtt, i, j, k
1925	IF ( i < nxl .OR. i > nxr &
1926	.OR. j < nys .OR. j > nyn ) CYCLE
1927	ta = (/ REAL(k), REAL(j), REAL(i) /)
1928
1929	DO isurfs = 1, nsurf
1930	IF ( .NOT. usm_facing(i, j, k, -1, &
1931	surf(ix, isurfs), surf(iy, isurfs), &
1932	surf(iz, isurfs), surf(id, isurfs)) ) THEN
1933	CYCLE
1934	ENDIF
1935
1936	sd = surf(id, isurfs)
1937	sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &
1938	REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &
1939	REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)
1940
1941	!-- unit vector source -> target
1942	uv = (/ (ta(1)-sa(1))dz, (ta(2)-sa(2))dy, (ta(3)-sa(3))*dx /)
1943	sqdist = SUM(uv(:)**2)
1944	uv = uv / SQRT(sqdist)
1945
1946	!-- irradiance factor - see svf. Here we consider that target face is always normal,
1947	!-- i.e. the second dot product equals 1
1948	rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) &
1949	/ (pi * sqdist) * facearea(sd)
1950
1951	!-- raytrace while not creating any canopy sink factors
1952	CALL usm_raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., &
1953	visible, transparency, win_lad)
1954	IF ( .NOT. visible ) CYCLE
1955
1956	!rsvf = rirrf * transparency
1957	WRITE(154, rec=imrtf, err=525) INT(imrtt, kind=4), &
1958	INT(surf(id, isurfs), kind=4), &
1959	INT(surf(iz, isurfs), kind=4), &
1960	INT(surf(iy, isurfs), kind=4), &
1961	INT(surf(ix, isurfs), kind=4), &
1962	REAL(rirrf, kind=8), REAL(transparency, kind=8)
1963	imrtf = imrtf + 1
1964
1965	ENDDO !< isurfs
1966	ENDDO !< MRT_TARGETS record
1967
1968	524 message_string = 'error reading file MRT_TARGETS'
1969	CALL message( 'usm_calc_svf', 'PA0524', 1, 2, 0, 6, 0 )
1970
1971	525 message_string = 'error writing file MRT_FACTORS'//myid_char
1972	CALL message( 'usm_calc_svf', 'PA0525', 1, 2, 0, 6, 0 )
1973
1974	526 CLOSE(153)
1975	CLOSE(154)
1976	ENDIF !< mrt_factors
1977
1978
1979	DO isurflt = 1, nsurfl
1980	!-- determine face centers
1981	td = surfl(id, isurflt)
1982	IF ( td >= isky .AND. .NOT. plant_canopy ) CYCLE
1983	ta = (/ REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td), &
1984	REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td), &
1985	REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td) /)
1986	DO isurfs = 1, nsurf
1987	IF ( .NOT. usm_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
1988	surfl(iz, isurflt), surfl(id, isurflt), &
1989	surf(ix, isurfs), surf(iy, isurfs), &
1990	surf(iz, isurfs), surf(id, isurfs)) ) THEN
1991	CYCLE
1992	ENDIF
1993
1994	sd = surf(id, isurfs)
1995	sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &
1996	REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &
1997	REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)
1998
1999	!-- unit vector source -> target
2000	uv = (/ (ta(1)-sa(1))dz, (ta(2)-sa(2))dy, (ta(3)-sa(3))*dx /)
2001	sqdist = SUM(uv(:)**2)
2002	uv = uv / SQRT(sqdist)
2003
2004	!-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
2005	rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
2006	* dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction
2007	/ (pi * sqdist) & ! square of distance between centers
2008	* facearea(sd)
2009
2010	!-- raytrace + process plant canopy sinks within
2011	CALL usm_raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
2012	visible, transparency, win_lad)
2013
2014	IF ( .NOT. visible ) CYCLE
2015	IF ( td >= isky ) CYCLE !< we calculated these only for raytracing
2016	!< to find plant canopy sinks, we don't need svf for them
2017	! rsvf = rirrf * transparency
2018
2019	!-- write to the svf array
2020	nsvfl = nsvfl + 1
2021	!-- check dimmension of asvf array and enlarge it if needed
2022	IF ( nsvfla < nsvfl ) THEN
2023	k = nsvfla * 2
2024	IF ( msvf == 0 ) THEN
2025	msvf = 1
2026	ALLOCATE( asvf1(k) )
2027	asvf => asvf1
2028	asvf1(1:nsvfla) = asvf2
2029	DEALLOCATE( asvf2 )
2030	ELSE
2031	msvf = 0
2032	ALLOCATE( asvf2(k) )
2033	asvf => asvf2
2034	asvf2(1:nsvfla) = asvf1
2035	DEALLOCATE( asvf1 )
2036	ENDIF
2037	nsvfla = k
2038	ENDIF
2039	!-- write svf values into the array
2040	asvf(nsvfl)%isurflt = isurflt
2041	asvf(nsvfl)%isurfs = isurfs
2042	asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency
2043	asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor
2044	ENDDO
2045	ENDDO
2046
2047	CALL location_message( ' waiting for completion of SVF and CSF calculation in all processes', .TRUE. )
2048	!-- deallocate temporary global arrays
2049	DEALLOCATE(nzterr)
2050
2051	IF ( plant_canopy ) THEN
2052	!-- finalize mpi_rma communication and deallocate temporary arrays
2053	#if defined( __parallel )
2054	IF ( usm_lad_rma ) THEN
2055	CALL MPI_Win_flush_all(win_lad, ierr)
2056	!-- unlock MPI window
2057	CALL MPI_Win_unlock_all(win_lad, ierr)
2058	!-- free MPI window
2059	CALL MPI_Win_free(win_lad, ierr)
2060
2061	!-- deallocate temporary arrays storing values for csf calculation during raytracing
2062	DEALLOCATE( lad_s_ray )
2063	!-- usm_lad is the pointer to lad_s_rma in case of usm_lad_rma
2064	!-- and must not be deallocated here
2065	ELSE
2066	DEALLOCATE(usm_lad)
2067	DEALLOCATE(usm_lad_g)
2068	ENDIF
2069	#else
2070	DEALLOCATE(usm_lad)
2071	#endif
2072	DEALLOCATE( boxes )
2073	DEALLOCATE( crlens )
2074	DEALLOCATE( plantt )
2075	ENDIF
2076
2077	CALL location_message( ' calculation of the complete SVF array', .TRUE. )
2078
2079	!-- sort svf ( a version of quicksort )
2080	CALL quicksort_svf(asvf,1,nsvfl)
2081
2082	ALLOCATE( svf(ndsvf,nsvfl) )
2083	ALLOCATE( svfsurf(idsvf,nsvfl) )
2084
2085	!< load svf from the structure array to plain arrays
2086	isurflt_prev = -1
2087	ksvf = 1
2088	svfsum = 0._wp
2089	DO isvf = 1, nsvfl
2090	!-- normalize svf per target face
2091	IF ( asvf(ksvf)%isurflt /= isurflt_prev ) THEN
2092	IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN
2093	!-- TODO detect and log when normalization differs too much from 1
2094	svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum
2095	ENDIF
2096	isurflt_prev = asvf(ksvf)%isurflt
2097	isvf_surflt = isvf
2098	svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
2099	ELSE
2100	svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
2101	ENDIF
2102
2103	svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
2104	svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
2105
2106	!-- next element
2107	ksvf = ksvf + 1
2108	ENDDO
2109
2110	IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN
2111	!-- TODO detect and log when normalization differs too much from 1
2112	svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum
2113	ENDIF
2114
2115	!-- deallocate temporary asvf array
2116	!-- DEALLOCATE(asvf) - ifort has a problem with deallocation of allocatable target
2117	!-- via pointing pointer - we need to test original targets
2118	IF ( ALLOCATED(asvf1) ) THEN
2119	DEALLOCATE(asvf1)
2120	ENDIF
2121	IF ( ALLOCATED(asvf2) ) THEN
2122	DEALLOCATE(asvf2)
2123	ENDIF
2124
2125	npcsfl = 0
2126	IF ( plant_canopy ) THEN
2127
2128	CALL location_message( ' calculation of the complete CSF array', .TRUE. )
2129
2130	!-- sort and merge csf for the last time, keeping the array size to minimum
2131	CALL usm_merge_and_grow_csf(-1)
2132
2133	!-- aggregate csb among processors
2134	!-- allocate necessary arrays
2135	ALLOCATE( csflt(ndcsf,max(ncsfl,ndcsf)) )
2136	ALLOCATE( kcsflt(kdcsf,max(ncsfl,kdcsf)) )
2137	ALLOCATE( icsflt(0:numprocs-1) )
2138	ALLOCATE( dcsflt(0:numprocs-1) )
2139	ALLOCATE( ipcsflt(0:numprocs-1) )
2140	ALLOCATE( dpcsflt(0:numprocs-1) )
2141
2142	!-- fill out arrays of csf values and
2143	!-- arrays of number of elements and displacements
2144	!-- for particular precessors
2145	icsflt = 0
2146	dcsflt = 0
2147	ip = -1
2148	j = -1
2149	d = 0
2150	DO kcsf = 1, ncsfl
2151	j = j+1
2152	IF ( acsf(kcsf)%ip /= ip ) THEN
2153	!-- new block of the processor
2154	!-- number of elements of previous block
2155	IF ( ip>=0) icsflt(ip) = j
2156	d = d+j
2157	!-- blank blocks
2158	DO jp = ip+1, acsf(kcsf)%ip-1
2159	!-- number of elements is zero, displacement is equal to previous
2160	icsflt(jp) = 0
2161	dcsflt(jp) = d
2162	ENDDO
2163	!-- the actual block
2164	ip = acsf(kcsf)%ip
2165	dcsflt(ip) = d
2166	j = 0
2167	ENDIF
2168	!-- fill out real values of rsvf, rtransp
2169	csflt(1,kcsf) = acsf(kcsf)%rsvf
2170	csflt(2,kcsf) = acsf(kcsf)%rtransp
2171	!-- fill out integer values of itz,ity,itx,isurfs
2172	kcsflt(1,kcsf) = acsf(kcsf)%itz
2173	kcsflt(2,kcsf) = acsf(kcsf)%ity
2174	kcsflt(3,kcsf) = acsf(kcsf)%itx
2175	kcsflt(4,kcsf) = acsf(kcsf)%isurfs
2176	ENDDO
2177	!-- last blank blocks at the end of array
2178	j = j+1
2179	IF ( ip>=0 ) icsflt(ip) = j
2180	d = d+j
2181	DO jp = ip+1, numprocs-1
2182	!-- number of elements is zero, displacement is equal to previous
2183	icsflt(jp) = 0
2184	dcsflt(jp) = d
2185	ENDDO
2186
2187	!-- deallocate temporary acsf array
2188	!-- DEALLOCATE(acsf) - ifort has a problem with deallocation of allocatable target
2189	!-- via pointing pointer - we need to test original targets
2190	IF ( ALLOCATED(acsf1) ) THEN
2191	DEALLOCATE(acsf1)
2192	ENDIF
2193	IF ( ALLOCATED(acsf2) ) THEN
2194	DEALLOCATE(acsf2)
2195	ENDIF
2196
2197	#if defined( __parallel )
2198	!-- scatter and gather the number of elements to and from all processor
2199	!-- and calculate displacements
2200	CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
2201
2202	npcsfl = SUM(ipcsflt)
2203	d = 0
2204	DO i = 0, numprocs-1
2205	dpcsflt(i) = d
2206	d = d + ipcsflt(i)
2207	ENDDO
2208
2209	!-- exchange csf fields between processors
2210	ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) )
2211	ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) )
2212	CALL MPI_AlltoAllv(csflt, ndcsficsflt, ndcsfdcsflt, MPI_REAL, &
2213	pcsflt, ndcsfipcsflt, ndcsfdpcsflt, MPI_REAL, comm2d, ierr)
2214	CALL MPI_AlltoAllv(kcsflt, kdcsficsflt, kdcsfdcsflt, MPI_INTEGER, &
2215	kpcsflt, kdcsfipcsflt, kdcsfdpcsflt, MPI_INTEGER, comm2d, ierr)
2216
2217	#else
2218	npcsfl = ncsfl
2219	ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) )
2220	ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) )
2221	pcsflt = csflt
2222	kpcsflt = kcsflt
2223	#endif
2224
2225	!-- deallocate temporary arrays
2226	DEALLOCATE( csflt )
2227	DEALLOCATE( kcsflt )
2228	DEALLOCATE( icsflt )
2229	DEALLOCATE( dcsflt )
2230	DEALLOCATE( ipcsflt )
2231	DEALLOCATE( dpcsflt )
2232
2233	!-- sort csf ( a version of quicksort )
2234	CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl)
2235
2236	!-- aggregate canopy sink factor records with identical box & source
2237	!-- againg across all values from all processors
2238	IF ( npcsfl > 0 ) THEN
2239	icsf = 1 !< reading index
2240	kcsf = 1 !< writing index
2241	DO while (icsf < npcsfl)
2242	!-- here kpcsf(kcsf) already has values from kpcsf(icsf)
2243	IF ( kpcsflt(3,icsf) == kpcsflt(3,icsf+1) .AND. &
2244	kpcsflt(2,icsf) == kpcsflt(2,icsf+1) .AND. &
2245	kpcsflt(1,icsf) == kpcsflt(1,icsf+1) .AND. &
2246	kpcsflt(4,icsf) == kpcsflt(4,icsf+1) ) THEN
2247	!-- We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
2248	!-- probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
2249	!-- might mean that the traced beam passes longer through the canopy box.
2250	IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) ) THEN
2251	pcsflt(2,kcsf) = pcsflt(2,icsf+1)
2252	ENDIF
2253	pcsflt(1,kcsf) = pcsflt(1,kcsf) + pcsflt(1,icsf+1)
2254
2255	!-- advance reading index, keep writing index
2256	icsf = icsf + 1
2257	ELSE
2258	!-- not identical, just advance and copy
2259	icsf = icsf + 1
2260	kcsf = kcsf + 1
2261	kpcsflt(:,kcsf) = kpcsflt(:,icsf)
2262	pcsflt(:,kcsf) = pcsflt(:,icsf)
2263	ENDIF
2264	ENDDO
2265	!-- last written item is now also the last item in valid part of array
2266	npcsfl = kcsf
2267	ENDIF
2268
2269	ncsfl = npcsfl
2270	IF ( ncsfl > 0 ) THEN
2271	ALLOCATE( csf(ndcsf,ncsfl) )
2272	ALLOCATE( csfsurf(idcsf,ncsfl) )
2273	DO icsf = 1, ncsfl
2274	csf(:,icsf) = pcsflt(:,icsf)
2275	csfsurf(1,icsf) = gridpcbl(kpcsflt(1,icsf),kpcsflt(2,icsf),kpcsflt(3,icsf))
2276	csfsurf(2,icsf) = kpcsflt(4,icsf)
2277	ENDDO
2278	ENDIF
2279
2280	!-- deallocation of temporary arrays
2281	DEALLOCATE( pcsflt )
2282	DEALLOCATE( kpcsflt )
2283
2284	ENDIF
2285
2286	RETURN
2287
2288	301 WRITE( message_string, * ) &
2289	'I/O error when processing shape view factors / ', &
2290	'plant canopy sink factors / direct irradiance factors.'
2291	CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
2292
2293	END SUBROUTINE usm_calc_svf
2294
2295
2296	!------------------------------------------------------------------------------!
2297	!
2298	! Description:
2299	! ------------
2300	!> Subroutine checks variables and assigns units.
2301	!> It is caaled out from subroutine check_parameters.
2302	!------------------------------------------------------------------------------!
2303	SUBROUTINE usm_check_data_output( variable, unit )
2304
2305	IMPLICIT NONE
2306
2307	CHARACTER (len=*),INTENT(IN) :: variable !:
2308	CHARACTER (len=*),INTENT(OUT) :: unit !:
2309
2310	CHARACTER (len=varnamelength) :: var
2311
2312	var = TRIM(variable)
2313	IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. &
2314	var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. &
2315	var(1:16) == 'usm_rad_inswdif_' .OR. var(1:16) == 'usm_rad_inswref_' .OR. &
2316	var(1:16) == 'usm_rad_inlwdif_' .OR. var(1:16) == 'usm_rad_inlwref_' .OR. &
2317	var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. &
2318	var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. &
2319	var(1:11) == 'usm_rad_hf_' .OR. &
2320	var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' ) THEN
2321	unit = 'W/m2'
2322	ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' ) THEN
2323	unit = 'K'
2324	ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. &
2325	var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. &
2326	var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis') THEN
2327	unit = '1'
2328	ELSE
2329	unit = 'illegal'
2330	ENDIF
2331
2332	END SUBROUTINE usm_check_data_output
2333
2334
2335	!------------------------------------------------------------------------------!
2336	! Description:
2337	! ------------
2338	!> Check parameters routine for urban surface model
2339	!------------------------------------------------------------------------------!
2340	SUBROUTINE usm_check_parameters
2341
2342	USE control_parameters, &
2343	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing, &
2344	lsf_surf, topography
2345
2346	!
2347	!-- Dirichlet boundary conditions are required as the surface fluxes are
2348	!-- calculated from the temperature/humidity gradients in the urban surface
2349	!-- model
2350	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
2351	message_string = 'urban surface model requires setting of '// &
2352	'bc_pt_b = "dirichlet" and '// &
2353	'bc_q_b = "dirichlet"'
2354	CALL message( 'check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
2355	ENDIF
2356
2357	IF ( .NOT. constant_flux_layer ) THEN
2358	message_string = 'urban surface model requires '// &
2359	'constant_flux_layer = .T.'
2360	CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
2361	ENDIF
2362	!
2363	!-- Surface forcing has to be disabled for LSF in case of enabled
2364	!-- urban surface module
2365	IF ( large_scale_forcing ) THEN
2366	lsf_surf = .FALSE.
2367	ENDIF
2368	!
2369	!-- Topography
2370	IF ( topography == 'flat' ) THEN
2371	message_string = 'topography /= "flat" is required '// &
2372	'when using the urban surface model'
2373	CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2374	ENDIF
2375
2376
2377	END SUBROUTINE usm_check_parameters
2378
2379
2380	!------------------------------------------------------------------------------!
2381	!
2382	! Description:
2383	! ------------
2384	!> Output of the 3D-arrays in netCDF and/or AVS format
2385	!> for variables of urban_surface model.
2386	!> It resorts the urban surface module output quantities from surf style
2387	!> indexing into temporary 3D array with indices (i,j,k).
2388	!> It is called from subroutine data_output_3d.
2389	!------------------------------------------------------------------------------!
2390	SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2391
2392	IMPLICIT NONE
2393
2394	INTEGER(iwp), INTENT(IN) :: av !<
2395	CHARACTER (len=*), INTENT(IN) :: variable !<
2396	INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0)
2397	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
2398	LOGICAL, INTENT(OUT) :: found !<
2399	REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) :: local_pf !< sp - it has to correspond to module data_output_3d
2400	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: temp_pf !< temp array for urban surface output procedure
2401
2402	CHARACTER (len=varnamelength) :: var, surfid
2403	INTEGER(iwp), PARAMETER :: nd = 5
2404	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
2405	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iroof, isouth, inorth, iwest, ieast /)
2406	INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart
2407	INTEGER(iwp), DIMENSION(0:nd-1) :: dirend
2408	INTEGER(iwp) :: ids,isurf,isvf,isurfs,isurflt
2409	INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat, l, m
2410	INTEGER(iwp) :: k_topo !< topography top index
2411
2412	dirstart = (/ startland, startwall, startwall, startwall, startwall /)
2413	dirend = (/ endland, endwall, endwall, endwall, endwall /)
2414
2415	found = .TRUE.
2416	temp_pf = -1._wp
2417
2418	ids = -1
2419	var = TRIM(variable)
2420	DO i = 0, nd-1
2421	k = len(TRIM(var))
2422	j = len(TRIM(dirname(i)))
2423	IF ( var(k-j+1:k) == dirname(i) ) THEN
2424	ids = i
2425	var = var(:k-j)
2426	EXIT
2427	ENDIF
2428	ENDDO
2429	IF ( ids == -1 ) THEN
2430	var = TRIM(variable)
2431	ENDIF
2432	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
2433	!-- wall layers
2434	READ(var(12:12), '(I1)', iostat=istat ) iwl
2435	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2436	var = var(1:10)
2437	ENDIF
2438	ENDIF
2439	IF ( (var(1:8) == 'usm_svf_' .OR. var(1:8) == 'usm_dif_') .AND. len(TRIM(var)) >= 13 ) THEN
2440	!-- svf values to particular surface
2441	surfid = var(9:)
2442	i = index(surfid,'_')
2443	j = index(surfid(i+1:),'_')
2444	READ(surfid(1:i-1),*, iostat=istat ) is
2445	IF ( istat == 0 ) THEN
2446	READ(surfid(i+1:i+j-1),*, iostat=istat ) js
2447	ENDIF
2448	IF ( istat == 0 ) THEN
2449	READ(surfid(i+j+1:),*, iostat=istat ) ks
2450	ENDIF
2451	IF ( istat == 0 ) THEN
2452	var = var(1:7)
2453	ENDIF
2454	ENDIF
2455
2456	SELECT CASE ( TRIM(var) )
2457
2458	CASE ( 'usm_surfz' )
2459	!-- array of lw radiation falling to local surface after i-th reflection
2460	DO m = 1, surf_usm_h%ns
2461	i = surf_usm_h%i(m)
2462	j = surf_usm_h%j(m)
2463	k = surf_usm_h%k(m)
2464	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
2465	ENDDO
2466	DO l = 0, 3
2467	DO m = 1, surf_usm_v(l)%ns
2468	i = surf_usm_v(l)%i(m)
2469	j = surf_usm_v(l)%j(m)
2470	k = surf_usm_v(l)%k(m)
2471	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
2472	ENDDO
2473	ENDDO
2474
2475	CASE ( 'usm_surfcat' )
2476	!-- surface category
2477	DO m = 1, surf_usm_h%ns
2478	i = surf_usm_h%i(m)
2479	j = surf_usm_h%j(m)
2480	k = surf_usm_h%k(m)
2481	temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2482	ENDDO
2483	DO l = 0, 3
2484	DO m = 1, surf_usm_v(l)%ns
2485	i = surf_usm_v(l)%i(m)
2486	j = surf_usm_v(l)%j(m)
2487	k = surf_usm_v(l)%k(m)
2488	temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
2489	ENDDO
2490	ENDDO
2491
2492	CASE ( 'usm_surfalb' )
2493	!-- surface albedo
2494	DO m = 1, surf_usm_h%ns
2495	i = surf_usm_h%i(m)
2496	j = surf_usm_h%j(m)
2497	k = surf_usm_h%k(m)
2498	temp_pf(k,j,i) = surf_usm_h%albedo_surf(m)
2499	ENDDO
2500	DO l = 0, 3
2501	DO m = 1, surf_usm_v(l)%ns
2502	i = surf_usm_v(l)%i(m)
2503	j = surf_usm_v(l)%j(m)
2504	k = surf_usm_v(l)%k(m)
2505	temp_pf(k,j,i) = surf_usm_v(l)%albedo_surf(m)
2506	ENDDO
2507	ENDDO
2508
2509	CASE ( 'usm_surfemis' )
2510	!-- surface albedo
2511	DO m = 1, surf_usm_h%ns
2512	i = surf_usm_h%i(m)
2513	j = surf_usm_h%j(m)
2514	k = surf_usm_h%k(m)
2515	temp_pf(k,j,i) = surf_usm_h%emiss_surf(m)
2516	ENDDO
2517	DO l = 0, 3
2518	DO m = 1, surf_usm_v(l)%ns
2519	i = surf_usm_v(l)%i(m)
2520	j = surf_usm_v(l)%j(m)
2521	k = surf_usm_v(l)%k(m)
2522	temp_pf(k,j,i) = surf_usm_v(l)%emiss_surf(m)
2523	ENDDO
2524	ENDDO
2525	!
2526	!-- Not adjusted so far
2527	CASE ( 'usm_svf', 'usm_dif' )
2528	!-- shape view factors or iradiance factors to selected surface
2529	IF ( TRIM(var)=='usm_svf' ) THEN
2530	k = 1
2531	ELSE
2532	k = 2
2533	ENDIF
2534	DO isvf = 1, nsvfl
2535	isurflt = svfsurf(1, isvf)
2536	isurfs = svfsurf(2, isvf)
2537
2538	IF ( surf(ix,isurfs) == is .AND. surf(iy,isurfs) == js .AND. &
2539	surf(iz,isurfs) == ks .AND. surf(id,isurfs) == ids ) THEN
2540	!-- correct source surface
2541	temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
2542	ENDIF
2543	ENDDO
2544
2545	CASE ( 'usm_rad_net' )
2546	!-- array of complete radiation balance
2547	IF ( av == 0 ) THEN
2548	DO m = 1, surf_usm_h%ns
2549	i = surf_usm_h%i(m)
2550	j = surf_usm_h%j(m)
2551	k = surf_usm_h%k(m)
2552	temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
2553	ENDDO
2554	DO l = 0, 3
2555	DO m = 1, surf_usm_v(l)%ns
2556	i = surf_usm_v(l)%i(m)
2557	j = surf_usm_v(l)%j(m)
2558	k = surf_usm_v(l)%k(m)
2559	temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
2560	ENDDO
2561	ENDDO
2562	ELSE
2563	DO m = 1, surf_usm_h%ns
2564	i = surf_usm_h%i(m)
2565	j = surf_usm_h%j(m)
2566	k = surf_usm_h%k(m)
2567	temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
2568	ENDDO
2569	DO l = 0, 3
2570	DO m = 1, surf_usm_v(l)%ns
2571	i = surf_usm_v(l)%i(m)
2572	j = surf_usm_v(l)%j(m)
2573	k = surf_usm_v(l)%k(m)
2574	temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
2575	ENDDO
2576	ENDDO
2577	ENDIF
2578
2579	CASE ( 'usm_rad_insw' )
2580	!-- array of sw radiation falling to surface after i-th reflection
2581	DO isurf = dirstart(ids), dirend(ids)
2582	IF ( surfl(id,isurf) == ids ) THEN
2583	IF ( av == 0 ) THEN
2584	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
2585	ELSE
2586	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
2587	ENDIF
2588	ENDIF
2589	ENDDO
2590
2591	CASE ( 'usm_rad_inlw' )
2592	!-- array of lw radiation falling to surface after i-th reflection
2593	DO isurf = dirstart(ids), dirend(ids)
2594	IF ( surfl(id,isurf) == ids ) THEN
2595	IF ( av == 0 ) THEN
2596	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
2597	ELSE
2598	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
2599	ENDIF
2600	ENDIF
2601	ENDDO
2602
2603	CASE ( 'usm_rad_inswdir' )
2604	!-- array of direct sw radiation falling to surface from sun
2605	DO isurf = dirstart(ids), dirend(ids)
2606	IF ( surfl(id,isurf) == ids ) THEN
2607	IF ( av == 0 ) THEN
2608	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
2609	ELSE
2610	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
2611	ENDIF
2612	ENDIF
2613	ENDDO
2614
2615	CASE ( 'usm_rad_inswdif' )
2616	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
2617	DO isurf = dirstart(ids), dirend(ids)
2618	IF ( surfl(id,isurf) == ids ) THEN
2619	IF ( av == 0 ) THEN
2620	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
2621	ELSE
2622	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
2623	ENDIF
2624	ENDIF
2625	ENDDO
2626
2627	CASE ( 'usm_rad_inswref' )
2628	!-- array of sw radiation falling to surface from reflections
2629	DO isurf = dirstart(ids), dirend(ids)
2630	IF ( surfl(id,isurf) == ids ) THEN
2631	IF ( av == 0 ) THEN
2632	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
2633	surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
2634	ELSE
2635	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
2636	ENDIF
2637	ENDIF
2638	ENDDO
2639
2640	CASE ( 'usm_rad_inlwref' )
2641	!-- array of lw radiation falling to surface from reflections
2642	DO isurf = dirstart(ids), dirend(ids)
2643	IF ( surfl(id,isurf) == ids ) THEN
2644	IF ( av == 0 ) THEN
2645	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
2646	ELSE
2647	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
2648	ENDIF
2649	ENDIF
2650	ENDDO
2651
2652	CASE ( 'usm_rad_outsw' )
2653	!-- array of sw radiation emitted from surface after i-th reflection
2654	DO isurf = dirstart(ids), dirend(ids)
2655	IF ( surfl(id,isurf) == ids ) THEN
2656	IF ( av == 0 ) THEN
2657	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
2658	ELSE
2659	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
2660	ENDIF
2661	ENDIF
2662	ENDDO
2663
2664	CASE ( 'usm_rad_outlw' )
2665	!-- array of lw radiation emitted from surface after i-th reflection
2666	DO isurf = dirstart(ids), dirend(ids)
2667	IF ( surfl(id,isurf) == ids ) THEN
2668	IF ( av == 0 ) THEN
2669	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
2670	ELSE
2671	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
2672	ENDIF
2673	ENDIF
2674	ENDDO
2675
2676	CASE ( 'usm_rad_ressw' )
2677	!-- average of array of residua of sw radiation absorbed in surface after last reflection
2678	DO isurf = dirstart(ids), dirend(ids)
2679	IF ( surfl(id,isurf) == ids ) THEN
2680	IF ( av == 0 ) THEN
2681	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
2682	ELSE
2683	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
2684	ENDIF
2685	ENDIF
2686	ENDDO
2687
2688	CASE ( 'usm_rad_reslw' )
2689	!-- average of array of residua of lw radiation absorbed in surface after last reflection
2690	DO isurf = dirstart(ids), dirend(ids)
2691	IF ( surfl(id,isurf) == ids ) THEN
2692	IF ( av == 0 ) THEN
2693	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
2694	ELSE
2695	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
2696	ENDIF
2697	ENDIF
2698	ENDDO
2699
2700	CASE ( 'usm_rad_hf' )
2701	!-- array of heat flux from radiation for surfaces after all reflections
2702	IF ( av == 0 ) THEN
2703	DO m = 1, surf_usm_h%ns
2704	i = surf_usm_h%i(m)
2705	j = surf_usm_h%j(m)
2706	k = surf_usm_h%k(m)
2707	temp_pf(k,j,i) = surf_usm_h%surfhf(m)
2708	ENDDO
2709	DO l = 0, 3
2710	DO m = 1, surf_usm_v(l)%ns
2711	i = surf_usm_v(l)%i(m)
2712	j = surf_usm_v(l)%j(m)
2713	k = surf_usm_v(l)%k(m)
2714	temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
2715	ENDDO
2716	ENDDO
2717	ELSE
2718	DO m = 1, surf_usm_h%ns
2719	i = surf_usm_h%i(m)
2720	j = surf_usm_h%j(m)
2721	k = surf_usm_h%k(m)
2722	temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
2723	ENDDO
2724	DO l = 0, 3
2725	DO m = 1, surf_usm_v(l)%ns
2726	i = surf_usm_v(l)%i(m)
2727	j = surf_usm_v(l)%j(m)
2728	k = surf_usm_v(l)%k(m)
2729	temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
2730	ENDDO
2731	ENDDO
2732	ENDIF
2733
2734	CASE ( 'usm_wshf' )
2735	!-- array of sensible heat flux from surfaces
2736	IF ( av == 0 ) THEN
2737	DO m = 1, surf_usm_h%ns
2738	i = surf_usm_h%i(m)
2739	j = surf_usm_h%j(m)
2740	k = surf_usm_h%k(m)
2741	temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2742	ENDDO
2743	DO l = 0, 3
2744	DO m = 1, surf_usm_v(l)%ns
2745	i = surf_usm_v(l)%i(m)
2746	j = surf_usm_v(l)%j(m)
2747	k = surf_usm_v(l)%k(m)
2748	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2749	ENDDO
2750	ENDDO
2751	ELSE
2752	DO m = 1, surf_usm_h%ns
2753	i = surf_usm_h%i(m)
2754	j = surf_usm_h%j(m)
2755	k = surf_usm_h%k(m)
2756	temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2757	ENDDO
2758	DO l = 0, 3
2759	DO m = 1, surf_usm_v(l)%ns
2760	i = surf_usm_v(l)%i(m)
2761	j = surf_usm_v(l)%j(m)
2762	k = surf_usm_v(l)%k(m)
2763	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2764	ENDDO
2765	ENDDO
2766	ENDIF
2767
2768
2769	CASE ( 'usm_wghf' )
2770	!-- array of heat flux from ground (land, wall, roof)
2771	IF ( av == 0 ) THEN
2772	DO m = 1, surf_usm_h%ns
2773	i = surf_usm_h%i(m)
2774	j = surf_usm_h%j(m)
2775	k = surf_usm_h%k(m)
2776	temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2777	ENDDO
2778	DO l = 0, 3
2779	DO m = 1, surf_usm_v(l)%ns
2780	i = surf_usm_v(l)%i(m)
2781	j = surf_usm_v(l)%j(m)
2782	k = surf_usm_v(l)%k(m)
2783	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2784	ENDDO
2785	ENDDO
2786	ELSE
2787	DO m = 1, surf_usm_h%ns
2788	i = surf_usm_h%i(m)
2789	j = surf_usm_h%j(m)
2790	k = surf_usm_h%k(m)
2791	temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2792	ENDDO
2793	DO l = 0, 3
2794	DO m = 1, surf_usm_v(l)%ns
2795	i = surf_usm_v(l)%i(m)
2796	j = surf_usm_v(l)%j(m)
2797	k = surf_usm_v(l)%k(m)
2798	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2799	ENDDO
2800	ENDDO
2801	ENDIF
2802
2803	CASE ( 'usm_t_surf' )
2804	!-- surface temperature for surfaces
2805	IF ( av == 0 ) THEN
2806	DO m = 1, surf_usm_h%ns
2807	i = surf_usm_h%i(m)
2808	j = surf_usm_h%j(m)
2809	k = surf_usm_h%k(m)
2810	temp_pf(k,j,i) = t_surf_h(m)
2811	ENDDO
2812	DO l = 0, 3
2813	DO m = 1, surf_usm_v(l)%ns
2814	i = surf_usm_v(l)%i(m)
2815	j = surf_usm_v(l)%j(m)
2816	k = surf_usm_v(l)%k(m)
2817	temp_pf(k,j,i) = t_surf_v(l)%t(m)
2818	ENDDO
2819	ENDDO
2820	ELSE
2821	DO m = 1, surf_usm_h%ns
2822	i = surf_usm_h%i(m)
2823	j = surf_usm_h%j(m)
2824	k = surf_usm_h%k(m)
2825	temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
2826	ENDDO
2827	DO l = 0, 3
2828	DO m = 1, surf_usm_v(l)%ns
2829	i = surf_usm_v(l)%i(m)
2830	j = surf_usm_v(l)%j(m)
2831	k = surf_usm_v(l)%k(m)
2832	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
2833	ENDDO
2834	ENDDO
2835	ENDIF
2836
2837	CASE ( 'usm_t_wall' )
2838	!-- wall temperature for iwl layer of walls and land
2839	IF ( av == 0 ) THEN
2840	DO m = 1, surf_usm_h%ns
2841	i = surf_usm_h%i(m)
2842	j = surf_usm_h%j(m)
2843	k = surf_usm_h%k(m)
2844	temp_pf(k,j,i) = t_wall_h(iwl,m)
2845	ENDDO
2846	DO l = 0, 3
2847	DO m = 1, surf_usm_v(l)%ns
2848	i = surf_usm_v(l)%i(m)
2849	j = surf_usm_v(l)%j(m)
2850	k = surf_usm_v(l)%k(m)
2851	temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
2852	ENDDO
2853	ENDDO
2854	ELSE
2855	DO m = 1, surf_usm_h%ns
2856	i = surf_usm_h%i(m)
2857	j = surf_usm_h%j(m)
2858	k = surf_usm_h%k(m)
2859	temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
2860	ENDDO
2861	DO l = 0, 3
2862	DO m = 1, surf_usm_v(l)%ns
2863	i = surf_usm_v(l)%i(m)
2864	j = surf_usm_v(l)%j(m)
2865	k = surf_usm_v(l)%k(m)
2866	temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
2867	ENDDO
2868	ENDDO
2869	ENDIF
2870
2871	CASE DEFAULT
2872	found = .FALSE.
2873
2874	END SELECT
2875
2876	!-- fill out array local_pf which is subsequently treated by data_output_3d
2877	CALL exchange_horiz( temp_pf, nbgp )
2878	!
2879	!-- To Do: why reversed loop order
2880	DO j = nysg,nyng
2881	DO i = nxlg,nxrg
2882	DO k = nzb_do, nzt_do
2883	local_pf(i,j,k) = temp_pf(k,j,i)
2884	ENDDO
2885	ENDDO
2886	ENDDO
2887
2888	END SUBROUTINE usm_data_output_3d
2889
2890
2891	!------------------------------------------------------------------------------!
2892	!
2893	! Description:
2894	! ------------
2895	!> Soubroutine defines appropriate grid for netcdf variables.
2896	!> It is called out from subroutine netcdf.
2897	!------------------------------------------------------------------------------!
2898	SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
2899
2900	IMPLICIT NONE
2901
2902	CHARACTER (len=*), INTENT(IN) :: variable !<
2903	LOGICAL, INTENT(OUT) :: found !<
2904	CHARACTER (len=*), INTENT(OUT) :: grid_x !<
2905	CHARACTER (len=*), INTENT(OUT) :: grid_y !<
2906	CHARACTER (len=*), INTENT(OUT) :: grid_z !<
2907
2908	CHARACTER (len=varnamelength) :: var
2909
2910	var = TRIM(variable)
2911	IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. &
2912	var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. &
2913	var(1:16) == 'usm_rad_inswdif_' .OR. var(1:16) == 'usm_rad_inswref_' .OR. &
2914	var(1:16) == 'usm_rad_inlwdif_' .OR. var(1:16) == 'usm_rad_inlwref_' .OR. &
2915	var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. &
2916	var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. &
2917	var(1:11) == 'usm_rad_hf_' .OR. &
2918	var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
2919	var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. &
2920	var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. &
2921	var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. &
2922	var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis' ) THEN
2923
2924	found = .TRUE.
2925	grid_x = 'x'
2926	grid_y = 'y'
2927	grid_z = 'zu'
2928	ELSE
2929	found = .FALSE.
2930	grid_x = 'none'
2931	grid_y = 'none'
2932	grid_z = 'none'
2933	ENDIF
2934
2935	END SUBROUTINE usm_define_netcdf_grid
2936
2937
2938	!------------------------------------------------------------------------------!
2939	!> Finds first model boundary crossed by a ray
2940	!------------------------------------------------------------------------------!
2941	PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross)
2942
2943	IMPLICIT NONE
2944
2945	INTEGER(iwp) :: d !<
2946	INTEGER(iwp) :: seldim !< found fist crossing index
2947
2948	INTEGER(iwp), DIMENSION(3) :: bdyd !< boundary direction
2949	INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross !< found boundary crossing (d, z, y, x)
2950
2951	REAL(wp) :: bdydim !<
2952	REAL(wp) :: dist !<
2953
2954	REAL(wp), DIMENSION(3) :: crossdist !< crossing distance
2955	REAL(wp), DIMENSION(3), INTENT(in) :: origin !< ray origin
2956	REAL(wp), DIMENSION(3), INTENT(in) :: uvect !< ray unit vector
2957
2958
2959	bdydim = nzut + .5_wp !< top boundary
2960	bdyd(1) = isky
2961	crossdist(1) = ( bdydim - origin(1) ) / uvect(1) !< subroutine called only when uvect(1)>0
2962
2963	IF ( uvect(2) == 0._wp ) THEN
2964	crossdist(2) = huge(1._wp)
2965	ELSE
2966	IF ( uvect(2) >= 0._wp ) THEN
2967	bdydim = ny + .5_wp !< north global boundary
2968	bdyd(2) = inorthb
2969	ELSE
2970	bdydim = -.5_wp !< south global boundary
2971	bdyd(2) = isouthb
2972	ENDIF
2973	crossdist(2) = ( bdydim - origin(2) ) / uvect(2)
2974	ENDIF
2975
2976	IF ( uvect(3) == 0._wp ) THEN
2977	crossdist(3) = huge(1._wp)
2978	ELSE
2979	IF ( uvect(3) >= 0._wp ) THEN
2980	bdydim = nx + .5_wp !< east global boundary
2981	bdyd(3) = ieastb
2982	ELSE
2983	bdydim = -.5_wp !< west global boundary
2984	bdyd(3) = iwestb
2985	ENDIF
2986	crossdist(3) = ( bdydim - origin(3) ) / uvect(3)
2987	ENDIF
2988
2989	seldim = minloc(crossdist, 1)
2990	dist = crossdist(seldim)
2991	d = bdyd(seldim)
2992
2993	bdycross(1) = d
2994	bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &
2995	+ .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
2996
2997	END SUBROUTINE
2998
2999
3000	!------------------------------------------------------------------------------!
3001	!> Determines whether two faces are oriented towards each other
3002	!------------------------------------------------------------------------------!
3003	PURE LOGICAL FUNCTION usm_facing(x, y, z, d, x2, y2, z2, d2)
3004	IMPLICIT NONE
3005	INTEGER(iwp), INTENT(in) :: x, y, z, d, x2, y2, z2, d2
3006
3007	usm_facing = .FALSE.
3008	IF ( d==iroof .AND. d2==iroof ) RETURN
3009	IF ( d==isky .AND. d2==isky ) RETURN
3010	IF ( (d==isouth .OR. d==inorthb) .AND. (d2==isouth.OR.d2==inorthb) ) RETURN
3011	IF ( (d==inorth .OR. d==isouthb) .AND. (d2==inorth.OR.d2==isouthb) ) RETURN
3012	IF ( (d==iwest .OR. d==ieastb) .AND. (d2==iwest.OR.d2==ieastb) ) RETURN
3013	IF ( (d==ieast .OR. d==iwestb) .AND. (d2==ieast.OR.d2==iwestb) ) RETURN
3014
3015	SELECT CASE (d)
3016	CASE (iroof) !< ground, roof
3017	IF ( z2 < z ) RETURN
3018	CASE (isky) !< sky
3019	IF ( z2 > z ) RETURN
3020	CASE (isouth, inorthb) !< south facing
3021	IF ( y2 > y ) RETURN
3022	CASE (inorth, isouthb) !< north facing
3023	IF ( y2 < y ) RETURN
3024	CASE (iwest, ieastb) !< west facing
3025	IF ( x2 > x ) RETURN
3026	CASE (ieast, iwestb) !< east facing
3027	IF ( x2 < x ) RETURN
3028	END SELECT
3029
3030	SELECT CASE (d2)
3031	CASE (iroof) !< ground, roof
3032	IF ( z < z2 ) RETURN
3033	CASE (isky) !< sky
3034	IF ( z > z2 ) RETURN
3035	CASE (isouth, inorthb) !< south facing
3036	IF ( y > y2 ) RETURN
3037	CASE (inorth, isouthb) !< north facing
3038	IF ( y < y2 ) RETURN
3039	CASE (iwest, ieastb) !< west facing
3040	IF ( x > x2 ) RETURN
3041	CASE (ieast, iwestb) !< east facing
3042	IF ( x < x2 ) RETURN
3043	CASE (-1)
3044	CONTINUE
3045	END SELECT
3046
3047	usm_facing = .TRUE.
3048
3049	END FUNCTION usm_facing
3050
3051
3052	!------------------------------------------------------------------------------!
3053	! Description:
3054	! ------------
3055	!> Initialization of the wall surface model
3056	!------------------------------------------------------------------------------!
3057	SUBROUTINE usm_init_material_model
3058
3059	IMPLICIT NONE
3060
3061	INTEGER(iwp) :: k, l, m !< running indices
3062
3063	CALL location_message( ' initialization of wall surface model', .TRUE. )
3064
3065	!-- Calculate wall grid spacings.
3066	!-- Temperature is defined at the center of the wall layers,
3067	!-- whereas gradients/fluxes are defined at the edges (_stag)
3068	DO k = nzb_wall, nzt_wall
3069	zwn(k) = zwn_default(k)
3070	ENDDO
3071	!
3072	!-- apply for all particular surface grids. First for horizontal surfaces
3073	DO m = 1, surf_usm_h%ns
3074	surf_usm_h%zw(:,m) = zwn(:) * &
3075	surf_usm_h%thickness_wall(m)
3076	surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
3077	DO k = nzb_wall+1, nzt_wall
3078	surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) - &
3079	surf_usm_h%zw(k-1,m)
3080	ENDDO
3081
3082	surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
3083
3084	DO k = nzb_wall, nzt_wall-1
3085	surf_usm_h%dz_wall_stag(k,m) = 0.5 * ( &
3086	surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
3087	ENDDO
3088	surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
3089	ENDDO
3090	surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall
3091	surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag
3092	!
3093	!-- For vertical surfaces
3094	DO l = 0, 3
3095	DO m = 1, surf_usm_v(l)%ns
3096	surf_usm_v(l)%zw(:,m) = zwn(:) * &
3097	surf_usm_v(l)%thickness_wall(m)
3098	surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
3099	DO k = nzb_wall+1, nzt_wall
3100	surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - &
3101	surf_usm_v(l)%zw(k-1,m)
3102	ENDDO
3103
3104	surf_usm_v(l)%dz_wall(nzt_wall+1,m) = &
3105	surf_usm_v(l)%dz_wall(nzt_wall,m)
3106
3107	DO k = nzb_wall, nzt_wall-1
3108	surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * ( &
3109	surf_usm_v(l)%dz_wall(k+1,m) + &
3110	surf_usm_v(l)%dz_wall(k,m) )
3111	ENDDO
3112	surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = &
3113	surf_usm_v(l)%dz_wall(nzt_wall,m)
3114	ENDDO
3115	surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall
3116	surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag
3117	ENDDO
3118
3119
3120	CALL location_message( ' wall structures filed out', .TRUE. )
3121
3122	CALL location_message( ' initialization of wall surface model finished', .TRUE. )
3123
3124	END SUBROUTINE usm_init_material_model
3125
3126
3127	!------------------------------------------------------------------------------!
3128	! Description:
3129	! ------------
3130	!> Initialization of the urban surface model
3131	!------------------------------------------------------------------------------!
3132	SUBROUTINE usm_init_urban_surface
3133
3134	IMPLICIT NONE
3135
3136	INTEGER(iwp) :: i, j, k, l, m !< running indices
3137	REAL(wp) :: c, d, tin, exn
3138
3139
3140	CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
3141	!-- surface forcing have to be disabled for LSF
3142	!-- in case of enabled urban surface module
3143	IF ( large_scale_forcing ) THEN
3144	lsf_surf = .FALSE.
3145	ENDIF
3146
3147	!-- init anthropogenic sources of heat
3148	CALL usm_allocate_urban_surface()
3149
3150	!-- read the surface_types array somewhere
3151	CALL usm_read_urban_surface_types()
3152
3153	!-- init material heat model
3154	CALL usm_init_material_model()
3155
3156	IF ( usm_anthropogenic_heat ) THEN
3157	!-- init anthropogenic sources of heat (from transportation for now)
3158	CALL usm_read_anthropogenic_heat()
3159	ENDIF
3160
3161	IF ( read_svf_on_init ) THEN
3162	!-- read svf, csf, svfsurf and csfsurf data from file
3163	CALL location_message( ' Start reading SVF from file', .TRUE. )
3164	CALL usm_read_svf_from_file()
3165	CALL location_message( ' Reading SVF from file has finished', .TRUE. )
3166	ELSE
3167	!-- calculate SFV and CSF
3168	CALL location_message( ' Start calculation of SVF', .TRUE. )
3169	CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'start' )
3170	CALL usm_calc_svf()
3171	CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'stop' )
3172	CALL location_message( ' Calculation of SVF has finished', .TRUE. )
3173	ENDIF
3174
3175	IF ( write_svf_on_init ) THEN
3176	!-- write svf, csf svfsurf and csfsurf data to file
3177	CALL location_message( ' Store SVF and CSF to file', .TRUE. )
3178	CALL usm_write_svf_to_file()
3179	ENDIF
3180
3181	IF ( plant_canopy ) THEN
3182	!-- gridpcbl was only necessary for initialization
3183	DEALLOCATE( gridpcbl )
3184	IF ( .NOT. ALLOCATED(pc_heating_rate) ) THEN
3185	!-- then pc_heating_rate is allocated in init_plant_canopy
3186	!-- in case of cthf /= 0 => we need to allocate it for our use here
3187	ALLOCATE( pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3188	ENDIF
3189	ENDIF
3190
3191	!-- Intitialization of the surface and wall/ground/roof temperature
3192
3193	!-- Initialization for restart runs
3194	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
3195
3196	!-- restore data from restart file
3197	CALL usm_read_restart_data()
3198	ELSE
3199
3200	!-- Calculate initial surface temperature from pt of adjacent gridbox
3201	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
3202
3203	!
3204	!-- At horizontal surfaces. Please note, t_surf_h is defined on a
3205	!-- different data type, but with the same dimension.
3206	DO m = 1, surf_usm_h%ns
3207	i = surf_usm_h%i(m)
3208	j = surf_usm_h%j(m)
3209	k = surf_usm_h%k(m)
3210
3211	t_surf_h(m) = pt(k,j,i) * exn
3212	ENDDO
3213	!
3214	!-- At vertical surfaces.
3215	DO l = 0, 3
3216	DO m = 1, surf_usm_v(l)%ns
3217	i = surf_usm_v(l)%i(m)
3218	j = surf_usm_v(l)%j(m)
3219	k = surf_usm_v(l)%k(m)
3220
3221	t_surf_v(l)%t(m) = pt(k,j,i) * exn
3222	ENDDO
3223	ENDDO
3224
3225
3226	!-- initial values for t_wall
3227	!-- outer value is set to surface temperature
3228	!-- inner value is set to wall_inner_temperature
3229	!-- and profile is logaritmic (linear in nz).
3230	!-- Horizontal surfaces
3231	DO m = 1, surf_usm_h%ns
3232	!
3233	!-- Roof
3234	IF ( surf_usm_h%isroof_surf(m) ) THEN
3235	tin = roof_inner_temperature
3236	!
3237	!-- Normal land surface
3238	ELSE
3239	tin = soil_inner_temperature
3240	ENDIF
3241
3242	DO k = nzb_wall, nzt_wall+1
3243	c = REAL( k - nzb_wall, wp ) / &
3244	REAL( nzt_wall + 1 - nzb_wall , wp )
3245
3246	t_wall_h(k,m) = ( 1.0_wp - c ) * t_surf_h(m) + c * tin
3247	ENDDO
3248	ENDDO
3249	!
3250	!-- Vertical surfaces
3251	DO l = 0, 3
3252	DO m = 1, surf_usm_v(l)%ns
3253	!
3254	!-- Inner wall
3255	tin = wall_inner_temperature
3256
3257	DO k = nzb_wall, nzt_wall+1
3258	c = REAL( k - nzb_wall, wp ) / &
3259	REAL( nzt_wall + 1 - nzb_wall , wp )
3260
3261	t_wall_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) + &
3262	c * tin
3263	ENDDO
3264	ENDDO
3265	ENDDO
3266
3267	ENDIF
3268
3269	!--
3270	!-- Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
3271	CALL user_init_urban_surface
3272
3273	!-- initialize prognostic values for the first timestep
3274	t_surf_h_p = t_surf_h
3275	t_surf_v_p = t_surf_v
3276
3277	t_wall_h_p = t_wall_h
3278	t_wall_v_p = t_wall_v
3279
3280	!-- Adjust radiative fluxes for urban surface at model start
3281	CALL usm_radiation
3282
3283	CALL cpu_log( log_point_s(78), 'usm_init', 'stop' )
3284
3285
3286	END SUBROUTINE usm_init_urban_surface
3287
3288
3289	!------------------------------------------------------------------------------!
3290	! Description:
3291	! ------------
3292	!
3293	!> Wall model as part of the urban surface model. The model predicts wall
3294	!> temperature.
3295	!------------------------------------------------------------------------------!
3296	SUBROUTINE usm_material_heat_model
3297
3298
3299	IMPLICIT NONE
3300
3301	INTEGER(iwp) :: i,j,k,l,kw, m !< running indices
3302
3303	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend !< tendency
3304
3305	!
3306	!-- For horizontal surfaces
3307	DO m = 1, surf_usm_h%ns
3308	!
3309	!-- Obtain indices
3310	i = surf_usm_h%i(m)
3311	j = surf_usm_h%j(m)
3312	k = surf_usm_h%k(m)
3313	!
3314	!-- prognostic equation for ground/roof temperature t_wall_h
3315	wtend(:) = 0.0_wp
3316	wtend(nzb_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m)) * &
3317	( surf_usm_h%lambda_h(nzb_wall,m) * &
3318	( t_wall_h(nzb_wall+1,m) &
3319	- t_wall_h(nzb_wall,m) ) * &
3320	surf_usm_h%ddz_wall(nzb_wall+1,m) &
3321	+ surf_usm_h%wghf_eb(m) ) * &
3322	surf_usm_h%ddz_wall_stag(nzb_wall,m)
3323
3324	DO kw = nzb_wall+1, nzt_wall
3325	wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m)) &
3326	* ( surf_usm_h%lambda_h(kw,m) &
3327	* ( t_wall_h(kw+1,m) - t_wall_h(kw,m) ) &
3328	* surf_usm_h%ddz_wall(kw+1,m) &
3329	- surf_usm_h%lambda_h(kw-1,m) &
3330	* ( t_wall_h(kw,m) - t_wall_h(kw-1,m) ) &
3331	* surf_usm_h%ddz_wall(kw,m) &
3332	) * surf_usm_h%ddz_wall_stag(kw,m)
3333	ENDDO
3334
3335	t_wall_h_p(nzb_wall:nzt_wall,m) = t_wall_h(nzb_wall:nzt_wall,m) &
3336	+ dt_3d * ( tsc(2) &
3337	* wtend(nzb_wall:nzt_wall) + tsc(3) &
3338	* surf_usm_h%tt_wall_m(nzb_wall:nzt_wall,m) )
3339
3340	!
3341	!-- calculate t_wall tendencies for the next Runge-Kutta step
3342	IF ( timestep_scheme(1:5) == 'runge' ) THEN
3343	IF ( intermediate_timestep_count == 1 ) THEN
3344	DO kw = nzb_wall, nzt_wall
3345	surf_usm_h%tt_wall_m(kw,m) = wtend(kw)
3346	ENDDO
3347	ELSEIF ( intermediate_timestep_count < &
3348	intermediate_timestep_count_max ) THEN
3349	DO kw = nzb_wall, nzt_wall
3350	surf_usm_h%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) + &
3351	5.3125_wp * surf_usm_h%tt_wall_m(kw,m)
3352	ENDDO
3353	ENDIF
3354	ENDIF
3355	ENDDO
3356	!
3357	!-- For vertical surfaces
3358	DO l = 0, 3
3359	DO m = 1, surf_usm_v(l)%ns
3360	!
3361	!-- Obtain indices
3362	i = surf_usm_v(l)%i(m)
3363	j = surf_usm_v(l)%j(m)
3364	k = surf_usm_v(l)%k(m)
3365	!
3366	!-- prognostic equation for wall temperature t_wall_v
3367	wtend(:) = 0.0_wp
3368	wtend(nzb_wall) = (1.0_wp / surf_usm_v(l)%rho_c_wall(nzb_wall,m)) * &
3369	( surf_usm_v(l)%lambda_h(nzb_wall,m) * &
3370	( t_wall_v(l)%t(nzb_wall+1,m) &
3371	- t_wall_v(l)%t(nzb_wall,m) ) * &
3372	surf_usm_v(l)%ddz_wall(nzb_wall+1,m) &
3373	+ surf_usm_v(l)%wghf_eb(m) ) * &
3374	surf_usm_v(l)%ddz_wall_stag(nzb_wall,m)
3375
3376	DO kw = nzb_wall+1, nzt_wall
3377	wtend(kw) = (1.0_wp / surf_usm_v(l)%rho_c_wall(kw,m)) &
3378	* ( surf_usm_v(l)%lambda_h(kw,m) &
3379	* ( t_wall_v(l)%t(kw+1,m) - t_wall_v(l)%t(kw,m) )&
3380	* surf_usm_v(l)%ddz_wall(kw+1,m) &
3381	- surf_usm_v(l)%lambda_h(kw-1,m) &
3382	* ( t_wall_v(l)%t(kw,m) - t_wall_v(l)%t(kw-1,m) )&
3383	* surf_usm_v(l)%ddz_wall(kw,m) &
3384	) * surf_usm_v(l)%ddz_wall_stag(kw,m)
3385	ENDDO
3386
3387	t_wall_v_p(l)%t(nzb_wall:nzt_wall,m) = &
3388	t_wall_v(l)%t(nzb_wall:nzt_wall,m) &
3389	+ dt_3d * ( tsc(2) &
3390	* wtend(nzb_wall:nzt_wall) + tsc(3) &
3391	* surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall,m) )
3392
3393	!
3394	!-- calculate t_wall tendencies for the next Runge-Kutta step
3395	IF ( timestep_scheme(1:5) == 'runge' ) THEN
3396	IF ( intermediate_timestep_count == 1 ) THEN
3397	DO kw = nzb_wall, nzt_wall
3398	surf_usm_v(l)%tt_wall_m(kw,m) = wtend(kw)
3399	ENDDO
3400	ELSEIF ( intermediate_timestep_count < &
3401	intermediate_timestep_count_max ) THEN
3402	DO kw = nzb_wall, nzt_wall
3403	surf_usm_v(l)%tt_wall_m(kw,m) = &
3404	- 9.5625_wp * wtend(kw) + &
3405	5.3125_wp * surf_usm_v(l)%tt_wall_m(kw,m)
3406	ENDDO
3407	ENDIF
3408	ENDIF
3409	ENDDO
3410	ENDDO
3411
3412	END SUBROUTINE usm_material_heat_model
3413
3414
3415	!------------------------------------------------------------------------------!
3416	! Description:
3417	! ------------
3418	!> Parin for &usm_par for urban surface model
3419	!------------------------------------------------------------------------------!
3420	SUBROUTINE usm_parin
3421
3422	IMPLICIT NONE
3423
3424	CHARACTER (LEN=80) :: line !< string containing current line of file PARIN
3425
3426	NAMELIST /urban_surface_par/ &
3427	land_category, &
3428	mrt_factors, &
3429	nrefsteps, &
3430	pedestrant_category, &
3431	ra_horiz_coef, &
3432	read_svf_on_init, &
3433	roof_category, &
3434	split_diffusion_radiation, &
3435	urban_surface, &
3436	usm_anthropogenic_heat, &
3437	usm_energy_balance_land, &
3438	usm_energy_balance_wall, &
3439	usm_material_model, &
3440	usm_lad_rma, &
3441	wall_category, &
3442	write_svf_on_init
3443
3444	line = ' '
3445
3446	!
3447	!-- Try to find urban surface model package
3448	REWIND ( 11 )
3449	line = ' '
3450	DO WHILE ( INDEX( line, '&urban_surface_par' ) == 0 )
3451	READ ( 11, '(A)', END=10 ) line
3452	ENDDO
3453	BACKSPACE ( 11 )
3454
3455	!
3456	!-- Read user-defined namelist
3457	READ ( 11, urban_surface_par )
3458	!
3459	!-- Set flag that indicates that the land surface model is switched on
3460	urban_surface = .TRUE.
3461
3462	10 CONTINUE
3463
3464	END SUBROUTINE usm_parin
3465
3466
3467	!------------------------------------------------------------------------------!
3468	! Description:
3469	! ------------
3470	!> This subroutine calculates interaction of the solar radiation
3471	!> with urban surface and updates surface, roofs and walls heatfluxes.
3472	!> It also updates rad_sw_out and rad_lw_out.
3473	!------------------------------------------------------------------------------!
3474	SUBROUTINE usm_radiation
3475
3476	IMPLICIT NONE
3477
3478	INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
3479	INTEGER(iwp) :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
3480	INTEGER(iwp), DIMENSION(4) :: bdycross
3481	REAL(wp), DIMENSION(3,3) :: mrot !< grid rotation matrix (xyz)
3482	REAL(wp), DIMENSION(3,0:9) :: vnorm !< face direction normal vectors (xyz)
3483	REAL(wp), DIMENSION(3) :: sunorig !< grid rotated solar direction unit vector (xyz)
3484	REAL(wp), DIMENSION(3) :: sunorig_grid !< grid squashed solar direction unit vector (zyx)
3485	REAL(wp), DIMENSION(0:9) :: costheta !< direct irradiance factor of solar angle
3486	REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temp tendency
3487	REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove)
3488	REAL(wp) :: rx, ry, rz
3489	REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff
3490	INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy
3491	INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /)
3492
3493
3494	IF ( plant_canopy ) THEN
3495	pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
3496	/ (cp * hyp(nzub:nzut) * dxdydz) !< equals to 1 / (rho * c_p * Vbox * T)
3497	ENDIF
3498
3499	sun_direction = .TRUE.
3500	CALL calc_zenith !< required also for diffusion radiation
3501
3502	!-- prepare rotated normal vectors and irradiance factor
3503	vnorm(1,:) = idir(:)
3504	vnorm(2,:) = jdir(:)
3505	vnorm(3,:) = kdir(:)
3506	mrot(1, :) = (/ cos(alpha), -sin(alpha), 0._wp /)
3507	mrot(2, :) = (/ sin(alpha), cos(alpha), 0._wp /)
3508	mrot(3, :) = (/ 0._wp, 0._wp, 1._wp /)
3509	sunorig = (/ sun_dir_lon, sun_dir_lat, zenith(0) /)
3510	sunorig = matmul(mrot, sunorig)
3511	DO d = 0, 9
3512	costheta(d) = dot_product(sunorig, vnorm(:,d))
3513	ENDDO
3514
3515	IF ( zenith(0) > 0 ) THEN
3516	!-- now we will "squash" the sunorig vector by grid box size in
3517	!-- each dimension, so that this new direction vector will allow us
3518	!-- to traverse the ray path within grid coordinates directly
3519	sunorig_grid = (/ sunorig(3)/dz, sunorig(2)/dy, sunorig(1)/dx /)
3520	!-- sunorig_grid = sunorig_grid / norm2(sunorig_grid)
3521	sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
3522
3523	IF ( plant_canopy ) THEN
3524	!-- precompute effective box depth with prototype Leaf Area Density
3525	pc_box_dimshift = maxloc(sunorig, 1) - 1
3526	CALL usm_box_absorb(cshift((/dx,dy,dz/), pc_box_dimshift), &
3527	60, prototype_lad, &
3528	cshift(sunorig, pc_box_dimshift), &
3529	pc_box_area, pc_abs_frac)
3530	pc_box_area = pc_box_area * sunorig(pc_box_dimshift+1) / sunorig(3)
3531	pc_abs_eff = log(1._wp - pc_abs_frac) / prototype_lad
3532	ENDIF
3533	ENDIF
3534
3535	!-- split diffusion and direct part of the solar downward radiation
3536	!-- comming from radiation model and store it in 2D arrays
3537	!-- rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff
3538	IF ( split_diffusion_radiation ) THEN
3539	CALL usm_calc_diffusion_radiation
3540	ELSE
3541	rad_sw_in_diff = 0.0_wp
3542	rad_sw_in_dir(:,:) = rad_sw_in(0,:,:)
3543	rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
3544	ENDIF
3545
3546	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3547	!-- First pass: direct + diffuse irradiance
3548	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3549	surfinswdir = 0._wp !nsurfl
3550	surfinswdif = 0._wp !nsurfl
3551	surfinlwdif = 0._wp !nsurfl
3552	surfins = 0._wp !nsurfl
3553	surfinl = 0._wp !nsurfl
3554	surfoutsl(:) = 0.0_wp !start-end
3555	surfoutll(:) = 0.0_wp !start-end
3556
3557	!-- Set up thermal radiation from surfaces
3558	!-- emiss_surf is defined only for surfaces for which energy balance is calculated
3559	!-- Workaround: reorder surface data type back on 1D array including all surfaces,
3560	!-- which implies to reorder horizontal and vertical surfaces
3561	!
3562	!-- Horizontal walls
3563	mm = 1
3564	DO i = nxl, nxr
3565	DO j = nys, nyn
3566
3567	DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
3568	surfoutll(mm) = surf_usm_h%emiss_surf(m) * sigma_sb &
3569	* t_surf_h(m)**4
3570	albedo_surf(mm) = surf_usm_h%albedo_surf(m)
3571	emiss_surf(mm) = surf_usm_h%emiss_surf(m)
3572	mm = mm + 1
3573	ENDDO
3574	ENDDO
3575	ENDDO
3576	!
3577	!-- Vertical walls
3578	DO i = nxl, nxr
3579	DO j = nys, nyn
3580	DO ll = 0, 3
3581	l = reorder(ll)
3582	DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
3583	surfoutll(mm) = surf_usm_v(l)%emiss_surf(m) * sigma_sb &
3584	* t_surf_v(l)%t(m)**4
3585	albedo_surf(mm) = surf_usm_v(l)%albedo_surf(m)
3586	emiss_surf(mm) = surf_usm_v(l)%emiss_surf(m)
3587	mm = mm + 1
3588	ENDDO
3589	ENDDO
3590	ENDDO
3591	ENDDO
3592
3593	#if defined( __parallel )
3594	!-- might be optimized and gather only values relevant for current processor
3595
3596	CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, &
3597	surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
3598	#else
3599	surfoutl(:) = surfoutll(:) !nsurf global
3600	#endif
3601
3602	isurf1 = -1 !< previous processed surface
3603	DO isvf = 1, nsvfl
3604	isurf = svfsurf(1, isvf)
3605	k = surfl(iz, isurf)
3606	j = surfl(iy, isurf)
3607	i = surfl(ix, isurf)
3608	isurfsrc = svfsurf(2, isvf)
3609	IF ( zenith(0) > 0 .AND. isurf /= isurf1 ) THEN
3610	!-- locate the virtual surface where the direct solar ray crosses domain boundary
3611	!-- (once per target surface)
3612	d = surfl(id, isurf)
3613	rz = REAL(k, wp) - 0.5_wp * kdir(d)
3614	ry = REAL(j, wp) - 0.5_wp * jdir(d)
3615	rx = REAL(i, wp) - 0.5_wp * idir(d)
3616
3617	CALL usm_find_boundary_face( (/ rz, ry, rx /), sunorig_grid, bdycross)
3618
3619	isurf1 = isurf
3620	ENDIF
3621
3622	IF ( surf(id, isurfsrc) >= isky ) THEN
3623	!-- diffuse rad from boundary surfaces. Since it is a simply
3624	!-- calculated value, it is not assigned to surfref(s/l),
3625	!-- instead it is used directly here
3626	!-- we consider the radiation from the radiation model falling on surface
3627	!-- as the radiation falling on the top of urban layer into the place of the source surface
3628	!-- we consider it as a very reasonable simplification which allow as avoid
3629	!-- necessity of other global range arrays and some all to all mpi communication
3630	surfinswdif(isurf) = surfinswdif(isurf) + rad_sw_in_diff(j,i) * svf(1,isvf) * svf(2,isvf)
3631	!< canopy shading is applied only to shortwave
3632	surfinlwdif(isurf) = surfinlwdif(isurf) + rad_lw_in_diff(j,i) * svf(1,isvf)
3633	ELSE
3634	!-- for surface-to-surface factors we calculate thermal radiation in 1st pass
3635	surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
3636	ENDIF
3637
3638	IF ( zenith(0) > 0 .AND. all( surf(1:4,isurfsrc) == bdycross ) ) THEN
3639	!-- found svf between model boundary and the face => face isn't shaded
3640	surfinswdir(isurf) = rad_sw_in_dir(j,i) &
3641	* costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0)
3642
3643	ENDIF
3644	ENDDO
3645
3646	IF ( plant_canopy ) THEN
3647
3648	pcbinsw(:) = 0._wp
3649	pcbinlw(:) = 0._wp !< will stay always 0 since we don't absorb lw anymore
3650	!
3651	!-- pcsf first pass
3652	isurf1 = -1 !< previous processed pcgb
3653	DO icsf = 1, ncsfl
3654	ipcgb = csfsurf(1, icsf)
3655	i = pcbl(ix,ipcgb)
3656	j = pcbl(iy,ipcgb)
3657	k = pcbl(iz,ipcgb)
3658	isurfsrc = csfsurf(2, icsf)
3659
3660	IF ( zenith(0) > 0 .AND. ipcgb /= isurf1 ) THEN
3661	!-- locate the virtual surface where the direct solar ray crosses domain boundary
3662	!-- (once per target PC gridbox)
3663	rz = REAL(k, wp)
3664	ry = REAL(j, wp)
3665	rx = REAL(i, wp)
3666	CALL usm_find_boundary_face( (/ rz, ry, rx /), &
3667	sunorig_grid, bdycross)
3668
3669	isurf1 = ipcgb
3670	ENDIF
3671
3672	IF ( surf(id, isurfsrc) >= isky ) THEN
3673	!-- Diffuse rad from boundary surfaces. See comments for svf above.
3674	pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * rad_sw_in_diff(j,i)
3675	!-- canopy shading is applied only to shortwave, therefore no absorbtion for lw
3676	!-- pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * rad_lw_in_diff(j,i)
3677	!ELSE
3678	!-- Thermal radiation in 1st pass
3679	!-- pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc)
3680	ENDIF
3681
3682	IF ( zenith(0) > 0 .AND. all( surf(1:4,isurfsrc) == bdycross ) ) THEN
3683	!-- found svf between model boundary and the pcgb => pcgb isn't shaded
3684	pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
3685	pcbinsw(ipcgb) = pcbinsw(ipcgb) &
3686	+ rad_sw_in_dir(j, i) * pc_box_area * csf(2,icsf) * pc_abs_frac
3687	ENDIF
3688	ENDDO
3689	ENDIF
3690
3691	surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy)
3692	surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy)
3693	surfinsw(:) = surfins(:)
3694	surfinlw(:) = surfinl(:)
3695	surfoutsw(:) = 0.0_wp
3696	surfoutlw(:) = surfoutll(:)
3697	! surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
3698	! - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
3699
3700	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3701	!-- Next passes - reflections
3702	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3703	DO refstep = 1, nrefsteps
3704
3705	surfoutsl(startenergy:endenergy) = albedo_surf(startenergy:endenergy) * surfins(startenergy:endenergy)
3706	!-- for non-transparent surfaces, longwave albedo is 1 - emissivity
3707	surfoutll(startenergy:endenergy) = (1._wp - emiss_surf(startenergy:endenergy)) * surfinl(startenergy:endenergy)
3708
3709	#if defined( __parallel )
3710	CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
3711	surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
3712	CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
3713	surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
3714	#else
3715	surfouts(:) = surfoutsl(:)
3716	surfoutl(:) = surfoutll(:)
3717	#endif
3718
3719	!-- reset for next pass input
3720	surfins(:) = 0._wp
3721	surfinl(:) = 0._wp
3722
3723	!-- reflected radiation
3724	DO isvf = 1, nsvfl
3725	isurf = svfsurf(1, isvf)
3726	isurfsrc = svfsurf(2, isvf)
3727
3728	!-- TODO: to remove if, use start+end for isvf
3729	IF ( surf(id, isurfsrc) < isky ) THEN
3730	surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
3731	surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
3732	ENDIF
3733	ENDDO
3734
3735	!-- radiation absorbed by plant canopy
3736	DO icsf = 1, ncsfl
3737	ipcgb = csfsurf(1, icsf)
3738	isurfsrc = csfsurf(2, icsf)
3739
3740	IF ( surf(id, isurfsrc) < isky ) THEN
3741	pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
3742	!-- pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc)
3743	ENDIF
3744	ENDDO
3745
3746	surfinsw(:) = surfinsw(:) + surfins(:)
3747	surfinlw(:) = surfinlw(:) + surfinl(:)
3748	surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy)
3749	surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy)
3750	! surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
3751	! - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
3752
3753	ENDDO
3754
3755	!-- push heat flux absorbed by plant canopy to respective 3D arrays
3756	IF ( plant_canopy ) THEN
3757	pc_heating_rate(:,:,:) = 0._wp
3758	DO ipcgb = 1, npcbl
3759	j = pcbl(iy, ipcgb)
3760	i = pcbl(ix, ipcgb)
3761	k = pcbl(iz, ipcgb)
3762	kk = k - MAXLOC( &
3763	MERGE( 1, 0, &
3764	BTEST( wall_flags_0(:,j,i), 12 ) &
3765	), DIM = 1 &
3766	) - 1 ! kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
3767	pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
3768	* pchf_prep(k) * pt(k, j, i) !-- = dT/dt
3769	ENDDO
3770	ENDIF
3771	!
3772	!-- Transfer radiation arrays required for energy balance to the respective data types
3773	DO i = startenergy, endenergy
3774	m = surfl(5,i)
3775	!
3776	!-- upward-facing
3777	IF ( surfl(1,i) == 0 ) THEN
3778	surf_usm_h%rad_in_sw(m) = surfinsw(i)
3779	surf_usm_h%rad_out_sw(m) = surfoutsw(i)
3780	surf_usm_h%rad_in_lw(m) = surfinlw(i)
3781	surf_usm_h%rad_out_lw(m) = surfoutlw(i)
3782	!
3783	!-- southward-facding
3784	ELSEIF ( surfl(1,i) == 1 ) THEN
3785	surf_usm_v(1)%rad_in_sw(m) = surfinsw(i)
3786	surf_usm_v(1)%rad_out_sw(m) = surfoutsw(i)
3787	surf_usm_v(1)%rad_in_lw(m) = surfinlw(i)
3788	surf_usm_v(1)%rad_out_lw(m) = surfoutlw(i)
3789	!
3790	!-- northward-facding
3791	ELSEIF ( surfl(1,i) == 2 ) THEN
3792	surf_usm_v(0)%rad_in_sw(m) = surfinsw(i)
3793	surf_usm_v(0)%rad_out_sw(m) = surfoutsw(i)
3794	surf_usm_v(0)%rad_in_lw(m) = surfinlw(i)
3795	surf_usm_v(0)%rad_out_lw(m) = surfoutlw(i)
3796	!
3797	!-- westward-facding
3798	ELSEIF ( surfl(1,i) == 3 ) THEN
3799	surf_usm_v(3)%rad_in_sw(m) = surfinsw(i)
3800	surf_usm_v(3)%rad_out_sw(m) = surfoutsw(i)
3801	surf_usm_v(3)%rad_in_lw(m) = surfinlw(i)
3802	surf_usm_v(3)%rad_out_lw(m) = surfoutlw(i)
3803	!
3804	!-- eastward-facing
3805	ELSEIF ( surfl(1,i) == 4 ) THEN
3806	surf_usm_v(2)%rad_in_sw(m) = surfinsw(i)
3807	surf_usm_v(2)%rad_out_sw(m) = surfoutsw(i)
3808	surf_usm_v(2)%rad_in_lw(m) = surfinlw(i)
3809	surf_usm_v(2)%rad_out_lw(m) = surfoutlw(i)
3810	ENDIF
3811
3812	ENDDO
3813
3814
3815	DO m = 1, surf_usm_h%ns
3816	surf_usm_h%surfhf(m) = surf_usm_h%rad_in_sw(m) + &
3817	surf_usm_h%rad_in_lw(m) - &
3818	surf_usm_h%rad_out_sw(m) - &
3819	surf_usm_h%rad_out_lw(m)
3820	ENDDO
3821
3822	DO l = 0, 3
3823	DO m = 1, surf_usm_v(l)%ns
3824	surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_in_sw(m) + &
3825	surf_usm_v(l)%rad_in_lw(m) - &
3826	surf_usm_v(l)%rad_out_sw(m) - &
3827	surf_usm_v(l)%rad_out_lw(m)
3828	ENDDO
3829	ENDDO
3830
3831	!-- return surface radiation to horizontal surfaces
3832	!-- to rad_sw_in, rad_lw_in and rad_net for outputs
3833	!!!!!!!!!!
3834	!-- we need the original radiation on urban top layer
3835	!-- for calculation of MRT so we can't do adjustment here for now
3836	!!!!!!!!!!
3837	!!!DO isurf = 1, nsurfl
3838	!!! i = surfl(ix,isurf)
3839	!!! j = surfl(iy,isurf)
3840	!!! k = surfl(iz,isurf)
3841	!!! d = surfl(id,isurf)
3842	!!! IF ( d==iroof ) THEN
3843	!!! rad_sw_in(:,j,i) = surfinsw(isurf)
3844	!!! rad_lw_in(:,j,i) = surfinlw(isurf)
3845	!!! 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)
3846	!!! ENDIF
3847	!!!ENDDO
3848
3849	END SUBROUTINE usm_radiation
3850
3851
3852	!------------------------------------------------------------------------------!
3853	! Description:
3854	! ------------
3855	!> Raytracing for detecting obstacles and calculating compound canopy sink
3856	!> factors. (A simple obstacle detection would only need to process faces in
3857	!> 3 dimensions without any ordering.)
3858	!> Assumtions:
3859	!> -----------
3860	!> 1. The ray always originates from a face midpoint (only one coordinate equals
3861	!> *.5, i.e. wall) and doesn't travel parallel to the surface (that would mean
3862	!> shape factor=0). Therefore, the ray may never travel exactly along a face
3863	!> or an edge.
3864	!> 2. From grid bottom to urban surface top the grid has to be equidistant
3865	!> within each of the dimensions, including vertical (but the resolution
3866	!> doesn't need to be the same in all three dimensions).
3867	!------------------------------------------------------------------------------!
3868	SUBROUTINE usm_raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency, win_lad)
3869	IMPLICIT NONE
3870
3871	REAL(wp), DIMENSION(3), INTENT(in) :: src, targ !< real coordinates z,y,x
3872	INTEGER(iwp), INTENT(in) :: isrc !< index of source face for csf
3873	REAL(wp), INTENT(in) :: rirrf !< irradiance factor for csf
3874	REAL(wp), INTENT(in) :: atarg !< target surface area for csf
3875	LOGICAL, INTENT(in) :: create_csf !< whether to generate new CSFs during raytracing
3876	LOGICAL, INTENT(out) :: visible
3877	REAL(wp), INTENT(out) :: transparency !< along whole path
3878	INTEGER(iwp), INTENT(in) :: win_lad
3879	INTEGER(iwp) :: i, j, k, d
3880	INTEGER(iwp) :: seldim !< dimension to be incremented
3881	INTEGER(iwp) :: ncsb !< no of written plant canopy sinkboxes
3882	INTEGER(iwp) :: maxboxes !< max no of gridboxes visited
3883	REAL(wp) :: distance !< euclidean along path
3884	REAL(wp) :: crlen !< length of gridbox crossing
3885	REAL(wp) :: lastdist !< beginning of current crossing
3886	REAL(wp) :: nextdist !< end of current crossing
3887	REAL(wp) :: realdist !< distance in meters per unit distance
3888	REAL(wp) :: crmid !< midpoint of crossing
3889	REAL(wp) :: cursink !< sink factor for current canopy box
3890	REAL(wp), DIMENSION(3) :: delta !< path vector
3891	REAL(wp), DIMENSION(3) :: uvect !< unit vector
3892	REAL(wp), DIMENSION(3) :: dimnextdist !< distance for each dimension increments
3893	INTEGER(iwp), DIMENSION(3) :: box !< gridbox being crossed
3894	INTEGER(iwp), DIMENSION(3) :: dimnext !< next dimension increments along path
3895	INTEGER(iwp), DIMENSION(3) :: dimdelta !< dimension direction = +- 1
3896	INTEGER(iwp) :: px, py !< number of processors in x and y dir before
3897	!< the processor in the question
3898	INTEGER(iwp) :: ip !< number of processor where gridbox reside
3899	INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array
3900	REAL(wp) :: lad_s_target !< recieved lad_s of particular grid box
3901	REAL(wp), PARAMETER :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
3902
3903	!
3904	!-- Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
3905	!-- the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
3906	maxboxes = SUM(ABS(NINT(targ) - NINT(src))) + 1
3907	IF ( plant_canopy .AND. ncsfl + maxboxes > ncsfla ) THEN
3908	!-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
3909	!-- k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
3910	!-- / log(grow_factor)), kind=wp))
3911	!-- or use this code to simply always keep some extra space after growing
3912	k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
3913
3914	CALL usm_merge_and_grow_csf(k)
3915	ENDIF
3916
3917	transparency = 1._wp
3918	ncsb = 0
3919
3920	delta(:) = targ(:) - src(:)
3921	distance = SQRT(SUM(delta(:)**2))
3922	IF ( distance == 0._wp ) THEN
3923	visible = .TRUE.
3924	RETURN
3925	ENDIF
3926	uvect(:) = delta(:) / distance
3927	realdist = SQRT(SUM( (uvect(:)(/dz,dy,dx/))*2 ))
3928
3929	lastdist = 0._wp
3930
3931	!-- Since all face coordinates have values *.5 and we'd like to use
3932	!-- integers, all these have .5 added
3933	DO d = 1, 3
3934	IF ( uvect(d) == 0._wp ) THEN
3935	dimnext(d) = 999999999
3936	dimdelta(d) = 999999999
3937	dimnextdist(d) = 1.0E20_wp
3938	ELSE IF ( uvect(d) > 0._wp ) THEN
3939	dimnext(d) = CEILING(src(d) + .5_wp)
3940	dimdelta(d) = 1
3941	dimnextdist(d) = (dimnext(d) - .5_wp - src(d)) / uvect(d)
3942	ELSE
3943	dimnext(d) = FLOOR(src(d) + .5_wp)
3944	dimdelta(d) = -1
3945	dimnextdist(d) = (dimnext(d) - .5_wp - src(d)) / uvect(d)
3946	ENDIF
3947	ENDDO
3948
3949	DO
3950	!-- along what dimension will the next wall crossing be?
3951	seldim = minloc(dimnextdist, 1)
3952	nextdist = dimnextdist(seldim)
3953	IF ( nextdist > distance ) nextdist = distance
3954
3955	crlen = nextdist - lastdist
3956	IF ( crlen > .001_wp ) THEN
3957	crmid = (lastdist + nextdist) * .5_wp
3958	box = NINT(src(:) + uvect(:) * crmid)
3959
3960	!-- calculate index of the grid with global indices (box(2),box(3))
3961	!-- in the array nzterr and plantt and id of the coresponding processor
3962	px = box(3)/nnx
3963	py = box(2)/nny
3964	ip = px*pdims(2)+py
3965	ig = ipnnxnny + (box(3)-pxnnx)nny + box(2)-py*nny
3966	IF ( box(1) <= nzterr(ig) ) THEN
3967	visible = .FALSE.
3968	RETURN
3969	ENDIF
3970
3971	IF ( plant_canopy ) THEN
3972	IF ( box(1) <= plantt(ig) ) THEN
3973	ncsb = ncsb + 1
3974	boxes(:,ncsb) = box
3975	crlens(ncsb) = crlen
3976	#if defined( __parallel )
3977	lad_ip(ncsb) = ip
3978	lad_disp(ncsb) = (box(3)-pxnnx)(nnynzu) + (box(2)-pynny)*nzu + box(1)-nzub
3979	#endif
3980	ENDIF
3981	ENDIF
3982	ENDIF
3983
3984	IF ( nextdist >= distance ) EXIT
3985	lastdist = nextdist
3986	dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
3987	dimnextdist(seldim) = (dimnext(seldim) - .5_wp - src(seldim)) / uvect(seldim)
3988	ENDDO
3989
3990	IF ( plant_canopy ) THEN
3991	#if defined( __parallel )
3992	IF ( usm_lad_rma ) THEN
3993	!-- send requests for lad_s to appropriate processor
3994	CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
3995	DO i = 1, ncsb
3996	CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
3997	1, MPI_REAL, win_lad, ierr)
3998	IF ( ierr /= 0 ) THEN
3999	WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
4000	CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
4001	ENDIF
4002	ENDDO
4003
4004	!-- wait for all pending local requests complete
4005	CALL MPI_Win_flush_local_all(win_lad, ierr)
4006	IF ( ierr /= 0 ) THEN
4007	WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
4008	CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
4009	ENDIF
4010	CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
4011
4012	ENDIF
4013	#endif
4014
4015	!-- calculate csf and transparency
4016	DO i = 1, ncsb
4017	#if defined( __parallel )
4018	IF ( usm_lad_rma ) THEN
4019	lad_s_target = lad_s_ray(i)
4020	ELSE
4021	lad_s_target = usm_lad_g(lad_ip(i)nnxnny*nzu + lad_disp(i))
4022	ENDIF
4023	#else
4024	lad_s_target = usm_lad(boxes(1,i),boxes(2,i),boxes(3,i))
4025	#endif
4026	cursink = 1._wp - exp(-ext_coef * lad_s_target * crlens(i)*realdist)
4027
4028	IF ( create_csf ) THEN
4029	!-- write svf values into the array
4030	ncsfl = ncsfl + 1
4031	acsf(ncsfl)%ip = lad_ip(i)
4032	acsf(ncsfl)%itx = boxes(3,i)
4033	acsf(ncsfl)%ity = boxes(2,i)
4034	acsf(ncsfl)%itz = boxes(1,i)
4035	acsf(ncsfl)%isurfs = isrc
4036	acsf(ncsfl)%rsvf = REAL(cursinkrirrfatarg, wp) !-- we postpone multiplication by transparency
4037	acsf(ncsfl)%rtransp = REAL(transparency, wp)
4038	ENDIF !< create_csf
4039
4040	transparency = transparency * (1._wp - cursink)
4041
4042	ENDDO
4043	ENDIF
4044
4045	visible = .TRUE.
4046
4047	END SUBROUTINE usm_raytrace
4048
4049
4050	!------------------------------------------------------------------------------!
4051	! Description:
4052	! ------------
4053	!
4054	!> This subroutine is part of the urban surface model.
4055	!> It reads daily heat produced by anthropogenic sources
4056	!> and the diurnal cycle of the heat.
4057	!------------------------------------------------------------------------------!
4058	SUBROUTINE usm_read_anthropogenic_heat
4059
4060	INTEGER(iwp) :: i,j,ii
4061	REAL(wp) :: heat
4062
4063	!-- allocation of array of sources of anthropogenic heat and their diural profile
4064	ALLOCATE( aheat(nys:nyn,nxl:nxr) )
4065	ALLOCATE( aheatprof(0:24) )
4066
4067	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4068	!-- read daily amount of heat and its daily cycle
4069	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4070	aheat = 0.0_wp
4071	DO ii = 0, io_blocks-1
4072	IF ( ii == io_group ) THEN
4073
4074	!-- open anthropogenic heat file
4075	OPEN( 151, file='ANTHROPOGENIC_HEAT'//TRIM(coupling_char), action='read', &
4076	status='old', form='formatted', err=11 )
4077	i = 0
4078	j = 0
4079	DO
4080	READ( 151, *, err=12, end=13 ) i, j, heat
4081	IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn ) THEN
4082	!-- write heat into the array
4083	aheat(j,i) = heat
4084	ENDIF
4085	CYCLE
4086	12 WRITE(message_string,'(a,2i4)') 'error in file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' after line ',i,j
4087	CALL message( 'usm_read_anthropogenic_heat', 'PA0515', 0, 1, 0, 6, 0 )
4088	ENDDO
4089	13 CLOSE(151)
4090	CYCLE
4091	11 message_string = 'file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' does not exist'
4092	CALL message( 'usm_read_anthropogenic_heat', 'PA0516', 1, 2, 0, 6, 0 )
4093	ENDIF
4094
4095	#if defined( __parallel ) && ! defined ( __check )
4096	CALL MPI_BARRIER( comm2d, ierr )
4097	#endif
4098	ENDDO
4099
4100	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4101	!-- read diurnal profiles of heat sources
4102	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4103	aheatprof = 0.0_wp
4104	DO ii = 0, io_blocks-1
4105	IF ( ii == io_group ) THEN
4106
4107	!-- open anthropogenic heat profile file
4108	OPEN( 151, file='ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char), action='read', &
4109	status='old', form='formatted', err=21 )
4110	i = 0
4111	DO
4112	READ( 151, *, err=22, end=23 ) i, heat
4113	IF ( i >= 0 .AND. i <= 24 ) THEN
4114	!-- write heat into the array
4115	aheatprof(i) = heat
4116	ENDIF
4117	CYCLE
4118	22 WRITE(message_string,'(a,i4)') 'error in file ANTHROPOGENIC_HEAT_PROFILE'// &
4119	TRIM(coupling_char)//' after line ',i
4120	CALL message( 'usm_read_anthropogenic_heat', 'PA0517', 0, 1, 0, 6, 0 )
4121	ENDDO
4122	aheatprof(24) = aheatprof(0)
4123	23 CLOSE(151)
4124	CYCLE
4125	21 message_string = 'file ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char)//' does not exist'
4126	CALL message( 'usm_read_anthropogenic_heat', 'PA0518', 1, 2, 0, 6, 0 )
4127	ENDIF
4128
4129	#if defined( __parallel ) && ! defined ( __check )
4130	CALL MPI_BARRIER( comm2d, ierr )
4131	#endif
4132	ENDDO
4133
4134	END SUBROUTINE usm_read_anthropogenic_heat
4135
4136
4137	!------------------------------------------------------------------------------!
4138	!
4139	! Description:
4140	! ------------
4141	!> Soubroutine reads t_surf and t_wall data from restart files
4142	!kanani: Renamed this routine according to corresponging routines in PALM
4143	!kanani: Modified the routine to match read_var_list, from where usm_read_restart_data
4144	! shall be called in the future. This part has not been tested yet. (see virtual_flight_mod)
4145	! Also, I had some trouble with the allocation of t_surf, since this is a pointer.
4146	! So, I added some directives here.
4147	!------------------------------------------------------------------------------!
4148	SUBROUTINE usm_read_restart_data
4149
4150
4151	IMPLICIT NONE
4152
4153	CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string
4154
4155	INTEGER(iwp) :: i !< running index
4156
4157
4158	DO i = 0, io_blocks-1
4159	IF ( i == io_group ) THEN
4160	READ ( 13 ) variable_chr
4161	DO WHILE ( TRIM( variable_chr ) /= '* end usm *' )
4162
4163	SELECT CASE ( TRIM( variable_chr ) )
4164
4165	CASE ( 't_surf_h' )
4166	#if defined( __nopointer )
4167	IF ( .NOT. ALLOCATED( t_surf_h ) ) &
4168	ALLOCATE( t_surf_h(1:surf_usm_h%ns) )
4169	READ ( 13 ) t_surf_h
4170	#else
4171	IF ( .NOT. ALLOCATED( t_surf_h_1 ) ) &
4172	ALLOCATE( t_surf_h_1(1:surf_usm_h%ns) )
4173	READ ( 13 ) t_surf_h_1
4174	#endif
4175	CASE ( 't_surf_v(0)' )
4176	#if defined( __nopointer )
4177	IF ( .NOT. ALLOCATED( t_surf_v(0)%t ) ) &
4178	ALLOCATE( t_surf_v(0)%t(1:surf_usm_v(0)%ns) )
4179	READ ( 13 ) t_surf_v(0)%t
4180	#else
4181	IF ( .NOT. ALLOCATED( t_surf_v_1(0)%t ) ) &
4182	ALLOCATE( t_surf_v_1(0)%t(1:surf_usm_v(0)%ns) )
4183	READ ( 13 ) t_surf_v_1(0)%t
4184	#endif
4185	CASE ( 't_surf_v(1)' )
4186	#if defined( __nopointer )
4187	IF ( .NOT. ALLOCATED( t_surf_v(1)%t ) ) &
4188	ALLOCATE( t_surf_v(1)%t(1:surf_usm_v(1)%ns) )
4189	READ ( 13 ) t_surf_v(1)%t
4190	#else
4191	IF ( .NOT. ALLOCATED( t_surf_v_1(1)%t ) ) &
4192	ALLOCATE( t_surf_v_1(1)%t(1:surf_usm_v(1)%ns) )
4193	READ ( 13 ) t_surf_v_1(1)%t
4194	#endif
4195	CASE ( 't_surf_v(2)' )
4196	#if defined( __nopointer )
4197	IF ( .NOT. ALLOCATED( t_surf_v(2)%t ) ) &
4198	ALLOCATE( t_surf_v(2)%t(1:surf_usm_v(2)%ns) )
4199	READ ( 13 ) t_surf_v(2)%t
4200	#else
4201	IF ( .NOT. ALLOCATED( t_surf_v_1(2)%t ) ) &
4202	ALLOCATE( t_surf_v_1(2)%t(1:surf_usm_v(2)%ns) )
4203	READ ( 13 ) t_surf_v_1(2)%t
4204	#endif
4205	CASE ( 't_surf_v(3)' )
4206	#if defined( __nopointer )
4207	IF ( .NOT. ALLOCATED( t_surf_v(3)%t ) ) &
4208	ALLOCATE( t_surf_v(3)%t(1:surf_usm_v(3)%ns) )
4209	READ ( 13 ) t_surf_v(3)%t
4210	#else
4211	IF ( .NOT. ALLOCATED( t_surf_v_1(3)%t ) ) &
4212	ALLOCATE( t_surf_v_1(3)%t(1:surf_usm_v(3)%ns) )
4213	READ ( 13 ) t_surf_v_1(3)%t
4214	#endif
4215	CASE ( 't_wall_h' )
4216	#if defined( __nopointer )
4217	IF ( .NOT. ALLOCATED( t_wall_h ) ) &
4218	ALLOCATE( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
4219	READ ( 13 ) t_wall_h
4220	#else
4221	IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) &
4222	ALLOCATE( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
4223	READ ( 13 ) t_wall_h_1
4224	#endif
4225	CASE ( 't_wall_v(0)' )
4226	#if defined( __nopointer )
4227	IF ( .NOT. ALLOCATED( t_wall_v(0)%t ) ) &
4228	ALLOCATE( t_wall_v(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
4229	READ ( 13 ) t_wall_v(0)%t
4230	#else
4231	IF ( .NOT. ALLOCATED( t_wall_v_1(0)%t ) ) &
4232	ALLOCATE( t_wall_v_1(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
4233	READ ( 13 ) t_wall_v_1(0)%t
4234	#endif
4235	CASE ( 't_wall_v(1)' )
4236	#if defined( __nopointer )
4237	IF ( .NOT. ALLOCATED( t_wall_v(1)%t ) ) &
4238	ALLOCATE( t_wall_v(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
4239	READ ( 13 ) t_wall_v(1)%t
4240	#else
4241	IF ( .NOT. ALLOCATED( t_wall_v_1(0)%t ) ) &
4242	ALLOCATE( t_wall_v_1(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
4243	READ ( 13 ) t_wall_v_1(1)%t
4244	#endif
4245	CASE ( 't_wall_v(2)' )
4246	#if defined( __nopointer )
4247	IF ( .NOT. ALLOCATED( t_wall_v(2)%t ) ) &
4248	ALLOCATE( t_wall_v(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
4249	READ ( 13 ) t_wall_v(2)%t
4250	#else
4251	IF ( .NOT. ALLOCATED( t_wall_v_1(2)%t ) ) &
4252	ALLOCATE( t_wall_v_1(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
4253	READ ( 13 ) t_wall_v_1(2)%t
4254	#endif
4255	CASE ( 't_wall_v(3)' )
4256	#if defined( __nopointer )
4257	IF ( .NOT. ALLOCATED( t_wall_v(3)%t ) ) &
4258	ALLOCATE( t_wall_v(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
4259	READ ( 13 ) t_wall_v(3)%t
4260	#else
4261	IF ( .NOT. ALLOCATED( t_wall_v_1(3)%t ) ) &
4262	ALLOCATE( t_wall_v_1(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
4263	READ ( 13 ) t_wall_v_1(3)%t
4264	#endif
4265
4266	CASE DEFAULT
4267	WRITE ( message_string, * ) 'unknown variable named "', &
4268	TRIM( variable_chr ), '" found in', &
4269	'&data from prior run on PE ', myid
4270	CALL message( 'user_read_restart_data', 'UI0012', 1, 2, 0, 6, 0 )
4271
4272	END SELECT
4273
4274	READ ( 13 ) variable_chr
4275
4276	ENDDO
4277	ENDIF
4278	#if defined( __parallel )
4279	CALL MPI_BARRIER( comm2d, ierr )
4280	#endif
4281	ENDDO
4282
4283	END SUBROUTINE usm_read_restart_data
4284
4285
4286	!------------------------------------------------------------------------------!
4287	!
4288	! Description:
4289	! ------------
4290	!> Soubroutine reads svf and svfsurf data from saved file
4291	!------------------------------------------------------------------------------!
4292	SUBROUTINE usm_read_svf_from_file
4293
4294	IMPLICIT NONE
4295	INTEGER(iwp) :: fsvf = 89
4296	INTEGER(iwp) :: i
4297	CHARACTER(usm_version_len) :: usm_version_field
4298	CHARACTER(svf_code_len) :: svf_code_field
4299
4300	DO i = 0, io_blocks-1
4301	IF ( i == io_group ) THEN
4302	OPEN ( fsvf, file=TRIM(svf_file_name)//TRIM(coupling_char)//myid_char, &
4303	form='unformatted', status='old' )
4304
4305	!-- read and check version
4306	READ ( fsvf ) usm_version_field
4307	IF ( TRIM(usm_version_field) /= TRIM(usm_version) ) THEN
4308	WRITE( message_string, * ) 'Version of binary SVF file "', &
4309	TRIM(usm_version_field), '" does not match ', &
4310	'the version of model "', TRIM(usm_version), '"'
4311	CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
4312	ENDIF
4313
4314	!-- read nsvfl, ncsfl
4315	READ ( fsvf ) nsvfl, ncsfl
4316	IF ( nsvfl <= 0 .OR. ncsfl < 0 ) THEN
4317	WRITE( message_string, * ) 'Wrong number of SVF or CSF'
4318	CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
4319	ELSE
4320	WRITE(message_string,*) ' Number of SVF and CSF to read', nsvfl, ncsfl
4321	CALL location_message( message_string, .TRUE. )
4322	ENDIF
4323
4324	ALLOCATE(svf(ndsvf,nsvfl))
4325	ALLOCATE(svfsurf(idsvf,nsvfl))
4326	READ(fsvf) svf
4327	READ(fsvf) svfsurf
4328	IF ( plant_canopy ) THEN
4329	ALLOCATE(csf(ndcsf,ncsfl))
4330	ALLOCATE(csfsurf(idcsf,ncsfl))
4331	READ(fsvf) csf
4332	READ(fsvf) csfsurf
4333	ENDIF
4334	READ ( fsvf ) svf_code_field
4335
4336	IF ( TRIM(svf_code_field) /= TRIM(svf_code) ) THEN
4337	WRITE( message_string, * ) 'Wrong structure of binary svf file'
4338	CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
4339	ENDIF
4340
4341	CLOSE (fsvf)
4342
4343	ENDIF
4344	#if defined( __parallel )
4345	CALL MPI_BARRIER( comm2d, ierr )
4346	#endif
4347	ENDDO
4348
4349	END SUBROUTINE usm_read_svf_from_file
4350
4351
4352	!------------------------------------------------------------------------------!
4353	! Description:
4354	! ------------
4355	!
4356	!> This subroutine reads walls, roofs and land categories and it parameters
4357	!> from input files.
4358	!------------------------------------------------------------------------------!
4359	SUBROUTINE usm_read_urban_surface_types
4360
4361	CHARACTER(12) :: wtn
4362	INTEGER(iwp) :: wtc
4363	REAL(wp), DIMENSION(n_surface_params) :: wtp
4364
4365	INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg) :: usm_par
4366	REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg) :: usm_val
4367	INTEGER(iwp) :: k, l, d, iw, jw, kw, it, ip, ii, ij, m
4368	INTEGER(iwp) :: i, j
4369	INTEGER(iwp) :: nz, roof, dirwe, dirsn
4370	INTEGER(iwp) :: category
4371	INTEGER(iwp) :: weheight1, wecat1, snheight1, sncat1
4372	INTEGER(iwp) :: weheight2, wecat2, snheight2, sncat2
4373	INTEGER(iwp) :: weheight3, wecat3, snheight3, sncat3
4374	REAL(wp) :: height, albedo, thick
4375	REAL(wp) :: wealbedo1, wethick1, snalbedo1, snthick1
4376	REAL(wp) :: wealbedo2, wethick2, snalbedo2, snthick2
4377	REAL(wp) :: wealbedo3, wethick3, snalbedo3, snthick3
4378
4379	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4380	!-- read categories of walls and their parameters
4381	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4382	DO ii = 0, io_blocks-1
4383	IF ( ii == io_group ) THEN
4384
4385	!-- open urban surface file
4386	OPEN( 151, file='SURFACE_PARAMETERS'//coupling_char, action='read', &
4387	status='old', form='formatted', err=15 )
4388	!-- first test and get n_surface_types
4389	k = 0
4390	l = 0
4391	DO
4392	l = l+1
4393	READ( 151, *, err=11, end=12 ) wtc, wtp, wtn
4394	k = k+1
4395	CYCLE
4396	11 CONTINUE
4397	ENDDO
4398	12 n_surface_types = k
4399	ALLOCATE( surface_type_names(n_surface_types) )
4400	ALLOCATE( surface_type_codes(n_surface_types) )
4401	ALLOCATE( surface_params(n_surface_params, n_surface_types) )
4402	!-- real reading
4403	rewind( 151 )
4404	k = 0
4405	DO
4406	READ( 151, *, err=13, end=14 ) wtc, wtp, wtn
4407	k = k+1
4408	surface_type_codes(k) = wtc
4409	surface_params(:,k) = wtp
4410	surface_type_names(k) = wtn
4411	CYCLE
4412	13 WRITE(6,'(i3,a,2i5)') myid, 'readparams2 error k=', k
4413	FLUSH(6)
4414	CONTINUE
4415	ENDDO
4416	14 CLOSE(151)
4417	CYCLE
4418	15 message_string = 'file SURFACE_PARAMETERS'//TRIM(coupling_char)//' does not exist'
4419	CALL message( 'usm_read_urban_surface_types', 'PA0513', 1, 2, 0, 6, 0 )
4420	ENDIF
4421	ENDDO
4422
4423	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4424	!-- read types of surfaces
4425	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4426	usm_par = 0
4427	DO ii = 0, io_blocks-1
4428	IF ( ii == io_group ) THEN
4429
4430	!
4431	!-- open csv urban surface file
4432	OPEN( 151, file='URBAN_SURFACE'//TRIM(coupling_char), action='read', &
4433	status='old', form='formatted', err=23 )
4434
4435	l = 0
4436	DO
4437	l = l+1
4438	!-- i, j, height, nz, roof, dirwe, dirsn, category, soilcat,
4439	!-- weheight1, wecat1, snheight1, sncat1, weheight2, wecat2, snheight2, sncat2,
4440	!-- weheight3, wecat3, snheight3, sncat3
4441	READ( 151, *, err=21, end=25 ) i, j, height, nz, roof, dirwe, dirsn, &
4442	category, albedo, thick, &
4443	weheight1, wecat1, wealbedo1, wethick1, &
4444	weheight2, wecat2, wealbedo2, wethick2, &
4445	weheight3, wecat3, wealbedo3, wethick3, &
4446	snheight1, sncat1, snalbedo1, snthick1, &
4447	snheight2, sncat2, snalbedo2, snthick2, &
4448	snheight3, sncat3, snalbedo3, snthick3
4449
4450	IF ( i >= nxlg .AND. i <= nxrg .AND. j >= nysg .AND. j <= nyng ) THEN
4451	!-- write integer variables into array
4452	usm_par(:,j,i) = (/1, nz, roof, dirwe, dirsn, category, &
4453	weheight1, wecat1, weheight2, wecat2, weheight3, wecat3, &
4454	snheight1, sncat1, snheight2, sncat2, snheight3, sncat3 /)
4455	!-- write real values into array
4456	usm_val(:,j,i) = (/ albedo, thick, &
4457	wealbedo1, wethick1, wealbedo2, wethick2, &
4458	wealbedo3, wethick3, snalbedo1, snthick1, &
4459	snalbedo2, snthick2, snalbedo3, snthick3 /)
4460	ENDIF
4461	CYCLE
4462	21 WRITE (message_string, "(A,I5)") 'errors in file URBAN_SURFACE'//TRIM(coupling_char)//' on line ', l
4463	CALL message( 'usm_read_urban_surface_types', 'PA0512', 0, 1, 0, 6, 0 )
4464	ENDDO
4465
4466	23 message_string = 'file URBAN_SURFACE'//TRIM(coupling_char)//' does not exist'
4467	CALL message( 'usm_read_urban_surface_types', 'PA0514', 1, 2, 0, 6, 0 )
4468
4469	25 CLOSE( 90 )
4470
4471	ENDIF
4472	#if defined( __parallel ) && ! defined ( __check )
4473	CALL MPI_BARRIER( comm2d, ierr )
4474	#endif
4475	ENDDO
4476
4477	!
4478	!-- check completeness and formal correctness of the data
4479	DO i = nxlg, nxrg
4480	DO j = nysg, nyng
4481	IF ( usm_par(0,j,i) /= 0 .AND. ( & !< incomplete data,supply default values later
4482	usm_par(1,j,i) < nzb .OR. &
4483	usm_par(1,j,i) > nzt .OR. & !< incorrect height (nz < nzb .OR. nz > nzt)
4484	usm_par(2,j,i) < 0 .OR. &
4485	usm_par(2,j,i) > 1 .OR. & !< incorrect roof sign
4486	usm_par(3,j,i) < nzb-nzt .OR. &
4487	usm_par(3,j,i) > nzt-nzb .OR. & !< incorrect west-east wall direction sign
4488	usm_par(4,j,i) < nzb-nzt .OR. &
4489	usm_par(4,j,i) > nzt-nzb .OR. & !< incorrect south-north wall direction sign
4490	usm_par(6,j,i) < nzb .OR. &
4491	usm_par(6,j,i) > nzt .OR. & !< incorrect pedestrian level height for west-east wall
4492	usm_par(8,j,i) > nzt .OR. &
4493	usm_par(10,j,i) > nzt .OR. & !< incorrect wall or roof level height for west-east wall
4494	usm_par(12,j,i) < nzb .OR. &
4495	usm_par(12,j,i) > nzt .OR. & !< incorrect pedestrian level height for south-north wall
4496	usm_par(14,j,i) > nzt .OR. &
4497	usm_par(16,j,i) > nzt & !< incorrect wall or roof level height for south-north wall
4498	) ) THEN
4499	!-- incorrect input data
4500	WRITE (message_string, "(A,2I5)") 'missing or incorrect data in file URBAN_SURFACE'// &
4501	TRIM(coupling_char)//' for i,j=', i,j
4502	CALL message( 'usm_read_urban_surface', 'PA0504', 1, 2, 0, 6, 0 )
4503	ENDIF
4504
4505	ENDDO
4506	ENDDO
4507	!
4508	!-- Assign the surface types to the respective data type.
4509	!-- First, for horizontal upward-facing surfaces.
4510	DO m = 1, surf_usm_h%ns
4511	iw = surf_usm_h%i(m)
4512	jw = surf_usm_h%j(m)
4513	kw = surf_usm_h%k(m)
4514
4515	IF ( usm_par(5,jw,iw) == 0 ) THEN
4516	IF ( zu(kw) >= roof_height_limit ) THEN
4517	surf_usm_h%isroof_surf(m) = .TRUE.
4518	surf_usm_h%surface_types(m) = roof_category !< default category for root surface
4519	ELSE
4520	surf_usm_h%isroof_surf(m) = .FALSE.
4521	surf_usm_h%surface_types(m) = land_category !< default category for land surface
4522	ENDIF
4523	surf_usm_h%albedo_surf(m) = -1.0_wp
4524	surf_usm_h%thickness_wall(m) = -1.0_wp
4525	ELSE
4526	IF ( usm_par(2,jw,iw)==0 ) THEN
4527	surf_usm_h%isroof_surf(m) = .FALSE.
4528	surf_usm_h%thickness_wall(m) = -1.0_wp
4529	ELSE
4530	surf_usm_h%isroof_surf(m) = .TRUE.
4531	surf_usm_h%thickness_wall(m) = usm_val(2,jw,iw)
4532	ENDIF
4533	surf_usm_h%surface_types(m) = usm_par(5,jw,iw)
4534	surf_usm_h%albedo_surf(m) = usm_val(1,jw,iw)
4535	ENDIF
4536	!
4537	!-- Find the type position
4538	it = surf_usm_h%surface_types(m)
4539	ip = -99999
4540	DO k = 1, n_surface_types
4541	IF ( surface_type_codes(k) == it ) THEN
4542	ip = k
4543	EXIT
4544	ENDIF
4545	ENDDO
4546	IF ( ip == -99999 ) THEN
4547	!-- wall category not found
4548	WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, &
4549	' not found for i,j,k=', iw,jw,kw
4550	CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
4551	ENDIF
4552	!
4553	!-- Albedo
4554	IF ( surf_usm_h%albedo_surf(m) < 0.0_wp ) THEN
4555	surf_usm_h%albedo_surf(m) = surface_params(ialbedo,ip)
4556	ENDIF
4557	!
4558	!-- emissivity of the wall
4559	surf_usm_h%emiss_surf(m) = surface_params(iemiss,ip)
4560	!
4561	!-- heat conductivity Î»S between air and wall ( W mâ2 Kâ1 )
4562	surf_usm_h%lambda_surf(m) = surface_params(ilambdas,ip)
4563	!
4564	!-- roughness relative to concrete
4565	surf_usm_h%roughness_wall(m) = surface_params(irough,ip)
4566	!
4567	!-- Surface skin layer heat capacity (J mâ2 Kâ1 )
4568	surf_usm_h%c_surface(m) = surface_params(icsurf,ip)
4569	!
4570	!-- wall material parameters:
4571	!-- thickness of the wall (m)
4572	!-- missing values are replaced by default value for category
4573	IF ( surf_usm_h%thickness_wall(m) <= 0.001_wp ) THEN
4574	surf_usm_h%thickness_wall(m) = surface_params(ithick,ip)
4575	ENDIF
4576	!
4577	!-- volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )
4578	surf_usm_h%rho_c_wall(:,m) = surface_params(irhoC,ip)
4579	!
4580	!-- thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
4581	surf_usm_h%lambda_h(:,m) = surface_params(ilambdah,ip)
4582
4583	ENDDO
4584	!
4585	!-- For vertical surface elements ( 0 -- northward-facing, 1 -- southward-facing,
4586	!-- 2 -- eastward-facing, 3 -- westward-facing )
4587	DO l = 0, 3
4588	DO m = 1, surf_usm_v(l)%ns
4589	i = surf_usm_v(l)%i(m)
4590	j = surf_usm_v(l)%j(m)
4591	kw = surf_usm_v(l)%k(m)
4592
4593	IF ( l == 3 ) THEN ! westward facing
4594	iw = i
4595	jw = j
4596	ii = 6
4597	ij = 3
4598	ELSEIF ( l == 2 ) THEN
4599	iw = i-1
4600	jw = j
4601	ii = 6
4602	ij = 3
4603	ELSEIF ( l == 1 ) THEN
4604	iw = i
4605	jw = j
4606	ii = 12
4607	ij = 9
4608	ELSEIF ( l == 0 ) THEN
4609	iw = i
4610	jw = j-1
4611	ii = 12
4612	ij = 9
4613	ENDIF
4614
4615	IF ( kw <= usm_par(ii,jw,iw) ) THEN
4616	!-- pedestrant zone
4617	IF ( usm_par(ii+1,jw,iw) == 0 ) THEN
4618	surf_usm_v(l)%surface_types(m) = pedestrant_category !< default category for wall surface in pedestrant zone
4619	surf_usm_v(l)%albedo_surf(m) = -1.0_wp
4620	surf_usm_v(l)%thickness_wall(m) = -1.0_wp
4621	ELSE
4622	surf_usm_v(l)%surface_types(m) = usm_par(ii+1,jw,iw)
4623	surf_usm_v(l)%albedo_surf(m) = usm_val(ij,jw,iw)
4624	surf_usm_v(l)%thickness_wall(m) = usm_val(ij+1,jw,iw)
4625	ENDIF
4626	ELSE IF ( kw <= usm_par(ii+2,jw,iw) ) THEN
4627	!-- wall zone
4628	IF ( usm_par(ii+3,jw,iw) == 0 ) THEN
4629	surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface
4630	surf_usm_v(l)%albedo_surf(m) = -1.0_wp
4631	surf_usm_v(l)%thickness_wall(m) = -1.0_wp
4632	ELSE
4633	surf_usm_v(l)%surface_types(m) = usm_par(ii+3,jw,iw)
4634	surf_usm_v(l)%albedo_surf(m) = usm_val(ij+2,jw,iw)
4635	surf_usm_v(l)%thickness_wall(m) = usm_val(ij+3,jw,iw)
4636	ENDIF
4637	ELSE IF ( kw <= usm_par(ii+4,jw,iw) ) THEN
4638	!-- roof zone
4639	IF ( usm_par(ii+5,jw,iw) == 0 ) THEN
4640	surf_usm_v(l)%surface_types(m) = roof_category !< default category for roof surface
4641	surf_usm_v(l)%albedo_surf(m) = -1.0_wp
4642	surf_usm_v(l)%thickness_wall(m) = -1.0_wp
4643	ELSE
4644	surf_usm_v(l)%surface_types(m) = usm_par(ii+5,jw,iw)
4645	surf_usm_v(l)%albedo_surf(m) = usm_val(ij+4,jw,iw)
4646	surf_usm_v(l)%thickness_wall(m) = usm_val(ij+5,jw,iw)
4647	ENDIF
4648	ELSE
4649	!-- something wrong
4650	CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 )
4651	ENDIF
4652
4653	!
4654	!-- Find the type position
4655	it = surf_usm_v(l)%surface_types(m)
4656	ip = -99999
4657	DO k = 1, n_surface_types
4658	IF ( surface_type_codes(k) == it ) THEN
4659	ip = k
4660	EXIT
4661	ENDIF
4662	ENDDO
4663	IF ( ip == -99999 ) THEN
4664	!-- wall category not found
4665	WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, &
4666	' not found for i,j,k=', iw,jw,kw
4667	CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
4668	ENDIF
4669	!
4670	!-- Albedo
4671	IF ( surf_usm_v(l)%albedo_surf(m) < 0.0_wp ) THEN
4672	surf_usm_v(l)%albedo_surf(m) = surface_params(ialbedo,ip)
4673	ENDIF
4674	!
4675	!-- emissivity of the wall
4676	surf_usm_v(l)%emiss_surf(m) = surface_params(iemiss,ip)
4677	!
4678	!-- heat conductivity Î»S between air and wall ( W mâ2 Kâ1 )
4679	surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip)
4680	!
4681	!-- roughness relative to concrete
4682	surf_usm_v(l)%roughness_wall(m) = surface_params(irough,ip)
4683	!
4684	!-- Surface skin layer heat capacity (J mâ2 Kâ1 )
4685	surf_usm_v(l)%c_surface(m) = surface_params(icsurf,ip)
4686	!
4687	!-- wall material parameters:
4688	!-- thickness of the wall (m)
4689	!-- missing values are replaced by default value for category
4690	IF ( surf_usm_v(l)%thickness_wall(m) <= 0.001_wp ) THEN
4691	surf_usm_v(l)%thickness_wall(m) = surface_params(ithick,ip)
4692	ENDIF
4693	!
4694	!-- volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )
4695	surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip)
4696	!
4697	!-- thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
4698	surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip)
4699
4700	ENDDO
4701	ENDDO
4702
4703	CALL location_message( ' types and parameters of urban surfaces read', .TRUE. )
4704
4705	END SUBROUTINE usm_read_urban_surface_types
4706
4707
4708	!------------------------------------------------------------------------------!
4709	! Description:
4710	! ------------
4711	!> Solver for the energy balance at the ground/roof/wall surface.
4712	!> It follows basic ideas and structure of lsm_energy_balance
4713	!> with many simplifications and adjustments.
4714	!> TODO better description
4715	!------------------------------------------------------------------------------!
4716	SUBROUTINE usm_surface_energy_balance
4717
4718	IMPLICIT NONE
4719
4720	INTEGER(iwp) :: i, j, k, l, d, m !< running indices
4721
4722	REAL(wp) :: pt1 !< temperature at first grid box adjacent to surface
4723	REAL(wp) :: u1,v1,w1 !< near wall u,v,w
4724	REAL(wp) :: stend !< surface tendency
4725	REAL(wp) :: coef_1 !< first coeficient for prognostic equation
4726	REAL(wp) :: coef_2 !< second coeficient for prognostic equation
4727	REAL(wp) :: rho_cp !< rho_wall_surface * cp
4728	REAL(wp) :: r_a !< aerodynamic resistance for horizontal and vertical surfaces
4729	REAL(wp) :: f_shf !< factor for shf_eb
4730	REAL(wp) :: lambda_surface !< current value of lambda_surface (heat conductivity between air and wall)
4731	REAL(wp) :: Ueff !< effective wind speed for calculation of heat transfer coefficients
4732	REAL(wp) :: httc !< heat transfer coefficient
4733	REAL(wp), DIMENSION(nzub:nzut) :: exn !< value of the Exner function in layers
4734
4735	REAL(wp), DIMENSION(0:4) :: dxdir !< surface normal direction gridbox length
4736	REAL(wp) :: dtime !< simulated time of day (in UTC)
4737	INTEGER(iwp) :: dhour !< simulated hour of day (in UTC)
4738	REAL(wp) :: acoef !< actual coefficient of diurnal profile of anthropogenic heat
4739
4740
4741	dxdir = (/dz,dy,dy,dx,dx/)
4742
4743	exn(:) = (hyp(nzub:nzut) / 100000.0_wp )**0.286_wp !< Exner function
4744	!
4745	!-- First, treat horizontal surface elements
4746
4747	DO m = 1, surf_usm_h%ns
4748	!
4749	!-- Get indices of respective grid point
4750	i = surf_usm_h%i(m)
4751	j = surf_usm_h%j(m)
4752	k = surf_usm_h%k(m)
4753	!
4754	!-- TODO - how to calculate lambda_surface for horizontal surfaces
4755	!-- (lambda_surface is set according to stratification in land surface model)
4756	!-- MS: ???
4757	IF ( surf_usm_h%ol(m) >= 0.0_wp ) THEN
4758	lambda_surface = surf_usm_h%lambda_surf(m)
4759	ELSE
4760	lambda_surface = surf_usm_h%lambda_surf(m)
4761	ENDIF
4762
4763	pt1 = pt(k,j,i)
4764	!
4765	!-- calculate rho * cp coefficient at surface layer
4766	rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) )
4767	!
4768	!-- Calculate aerodyamic resistance.
4769	!-- Calculation for horizontal surfaces follows LSM formulation
4770	!-- pt, us, ts are not available for the prognostic time step,
4771	!-- data from the last time step is used here.
4772
4773	r_a = ( pt1 - t_surf_h(m) / exn(k) ) / &
4774	( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp )
4775
4776	!-- make sure that the resistance does not drop to zero
4777	IF ( ABS(r_a) < 1.0E-10_wp ) r_a = 1.0E-10_wp
4778
4779	!-- the parameterization is developed originally for larger scales
4780	!-- (compare with remark in TUF-3D)
4781	!-- our first experiences show that the parameterization underestimates
4782	!-- r_a in meter resolution.
4783	!-- temporary solution - multiplication by magic constant :-(.
4784	r_a = r_a * ra_horiz_coef
4785
4786	!-- factor for shf_eb
4787	f_shf = rho_cp / r_a
4788
4789	!-- add LW up so that it can be removed in prognostic equation
4790	surf_usm_h%rad_net_l(m) = surf_usm_h%rad_in_sw(m) - &
4791	surf_usm_h%rad_out_sw(m) + &
4792	surf_usm_h%rad_in_lw(m) - &
4793	surf_usm_h%rad_out_lw(m)
4794
4795	!-- numerator of the prognostic equation
4796	coef_1 = surf_usm_h%rad_net_l(m) + &
4797	( 3.0_wp + 1.0_wp ) * surf_usm_h%emiss_surf(m) * sigma_sb * &
4798	t_surf_h(m) ** 4 + &
4799	f_shf * pt1 + &
4800	lambda_surface * t_wall_h(nzb_wall,m)
4801
4802	!-- denominator of the prognostic equation
4803	coef_2 = 4.0_wp * surf_usm_h%emiss_surf(m) * sigma_sb * &
4804	t_surf_h(m) ** 3 &
4805	+ lambda_surface + f_shf / exn(k)
4806
4807	!-- implicit solution when the surface layer has no heat capacity,
4808	!-- otherwise use RK3 scheme.
4809	t_surf_h_p(m) = ( coef_1 * dt_3d * tsc(2) + &
4810	surf_usm_h%c_surface(m) * t_surf_h(m) ) / &
4811	( surf_usm_h%c_surface(m) + coef_2 * dt_3d * tsc(2) )
4812
4813	!-- add RK3 term
4814	t_surf_h_p(m) = t_surf_h_p(m) + dt_3d * tsc(3) * &
4815	surf_usm_h%tt_surface_m(m)
4816
4817	!-- calculate true tendency
4818	stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) * &
4819	surf_usm_h%tt_surface_m(m)) / ( dt_3d * tsc(2) )
4820
4821	!-- calculate t_surf tendencies for the next Runge-Kutta step
4822	IF ( timestep_scheme(1:5) == 'runge' ) THEN
4823	IF ( intermediate_timestep_count == 1 ) THEN
4824	surf_usm_h%tt_surface_m(m) = stend
4825	ELSEIF ( intermediate_timestep_count < &
4826	intermediate_timestep_count_max ) THEN
4827	surf_usm_h%tt_surface_m(m) = -9.5625_wp * stend + &
4828	5.3125_wp * surf_usm_h%tt_surface_m(m)
4829	ENDIF
4830	ENDIF
4831
4832	!-- in case of fast changes in the skin temperature, it is required to
4833	!-- update the radiative fluxes in order to keep the solution stable
4834	IF ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp ) THEN
4835	force_radiation_call_l = .TRUE.
4836	ENDIF
4837	!
4838	!-- for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf.
4839	!-- there is no equivalent surface gridpoint for vertical surfaces.
4840	!-- pt(k,j,i) is calculated for all directions in diffusion_s
4841	!-- using surface and wall heat fluxes
4842	pt(k-1,j,i) = t_surf_h_p(m) / exn(k) ! not for vertical surfaces
4843
4844	!-- calculate fluxes
4845	!-- rad_net_l is never used!
4846	surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) + &
4847	3.0_wp * sigma_sb * &
4848	t_surf_h(m)*4 - 4.0_wp sigma_sb * &
4849	t_surf_h(m)*3 t_surf_h_p(m)
4850	surf_usm_h%wghf_eb(m) = lambda_surface * &
4851	( t_surf_h_p(m) - t_wall_h(nzb_wall,m) )
4852	!
4853	!-- ground/wall/roof surface heat flux
4854	surf_usm_h%wshf_eb(m) = - f_shf * ( pt1 - t_surf_h_p(m) )
4855	!
4856	!-- store kinematic surface heat fluxes for utilization in other processes
4857	!-- diffusion_s, surface_layer_fluxes,...
4858	surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / rho_cp
4859
4860	ENDDO
4861	!
4862	!-- Now, treat vertical surface elements
4863	DO l = 0, 3
4864	DO m = 1, surf_usm_v(l)%ns
4865	!
4866	!-- Get indices of respective grid point
4867	i = surf_usm_v(l)%i(m)
4868	j = surf_usm_v(l)%j(m)
4869	k = surf_usm_v(l)%k(m)
4870	!
4871	!-- TODO - how to calculate lambda_surface for horizontal (??? do you mean verical ???) surfaces
4872	!-- (lambda_surface is set according to stratification in land surface model).
4873	!-- Please note, for vertical surfaces no ol is defined, since
4874	!-- stratification is not considered in this case.
4875	lambda_surface = surf_usm_v(l)%lambda_surf(m)
4876
4877	pt1 = pt(k,j,i)
4878	!
4879	!-- calculate rho * cp coefficient at surface layer
4880	rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) )
4881
4882	!-- Calculation of r_a for vertical surfaces
4883	!--
4884	!-- heat transfer coefficient for forced convection along vertical walls
4885	!-- follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
4886	!--
4887	!-- H = httc (Tsfc - Tair)
4888	!-- httc = rw * (11.8 + 4.2 * Ueff) - 4.0
4889	!--
4890	!-- rw: wall patch roughness relative to 1.0 for concrete
4891	!-- Ueff: effective wind speed
4892	!-- - 4.0 is a reduction of Rowley et al (1930) formulation based on
4893	!-- Cole and Sturrock (1977)
4894	!--
4895	!-- Ucan: Canyon wind speed
4896	!-- wstar: convective velocity
4897	!-- Qs: surface heat flux
4898	!-- zH: height of the convective layer
4899	!-- wstar = (g/TcanQszH)**(1./3.)
4900
4901	!-- Effective velocity components must always
4902	!-- be defined at scalar grid point. The wall normal component is
4903	!-- obtained by simple linear interpolation. ( An alternative would
4904	!-- be an logarithmic interpolation. )
4905	u1 = ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp
4906	v1 = ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp
4907	w1 = ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp
4908
4909	Ueff = SQRT( u12 + v12 + w1**2 )
4910	httc = surf_usm_v(l)%roughness_wall(m) * &
4911	( 11.8_wp + 4.2_wp * Ueff ) - 4.0_wp
4912	f_shf = httc
4913
4914	!-- add LW up so that it can be removed in prognostic equation
4915	surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_in_sw(m) - &
4916	surf_usm_v(l)%rad_out_sw(m) + &
4917	surf_usm_v(l)%rad_in_lw(m) - &
4918	surf_usm_v(l)%rad_out_lw(m)
4919
4920	!-- numerator of the prognostic equation
4921	coef_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout included in calculation of radnet_l
4922	( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emiss_surf(m) * sigma_sb * &
4923	t_surf_v(l)%t(m) ** 4 + &
4924	f_shf * pt1 + &
4925	lambda_surface * t_wall_v(l)%t(nzb_wall,m)
4926
4927	!-- denominator of the prognostic equation
4928	coef_2 = 4.0_wp * surf_usm_v(l)%emiss_surf(m) * sigma_sb * &
4929	t_surf_v(l)%t(m) ** 3 &
4930	+ lambda_surface + f_shf / exn(k)
4931
4932	!-- implicit solution when the surface layer has no heat capacity,
4933	!-- otherwise use RK3 scheme.
4934	t_surf_v_p(l)%t(m) = ( coef_1 * dt_3d * tsc(2) + &
4935	surf_usm_v(l)%c_surface(m) * t_surf_v(l)%t(m) ) / &
4936	( surf_usm_v(l)%c_surface(m) + coef_2 * dt_3d * tsc(2) )
4937
4938	!-- add RK3 term
4939	t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) * &
4940	surf_usm_v(l)%tt_surface_m(m)
4941
4942	!-- calculate true tendency
4943	stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *&
4944	surf_usm_v(l)%tt_surface_m(m) ) / ( dt_3d * tsc(2) )
4945
4946	!-- calculate t_surf tendencies for the next Runge-Kutta step
4947	IF ( timestep_scheme(1:5) == 'runge' ) THEN
4948	IF ( intermediate_timestep_count == 1 ) THEN
4949	surf_usm_v(l)%tt_surface_m(m) = stend
4950	ELSEIF ( intermediate_timestep_count < &
4951	intermediate_timestep_count_max ) THEN
4952	surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend + &
4953	5.3125_wp * surf_usm_h%tt_surface_m(m)
4954	ENDIF
4955	ENDIF
4956
4957	!-- in case of fast changes in the skin temperature, it is required to
4958	!-- update the radiative fluxes in order to keep the solution stable
4959	IF ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp ) THEN
4960	force_radiation_call_l = .TRUE.
4961	ENDIF
4962
4963	!-- calculate fluxes
4964	!-- rad_net_l is never used!
4965	surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_net_l(m) + &
4966	3.0_wp * sigma_sb * &
4967	t_surf_v(l)%t(m)*4 - 4.0_wp sigma_sb * &
4968	t_surf_v(l)%t(m)*3 t_surf_v_p(l)%t(m)
4969
4970	surf_usm_v(l)%wghf_eb(m) = lambda_surface * &
4971	( t_surf_v_p(l)%t(m) - t_wall_v(l)%t(nzb_wall,m) )
4972
4973	!-- ground/wall/roof surface heat flux
4974	surf_usm_v(l)%wshf_eb(m) = - f_shf * ( pt1 - t_surf_v_p(l)%t(m) )
4975
4976	!
4977	!-- store kinematic surface heat fluxes for utilization in other processes
4978	!-- diffusion_s, surface_layer_fluxes,...
4979	surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / rho_cp
4980
4981	ENDDO
4982
4983	ENDDO
4984	!
4985	!-- Add-up anthropogenic heat, for now only at upward-facing surfaces
4986	IF ( usm_anthropogenic_heat .AND. &
4987	intermediate_timestep_count == intermediate_timestep_count_max ) THEN
4988	!-- application of the additional anthropogenic heat sources
4989	!-- we considere the traffic for now so all heat is absorbed
4990	!-- to the first layer, generalization would be worth.
4991
4992	!-- calculation of actual profile coefficient
4993	!-- ??? check time_since_reference_point ???
4994	dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
4995	dhour = INT(dtime/3600.0_wp)
4996	!-- linear interpolation of coeficient
4997	acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))aheatprof(dhour+1)
4998
4999	DO m = 1, surf_usm_h%ns
5000	!
5001	!-- Get indices of respective grid point
5002	i = surf_usm_h%i(m)
5003	j = surf_usm_h%j(m)
5004	k = surf_usm_h%k(m)
5005
5006	IF ( aheat(j,i) > 0.0_wp ) THEN
5007	!-- TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d
5008	!-- given to anthropogenic heat aheatacoef (Wm-2)
5009	!-- k = nzb_s_inner(j,i)+1
5010	!-- pt(k,j,i) = pt(k,j,i) + aheat(j,i)acoefdt_3d/(exn(k)rho_cpdz)
5011	!-- Instead of this, we can adjust shf in case AH only at surface
5012	surf_usm_h%shf(m) = surf_usm_h%shf(m) + &
5013	aheat(j,i) * acoef * ddx * ddy / rho_cp
5014	ENDIF
5015	ENDDO
5016
5017	ENDIF
5018
5019	!-- pt and shf are defined on nxlg:nxrg,nysg:nyng
5020	!-- get the borders from neighbours
5021	CALL exchange_horiz( pt, nbgp )
5022
5023
5024	!-- calculation of force_radiation_call:
5025	!-- Make logical OR for all processes.
5026	!-- Force radiation call if at least one processor forces it.
5027	IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )&
5028	THEN
5029	#if defined( __parallel )
5030	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
5031	CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, &
5032	1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
5033	#else
5034	force_radiation_call = force_radiation_call_l
5035	#endif
5036	force_radiation_call_l = .FALSE.
5037	ENDIF
5038
5039	END SUBROUTINE usm_surface_energy_balance
5040
5041
5042	!------------------------------------------------------------------------------!
5043	! Description:
5044	! ------------
5045	!> Swapping of timelevels for t_surf and t_wall
5046	!> called out from subroutine swap_timelevel
5047	!------------------------------------------------------------------------------!
5048	SUBROUTINE usm_swap_timelevel ( mod_count )
5049
5050	IMPLICIT NONE
5051
5052	INTEGER(iwp), INTENT(IN) :: mod_count
5053	INTEGER(iwp) :: i
5054
5055	#if defined( __nopointer )
5056	t_surf_h = t_surf_h_p
5057	t_wall_h = t_wall_h_p
5058	t_surf_v = t_surf_v_p
5059	t_wall_v = t_wall_v_p
5060	#else
5061	SELECT CASE ( mod_count )
5062	CASE ( 0 )
5063	!
5064	!-- Horizontal surfaces
5065	t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2
5066	t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2
5067	!
5068	!-- Vertical surfaces
5069	t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
5070	t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2
5071	CASE ( 1 )
5072	!
5073	!-- Horizontal surfaces
5074	t_surf_h => t_surf_h_2; t_surf_h_p => t_surf_h_1
5075	t_wall_h => t_wall_h_2; t_wall_h_p => t_wall_h_1
5076	!
5077	!-- Vertical surfaces
5078	t_surf_v => t_surf_v_2; t_surf_v_p => t_surf_v_1
5079	t_wall_v => t_wall_v_2; t_wall_v_p => t_wall_v_1
5080	END SELECT
5081	#endif
5082
5083	END SUBROUTINE usm_swap_timelevel
5084
5085	!------------------------------------------------------------------------------!
5086	!
5087	! Description:
5088	! ------------
5089	!> Subroutine writes t_surf and t_wall data into restart files
5090	!kanani: Renamed this routine according to corresponging routines in PALM
5091	!kanani: Modified the routine to match write_var_list, from where usm_write_restart_data
5092	! shall be called in the future. This part has not been tested yet. (see virtual_flight_mod)
5093	! Also, I had some trouble with the allocation of t_surf, since this is a pointer.
5094	! So, I added some directives here.
5095	!------------------------------------------------------------------------------!
5096	SUBROUTINE usm_write_restart_data
5097
5098	IMPLICIT NONE
5099
5100	INTEGER(iwp) :: i
5101
5102	DO i = 0, io_blocks-1
5103	IF ( i == io_group ) THEN
5104
5105	WRITE ( 14 ) 't_surf_h '
5106	#if defined( __nopointer )
5107	WRITE ( 14 ) t_surf_h
5108	#else
5109	WRITE ( 14 ) t_surf_h_1
5110	#endif
5111	WRITE ( 14 ) 't_surf_v(0) '
5112	#if defined( __nopointer )
5113	WRITE ( 14 ) t_surf_v(0)%t
5114	#else
5115	WRITE ( 14 ) t_surf_v_1(0)%t
5116	#endif
5117	WRITE ( 14 ) 't_surf_v(1) '
5118	#if defined( __nopointer )
5119	WRITE ( 14 ) t_surf_v(1)%t
5120	#else
5121	WRITE ( 14 ) t_surf_v_1(1)%t
5122	#endif
5123	WRITE ( 14 ) 't_surf_v(2) '
5124	#if defined( __nopointer )
5125	WRITE ( 14 ) t_surf_v(2)%t
5126	#else
5127	WRITE ( 14 ) t_surf_v_1(2)%t
5128	#endif
5129	WRITE ( 14 ) 't_surf_v(3) '
5130	#if defined( __nopointer )
5131	WRITE ( 14 ) t_surf_v(3)%t
5132	#else
5133	WRITE ( 14 ) t_surf_v_1(3)%t
5134	#endif
5135	WRITE ( 14 ) 't_wall_h '
5136	#if defined( __nopointer )
5137	WRITE ( 14 ) t_wall_h
5138	#else
5139	WRITE ( 14 ) t_wall_h_1
5140	#endif
5141	WRITE ( 14 ) 't_wall_v(0) '
5142	#if defined( __nopointer )
5143	WRITE ( 14 ) t_wall_v(0)%t
5144	#else
5145	WRITE ( 14 ) t_wall_v_1(0)%t
5146	#endif
5147	WRITE ( 14 ) 't_wall_v(1) '
5148	#if defined( __nopointer )
5149	WRITE ( 14 ) t_wall_v(1)%t
5150	#else
5151	WRITE ( 14 ) t_wall_v_1(1)%t
5152	#endif
5153	WRITE ( 14 ) 't_wall_v(2) '
5154	#if defined( __nopointer )
5155	WRITE ( 14 ) t_wall_v(2)%t
5156	#else
5157	WRITE ( 14 ) t_wall_v_1(2)%t
5158	#endif
5159	WRITE ( 14 ) 't_wall_v(3) '
5160	#if defined( __nopointer )
5161	WRITE ( 14 ) t_wall_v(3)%t
5162	#else
5163	WRITE ( 14 ) t_wall_v_1(3)%t
5164	#endif
5165	WRITE ( 14 ) '* end usm * '
5166	ENDIF
5167	#if defined( __parallel )
5168	CALL MPI_BARRIER( comm2d, ierr )
5169	#endif
5170	ENDDO
5171
5172
5173	END SUBROUTINE usm_write_restart_data
5174
5175
5176	!------------------------------------------------------------------------------!
5177	!
5178	! Description:
5179	! ------------
5180	!> Subroutine stores svf, svfsurf, csf and csfsurf data to a file.
5181	!------------------------------------------------------------------------------!
5182	SUBROUTINE usm_write_svf_to_file
5183
5184	IMPLICIT NONE
5185	INTEGER(iwp) :: fsvf = 89
5186	INTEGER(iwp) :: i
5187
5188	DO i = 0, io_blocks-1
5189	IF ( i == io_group ) THEN
5190	OPEN ( fsvf, file=TRIM(svf_file_name)//TRIM(coupling_char)//myid_char, &
5191	form='unformatted', status='new' )
5192
5193	WRITE ( fsvf ) usm_version
5194	WRITE ( fsvf ) nsvfl, ncsfl
5195	WRITE ( fsvf ) svf
5196	WRITE ( fsvf ) svfsurf
5197	IF ( plant_canopy ) THEN
5198	WRITE ( fsvf ) csf
5199	WRITE ( fsvf ) csfsurf
5200	ENDIF
5201	WRITE ( fsvf ) TRIM(svf_code)
5202
5203	CLOSE (fsvf)
5204	#if defined( __parallel )
5205	CALL MPI_BARRIER( comm2d, ierr )
5206	#endif
5207	ENDIF
5208	ENDDO
5209	END SUBROUTINE usm_write_svf_to_file
5210
5211
5212	!------------------------------------------------------------------------------!
5213	!
5214	! Description:
5215	! ------------
5216	!> Block of auxiliary subroutines:
5217	!> 1. quicksort and corresponding comparison
5218	!> 2. usm_merge_and_grow_csf for implementation of "dynamical growing"
5219	!> array for csf
5220	!------------------------------------------------------------------------------!
5221	PURE FUNCTION svf_lt(svf1,svf2) result (res)
5222	TYPE (t_svf), INTENT(in) :: svf1,svf2
5223	LOGICAL :: res
5224	IF ( svf1%isurflt < svf2%isurflt .OR. &
5225	(svf1%isurflt == svf2%isurflt .AND. svf1%isurfs < svf2%isurfs) ) THEN
5226	res = .TRUE.
5227	ELSE
5228	res = .FALSE.
5229	ENDIF
5230	END FUNCTION svf_lt
5231
5232
5233	!-- quicksort.f --f90--
5234	!-- Author: t-nissie, adaptation J.Resler
5235	!-- License: GPLv3
5236	!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
5237	RECURSIVE SUBROUTINE quicksort_svf(svfl, first, last)
5238	IMPLICIT NONE
5239	TYPE(t_svf), DIMENSION(:), INTENT(INOUT) :: svfl
5240	INTEGER(iwp), INTENT(IN) :: first, last
5241	TYPE(t_svf) :: x, t
5242	INTEGER(iwp) :: i, j
5243
5244	IF ( first>=last ) RETURN
5245	x = svfl( (first+last) / 2 )
5246	i = first
5247	j = last
5248	DO
5249	DO while ( svf_lt(svfl(i),x) )
5250	i=i+1
5251	ENDDO
5252	DO while ( svf_lt(x,svfl(j)) )
5253	j=j-1
5254	ENDDO
5255	IF ( i >= j ) EXIT
5256	t = svfl(i); svfl(i) = svfl(j); svfl(j) = t
5257	i=i+1
5258	j=j-1
5259	ENDDO
5260	IF ( first < i-1 ) CALL quicksort_svf(svfl, first, i-1)
5261	IF ( j+1 < last ) CALL quicksort_svf(svfl, j+1, last)
5262	END SUBROUTINE quicksort_svf
5263
5264
5265	PURE FUNCTION csf_lt(csf1,csf2) result (res)
5266	TYPE (t_csf), INTENT(in) :: csf1,csf2
5267	LOGICAL :: res
5268	IF ( csf1%ip < csf2%ip .OR. &
5269	(csf1%ip == csf2%ip .AND. csf1%itx < csf2%itx) .OR. &
5270	(csf1%ip == csf2%ip .AND. csf1%itx == csf2%itx .AND. csf1%ity < csf2%ity) .OR. &
5271	(csf1%ip == csf2%ip .AND. csf1%itx == csf2%itx .AND. csf1%ity == csf2%ity .AND. &
5272	csf1%itz < csf2%itz) .OR. &
5273	(csf1%ip == csf2%ip .AND. csf1%itx == csf2%itx .AND. csf1%ity == csf2%ity .AND. &
5274	csf1%itz == csf2%itz .AND. csf1%isurfs < csf2%isurfs) ) THEN
5275	res = .TRUE.
5276	ELSE
5277	res = .FALSE.
5278	ENDIF
5279	END FUNCTION csf_lt
5280
5281
5282	!-- quicksort.f --f90--
5283	!-- Author: t-nissie, adaptation J.Resler
5284	!-- License: GPLv3
5285	!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
5286	RECURSIVE SUBROUTINE quicksort_csf(csfl, first, last)
5287	IMPLICIT NONE
5288	TYPE(t_csf), DIMENSION(:), INTENT(INOUT) :: csfl
5289	INTEGER(iwp), INTENT(IN) :: first, last
5290	TYPE(t_csf) :: x, t
5291	INTEGER(iwp) :: i, j
5292
5293	IF ( first>=last ) RETURN
5294	x = csfl( (first+last)/2 )
5295	i = first
5296	j = last
5297	DO
5298	DO while ( csf_lt(csfl(i),x) )
5299	i=i+1
5300	ENDDO
5301	DO while ( csf_lt(x,csfl(j)) )
5302	j=j-1
5303	ENDDO
5304	IF ( i >= j ) EXIT
5305	t = csfl(i); csfl(i) = csfl(j); csfl(j) = t
5306	i=i+1
5307	j=j-1
5308	ENDDO
5309	IF ( first < i-1 ) CALL quicksort_csf(csfl, first, i-1)
5310	IF ( j+1 < last ) CALL quicksort_csf(csfl, j+1, last)
5311	END SUBROUTINE quicksort_csf
5312
5313
5314	SUBROUTINE usm_merge_and_grow_csf(newsize)
5315	INTEGER(iwp), INTENT(in) :: newsize !< new array size after grow, must be >= ncsfl
5316	!< or -1 to shrink to minimum
5317	INTEGER(iwp) :: iread, iwrite
5318	TYPE(t_csf), DIMENSION(:), POINTER :: acsfnew
5319
5320	IF ( newsize == -1 ) THEN
5321	!-- merge in-place
5322	acsfnew => acsf
5323	ELSE
5324	!-- allocate new array
5325	IF ( mcsf == 0 ) THEN
5326	ALLOCATE( acsf1(newsize) )
5327	acsfnew => acsf1
5328	ELSE
5329	ALLOCATE( acsf2(newsize) )
5330	acsfnew => acsf2
5331	ENDIF
5332	ENDIF
5333
5334	IF ( ncsfl >= 1 ) THEN
5335	!-- sort csf in place (quicksort)
5336	CALL quicksort_csf(acsf,1,ncsfl)
5337
5338	!-- while moving to a new array, aggregate canopy sink factor records with identical box & source
5339	acsfnew(1) = acsf(1)
5340	iwrite = 1
5341	DO iread = 2, ncsfl
5342	!-- here acsf(kcsf) already has values from acsf(icsf)
5343	IF ( acsfnew(iwrite)%itx == acsf(iread)%itx &
5344	.AND. acsfnew(iwrite)%ity == acsf(iread)%ity &
5345	.AND. acsfnew(iwrite)%itz == acsf(iread)%itz &
5346	.AND. acsfnew(iwrite)%isurfs == acsf(iread)%isurfs ) THEN
5347	!-- We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
5348	!-- probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
5349	!-- might mean that the traced beam passes longer through the canopy box.
5350	IF ( acsfnew(iwrite)%rsvf < acsf(iread)%rsvf ) THEN
5351	acsfnew(iwrite)%rtransp = acsf(iread)%rtransp
5352	ENDIF
5353	acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf
5354	!-- advance reading index, keep writing index
5355	ELSE
5356	!-- not identical, just advance and copy
5357	iwrite = iwrite + 1
5358	acsfnew(iwrite) = acsf(iread)
5359	ENDIF
5360	ENDDO
5361	ncsfl = iwrite
5362	ENDIF
5363
5364	IF ( newsize == -1 ) THEN
5365	!-- allocate new array and copy shrinked data
5366	IF ( mcsf == 0 ) THEN
5367	ALLOCATE( acsf1(ncsfl) )
5368	acsf1(1:ncsfl) = acsf2(1:ncsfl)
5369	ELSE
5370	ALLOCATE( acsf2(ncsfl) )
5371	acsf2(1:ncsfl) = acsf1(1:ncsfl)
5372	ENDIF
5373	ENDIF
5374
5375	!-- deallocate old array
5376	IF ( mcsf == 0 ) THEN
5377	mcsf = 1
5378	acsf => acsf1
5379	DEALLOCATE( acsf2 )
5380	ELSE
5381	mcsf = 0
5382	acsf => acsf2
5383	DEALLOCATE( acsf1 )
5384	ENDIF
5385	ncsfla = newsize
5386	END SUBROUTINE usm_merge_and_grow_csf
5387
5388
5389	!-- quicksort.f --f90--
5390	!-- Author: t-nissie, adaptation J.Resler
5391	!-- License: GPLv3
5392	!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
5393	RECURSIVE SUBROUTINE quicksort_csf2(kpcsflt, pcsflt, first, last)
5394	IMPLICIT NONE
5395	INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT) :: kpcsflt
5396	REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: pcsflt
5397	INTEGER(iwp), INTENT(IN) :: first, last
5398	REAL(wp), DIMENSION(ndcsf) :: t2
5399	INTEGER(iwp), DIMENSION(kdcsf) :: x, t1
5400	INTEGER(iwp) :: i, j
5401
5402	IF ( first>=last ) RETURN
5403	x = kpcsflt(:, (first+last)/2 )
5404	i = first
5405	j = last
5406	DO
5407	DO while ( csf_lt2(kpcsflt(:,i),x) )
5408	i=i+1
5409	ENDDO
5410	DO while ( csf_lt2(x,kpcsflt(:,j)) )
5411	j=j-1
5412	ENDDO
5413	IF ( i >= j ) EXIT
5414	t1 = kpcsflt(:,i); kpcsflt(:,i) = kpcsflt(:,j); kpcsflt(:,j) = t1
5415	t2 = pcsflt(:,i); pcsflt(:,i) = pcsflt(:,j); pcsflt(:,j) = t2
5416	i=i+1
5417	j=j-1
5418	ENDDO
5419	IF ( first < i-1 ) CALL quicksort_csf2(kpcsflt, pcsflt, first, i-1)
5420	IF ( j+1 < last ) CALL quicksort_csf2(kpcsflt, pcsflt, j+1, last)
5421	END SUBROUTINE quicksort_csf2
5422
5423
5424	PURE FUNCTION csf_lt2(item1, item2) result(res)
5425	INTEGER(iwp), DIMENSION(kdcsf), INTENT(in) :: item1, item2
5426	LOGICAL :: res
5427	res = ( (item1(3) < item2(3)) &
5428	.OR. (item1(3) == item2(3) .AND. item1(2) < item2(2)) &
5429	.OR. (item1(3) == item2(3) .AND. item1(2) == item2(2) .AND. item1(1) < item2(1)) &
5430	.OR. (item1(3) == item2(3) .AND. item1(2) == item2(2) .AND. item1(1) == item2(1) &
5431	.AND. item1(4) < item2(4)) )
5432	END FUNCTION csf_lt2
5433
5434
5435	END MODULE urban_surface_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |