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

Last change on this file since 2114 was 2114, checked in by kanani, 5 years ago

last commit documented

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