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

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

last commit documented

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