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

Last change on this file since 2287 was 2287, checked in by suehring, 4 years ago

Bugfix in determination topography-top index

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