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

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

Enable restarts with USM with different number of PEs; some bugfixes in new surface structure in USM; formatting adjustments and descriptions in surface_mod

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