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

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

changes in the course of urban surface model implementation

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