Home

Context Navigation

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

Last change on this file since 2296 was 2296, checked in by maronga, 8 years ago
added new spinup mechanism for surface/radiation models
Property svn:keywords set to `Id`
File size: 265.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |