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

Last change on this file since 2723 was 2723, checked in by maronga, 4 years ago

bugfixes for spinup mechanism to work with lsm+usm+radiation

  • Property svn:keywords set to Id
File size: 423.3 KB
Line 
1!> @file urban_surface_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2015-2018 Czech Technical University in Prague
18! Copyright 1997-2018 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: urban_surface_mod.f90 2723 2018-01-05 09:27:03Z maronga $
28! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
29!
30! 2720 2018-01-02 16:27:15Z kanani
31! Correction of comment
32!
33! 2718 2018-01-02 08:49:38Z maronga
34! Corrected "Former revisions" section
35!
36! 2705 2017-12-18 11:26:23Z maronga
37! Changes from last commit documented
38!
39! 2703 2017-12-15 20:12:38Z maronga
40! Workaround for calculation of r_a
41!
42! 2696 2017-12-14 17:12:51Z kanani
43! - Change in file header (GPL part)
44! - Bugfix in calculation of pt_surface and related fluxes. (BM)
45! - Do not write surface temperatures onto pt array as this might cause
46!   problems with nesting. (MS)
47! - Revised calculation of pt1 (now done in surface_layer_fluxes).
48!   Bugfix, f_shf_window and f_shf_green were not set at vertical surface
49!   elements. (MS)
50! - merged with branch ebsolver
51!   green building surfaces do not evaporate yet
52!   properties of green wall layers and window layers are taken from wall layers
53!   this input data is missing. (RvT)
54! - Merged with branch radiation (developed by Mohamed Salim)
55! - Revised initialization. (MS)
56! - Rename emiss_surf into emissivity, roughness_wall into z0, albedo_surf into
57!   albedo. (MS)
58! - Move first call of usm_radiatin from usm_init to init_3d_model
59! - fixed problem with near surface temperature
60! - added near surface temperature t_surf_10cm_h(m), t_surf_10cm_v(l)%t(m)
61! - does not work with temp profile including stability, ol
62!   t_surf_10cm = pt1 now
63! - merged with 2357 bugfix, error message for nopointer version
64! - added indoor model coupling with wall heat flux
65! - added green substrate/ dry vegetation layer for buildings
66! - merged with 2232 new surface-type structure
67! - added transmissivity of window tiles
68! - added MOSAIK tile approach for 3 different surfaces (RvT)
69!
70! 2583 2017-10-26 13:58:38Z knoop
71! Bugfix: reverted MPI_Win_allocate_cptr introduction in last commit
72!
73! 2582 2017-10-26 13:19:46Z hellstea
74! Workaround for gnufortran compiler added in usm_calc_svf. CALL MPI_Win_allocate is
75! replaced by CALL MPI_Win_allocate_cptr if defined ( __gnufortran ).
76!
77! 2544 2017-10-13 18:09:32Z maronga
78! Date and time quantities are now read from date_and_time_mod. Solar constant is
79! read from radiation_model_mod
80!
81! 2516 2017-10-04 11:03:04Z suehring
82! Remove tabs
83!
84! 2514 2017-10-04 09:52:37Z suehring
85! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
86! no output of ghost layer data
87!
88! 2350 2017-08-15 11:48:26Z kanani
89! Bugfix and error message for nopointer version.
90! Additional "! defined(__nopointer)" as workaround to enable compilation of
91! nopointer version.
92!
93! 2318 2017-07-20 17:27:44Z suehring
94! Get topography top index via Function call
95!
96! 2317 2017-07-20 17:27:19Z suehring
97! Bugfix: adjust output of shf. Added support for spinups
98!
99! 2287 2017-06-15 16:46:30Z suehring
100! Bugfix in determination topography-top index
101!
102! 2269 2017-06-09 11:57:32Z suehring
103! Enable restart runs with different number of PEs
104! Bugfixes nopointer branch
105!
106! 2258 2017-06-08 07:55:13Z suehring
107! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
108!
109! 2233 2017-05-30 18:08:54Z suehring
110!
111! 2232 2017-05-30 17:47:52Z suehring
112! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
113! insteat, heat fluxes are directly applied in diffusion_s.
114!
115! 2213 2017-04-24 15:10:35Z kanani
116! Removal of output quantities usm_lad and usm_canopy_hr
117!
118! 2209 2017-04-19 09:34:46Z kanani
119! cpp switch __mpi3 removed,
120! minor formatting,
121! small bugfix for division by zero (Krc)
122!
123! 2113 2017-01-12 13:40:46Z kanani
124! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
125!
126! 2071 2016-11-17 11:22:14Z maronga
127! Small bugfix (Resler)
128!
129! 2031 2016-10-21 15:11:58Z knoop
130! renamed variable rho to rho_ocean
131!
132! 2024 2016-10-12 16:42:37Z kanani
133! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
134! optimization of MPI-RMA operations,
135! declaration of pcbl as integer,
136! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
137! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
138! use of new control parameter varnamelength,
139! added output variables usm_rad_ressw, usm_rad_reslw,
140! minor formatting changes,
141! minor optimizations.
142!
143! 2011 2016-09-19 17:29:57Z kanani
144! Major reformatting according to PALM coding standard (comments, blanks,
145! alphabetical ordering, etc.),
146! removed debug_prints,
147! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
148! defined in MODULE control_parameters (modules.f90) to avoid circular
149! dependencies,
150! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
151!
152! 2007 2016-08-24 15:47:17Z kanani
153! Initial revision
154!
155!
156! Description:
157! ------------
158! 2016/6/9 - Initial version of the USM (Urban Surface Model)
159!            authors: Jaroslav Resler, Pavel Krc
160!                     (Czech Technical University in Prague and Institute of
161!                      Computer Science of the Czech Academy of Sciences, Prague)
162!            with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
163!            partly inspired by PALM LSM (B. Maronga)
164!            parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
165!> Module for Urban Surface Model (USM)
166!> The module includes:
167!>    1. radiation model with direct/diffuse radiation, shading, reflections
168!>       and integration with plant canopy
169!>    2. wall and wall surface model
170!>    3. surface layer energy balance
171!>    4. anthropogenic heat (only from transportation so far)
172!>    5. necessary auxiliary subroutines (reading inputs, writing outputs,
173!>       restart simulations, ...)
174!> It also make use of standard radiation and integrates it into
175!> urban surface model.
176!>
177!> Further work:
178!> -------------
179!> 1. Reduce number of shape view factors by merging factors for distant surfaces
180!>    under shallow angles. Idea: Iteratively select the smallest shape view
181!>    factor by value (among all sources and targets) which has a similarly
182!>    oriented source neighbor (or near enough) SVF and merge them by adding
183!>    value of the smaller SVF to the larger one and deleting the smaller one.
184!>    This will allow for better scaling at higher resolutions.
185!>
186!> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity
187!>    from surfaces that are visible from local surfaces (i.e. there is a SVF
188!>    where target is local). To do that, radiosity will be exchanged after each
189!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
190!>
191!> 3. Temporarily large values of surface heat flux can be observed, up to
192!>    1.2 Km/s, which seem to be not realistic.
193!>
194!> @todo Output of _av variables in case of restarts
195!> @todo Revise flux conversion in energy-balance solver
196!> @todo Bugfixing in nopointer branch
197!> @todo Check optimizations for RMA operations
198!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
199!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
200!>       factor 3 between min and max time
201!> @todo Move setting of flag indoor_model to indoor_model_mod once available
202!> @todo Check divisions in wtend (etc.) calculations for possible division
203!>       by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
204!> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
205!------------------------------------------------------------------------------!
206 MODULE urban_surface_mod
207
208#if ! defined( __nopointer )
209    USE arrays_3d,                                                             &
210        ONLY:  zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
211#endif
212
213    USE cloud_parameters,                                                      &
214        ONLY:  cp, r_d
215
216    USE constants,                                                             &
217        ONLY:  pi
218   
219    USE control_parameters,                                                    &
220        ONLY:  coupling_start_time, dz, topography, dt_3d,                     &
221               intermediate_timestep_count, initializing_actions,              &
222               intermediate_timestep_count_max, simulated_time, end_time,      &
223               timestep_scheme, tsc, coupling_char, io_blocks, io_group,       &
224               message_string, time_since_reference_point, surface_pressure,   &
225               g, pt_surface, large_scale_forcing, lsf_surf, spinup,           &
226               spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
227               average_count_3d, varnamelength, urban_surface, kappa
228
229    USE cpulog,                                                                &
230        ONLY:  cpu_log, log_point, log_point_s
231
232    USE date_and_time_mod,                                                     &
233        ONLY:  d_seconds_year, day_of_year_init, time_utc_init
234
235    USE grid_variables,                                                        &
236        ONLY:  dx, dy, ddx, ddy, ddx2, ddy2
237
238    USE indices,                                                               &
239        ONLY:  nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys,    &
240               nysg, nzb, nzt, nbgp, wall_flags_0
241
242    USE, INTRINSIC :: iso_c_binding
243
244    USE kinds
245             
246    USE pegrid
247   
248    USE plant_canopy_model_mod,                                                &
249        ONLY:  pc_heating_rate, plant_canopy, usm_lad_rma
250   
251    USE radiation_model_mod,                                                   &
252        ONLY:  albedo_type, radiation, calc_zenith, zenith,                    &
253               rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,                   &
254               sigma_sb, solar_constant, sun_direction, sun_dir_lat,           &
255               sun_dir_lon,                                                    &
256               force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
257               surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
258               surfinl, surfinlwdif, energy_balance_surf_h,                    &
259               energy_balance_surf_v, rad_sw_in_dir, rad_sw_in_diff,           &
260               rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
261               surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, pcbl, npcbl, startenergy,      &
262               endenergy, nenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,  &
263               inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,id,      &
264               iz, iy, ix, idir, jdir, kdir, startborder, endborder, nsurf_type, nzub, nzut,   &
265               isky, inorth_b,idown_a, isouth_b, ieast_b, iwest_b, nzu, pch, nsurf,  &
266               iup_a, inorth_a, isouth_a, ieast_a, iwest_a, idsvf, ndsvf,      &
267               idcsf, ndcsf, kdcsf, pct, startland, endland, startwall, endwall
268
269    USE statistics,                                                            &
270        ONLY:  hom, statistic_regions
271
272    USE surface_mod
273
274
275    IMPLICIT NONE
276
277
278!-- configuration parameters (they can be setup in PALM config)
279    LOGICAL                                        ::  usm_material_model = .TRUE.        !< flag parameter indicating wheather the  model of heat in materials is used
280    LOGICAL                                        ::  usm_anthropogenic_heat = .FALSE.   !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used
281    LOGICAL                                        ::  force_radiation_call_l = .FALSE.   !< flag parameter for unscheduled radiation model calls
282    LOGICAL                                        ::  indoor_model = .FALSE.              !< whether to use the indoor model
283
284   
285    INTEGER(iwp)                                   ::  building_type = 1                  !< default building type (preleminary setting)
286    INTEGER(iwp)                                   ::  land_category = 2                  !< default category for land surface
287    INTEGER(iwp)                                   ::  wall_category = 2                  !< default category for wall surface over pedestrian zone
288    INTEGER(iwp)                                   ::  pedestrant_category = 2            !< default category for wall surface in pedestrian zone
289    INTEGER(iwp)                                   ::  roof_category = 2                  !< default category for root surface
290
291!
292!-- Indices of input attributes for (above) ground floor level
293    INTEGER(iwp) ::  ind_alb_wall          = 38 !< index in input list for albedo_type of wall fraction
294    INTEGER(iwp) ::  ind_alb_green         = 39 !< index in input list for albedo_type of green fraction
295    INTEGER(iwp) ::  ind_alb_win           = 40 !< index in input list for albedo_type of window fraction
296    INTEGER(iwp) ::  ind_emis_wall_agfl    = 14 !< index in input list for wall emissivity, above ground floor level
297    INTEGER(iwp) ::  ind_emis_wall_gfl     = 32 !< index in input list for wall emissivity, ground floor level
298    INTEGER(iwp) ::  ind_emis_green_agfl   = 15 !< index in input list for green emissivity, above ground floor level
299    INTEGER(iwp) ::  ind_emis_green_gfl    = 33 !< index in input list for green emissivity, ground floor level
300    INTEGER(iwp) ::  ind_emis_win_agfl     = 16 !< index in input list for window emissivity, above ground floor level
301    INTEGER(iwp) ::  ind_emis_win_gfl      = 34 !< index in input list for window emissivity, ground floor level
302    INTEGER(iwp) ::  ind_green_frac_w_agfl = 2  !< index in input list for green fraction on wall, above ground floor level
303    INTEGER(iwp) ::  ind_green_frac_w_gfl  = 23 !< index in input list for green fraction on wall, ground floor level
304    INTEGER(iwp) ::  ind_green_frac_r_agfl = 3  !< index in input list for green fraction on roof, above ground floor level
305    INTEGER(iwp) ::  ind_green_frac_r_gfl  = 24 !< index in input list for green fraction on roof, ground floor level
306    INTEGER(iwp) ::  ind_hc1_agfl          =  6 !< index in input list for heat capacity at first wall layer, above ground floor level
307    INTEGER(iwp) ::  ind_hc1_gfl           = 26 !< index in input list for heat capacity at first wall layer, ground floor level
308    INTEGER(iwp) ::  ind_hc2_agfl          = 7  !< index in input list for heat capacity at second wall layer, above ground floor level
309    INTEGER(iwp) ::  ind_hc2_gfl           = 27 !< index in input list for heat capacity at second wall layer, ground floor level
310    INTEGER(iwp) ::  ind_hc3_agfl          = 8  !< index in input list for heat capacity at third wall layer, above ground floor level
311    INTEGER(iwp) ::  ind_hc3_gfl           = 28 !< index in input list for heat capacity at third wall layer, ground floor level
312    INTEGER(iwp) ::  ind_gflh              = 20 !< index in input list for ground floor level height
313    INTEGER(iwp) ::  ind_lai_r_agfl        = 4  !< index in input list for LAI on roof, above ground floor level
314    INTEGER(iwp) ::  ind_lai_r_gfl         = 4  !< index in input list for LAI on roof, ground floor level
315    INTEGER(iwp) ::  ind_lai_w_agfl        = 5  !< index in input list for LAI on wall, above ground floor level
316    INTEGER(iwp) ::  ind_lai_w_gfl         = 25 !< index in input list for LAI on wall, ground floor level
317    INTEGER(iwp) ::  ind_tc1_agfl          = 9  !< index in input list for thermal conductivity at first wall layer, above ground floor level
318    INTEGER(iwp) ::  ind_tc1_gfl           = 29 !< index in input list for thermal conductivity at first wall layer, ground floor level
319    INTEGER(iwp) ::  ind_tc2_agfl          = 10 !< index in input list for thermal conductivity at second wall layer, above ground floor level
320    INTEGER(iwp) ::  ind_tc2_gfl           = 30 !< index in input list for thermal conductivity at second wall layer, ground floor level
321    INTEGER(iwp) ::  ind_tc3_agfl          = 11 !< index in input list for thermal conductivity at third wall layer, above ground floor level
322    INTEGER(iwp) ::  ind_tc3_gfl           = 31 !< index in input list for thermal conductivity at third wall layer, ground floor level
323    INTEGER(iwp) ::  ind_thick_1           = 41 !< index for wall layer thickness - 1st layer
324    INTEGER(iwp) ::  ind_thick_2           = 42 !< index for wall layer thickness - 2nd layer
325    INTEGER(iwp) ::  ind_thick_3           = 43 !< index for wall layer thickness - 3rd layer
326    INTEGER(iwp) ::  ind_thick_4           = 44 !< index for wall layer thickness - 4th layer
327    INTEGER(iwp) ::  ind_trans_agfl        = 17 !< index in input list for window transmissivity, above ground floor level
328    INTEGER(iwp) ::  ind_trans_gfl         = 35 !< index in input list for window transmissivity, ground floor level
329    INTEGER(iwp) ::  ind_wall_frac_agfl    = 0  !< index in input list for wall fraction, above ground floor level
330    INTEGER(iwp) ::  ind_wall_frac_gfl     = 21 !< index in input list for wall fraction, ground floor level
331    INTEGER(iwp) ::  ind_win_frac_agfl     = 1  !< index in input list for window fraction, above ground floor level
332    INTEGER(iwp) ::  ind_win_frac_gfl      = 22 !< index in input list for window fraction, ground floor level
333    INTEGER(iwp) ::  ind_z0_agfl           = 18 !< index in input list for z0, above ground floor level
334    INTEGER(iwp) ::  ind_z0_gfl            = 36 !< index in input list for z0, ground floor level
335    INTEGER(iwp) ::  ind_z0qh_agfl         = 19 !< index in input list for z0h / z0q, above ground floor level
336    INTEGER(iwp) ::  ind_z0qh_gfl          = 37 !< index in input list for z0h / z0q, ground floor level
337
338
339    REAL(wp)  ::  roof_height_limit = 4._wp          !< height for distinguish between land surfaces and roofs
340    REAL(wp)  ::  ground_floor_level = 4.0_wp        !< default ground floor level
341    REAL(wp)  ::  ra_horiz_coef = 5.0_wp             !< mysterious coefficient for correction of overestimation
342                                                                                          !< of r_a for horizontal surfaces -> TODO
343
344    CHARACTER(37), DIMENSION(0:6), PARAMETER :: building_type_name = (/     &
345                                   'user-defined                         ', & !  0
346                                   'residential - 1950                   ', & !  1
347                                   'residential 1951 - 2000              ', & !  2
348                                   'residential 2001 -                   ', & !  3
349                                   'office - 1950                        ', & !  4
350                                   'office 1951 - 2000                   ', & !  5
351                                   'office 2001 -                        '  & !  6
352                                                                     /)
353!
354!-- building parameters, 4 different types
355!-- 0 - wall fraction, 1- window fraction, 2 - green fraction on wall, 3- green fraction
356!-- at roof, 4 - lai of green fraction at roof,  5 - lai of green fraction at wall,
357!-- 6 - heat capacity of wall layer 1, 7 - heat capacity of wall layer 2,
358!-- 8 - heat capacity of wall layer 3, 9 - thermal conductivity of wall layer 1,
359!-- 10 - thermal conductivity of wall layer 2, 11 - thermal conductivity of wall layer 3, 
360!-- 12 - indoor target summer temperature ( K ), 13 - indoor target winter temperature (K),
361!-- 14 - emissivity of wall fraction, 15 - emissivity of green fraction, 16 - emissivity of window fraction,
362!-- 17 - transmissivity of window fraction, 18 - z0, 19 - z0h/z0q, 20 - ground floor height,
363!-- 21 - ground floor wall fraction, 22 - ground floor window fraction, 23 ground floor green fraction,
364!-- 24 - ground floor green fraction on roof, 25 - ground floor lai of green fraction,
365!-- 26 - ground floor heat capacity of wall layer 1, 27 - ground floor heat capacity of wall layer 1,
366!-- 28 - ground floor heat capacity of wall layer 3, 29 - ground floor thermal conductivity of wall layer 1,
367!-- 30 - ground floor thermal conductivity of wall layer 2, 31 - ground floor thermal conductivity of wall layer 3,
368!-- 32 - ground floor emissivity of wall fraction, 33 - ground floor emissivity of green fraction,
369!-- 34 - ground floor emissivity of window fraction, 35 - ground floor transmissivity of window fraction,
370!-- 36 - ground floor z0, 37 - ground floor z0h/z0q, 38 - albedo type wall fraction
371!-- 39 - albedo type green fraction, 40 - albedo type window fraction
372!-- 41 - wall layer thickness - 1st layer, 42 - wall layer thickness - 2nd layer,
373!-- 43 - wall layer thickness - 3rd layer, 44 - wall layer thickness - 4th layer,
374!-- 45 - heat capacity of the wall surface, 46 - heat conductivity
375!-- Please note, only preleminary dummy values so far!
376    REAL(wp), DIMENSION(0:46,1:6), PARAMETER :: building_pars = RESHAPE( (/    &
377        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
378        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
379        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
380        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
381        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
382        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
383        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
384        0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
385        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
386        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
387        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 1
388        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
389        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
390        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
391        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
392        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
393        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
394        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
395        0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
396        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
397        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
398        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 2
399        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
400        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
401        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
402        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
403        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
404        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
405        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
406        0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
407        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
408        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
409        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 3
410        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
411        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
412        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
413        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
414        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
415        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
416        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
417        0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
418        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
419        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
420        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 4
421        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
422        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
423        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
424        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
425        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
426        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
427        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
428        0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
429        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
430        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
431        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 5
432        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
433        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
434        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
435        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
436        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
437        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
438        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
439        0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
440        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
441        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
442        20000.0_wp, 10.0_wp                                                    & !parameter 45-46 - end of type 6
443                                                                          /),  &
444                                                               (/47, 6/) )
445
446!
447!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
448    TYPE t_surf_vertical
449       REAL(wp), DIMENSION(:), ALLOCATABLE         :: t
450    END TYPE t_surf_vertical
451!
452!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
453    TYPE t_wall_vertical
454       REAL(wp), DIMENSION(:,:), ALLOCATABLE       :: t
455    END TYPE t_wall_vertical
456
457
458!-- arrays for time averages
459!-- Attention: the variable rad_net_av is also used in the 3d field variable in radiation_model_mod.f90. It may be better to rename it
460    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rad_net_av       !< average of rad_net_l
461    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
462    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
463    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
464    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
465    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
466    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
467    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
468    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
469    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
470    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
471    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
472    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface 
473    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_av       !< average of wghf_eb
474    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wshf_eb_av       !< average of wshf_eb
475    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  t_wall_av        !< Average of t_wall
476    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_green_av !< average of wghf_eb_green
477    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  t_green_av       !< Average of t_green
478    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_window_av !< average of wghf_eb_window
479    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  t_window_av      !< Average of t_window   
480   
481
482!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483!-- anthropogenic heat sources
484!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheat             !< daily average of anthropogenic heat (W/m2)
486    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  aheatprof         !< diurnal profile of anthropogenic heat
487
488!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
489!-- wall surface model
490!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
491!-- wall surface model constants
492    INTEGER(iwp), PARAMETER                        :: nzb_wall = 0       !< inner side of the wall model (to be switched)
493    INTEGER(iwp), PARAMETER                        :: nzt_wall = 3       !< outer side of the wall model (to be switched)
494    INTEGER(iwp), PARAMETER                        :: nzw = 4            !< number of wall layers (fixed for now)
495
496    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
497                                                                         !< normalized soil, wall and roof layer depths (m/m)
498!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default = (/0.33_wp, 0.66_wp, 1.0_wp /)
499    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
500!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.33_wp, 0.66_wp, 1.0_wp /)
501!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
502                                                                         !< normalized window layer depths (m/m)
503!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
504                                                                         !< normalized green layer depths (m/m)
505    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
506!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /)
507
508                                                                       
509    REAL(wp)                                       ::   wall_inner_temperature = 296.0_wp    !< temperature of the inner wall surface (~23 degrees C) (K)
510    REAL(wp)                                       ::   roof_inner_temperature = 296.0_wp    !< temperature of the inner roof surface (~23 degrees C) (K)
511    REAL(wp)                                       ::   soil_inner_temperature = 283.0_wp    !< temperature of the deep soil (~10 degrees C) (K)
512    REAL(wp)                                       ::   window_inner_temperature = 296.0_wp  !< temperature of the inner window surface (~23 degrees C) (K)
513
514!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
515!-- surface and material model variables for walls, ground, roofs
516!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
517    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn                !< normalized wall layer depths (m)
518    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn_window         !< normalized window layer depths (m)
519    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn_green          !< normalized green layer depths (m)
520
521#if defined( __nopointer )
522    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h           !< wall surface temperature (K) at horizontal walls
523    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_p         !< progn. wall surface temperature (K) at horizontal walls
524    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h    !< window surface temperature (K) at horizontal walls
525    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_p  !< progn. window surface temperature (K) at horizontal walls
526    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h     !< green surface temperature (K) at horizontal walls
527    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_p   !< progn. green surface temperature (K) at horizontal walls
528    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_whole_h     !< whole surface temperature (K) at horizontal walls
529    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_whole_h_p   !< progn. whole surface temperature (K) at horizontal walls
530    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h      !< near surface temperature (10cm) (K) at horizontal walls
531    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h_p    !< progn. near surface temperature (10cm) (K) at horizontal walls
532    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v
533    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v_p
534    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_window_v
535    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_window_v_p
536    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_green_v
537    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_green_v_p
538    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_whole_v
539    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_whole_v_p
540    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_10cm_v
541    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_10cm_v_p
542#else
543    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h
544    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h_p
545    REAL(wp), DIMENSION(:), POINTER                :: t_surf_window_h
546    REAL(wp), DIMENSION(:), POINTER                :: t_surf_window_h_p
547    REAL(wp), DIMENSION(:), POINTER                :: t_surf_green_h
548    REAL(wp), DIMENSION(:), POINTER                :: t_surf_green_h_p
549    REAL(wp), DIMENSION(:), POINTER                :: t_surf_whole_h
550    REAL(wp), DIMENSION(:), POINTER                :: t_surf_whole_h_p
551    REAL(wp), DIMENSION(:), POINTER                :: t_surf_10cm_h
552    REAL(wp), DIMENSION(:), POINTER                :: t_surf_10cm_h_p
553
554    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_1
555    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_2
556    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_1
557    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_2
558    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_1
559    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_2
560    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_whole_h_1
561    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_whole_h_2
562    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h_1
563    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h_2
564
565    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v
566    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v_p
567    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_window_v
568    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_window_v_p
569    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_green_v
570    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_green_v_p
571    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_whole_v
572    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_whole_v_p
573    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_10cm_v
574    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_10cm_v_p
575
576    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_1
577    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_2
578    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_window_v_1
579    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_window_v_2
580    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_green_v_1
581    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_green_v_2
582    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_whole_v_1
583    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_whole_v_2
584    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_10cm_v_1
585    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_10cm_v_2
586   
587#endif
588    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_av          !< average of wall surface temperature (K)
589    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_av   !< average of window surface temperature (K)
590    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_av    !< average of green wall surface temperature (K)
591    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_whole_av    !< average of whole wall surface temperature (K)
592    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_av    !< average of whole wall surface temperature (K)
593
594!-- Temporal tendencies for time stepping           
595    REAL(wp), DIMENSION(:), ALLOCATABLE            :: tt_surface_m       !< surface temperature tendency of wall (K)
596    REAL(wp), DIMENSION(:), ALLOCATABLE            :: tt_surface_window_m !< surface temperature tendency of window (K)
597    REAL(wp), DIMENSION(:), ALLOCATABLE            :: tt_surface_green_m !< surface temperature tendency of green wall (K)
598
599!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
600!-- Energy balance variables
601!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
602!-- parameters of the land, roof and wall surfaces
603
604#if defined( __nopointer )
605    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h             !< Wall temperature (K)
606    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_av          !< Average of t_wall
607    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_p           !< Prog. wall temperature (K)
608    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h           !< Window temperature (K)
609    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_av        !< Average of t_window
610    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_p         !< Prog. window temperature (K)
611    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h            !< Green temperature (K)
612    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_av         !< Average of t_green
613    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_p          !< Prog. green temperature (K)
614
615    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v             !< Wall temperature (K)
616    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av          !< Average of t_wall
617    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_p           !< Prog. wall temperature (K)
618    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v           !< Window temperature (K)
619    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_av        !< Average of t_window
620    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_p         !< Prog. window temperature (K)
621    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v            !< Green temperature (K)
622    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_av         !< Average of t_green
623    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_p          !< Prog. green temperature (K)
624#else
625    REAL(wp), DIMENSION(:,:), POINTER                :: t_wall_h, t_wall_h_p
626    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_av, t_wall_h_1, t_wall_h_2
627    REAL(wp), DIMENSION(:,:), POINTER                :: t_window_h, t_window_h_p
628    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_av, t_window_h_1, t_window_h_2
629    REAL(wp), DIMENSION(:,:), POINTER                :: t_green_h, t_green_h_p
630    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_av, t_green_h_1, t_green_h_2
631
632    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_wall_v, t_wall_v_p
633    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av, t_wall_v_1, t_wall_v_2
634    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_window_v, t_window_v_p
635    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_av, t_window_v_1, t_window_v_2
636    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_green_v, t_green_v_p
637    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_av, t_green_v_1, t_green_v_2
638#endif
639
640!-- Wall temporal tendencies for time stepping
641    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: tt_wall_m          !< t_wall prognostic array
642    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: tt_window_m        !< t_window prognostic array
643    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: tt_green_m         !< t_green prognostic array
644
645!-- Surface and material parameters classes (surface_type)
646!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
647    INTEGER(iwp)                                   :: n_surface_types      !< number of the wall type categories
648    INTEGER(iwp), PARAMETER                        :: n_surface_params = 8 !< number of parameters for each type of the wall
649    INTEGER(iwp), PARAMETER                        :: ialbedo  = 1         !< albedo of the surface
650    INTEGER(iwp), PARAMETER                        :: iemiss   = 2         !< emissivity of the surface
651    INTEGER(iwp), PARAMETER                        :: ilambdas = 3         !< heat conductivity λS between air and surface ( W m−2 K−1 )
652    INTEGER(iwp), PARAMETER                        :: irough   = 4         !< roughness relative to concrete
653    INTEGER(iwp), PARAMETER                        :: icsurf   = 5         !< Surface skin layer heat capacity (J m−2 K−1 )
654    INTEGER(iwp), PARAMETER                        :: ithick   = 6         !< thickness of the surface (wall, roof, land)  ( m )
655    INTEGER(iwp), PARAMETER                        :: irhoC    = 7         !< volumetric heat capacity rho*C of the material ( J m−3 K−1 )
656    INTEGER(iwp), PARAMETER                        :: ilambdah = 8         !< thermal conductivity λH of the wall (W m−1 K−1 )
657    CHARACTER(12), DIMENSION(:), ALLOCATABLE       :: surface_type_names   !< names of wall types (used only for reports)
658    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        :: surface_type_codes   !< codes of wall types
659    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: surface_params       !< parameters of wall types
660
661   
662!-- interfaces of subroutines accessed from outside of this module
663    INTERFACE usm_boundary_condition
664       MODULE PROCEDURE usm_boundary_condition
665    END INTERFACE usm_boundary_condition
666
667    INTERFACE usm_check_data_output
668       MODULE PROCEDURE usm_check_data_output
669    END INTERFACE usm_check_data_output
670   
671    INTERFACE usm_check_parameters
672       MODULE PROCEDURE usm_check_parameters
673    END INTERFACE usm_check_parameters
674   
675    INTERFACE usm_data_output_3d
676       MODULE PROCEDURE usm_data_output_3d
677    END INTERFACE usm_data_output_3d
678   
679    INTERFACE usm_define_netcdf_grid
680       MODULE PROCEDURE usm_define_netcdf_grid
681    END INTERFACE usm_define_netcdf_grid
682
683    INTERFACE usm_init_urban_surface
684       MODULE PROCEDURE usm_init_urban_surface
685    END INTERFACE usm_init_urban_surface
686
687    INTERFACE usm_material_heat_model
688       MODULE PROCEDURE usm_material_heat_model
689    END INTERFACE usm_material_heat_model
690   
691    INTERFACE usm_green_heat_model
692       MODULE PROCEDURE usm_green_heat_model
693    END INTERFACE usm_green_heat_model
694   
695    INTERFACE usm_parin
696       MODULE PROCEDURE usm_parin
697    END INTERFACE usm_parin
698   
699    INTERFACE usm_temperature_near_surface
700       MODULE PROCEDURE usm_temperature_near_surface
701    END INTERFACE usm_temperature_near_surface
702
703    INTERFACE usm_read_restart_data
704       MODULE PROCEDURE usm_read_restart_data
705    END INTERFACE usm_read_restart_data
706
707    INTERFACE usm_surface_energy_balance
708       MODULE PROCEDURE usm_surface_energy_balance
709    END INTERFACE usm_surface_energy_balance
710   
711    INTERFACE usm_swap_timelevel
712       MODULE PROCEDURE usm_swap_timelevel
713    END INTERFACE usm_swap_timelevel
714       
715    INTERFACE usm_write_restart_data
716       MODULE PROCEDURE usm_write_restart_data
717    END INTERFACE usm_write_restart_data
718
719    INTERFACE usm_allocate_surface
720       MODULE PROCEDURE usm_allocate_surface
721    END INTERFACE usm_allocate_surface
722
723    INTERFACE usm_average_3d_data
724       MODULE PROCEDURE usm_average_3d_data
725    END INTERFACE usm_average_3d_data
726
727   
728    SAVE
729
730    PRIVATE 
731   
732!-- Public functions
733    PUBLIC usm_boundary_condition, usm_check_parameters, usm_init_urban_surface,&
734           usm_read_restart_data,                                               & 
735           usm_surface_energy_balance, usm_material_heat_model,                 &
736           usm_swap_timelevel, usm_check_data_output, usm_average_3d_data,      &
737           usm_data_output_3d, usm_define_netcdf_grid, usm_parin,               &
738           usm_write_restart_data, usm_allocate_surface
739
740!-- Public parameters, constants and initial values
741    PUBLIC usm_anthropogenic_heat, usm_material_model, ra_horiz_coef,           &
742           usm_green_heat_model, usm_temperature_near_surface
743
744
745
746 CONTAINS
747
748!------------------------------------------------------------------------------!
749! Description:
750! ------------
751!> This subroutine creates the necessary indices of the urban surfaces
752!> and plant canopy and it allocates the needed arrays for USM
753!------------------------------------------------------------------------------!
754    SUBROUTINE usm_allocate_surface
755   
756        IMPLICIT NONE
757       
758        INTEGER(iwp) ::  l
759
760!
761!--     Allocate radiation arrays which are part of the new data type.
762!--     For horizontal surfaces.
763        ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns)    )
764        ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
765!
766!--     For vertical surfaces
767        DO  l = 0, 3
768           ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns)    )
769           ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
770        ENDDO
771
772!--     Wall surface model
773!--     allocate arrays for wall surface model and define pointers
774       
775!--     allocate array of wall types and wall parameters
776        ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
777        DO  l = 0, 3
778           ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
779        ENDDO
780!
781!--     Allocate albedo_type and albedo. Each surface element
782!--     has 3 values, 0: wall fraction, 1: green fraction, 2: window fraction.
783        ALLOCATE( surf_usm_h%albedo_type(0:2,1:surf_usm_h%ns) )
784        ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)      )
785        surf_usm_h%albedo_type = albedo_type
786        DO  l = 0, 3
787           ALLOCATE( surf_usm_v(l)%albedo_type(0:2,1:surf_usm_v(l)%ns) )
788           ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns)      )
789           surf_usm_v(l)%albedo_type = albedo_type
790        ENDDO       
791
792
793!
794!--     Allocate indoor target temperature for summer and winter
795        ALLOCATE( surf_usm_h%target_temp_summer(1:surf_usm_h%ns) )
796        ALLOCATE( surf_usm_h%target_temp_winter(1:surf_usm_h%ns) )
797        DO  l = 0, 3
798           ALLOCATE( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) )
799           ALLOCATE( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) )
800        ENDDO   
801!
802!--     Allocate flag indicating ground floor level surface elements
803        ALLOCATE ( surf_usm_h%ground_level(1:surf_usm_h%ns) ) 
804        DO  l = 0, 3
805           ALLOCATE( surf_usm_v(l)%ground_level(1:surf_usm_v(l)%ns) )
806        ENDDO   
807!
808!--      Allocate arrays for relative surface fraction.
809!--      0 - wall fraction, 1 - green fraction, 2 - window fraction
810         ALLOCATE( surf_usm_h%frac(0:2,1:surf_usm_h%ns) )
811         surf_usm_h%frac = 0.0_wp
812         DO  l = 0, 3
813            ALLOCATE( surf_usm_v(l)%frac(0:2,1:surf_usm_v(l)%ns) )
814            surf_usm_v(l)%frac = 0.0_wp
815         ENDDO
816       
817!--     wall and roof surface parameters. First for horizontal surfaces
818        ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns)     )
819        ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns)     )
820        ALLOCATE ( surf_usm_h%lambda_surf_window(1:surf_usm_h%ns) )
821        ALLOCATE ( surf_usm_h%lambda_surf_green(1:surf_usm_h%ns)  )
822        ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns)       )
823        ALLOCATE ( surf_usm_h%c_surface_window(1:surf_usm_h%ns)   )
824        ALLOCATE ( surf_usm_h%c_surface_green(1:surf_usm_h%ns)    )
825        ALLOCATE ( surf_usm_h%transmissivity(1:surf_usm_h%ns)  )
826        ALLOCATE ( surf_usm_h%lai(1:surf_usm_h%ns)             )
827        ALLOCATE ( surf_usm_h%emissivity(0:2,1:surf_usm_h%ns)  )
828
829!
830!--     For vertical surfaces.
831        DO  l = 0, 3
832           ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns)     )
833           ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns)       )
834           ALLOCATE ( surf_usm_v(l)%lambda_surf_window(1:surf_usm_v(l)%ns) )
835           ALLOCATE ( surf_usm_v(l)%c_surface_window(1:surf_usm_v(l)%ns)   )
836           ALLOCATE ( surf_usm_v(l)%lambda_surf_green(1:surf_usm_v(l)%ns)  )
837           ALLOCATE ( surf_usm_v(l)%c_surface_green(1:surf_usm_v(l)%ns)    )
838           ALLOCATE ( surf_usm_v(l)%transmissivity(1:surf_usm_v(l)%ns)  )
839           ALLOCATE ( surf_usm_v(l)%lai(1:surf_usm_v(l)%ns)             )
840
841           ALLOCATE ( surf_usm_v(l)%emissivity(0:2,1:surf_usm_v(l)%ns)  )
842        ENDDO
843
844!       
845!--     allocate wall and roof material parameters. First for horizontal surfaces
846        ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns)               )
847        ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns)                  )
848        ALLOCATE ( surf_usm_h%thickness_green(1:surf_usm_h%ns)                   )
849        ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns)   )
850        ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
851        ALLOCATE ( surf_usm_h%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
852        ALLOCATE ( surf_usm_h%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_h%ns)    )
853        ALLOCATE ( surf_usm_h%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
854        ALLOCATE ( surf_usm_h%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
855
856!
857!--     For vertical surfaces.
858        DO  l = 0, 3
859           ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns)               )
860           ALLOCATE ( surf_usm_v(l)%thickness_window(1:surf_usm_v(l)%ns)                  )
861           ALLOCATE ( surf_usm_v(l)%thickness_green(1:surf_usm_v(l)%ns)                   )
862           ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)   )
863           ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
864           ALLOCATE ( surf_usm_v(l)%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
865           ALLOCATE ( surf_usm_v(l)%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)    )
866           ALLOCATE ( surf_usm_v(l)%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
867           ALLOCATE ( surf_usm_v(l)%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)     )
868        ENDDO
869
870!--     allocate wall and roof layers sizes. For horizontal surfaces.
871        ALLOCATE ( zwn(nzb_wall:nzt_wall) )
872        ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
873        ALLOCATE ( zwn_window(nzb_wall:nzt_wall) )
874        ALLOCATE ( surf_usm_h%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
875        ALLOCATE ( zwn_green(nzb_wall:nzt_wall) )
876        ALLOCATE ( surf_usm_h%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)      )
877        ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)    )
878        ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
879        ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
880        ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns)            )
881        ALLOCATE ( surf_usm_h%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)    )
882        ALLOCATE ( surf_usm_h%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
883        ALLOCATE ( surf_usm_h%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
884        ALLOCATE ( surf_usm_h%zw_window(nzb_wall:nzt_wall,1:surf_usm_h%ns)       )
885        ALLOCATE ( surf_usm_h%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
886        ALLOCATE ( surf_usm_h%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)   )
887        ALLOCATE ( surf_usm_h%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
888        ALLOCATE ( surf_usm_h%zw_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)        )
889!
890!--     For vertical surfaces.
891        DO  l = 0, 3
892           ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
893           ALLOCATE ( surf_usm_v(l)%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
894           ALLOCATE ( surf_usm_v(l)%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)      )
895           ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)    )
896           ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
897           ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
898           ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)            )
899           ALLOCATE ( surf_usm_v(l)%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)    )
900           ALLOCATE ( surf_usm_v(l)%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
901           ALLOCATE ( surf_usm_v(l)%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
902           ALLOCATE ( surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)       )
903           ALLOCATE ( surf_usm_v(l)%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
904           ALLOCATE ( surf_usm_v(l)%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)   )
905           ALLOCATE ( surf_usm_v(l)%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
906           ALLOCATE ( surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)        )
907        ENDDO
908
909!--     allocate wall and roof temperature arrays, for horizontal walls
910#if defined( __nopointer )
911        IF ( .NOT. ALLOCATED( t_surf_h ) )                                     &
912           ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
913        IF ( .NOT. ALLOCATED( t_surf_h_p ) )                                   &
914           ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
915        IF ( .NOT. ALLOCATED( t_wall_h ) )                                     &           
916           ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
917        IF ( .NOT. ALLOCATED( t_wall_h_p ) )                                   &           
918           ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
919        IF ( .NOT. ALLOCATED( t_surf_window_h ) )                              &
920           ALLOCATE ( t_surf_window_h(1:surf_usm_h%ns) )
921        IF ( .NOT. ALLOCATED( t_surf_window_h_p ) )                            &
922           ALLOCATE ( t_surf_window_h_p(1:surf_usm_h%ns) )
923        IF ( .NOT. ALLOCATED( t_window_h ) )                                   &           
924           ALLOCATE ( t_window_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
925        IF ( .NOT. ALLOCATED( t_window_h_p ) )                                 &           
926           ALLOCATE ( t_window_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
927        IF ( .NOT. ALLOCATED( t_surf_green_h ) )                               &
928           ALLOCATE ( t_surf_green_h(1:surf_usm_h%ns) )
929        IF ( .NOT. ALLOCATED( t_surf_green_h_p ) )                             &
930           ALLOCATE ( t_surf_green_h_p(1:surf_usm_h%ns) )           
931        IF ( .NOT. ALLOCATED( t_green_h ) )                                    &           
932           ALLOCATE ( t_green_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
933        IF ( .NOT. ALLOCATED( t_green_h_p ) )                                  &           
934           ALLOCATE ( t_green_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
935        IF ( .NOT. ALLOCATED( t_surf_whole_h ) )                               &
936           ALLOCATE ( t_surf_whole_h(1:surf_usm_h%ns) )
937        IF ( .NOT. ALLOCATED( t_surf_whole_h_p ) )                             &
938           ALLOCATE ( t_surf_whole_h_p(1:surf_usm_h%ns) )           
939        IF ( .NOT. ALLOCATED( t_surf_10cm_h ) )                                &
940           ALLOCATE ( t_surf_10cm_h(1:surf_usm_h%ns) )
941        IF ( .NOT. ALLOCATED( t_surf_10cm_h_p ) )                              &
942           ALLOCATE ( t_surf_10cm_h_p(1:surf_usm_h%ns) )
943#else
944!
945!--     Allocate if required. Note, in case of restarts, some of these arrays
946!--     might be already allocated.
947        IF ( .NOT. ALLOCATED( t_surf_h_1 ) )                                   &
948           ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) )
949        IF ( .NOT. ALLOCATED( t_surf_h_2 ) )                                   &
950           ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) )
951        IF ( .NOT. ALLOCATED( t_wall_h_1 ) )                                   &           
952           ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
953        IF ( .NOT. ALLOCATED( t_wall_h_2 ) )                                   &           
954           ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )         
955        IF ( .NOT. ALLOCATED( t_surf_window_h_1 ) )                            &
956           ALLOCATE ( t_surf_window_h_1(1:surf_usm_h%ns) )
957        IF ( .NOT. ALLOCATED( t_surf_window_h_2 ) )                            &
958           ALLOCATE ( t_surf_window_h_2(1:surf_usm_h%ns) )
959        IF ( .NOT. ALLOCATED( t_window_h_1 ) )                                 &           
960           ALLOCATE ( t_window_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
961        IF ( .NOT. ALLOCATED( t_window_h_2 ) )                                 &           
962           ALLOCATE ( t_window_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )         
963        IF ( .NOT. ALLOCATED( t_surf_green_h_1 ) )                             &
964           ALLOCATE ( t_surf_green_h_1(1:surf_usm_h%ns) )
965        IF ( .NOT. ALLOCATED( t_surf_green_h_2 ) )                             &
966           ALLOCATE ( t_surf_green_h_2(1:surf_usm_h%ns) )
967        IF ( .NOT. ALLOCATED( t_green_h_1 ) )                                  &           
968           ALLOCATE ( t_green_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
969        IF ( .NOT. ALLOCATED( t_green_h_2 ) )                                  &           
970           ALLOCATE ( t_green_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )         
971        IF ( .NOT. ALLOCATED( t_surf_whole_h_1 ) )                             &
972           ALLOCATE ( t_surf_whole_h_1(1:surf_usm_h%ns) )
973        IF ( .NOT. ALLOCATED( t_surf_whole_h_2 ) )                             &
974           ALLOCATE ( t_surf_whole_h_2(1:surf_usm_h%ns) )
975        IF ( .NOT. ALLOCATED( t_surf_10cm_h_1 ) )                              &
976           ALLOCATE ( t_surf_10cm_h_1(1:surf_usm_h%ns) )
977        IF ( .NOT. ALLOCATED( t_surf_10cm_h_2 ) )                              &
978           ALLOCATE ( t_surf_10cm_h_2(1:surf_usm_h%ns) )
979!           
980!--     initial assignment of the pointers
981        t_wall_h    => t_wall_h_1;    t_wall_h_p    => t_wall_h_2
982        t_window_h    => t_window_h_1;    t_window_h_p    => t_window_h_2
983        t_green_h    => t_green_h_1;    t_green_h_p    => t_green_h_2
984        t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2           
985        t_surf_window_h => t_surf_window_h_1; t_surf_window_h_p => t_surf_window_h_2 
986        t_surf_green_h => t_surf_green_h_1; t_surf_green_h_p => t_surf_green_h_2           
987        t_surf_whole_h => t_surf_whole_h_1; t_surf_whole_h_p => t_surf_whole_h_2
988        t_surf_10cm_h => t_surf_10cm_h_1; t_surf_10cm_h_p => t_surf_10cm_h_2 
989 
990#endif
991
992!--     allocate wall and roof temperature arrays, for vertical walls if required
993#if defined( __nopointer )
994        DO  l = 0, 3
995           IF ( .NOT. ALLOCATED( t_surf_v(l)%t ) )                             &
996              ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) )
997           IF ( .NOT. ALLOCATED( t_surf_v_p(l)%t ) )                           &
998              ALLOCATE ( t_surf_v_p(l)%t(1:surf_usm_v(l)%ns) )
999           IF ( .NOT. ALLOCATED( t_wall_v(l)%t ) )                             &
1000              ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1001           IF ( .NOT. ALLOCATED( t_wall_v_p(l)%t ) )                           &                 
1002              ALLOCATE ( t_wall_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1003           IF ( .NOT. ALLOCATED( t_surf_window_v(l)%t ) )                      &
1004              ALLOCATE ( t_surf_window_v(l)%t(1:surf_usm_v(l)%ns) )
1005           IF ( .NOT. ALLOCATED( t_surf_window_v_p(l)%t ) )                    &
1006              ALLOCATE ( t_surf_window_v_p(l)%t(1:surf_usm_v(l)%ns) )
1007           IF ( .NOT. ALLOCATED( t_window_v(l)%t ) )                           &
1008              ALLOCATE ( t_window_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1009           IF ( .NOT. ALLOCATED( t_window_v_p(l)%t ) )                         &                 
1010              ALLOCATE ( t_window_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1011           IF ( .NOT. ALLOCATED( t_green_v(l)%t ) )                            &
1012              ALLOCATE ( t_green_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1013           IF ( .NOT. ALLOCATED( t_green_v_p(l)%t ) )                          &                 
1014              ALLOCATE ( t_green_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1015           IF ( .NOT. ALLOCATED( t_surf_green_v(l)%t ) )                       &
1016              ALLOCATE ( t_surf_green_v(l)%t(1:surf_usm_v(l)%ns) )
1017           IF ( .NOT. ALLOCATED( t_surf_green_v_p(l)%t ) )                     &
1018              ALLOCATE ( t_surf_green_v_p(l)%t(1:surf_usm_v(l)%ns) )
1019           IF ( .NOT. ALLOCATED( t_surf_whole_v(l)%t ) )                       &
1020              ALLOCATE ( t_surf_whole_v(l)%t(1:surf_usm_v(l)%ns) )
1021           IF ( .NOT. ALLOCATED( t_surf_whole_v_p(l)%t ) )                     &
1022              ALLOCATE ( t_surf_whole_v_p(l)%t(1:surf_usm_v(l)%ns) )
1023           IF ( .NOT. ALLOCATED( t_surf_10cm_v(l)%t ) )                        &
1024              ALLOCATE ( t_surf_10cm_v(l)%t(1:surf_usm_v(l)%ns) )
1025           IF ( .NOT. ALLOCATED( t_surf_10cm_v_p(l)%t ) )                        &
1026              ALLOCATE ( t_surf_10cm_v_p(l)%t(1:surf_usm_v(l)%ns) )
1027        ENDDO
1028#else
1029!
1030!--     Allocate if required. Note, in case of restarts, some of these arrays
1031!--     might be already allocated.
1032        DO  l = 0, 3
1033           IF ( .NOT. ALLOCATED( t_surf_v_1(l)%t ) )                           &
1034              ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) )
1035           IF ( .NOT. ALLOCATED( t_surf_v_2(l)%t ) )                           &
1036              ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) )
1037           IF ( .NOT. ALLOCATED( t_wall_v_1(l)%t ) )                           &           
1038              ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
1039           IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) )                           &           
1040              ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
1041           IF ( .NOT. ALLOCATED( t_surf_window_v_1(l)%t ) )                    &
1042              ALLOCATE ( t_surf_window_v_1(l)%t(1:surf_usm_v(l)%ns) )
1043           IF ( .NOT. ALLOCATED( t_surf_window_v_2(l)%t ) )                    &
1044              ALLOCATE ( t_surf_window_v_2(l)%t(1:surf_usm_v(l)%ns) )
1045           IF ( .NOT. ALLOCATED( t_window_v_1(l)%t ) )                         &           
1046              ALLOCATE ( t_window_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
1047           IF ( .NOT. ALLOCATED( t_window_v_2(l)%t ) )                         &           
1048              ALLOCATE ( t_window_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
1049           IF ( .NOT. ALLOCATED( t_surf_green_v_1(l)%t ) )                     &
1050              ALLOCATE ( t_surf_green_v_1(l)%t(1:surf_usm_v(l)%ns) )
1051           IF ( .NOT. ALLOCATED( t_surf_green_v_2(l)%t ) )                     &
1052              ALLOCATE ( t_surf_green_v_2(l)%t(1:surf_usm_v(l)%ns) )
1053           IF ( .NOT. ALLOCATED( t_green_v_1(l)%t ) )                          &           
1054              ALLOCATE ( t_green_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
1055           IF ( .NOT. ALLOCATED( t_green_v_2(l)%t ) )                          &           
1056              ALLOCATE ( t_green_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
1057           IF ( .NOT. ALLOCATED( t_surf_whole_v_1(l)%t ) )                     &
1058              ALLOCATE ( t_surf_whole_v_1(l)%t(1:surf_usm_v(l)%ns) )
1059           IF ( .NOT. ALLOCATED( t_surf_whole_v_2(l)%t ) )                     &
1060              ALLOCATE ( t_surf_whole_v_2(l)%t(1:surf_usm_v(l)%ns) )
1061           IF ( .NOT. ALLOCATED( t_surf_10cm_v_1(l)%t ) )                     &
1062              ALLOCATE ( t_surf_10cm_v_1(l)%t(1:surf_usm_v(l)%ns) )
1063           IF ( .NOT. ALLOCATED( t_surf_10cm_v_2(l)%t ) )                     &
1064              ALLOCATE ( t_surf_10cm_v_2(l)%t(1:surf_usm_v(l)%ns) )
1065        ENDDO
1066!
1067!--     initial assignment of the pointers
1068        t_wall_v    => t_wall_v_1;    t_wall_v_p    => t_wall_v_2
1069        t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
1070        t_window_v    => t_window_v_1;    t_window_v_p    => t_window_v_2
1071        t_green_v    => t_green_v_1;    t_green_v_p    => t_green_v_2
1072        t_surf_window_v => t_surf_window_v_1; t_surf_window_v_p => t_surf_window_v_2
1073        t_surf_green_v => t_surf_green_v_1; t_surf_green_v_p => t_surf_green_v_2
1074        t_surf_whole_v => t_surf_whole_v_1; t_surf_whole_v_p => t_surf_whole_v_2
1075        t_surf_10cm_v => t_surf_10cm_v_1; t_surf_10cm_v_p => t_surf_10cm_v_2
1076
1077#endif
1078!
1079!--     Allocate intermediate timestep arrays. For horizontal surfaces.
1080        ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns)                  )
1081        ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1082        ALLOCATE ( surf_usm_h%tt_surface_window_m(1:surf_usm_h%ns)             )
1083        ALLOCATE ( surf_usm_h%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1084        ALLOCATE ( surf_usm_h%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)  )
1085        ALLOCATE ( surf_usm_h%tt_surface_green_m(1:surf_usm_h%ns)              )
1086
1087!
1088!--     Set inital values for prognostic quantities
1089        IF ( ALLOCATED( surf_usm_h%tt_surface_m ) )  surf_usm_h%tt_surface_m = 0.0_wp
1090        IF ( ALLOCATED( surf_usm_h%tt_wall_m    ) )  surf_usm_h%tt_wall_m    = 0.0_wp
1091        IF ( ALLOCATED( surf_usm_h%tt_surface_window_m ) )  surf_usm_h%tt_surface_window_m = 0.0_wp
1092        IF ( ALLOCATED( surf_usm_h%tt_window_m    )      )  surf_usm_h%tt_window_m         = 0.0_wp
1093        IF ( ALLOCATED( surf_usm_h%tt_green_m    )       )  surf_usm_h%tt_green_m          = 0.0_wp
1094        IF ( ALLOCATED( surf_usm_h%tt_surface_green_m )  )  surf_usm_h%tt_surface_green_m  = 0.0_wp
1095!
1096!--     Now, for vertical surfaces
1097        DO  l = 0, 3
1098           ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns)                  )
1099           ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1100           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) )  surf_usm_v(l)%tt_surface_m = 0.0_wp
1101           IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m    ) )  surf_usm_v(l)%tt_wall_m    = 0.0_wp
1102           ALLOCATE ( surf_usm_v(l)%tt_surface_window_m(1:surf_usm_v(l)%ns)             )
1103           ALLOCATE ( surf_usm_v(l)%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1104           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_window_m ) )  surf_usm_v(l)%tt_surface_window_m = 0.0_wp
1105           IF ( ALLOCATED( surf_usm_v(l)%tt_window_m  ) )  surf_usm_v(l)%tt_window_m    = 0.0_wp
1106           ALLOCATE ( surf_usm_v(l)%tt_surface_green_m(1:surf_usm_v(l)%ns)              )
1107           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_green_m ) )  surf_usm_v(l)%tt_surface_green_m = 0.0_wp
1108           ALLOCATE ( surf_usm_v(l)%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)  )
1109           IF ( ALLOCATED( surf_usm_v(l)%tt_green_m   ) )  surf_usm_v(l)%tt_green_m    = 0.0_wp
1110        ENDDO
1111
1112!--     allocate wall heat flux output array and set initial values. For horizontal surfaces
1113!         ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns)    )  !can be removed
1114        ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
1115        ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
1116        ALLOCATE ( surf_usm_h%wghf_eb_window(1:surf_usm_h%ns) )
1117        ALLOCATE ( surf_usm_h%wghf_eb_green(1:surf_usm_h%ns) )
1118        ALLOCATE ( surf_usm_h%iwghf_eb(1:surf_usm_h%ns) )
1119        ALLOCATE ( surf_usm_h%iwghf_eb_window(1:surf_usm_h%ns) )
1120        IF ( ALLOCATED( surf_usm_h%wshf    ) )  surf_usm_h%wshf    = 0.0_wp
1121        IF ( ALLOCATED( surf_usm_h%wshf_eb ) )  surf_usm_h%wshf_eb = 0.0_wp
1122        IF ( ALLOCATED( surf_usm_h%wghf_eb ) )  surf_usm_h%wghf_eb = 0.0_wp
1123        IF ( ALLOCATED( surf_usm_h%wghf_eb_window ) )  surf_usm_h%wghf_eb_window = 0.0_wp
1124        IF ( ALLOCATED( surf_usm_h%wghf_eb_green ) )  surf_usm_h%wghf_eb_green = 0.0_wp
1125        IF ( ALLOCATED( surf_usm_h%iwghf_eb ) )  surf_usm_h%iwghf_eb = 0.0_wp
1126        IF ( ALLOCATED( surf_usm_h%iwghf_eb_window ) )  surf_usm_h%iwghf_eb_window = 0.0_wp
1127!
1128!--     Now, for vertical surfaces
1129        DO  l = 0, 3
1130!            ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns)    )    ! can be removed
1131           ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
1132           ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
1133           ALLOCATE ( surf_usm_v(l)%wghf_eb_window(1:surf_usm_v(l)%ns) )
1134           ALLOCATE ( surf_usm_v(l)%wghf_eb_green(1:surf_usm_v(l)%ns) )
1135           ALLOCATE ( surf_usm_v(l)%iwghf_eb(1:surf_usm_v(l)%ns) )
1136           ALLOCATE ( surf_usm_v(l)%iwghf_eb_window(1:surf_usm_v(l)%ns) )
1137           IF ( ALLOCATED( surf_usm_v(l)%wshf    ) )  surf_usm_v(l)%wshf    = 0.0_wp
1138           IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) )  surf_usm_v(l)%wshf_eb = 0.0_wp
1139           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) )  surf_usm_v(l)%wghf_eb = 0.0_wp
1140           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_window ) )  surf_usm_v(l)%wghf_eb_window = 0.0_wp
1141           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_green ) )  surf_usm_v(l)%wghf_eb_green = 0.0_wp
1142           IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb ) )  surf_usm_v(l)%iwghf_eb = 0.0_wp
1143           IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) )  surf_usm_v(l)%iwghf_eb_window = 0.0_wp
1144        ENDDO
1145       
1146    END SUBROUTINE usm_allocate_surface
1147
1148
1149!------------------------------------------------------------------------------!
1150! Description:
1151! ------------
1152!> Sum up and time-average urban surface output quantities as well as allocate
1153!> the array necessary for storing the average.
1154!------------------------------------------------------------------------------!
1155    SUBROUTINE usm_average_3d_data( mode, variable )
1156
1157        IMPLICIT NONE
1158
1159        CHARACTER (len=*), INTENT(IN) ::  mode
1160        CHARACTER (len=*), INTENT(IN) :: variable
1161 
1162        INTEGER(iwp)                                       :: i, j, k, l, m, ids, iwl,istat
1163        CHARACTER (len=varnamelength)                      :: var, surfid
1164        INTEGER(iwp), PARAMETER                            :: nd = 5
1165        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER     :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
1166
1167!--     find the real name of the variable
1168        var = TRIM(variable)
1169        DO i = 0, nd-1
1170            k = len(TRIM(var))
1171            j = len(TRIM(dirname(i)))
1172            IF ( var(k-j+1:k) == dirname(i) )  THEN
1173                ids = i
1174                var = var(:k-j)
1175                EXIT
1176            ENDIF
1177        ENDDO
1178        IF ( ids == -1 )  THEN
1179            var = TRIM(variable)
1180        ENDIF
1181        IF ( var(1:11) == 'usm_t_wall_'  .AND.  len(TRIM(var)) >= 12 )  THEN
1182!--          wall layers
1183            READ(var(12:12), '(I1)', iostat=istat ) iwl
1184            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
1185                var = var(1:10)
1186            ELSE
1187!--             wrong wall layer index
1188                RETURN
1189            ENDIF
1190        ENDIF
1191        IF ( var(1:13) == 'usm_t_window_'  .AND.  len(TRIM(var)) >= 14 )  THEN
1192!--          wall layers
1193            READ(var(14:14), '(I1)', iostat=istat ) iwl
1194            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
1195                var = var(1:12)
1196            ELSE
1197!--             wrong window layer index
1198                RETURN
1199            ENDIF
1200        ENDIF
1201        IF ( var(1:12) == 'usm_t_green_'  .AND.  len(TRIM(var)) >= 13 )  THEN
1202!--          wall layers
1203            READ(var(13:13), '(I1)', iostat=istat ) iwl
1204            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
1205                var = var(1:11)
1206            ELSE
1207!--             wrong green layer index
1208                RETURN
1209            ENDIF
1210        ENDIF
1211
1212        IF ( mode == 'allocate' )  THEN
1213           
1214           SELECT CASE ( TRIM( var ) )
1215               
1216                CASE ( 'usm_rad_net' )
1217!--                 array of complete radiation balance
1218                    IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
1219                        ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
1220                        surf_usm_h%rad_net_av = 0.0_wp
1221                    ENDIF
1222                    DO  l = 0, 3
1223                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
1224                           ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
1225                           surf_usm_v(l)%rad_net_av = 0.0_wp
1226                       ENDIF
1227                    ENDDO
1228                   
1229                CASE ( 'usm_rad_insw' )
1230!--                 array of sw radiation falling to surface after i-th reflection
1231                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfinsw_av) )  THEN
1232                        ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
1233                        surf_usm_h%surfinsw_av = 0.0_wp
1234                    ENDIF
1235                    DO  l = 0, 3
1236                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinsw_av) )  THEN
1237                           ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
1238                           surf_usm_v(l)%surfinsw_av = 0.0_wp
1239                       ENDIF
1240                    ENDDO
1241                                   
1242                CASE ( 'usm_rad_inlw' )
1243!--                 array of lw radiation falling to surface after i-th reflection
1244                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfinlw_av) )  THEN
1245                        ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
1246                        surf_usm_h%surfinlw_av = 0.0_wp
1247                    ENDIF
1248                    DO  l = 0, 3
1249                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinlw_av) )  THEN
1250                           ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
1251                           surf_usm_v(l)%surfinlw_av = 0.0_wp
1252                       ENDIF
1253                    ENDDO
1254
1255                CASE ( 'usm_rad_inswdir' )
1256!--                 array of direct sw radiation falling to surface from sun
1257                    IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
1258                        ALLOCATE( surfinswdir_av(startenergy:endenergy) )
1259                        surfinswdir_av = 0.0_wp
1260                    ENDIF
1261
1262                CASE ( 'usm_rad_inswdif' )
1263!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
1264                    IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
1265                        ALLOCATE( surfinswdif_av(startenergy:endenergy) )
1266                        surfinswdif_av = 0.0_wp
1267                    ENDIF
1268
1269                CASE ( 'usm_rad_inswref' )
1270!--                 array of sw radiation falling to surface from reflections
1271                    IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
1272                        ALLOCATE( surfinswref_av(startenergy:endenergy) )
1273                        surfinswref_av = 0.0_wp
1274                    ENDIF
1275
1276                CASE ( 'usm_rad_inlwdif' )
1277!--                 array of sw radiation falling to surface after i-th reflection
1278                   IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
1279                        ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
1280                        surfinlwdif_av = 0.0_wp
1281                    ENDIF
1282
1283                CASE ( 'usm_rad_inlwref' )
1284!--                 array of lw radiation falling to surface from reflections
1285                    IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
1286                        ALLOCATE( surfinlwref_av(startenergy:endenergy) )
1287                        surfinlwref_av = 0.0_wp
1288                    ENDIF
1289
1290                CASE ( 'usm_rad_outsw' )
1291!--                 array of sw radiation emitted from surface after i-th reflection
1292                    IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
1293                        ALLOCATE( surfoutsw_av(startenergy:endenergy) )
1294                        surfoutsw_av = 0.0_wp
1295                    ENDIF
1296
1297                CASE ( 'usm_rad_outlw' )
1298!--                 array of lw radiation emitted from surface after i-th reflection
1299                    IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
1300                        ALLOCATE( surfoutlw_av(startenergy:endenergy) )
1301                        surfoutlw_av = 0.0_wp
1302                    ENDIF
1303                CASE ( 'usm_rad_ressw' )
1304!--                 array of residua of sw radiation absorbed in surface after last reflection
1305                    IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
1306                        ALLOCATE( surfins_av(startenergy:endenergy) )
1307                        surfins_av = 0.0_wp
1308                    ENDIF
1309                                   
1310                CASE ( 'usm_rad_reslw' )
1311!--                 array of residua of lw radiation absorbed in surface after last reflection
1312                    IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
1313                        ALLOCATE( surfinl_av(startenergy:endenergy) )
1314                        surfinl_av = 0.0_wp
1315                    ENDIF
1316                                   
1317                CASE ( 'usm_rad_hf' )
1318!--                 array of heat flux from radiation for surfaces after i-th reflection
1319                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
1320                        ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
1321                        surf_usm_h%surfhf_av = 0.0_wp
1322                    ENDIF
1323                    DO  l = 0, 3
1324                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
1325                           ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
1326                           surf_usm_v(l)%surfhf_av = 0.0_wp
1327                       ENDIF
1328                    ENDDO
1329
1330                CASE ( 'usm_wshf' )
1331!--                 array of sensible heat flux from surfaces
1332!--                 land surfaces
1333                    IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
1334                        ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1335                        surf_usm_h%wshf_eb_av = 0.0_wp
1336                    ENDIF
1337                    DO  l = 0, 3
1338                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
1339                           ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1340                           surf_usm_v(l)%wshf_eb_av = 0.0_wp
1341                       ENDIF
1342                    ENDDO
1343
1344                CASE ( 'usm_wghf' )
1345!--                 array of heat flux from ground (wall, roof, land)
1346                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
1347                        ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1348                        surf_usm_h%wghf_eb_av = 0.0_wp
1349                    ENDIF
1350                    DO  l = 0, 3
1351                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
1352                           ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1353                           surf_usm_v(l)%wghf_eb_av = 0.0_wp
1354                       ENDIF
1355                    ENDDO
1356
1357                CASE ( 'usm_wghf_window' )
1358!--                 array of heat flux from window ground (wall, roof, land)
1359                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
1360                        ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
1361                        surf_usm_h%wghf_eb_window_av = 0.0_wp
1362                    ENDIF
1363                    DO  l = 0, 3
1364                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
1365                           ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
1366                           surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
1367                       ENDIF
1368                    ENDDO
1369
1370                CASE ( 'usm_wghf_green' )
1371!--                 array of heat flux from green ground (wall, roof, land)
1372                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
1373                        ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
1374                        surf_usm_h%wghf_eb_green_av = 0.0_wp
1375                    ENDIF
1376                    DO  l = 0, 3
1377                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
1378                           ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
1379                           surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
1380                       ENDIF
1381                    ENDDO
1382
1383                CASE ( 'usm_iwghf' )
1384!--                 array of heat flux from indoor ground (wall, roof, land)
1385                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
1386                        ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
1387                        surf_usm_h%iwghf_eb_av = 0.0_wp
1388                    ENDIF
1389                    DO  l = 0, 3
1390                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
1391                           ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
1392                           surf_usm_v(l)%iwghf_eb_av = 0.0_wp
1393                       ENDIF
1394                    ENDDO
1395
1396                CASE ( 'usm_iwghf_window' )
1397!--                 array of heat flux from indoor window ground (wall, roof, land)
1398                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
1399                        ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
1400                        surf_usm_h%iwghf_eb_window_av = 0.0_wp
1401                    ENDIF
1402                    DO  l = 0, 3
1403                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
1404                           ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
1405                           surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
1406                       ENDIF
1407                    ENDDO
1408                   
1409                CASE ( 'usm_t_surf' )
1410!--                 surface temperature for surfaces
1411                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
1412                        ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
1413                        surf_usm_h%t_surf_av = 0.0_wp
1414                    ENDIF
1415                    DO  l = 0, 3
1416                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_av) )  THEN
1417                           ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
1418                           surf_usm_v(l)%t_surf_av = 0.0_wp
1419                       ENDIF
1420                    ENDDO
1421
1422                CASE ( 'usm_t_surf_window' )
1423!--                 surface temperature for window surfaces
1424                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
1425                        ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
1426                        surf_usm_h%t_surf_window_av = 0.0_wp
1427                    ENDIF
1428                    DO  l = 0, 3
1429                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
1430                           ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
1431                           surf_usm_v(l)%t_surf_window_av = 0.0_wp
1432                       ENDIF
1433                    ENDDO
1434                   
1435                CASE ( 'usm_t_surf_green' )
1436!--                 surface temperature for green surfaces
1437                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
1438                        ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
1439                        surf_usm_h%t_surf_green_av = 0.0_wp
1440                    ENDIF
1441                    DO  l = 0, 3
1442                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
1443                           ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
1444                           surf_usm_v(l)%t_surf_green_av = 0.0_wp
1445                       ENDIF
1446                    ENDDO
1447               
1448                CASE ( 'usm_t_surf_whole' )
1449!--                 surface temperature for whole surfaces
1450                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_whole_av) )  THEN
1451                        ALLOCATE( surf_usm_h%t_surf_whole_av(1:surf_usm_h%ns) )
1452                        surf_usm_h%t_surf_whole_av = 0.0_wp
1453                    ENDIF
1454                    DO  l = 0, 3
1455                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_whole_av) )  THEN
1456                           ALLOCATE( surf_usm_v(l)%t_surf_whole_av(1:surf_usm_v(l)%ns) )
1457                           surf_usm_v(l)%t_surf_whole_av = 0.0_wp
1458                       ENDIF
1459                    ENDDO
1460                   
1461                CASE ( 'usm_t_surf_10cm' )
1462!--                 near surface temperature for whole surfaces
1463                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
1464                        ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
1465                        surf_usm_h%t_surf_10cm_av = 0.0_wp
1466                    ENDIF
1467                    DO  l = 0, 3
1468                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) )  THEN
1469                           ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) )
1470                           surf_usm_v(l)%t_surf_10cm_av = 0.0_wp
1471                       ENDIF
1472                    ENDDO
1473
1474                CASE ( 'usm_t_wall' )
1475!--                 wall temperature for iwl layer of walls and land
1476                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
1477                        ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1478                        surf_usm_h%t_wall_av = 0.0_wp
1479                    ENDIF
1480                    DO  l = 0, 3
1481                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
1482                           ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1483                           surf_usm_v(l)%t_wall_av = 0.0_wp
1484                       ENDIF
1485                    ENDDO
1486
1487                CASE ( 'usm_t_window' )
1488!--                 window temperature for iwl layer of walls and land
1489                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
1490                        ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1491                        surf_usm_h%t_window_av = 0.0_wp
1492                    ENDIF
1493                    DO  l = 0, 3
1494                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
1495                           ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1496                           surf_usm_v(l)%t_window_av = 0.0_wp
1497                       ENDIF
1498                    ENDDO
1499
1500                CASE ( 'usm_t_green' )
1501!--                 green temperature for iwl layer of walls and land
1502                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
1503                        ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1504                        surf_usm_h%t_green_av = 0.0_wp
1505                    ENDIF
1506                    DO  l = 0, 3
1507                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
1508                           ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1509                           surf_usm_v(l)%t_green_av = 0.0_wp
1510                       ENDIF
1511                    ENDDO
1512
1513               CASE DEFAULT
1514                   CONTINUE
1515
1516           END SELECT
1517
1518        ELSEIF ( mode == 'sum' )  THEN
1519           
1520           SELECT CASE ( TRIM( var ) )
1521               
1522                CASE ( 'usm_rad_net' )
1523!--                 array of complete radiation balance
1524                    DO  m = 1, surf_usm_h%ns
1525                       surf_usm_h%rad_net_av(m) =                              &
1526                                          surf_usm_h%rad_net_av(m) +           &
1527                                          surf_usm_h%rad_net_l(m)
1528                    ENDDO
1529                    DO  l = 0, 3
1530                       DO  m = 1, surf_usm_v(l)%ns
1531                          surf_usm_v(l)%rad_net_av(m) =                        &
1532                                          surf_usm_v(l)%rad_net_av(m) +        &
1533                                          surf_usm_v(l)%rad_net_l(m)
1534                       ENDDO
1535                    ENDDO
1536                   
1537                CASE ( 'usm_rad_insw' )
1538!--                 array of sw radiation falling to surface after i-th reflection
1539                    DO l = startenergy, endenergy
1540                        IF ( surfl(id,l) == ids )  THEN
1541                            surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
1542                        ENDIF
1543                    ENDDO
1544                             
1545                CASE ( 'usm_rad_inlw' )
1546!--                 array of lw radiation falling to surface after i-th reflection
1547                    DO l = startenergy, endenergy
1548                        IF ( surfl(id,l) == ids )  THEN
1549                            surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
1550                        ENDIF
1551                    ENDDO
1552                   
1553                CASE ( 'usm_rad_inswdir' )
1554!--                 array of direct sw radiation falling to surface from sun
1555                    DO l = startenergy, endenergy
1556                        IF ( surfl(id,l) == ids )  THEN
1557                            surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
1558                        ENDIF
1559                    ENDDO
1560                   
1561                CASE ( 'usm_rad_inswdif' )
1562!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
1563                    DO l = startenergy, endenergy
1564                        IF ( surfl(id,l) == ids )  THEN
1565                            surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
1566                        ENDIF
1567                    ENDDO
1568                   
1569                CASE ( 'usm_rad_inswref' )
1570!--                 array of sw radiation falling to surface from reflections
1571                    DO l = startenergy, endenergy
1572                        IF ( surfl(id,l) == ids )  THEN
1573                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1574                                                surfinswdir(l) - surfinswdif(l)
1575                        ENDIF
1576                    ENDDO
1577
1578                   
1579                CASE ( 'usm_rad_inlwdif' )
1580!--                 array of sw radiation falling to surface after i-th reflection
1581                    DO l = startenergy, endenergy
1582                        IF ( surfl(id,l) == ids )  THEN
1583                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1584                                                surfinswdir(l) - surfinswdif(l)
1585                        ENDIF
1586                    ENDDO
1587!                     
1588                CASE ( 'usm_rad_inlwref' )
1589!--                 array of lw radiation falling to surface from reflections
1590                    DO l = startenergy, endenergy
1591                        IF ( surfl(id,l) == ids )  THEN
1592                            surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
1593                        ENDIF
1594                    ENDDO
1595                   
1596                CASE ( 'usm_rad_outsw' )
1597!--                 array of sw radiation emitted from surface after i-th reflection
1598                    DO l = startenergy, endenergy
1599                        IF ( surfl(id,l) == ids )  THEN
1600                            surfinlwref_av(l) = surfinlwref_av(l) + &
1601                                                surfinlw(l) - surfinlwdif(l)
1602                        ENDIF
1603                    ENDDO
1604                   
1605                CASE ( 'usm_rad_outlw' )
1606!--                 array of lw radiation emitted from surface after i-th reflection
1607                    DO l = startenergy, endenergy
1608                        IF ( surfl(id,l) == ids )  THEN
1609                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
1610                        ENDIF
1611                    ENDDO
1612                   
1613                CASE ( 'usm_rad_ressw' )
1614!--                 array of residua of sw radiation absorbed in surface after last reflection
1615                    DO l = startenergy, endenergy
1616                        IF ( surfl(id,l) == ids )  THEN
1617                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
1618                        ENDIF
1619                    ENDDO
1620                                   
1621                CASE ( 'usm_rad_reslw' )
1622!--                 array of residua of lw radiation absorbed in surface after last reflection
1623                    DO l = startenergy, endenergy
1624                        IF ( surfl(id,l) == ids )  THEN
1625                            surfins_av(l) = surfins_av(l) + surfins(l)
1626                        ENDIF
1627                    ENDDO
1628                   
1629                CASE ( 'usm_rad_hf' )
1630!--                 array of heat flux from radiation for surfaces after i-th reflection
1631                    DO  m = 1, surf_usm_h%ns
1632                       surf_usm_h%surfhf_av(m) =                               &
1633                                          surf_usm_h%surfhf_av(m) +            &
1634                                          surf_usm_h%surfhf(m)
1635                    ENDDO
1636                    DO  l = 0, 3
1637                       DO  m = 1, surf_usm_v(l)%ns
1638                          surf_usm_v(l)%surfhf_av(m) =                         &
1639                                          surf_usm_v(l)%surfhf_av(m) +         &
1640                                          surf_usm_v(l)%surfhf(m)
1641                       ENDDO
1642                    ENDDO
1643                   
1644                CASE ( 'usm_wshf' )
1645!--                 array of sensible heat flux from surfaces (land, roof, wall)
1646                    DO  m = 1, surf_usm_h%ns
1647                       surf_usm_h%wshf_eb_av(m) =                              &
1648                                          surf_usm_h%wshf_eb_av(m) +           &
1649                                          surf_usm_h%wshf_eb(m)
1650                    ENDDO
1651                    DO  l = 0, 3
1652                       DO  m = 1, surf_usm_v(l)%ns
1653                          surf_usm_v(l)%wshf_eb_av(m) =                        &
1654                                          surf_usm_v(l)%wshf_eb_av(m) +        &
1655                                          surf_usm_v(l)%wshf_eb(m)
1656                       ENDDO
1657                    ENDDO
1658                   
1659                CASE ( 'usm_wghf' )
1660!--                 array of heat flux from ground (wall, roof, land)
1661                    DO  m = 1, surf_usm_h%ns
1662                       surf_usm_h%wghf_eb_av(m) =                              &
1663                                          surf_usm_h%wghf_eb_av(m) +           &
1664                                          surf_usm_h%wghf_eb(m)
1665                    ENDDO
1666                    DO  l = 0, 3
1667                       DO  m = 1, surf_usm_v(l)%ns
1668                          surf_usm_v(l)%wghf_eb_av(m) =                        &
1669                                          surf_usm_v(l)%wghf_eb_av(m) +        &
1670                                          surf_usm_v(l)%wghf_eb(m)
1671                       ENDDO
1672                    ENDDO
1673                   
1674                CASE ( 'usm_wghf_window' )
1675!--                 array of heat flux from window ground (wall, roof, land)
1676                    DO  m = 1, surf_usm_h%ns
1677                       surf_usm_h%wghf_eb_window_av(m) =                              &
1678                                          surf_usm_h%wghf_eb_window_av(m) +           &
1679                                          surf_usm_h%wghf_eb_window(m)
1680                    ENDDO
1681                    DO  l = 0, 3
1682                       DO  m = 1, surf_usm_v(l)%ns
1683                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
1684                                          surf_usm_v(l)%wghf_eb_window_av(m) +        &
1685                                          surf_usm_v(l)%wghf_eb_window(m)
1686                       ENDDO
1687                    ENDDO
1688
1689                CASE ( 'usm_wghf_green' )
1690!--                 array of heat flux from green ground (wall, roof, land)
1691                    DO  m = 1, surf_usm_h%ns
1692                       surf_usm_h%wghf_eb_green_av(m) =                              &
1693                                          surf_usm_h%wghf_eb_green_av(m) +           &
1694                                          surf_usm_h%wghf_eb_green(m)
1695                    ENDDO
1696                    DO  l = 0, 3
1697                       DO  m = 1, surf_usm_v(l)%ns
1698                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
1699                                          surf_usm_v(l)%wghf_eb_green_av(m) +        &
1700                                          surf_usm_v(l)%wghf_eb_green(m)
1701                       ENDDO
1702                    ENDDO
1703                   
1704                CASE ( 'usm_iwghf' )
1705!--                 array of heat flux from indoor ground (wall, roof, land)
1706                    DO  m = 1, surf_usm_h%ns
1707                       surf_usm_h%iwghf_eb_av(m) =                              &
1708                                          surf_usm_h%iwghf_eb_av(m) +           &
1709                                          surf_usm_h%iwghf_eb(m)
1710                    ENDDO
1711                    DO  l = 0, 3
1712                       DO  m = 1, surf_usm_v(l)%ns
1713                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
1714                                          surf_usm_v(l)%iwghf_eb_av(m) +        &
1715                                          surf_usm_v(l)%iwghf_eb(m)
1716                       ENDDO
1717                    ENDDO
1718                   
1719                CASE ( 'usm_iwghf_window' )
1720!--                 array of heat flux from indoor window ground (wall, roof, land)
1721                    DO  m = 1, surf_usm_h%ns
1722                       surf_usm_h%iwghf_eb_window_av(m) =                              &
1723                                          surf_usm_h%iwghf_eb_window_av(m) +           &
1724                                          surf_usm_h%iwghf_eb_window(m)
1725                    ENDDO
1726                    DO  l = 0, 3
1727                       DO  m = 1, surf_usm_v(l)%ns
1728                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
1729                                          surf_usm_v(l)%iwghf_eb_window_av(m) +        &
1730                                          surf_usm_v(l)%iwghf_eb_window(m)
1731                       ENDDO
1732                    ENDDO
1733                   
1734                CASE ( 'usm_t_surf' )
1735!--                 surface temperature for surfaces
1736                    DO  m = 1, surf_usm_h%ns
1737                       surf_usm_h%t_surf_av(m) =                               & 
1738                                          surf_usm_h%t_surf_av(m) +            &
1739                                          t_surf_h(m)
1740                    ENDDO
1741                    DO  l = 0, 3
1742                       DO  m = 1, surf_usm_v(l)%ns
1743                          surf_usm_v(l)%t_surf_av(m) =                         &
1744                                          surf_usm_v(l)%t_surf_av(m) +         &
1745                                          t_surf_v(l)%t(m)
1746                       ENDDO
1747                    ENDDO
1748                   
1749                CASE ( 'usm_t_surf_window' )
1750!--                 surface temperature for window surfaces
1751                    DO  m = 1, surf_usm_h%ns
1752                       surf_usm_h%t_surf_window_av(m) =                               & 
1753                                          surf_usm_h%t_surf_window_av(m) +            &
1754                                          t_surf_window_h(m)
1755                    ENDDO
1756                    DO  l = 0, 3
1757                       DO  m = 1, surf_usm_v(l)%ns
1758                          surf_usm_v(l)%t_surf_window_av(m) =                         &
1759                                          surf_usm_v(l)%t_surf_window_av(m) +         &
1760                                          t_surf_window_v(l)%t(m)
1761                       ENDDO
1762                    ENDDO
1763                   
1764                CASE ( 'usm_t_surf_green' )
1765!--                 surface temperature for green surfaces
1766                    DO  m = 1, surf_usm_h%ns
1767                       surf_usm_h%t_surf_green_av(m) =                               & 
1768                                          surf_usm_h%t_surf_green_av(m) +            &
1769                                          t_surf_green_h(m)
1770                    ENDDO
1771                    DO  l = 0, 3
1772                       DO  m = 1, surf_usm_v(l)%ns
1773                          surf_usm_v(l)%t_surf_green_av(m) =                         &
1774                                          surf_usm_v(l)%t_surf_green_av(m) +         &
1775                                          t_surf_green_v(l)%t(m)
1776                       ENDDO
1777                    ENDDO
1778               
1779                CASE ( 'usm_t_surf_whole' )
1780!--                 surface temperature for whole surfaces
1781                    DO  m = 1, surf_usm_h%ns
1782                       surf_usm_h%t_surf_whole_av(m) =                               & 
1783                                          surf_usm_h%t_surf_whole_av(m) +            &
1784                                          t_surf_whole_h(m)
1785                    ENDDO
1786                    DO  l = 0, 3
1787                       DO  m = 1, surf_usm_v(l)%ns
1788                          surf_usm_v(l)%t_surf_whole_av(m) =                         &
1789                                          surf_usm_v(l)%t_surf_whole_av(m) +         &
1790                                          t_surf_whole_v(l)%t(m)
1791                       ENDDO
1792                    ENDDO
1793                   
1794                CASE ( 'usm_t_surf_10cm' )
1795!--                 near surface temperature for whole surfaces
1796                    DO  m = 1, surf_usm_h%ns
1797                       surf_usm_h%t_surf_10cm_av(m) =                               & 
1798                                          surf_usm_h%t_surf_10cm_av(m) +            &
1799                                          t_surf_10cm_h(m)
1800                    ENDDO
1801                    DO  l = 0, 3
1802                       DO  m = 1, surf_usm_v(l)%ns
1803                          surf_usm_v(l)%t_surf_10cm_av(m) =                         &
1804                                          surf_usm_v(l)%t_surf_10cm_av(m) +         &
1805                                          t_surf_10cm_v(l)%t(m)
1806                       ENDDO
1807                    ENDDO
1808
1809                   
1810                CASE ( 'usm_t_wall' )
1811!--                 wall temperature for  iwl layer of walls and land
1812                    DO  m = 1, surf_usm_h%ns
1813                       surf_usm_h%t_wall_av(iwl,m) =                           &
1814                                          surf_usm_h%t_wall_av(iwl,m) +        &
1815                                          t_wall_h(iwl,m)
1816                    ENDDO
1817                    DO  l = 0, 3
1818                       DO  m = 1, surf_usm_v(l)%ns
1819                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
1820                                          surf_usm_v(l)%t_wall_av(iwl,m) +     &
1821                                          t_wall_v(l)%t(iwl,m)
1822                       ENDDO
1823                    ENDDO
1824                   
1825                CASE ( 'usm_t_window' )
1826!--                 window temperature for  iwl layer of walls and land
1827                    DO  m = 1, surf_usm_h%ns
1828                       surf_usm_h%t_window_av(iwl,m) =                           &
1829                                          surf_usm_h%t_window_av(iwl,m) +        &
1830                                          t_window_h(iwl,m)
1831                    ENDDO
1832                    DO  l = 0, 3
1833                       DO  m = 1, surf_usm_v(l)%ns
1834                          surf_usm_v(l)%t_window_av(iwl,m) =                     &
1835                                          surf_usm_v(l)%t_window_av(iwl,m) +     &
1836                                          t_window_v(l)%t(iwl,m)
1837                       ENDDO
1838                    ENDDO
1839
1840                CASE ( 'usm_t_green' )
1841!--                 green temperature for  iwl layer of walls and land
1842                    DO  m = 1, surf_usm_h%ns
1843                       surf_usm_h%t_green_av(iwl,m) =                           &
1844                                          surf_usm_h%t_green_av(iwl,m) +        &
1845                                          t_green_h(iwl,m)
1846                    ENDDO
1847                    DO  l = 0, 3
1848                       DO  m = 1, surf_usm_v(l)%ns
1849                          surf_usm_v(l)%t_green_av(iwl,m) =                     &
1850                                          surf_usm_v(l)%t_green_av(iwl,m) +     &
1851                                          t_green_v(l)%t(iwl,m)
1852                       ENDDO
1853                    ENDDO
1854
1855                CASE DEFAULT
1856                    CONTINUE
1857
1858           END SELECT
1859
1860        ELSEIF ( mode == 'average' )  THEN
1861           
1862           SELECT CASE ( TRIM( var ) )
1863               
1864                CASE ( 'usm_rad_net' )
1865!--                 array of complete radiation balance
1866                    DO  m = 1, surf_usm_h%ns
1867                       surf_usm_h%rad_net_av(m) =                              &
1868                                          surf_usm_h%rad_net_av(m) /           &
1869                                          REAL( average_count_3d, kind=wp )
1870                    ENDDO
1871                    DO  l = 0, 3
1872                       DO  m = 1, surf_usm_v(l)%ns
1873                          surf_usm_v(l)%rad_net_av(m) =                        &
1874                                          surf_usm_v(l)%rad_net_av(m) /        &
1875                                          REAL( average_count_3d, kind=wp )
1876                       ENDDO
1877                    ENDDO
1878                   
1879                CASE ( 'usm_rad_insw' )
1880!--                 array of sw radiation falling to surface after i-th reflection
1881                    DO l = startenergy, endenergy
1882                        IF ( surfl(id,l) == ids )  THEN
1883                            surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
1884                        ENDIF
1885                    ENDDO
1886                             
1887                CASE ( 'usm_rad_inlw' )
1888!--                 array of lw radiation falling to surface after i-th reflection
1889                    DO l = startenergy, endenergy
1890                        IF ( surfl(id,l) == ids )  THEN
1891                            surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
1892                        ENDIF
1893                    ENDDO
1894                   
1895                CASE ( 'usm_rad_inswdir' )
1896!--                 array of direct sw radiation falling to surface from sun
1897                    DO l = startenergy, endenergy
1898                        IF ( surfl(id,l) == ids )  THEN
1899                            surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
1900                        ENDIF
1901                    ENDDO
1902                   
1903                CASE ( 'usm_rad_inswdif' )
1904!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
1905                    DO l = startenergy, endenergy
1906                        IF ( surfl(id,l) == ids )  THEN
1907                            surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
1908                        ENDIF
1909                    ENDDO
1910                   
1911                CASE ( 'usm_rad_inswref' )
1912!--                 array of sw radiation falling to surface from reflections
1913                    DO l = startenergy, endenergy
1914                        IF ( surfl(id,l) == ids )  THEN
1915                            surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
1916                        ENDIF
1917                    ENDDO
1918                   
1919                CASE ( 'usm_rad_inlwdif' )
1920!--                 array of sw radiation falling to surface after i-th reflection
1921                    DO l = startenergy, endenergy
1922                        IF ( surfl(id,l) == ids )  THEN
1923                            surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
1924                        ENDIF
1925                    ENDDO
1926                   
1927                CASE ( 'usm_rad_inlwref' )
1928!--                 array of lw radiation falling to surface from reflections
1929                    DO l = startenergy, endenergy
1930                        IF ( surfl(id,l) == ids )  THEN
1931                            surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
1932                        ENDIF
1933                    ENDDO
1934                   
1935                CASE ( 'usm_rad_outsw' )
1936!--                 array of sw radiation emitted from surface after i-th reflection
1937                    DO l = startenergy, endenergy
1938                        IF ( surfl(id,l) == ids )  THEN
1939                            surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
1940                        ENDIF
1941                    ENDDO
1942                   
1943                CASE ( 'usm_rad_outlw' )
1944!--                 array of lw radiation emitted from surface after i-th reflection
1945                    DO l = startenergy, endenergy
1946                        IF ( surfl(id,l) == ids )  THEN
1947                            surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
1948                        ENDIF
1949                    ENDDO
1950                   
1951                CASE ( 'usm_rad_ressw' )
1952!--                 array of residua of sw radiation absorbed in surface after last reflection
1953                    DO l = startenergy, endenergy
1954                        IF ( surfl(id,l) == ids )  THEN
1955                            surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
1956                        ENDIF
1957                    ENDDO
1958                                   
1959                CASE ( 'usm_rad_reslw' )
1960!--                 array of residua of lw radiation absorbed in surface after last reflection
1961                    DO l = startenergy, endenergy
1962                        IF ( surfl(id,l) == ids )  THEN
1963                            surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
1964                        ENDIF
1965                    ENDDO
1966                   
1967                CASE ( 'usm_rad_hf' )
1968!--                 array of heat flux from radiation for surfaces after i-th reflection
1969                    DO  m = 1, surf_usm_h%ns
1970                       surf_usm_h%surfhf_av(m) =                               &
1971                                          surf_usm_h%surfhf_av(m) /            &
1972                                          REAL( average_count_3d, kind=wp )
1973                    ENDDO
1974                    DO  l = 0, 3
1975                       DO  m = 1, surf_usm_v(l)%ns
1976                          surf_usm_v(l)%surfhf_av(m) =                         &
1977                                          surf_usm_v(l)%surfhf_av(m) /         &
1978                                          REAL( average_count_3d, kind=wp )
1979                       ENDDO
1980                    ENDDO
1981                   
1982                CASE ( 'usm_wshf' )
1983!--                 array of sensible heat flux from surfaces (land, roof, wall)
1984                    DO  m = 1, surf_usm_h%ns
1985                       surf_usm_h%wshf_eb_av(m) =                              &
1986                                          surf_usm_h%wshf_eb_av(m) /           &
1987                                          REAL( average_count_3d, kind=wp )
1988                    ENDDO
1989                    DO  l = 0, 3
1990                       DO  m = 1, surf_usm_v(l)%ns
1991                          surf_usm_v(l)%wshf_eb_av(m) =                        &
1992                                          surf_usm_v(l)%wshf_eb_av(m) /        &
1993                                          REAL( average_count_3d, kind=wp )
1994                       ENDDO
1995                    ENDDO
1996                   
1997                CASE ( 'usm_wghf' )
1998!--                 array of heat flux from ground (wall, roof, land)
1999                    DO  m = 1, surf_usm_h%ns
2000                       surf_usm_h%wghf_eb_av(m) =                              &
2001                                          surf_usm_h%wghf_eb_av(m) /           &
2002                                          REAL( average_count_3d, kind=wp )
2003                    ENDDO
2004                    DO  l = 0, 3
2005                       DO  m = 1, surf_usm_v(l)%ns
2006                          surf_usm_v(l)%wghf_eb_av(m) =                        &
2007                                          surf_usm_v(l)%wghf_eb_av(m) /        &
2008                                          REAL( average_count_3d, kind=wp )
2009                       ENDDO
2010                    ENDDO
2011                   
2012                CASE ( 'usm_wghf_window' )
2013!--                 array of heat flux from window ground (wall, roof, land)
2014                    DO  m = 1, surf_usm_h%ns
2015                       surf_usm_h%wghf_eb_window_av(m) =                              &
2016                                          surf_usm_h%wghf_eb_window_av(m) /           &
2017                                          REAL( average_count_3d, kind=wp )
2018                    ENDDO
2019                    DO  l = 0, 3
2020                       DO  m = 1, surf_usm_v(l)%ns
2021                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
2022                                          surf_usm_v(l)%wghf_eb_window_av(m) /        &
2023                                          REAL( average_count_3d, kind=wp )
2024                       ENDDO
2025                    ENDDO
2026
2027                CASE ( 'usm_wghf_green' )
2028!--                 array of heat flux from green ground (wall, roof, land)
2029                    DO  m = 1, surf_usm_h%ns
2030                       surf_usm_h%wghf_eb_green_av(m) =                              &
2031                                          surf_usm_h%wghf_eb_green_av(m) /           &
2032                                          REAL( average_count_3d, kind=wp )
2033                    ENDDO
2034                    DO  l = 0, 3
2035                       DO  m = 1, surf_usm_v(l)%ns
2036                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
2037                                          surf_usm_v(l)%wghf_eb_green_av(m) /        &
2038                                          REAL( average_count_3d, kind=wp )
2039                       ENDDO
2040                    ENDDO
2041
2042                CASE ( 'usm_iwghf' )
2043!--                 array of heat flux from indoor ground (wall, roof, land)
2044                    DO  m = 1, surf_usm_h%ns
2045                       surf_usm_h%iwghf_eb_av(m) =                              &
2046                                          surf_usm_h%iwghf_eb_av(m) /           &
2047                                          REAL( average_count_3d, kind=wp )
2048                    ENDDO
2049                    DO  l = 0, 3
2050                       DO  m = 1, surf_usm_v(l)%ns
2051                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
2052                                          surf_usm_v(l)%iwghf_eb_av(m) /        &
2053                                          REAL( average_count_3d, kind=wp )
2054                       ENDDO
2055                    ENDDO
2056                   
2057                CASE ( 'usm_iwghf_window' )
2058!--                 array of heat flux from indoor window ground (wall, roof, land)
2059                    DO  m = 1, surf_usm_h%ns
2060                       surf_usm_h%iwghf_eb_window_av(m) =                              &
2061                                          surf_usm_h%iwghf_eb_window_av(m) /           &
2062                                          REAL( average_count_3d, kind=wp )
2063                    ENDDO
2064                    DO  l = 0, 3
2065                       DO  m = 1, surf_usm_v(l)%ns
2066                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
2067                                          surf_usm_v(l)%iwghf_eb_window_av(m) /        &
2068                                          REAL( average_count_3d, kind=wp )
2069                       ENDDO
2070                    ENDDO
2071                   
2072                CASE ( 'usm_t_surf' )
2073!--                 surface temperature for surfaces
2074                    DO  m = 1, surf_usm_h%ns
2075                       surf_usm_h%t_surf_av(m) =                               & 
2076                                          surf_usm_h%t_surf_av(m) /            &
2077                                          REAL( average_count_3d, kind=wp )
2078                    ENDDO
2079                    DO  l = 0, 3
2080                       DO  m = 1, surf_usm_v(l)%ns
2081                          surf_usm_v(l)%t_surf_av(m) =                         &
2082                                          surf_usm_v(l)%t_surf_av(m) /         &
2083                                          REAL( average_count_3d, kind=wp )
2084                       ENDDO
2085                    ENDDO
2086                   
2087                CASE ( 'usm_t_surf_window' )
2088!--                 surface temperature for window surfaces
2089                    DO  m = 1, surf_usm_h%ns
2090                       surf_usm_h%t_surf_window_av(m) =                               & 
2091                                          surf_usm_h%t_surf_window_av(m) /            &
2092                                          REAL( average_count_3d, kind=wp )
2093                    ENDDO
2094                    DO  l = 0, 3
2095                       DO  m = 1, surf_usm_v(l)%ns
2096                          surf_usm_v(l)%t_surf_window_av(m) =                         &
2097                                          surf_usm_v(l)%t_surf_window_av(m) /         &
2098                                          REAL( average_count_3d, kind=wp )
2099                       ENDDO
2100                    ENDDO
2101                   
2102                CASE ( 'usm_t_surf_green' )
2103!--                 surface temperature for green surfaces
2104                    DO  m = 1, surf_usm_h%ns
2105                       surf_usm_h%t_surf_green_av(m) =                               & 
2106                                          surf_usm_h%t_surf_green_av(m) /            &
2107                                          REAL( average_count_3d, kind=wp )
2108                    ENDDO
2109                    DO  l = 0, 3
2110                       DO  m = 1, surf_usm_v(l)%ns
2111                          surf_usm_v(l)%t_surf_green_av(m) =                         &
2112                                          surf_usm_v(l)%t_surf_green_av(m) /         &
2113                                          REAL( average_count_3d, kind=wp )
2114                       ENDDO
2115                    ENDDO
2116                   
2117                CASE ( 'usm_t_surf_whole' )
2118!--                 surface temperature for whole surfaces
2119                    DO  m = 1, surf_usm_h%ns
2120                       surf_usm_h%t_surf_whole_av(m) =                               & 
2121                                          surf_usm_h%t_surf_whole_av(m) /            &
2122                                          REAL( average_count_3d, kind=wp )
2123                    ENDDO
2124                    DO  l = 0, 3
2125                       DO  m = 1, surf_usm_v(l)%ns
2126                          surf_usm_v(l)%t_surf_whole_av(m) =                         &
2127                                          surf_usm_v(l)%t_surf_whole_av(m) /         &
2128                                          REAL( average_count_3d, kind=wp )
2129                       ENDDO
2130                    ENDDO
2131                   
2132                CASE ( 'usm_t_surf_10cm' )
2133!--                 near surface temperature for whole surfaces
2134                    DO  m = 1, surf_usm_h%ns
2135                       surf_usm_h%t_surf_10cm_av(m) =                               & 
2136                                          surf_usm_h%t_surf_10cm_av(m) /            &
2137                                          REAL( average_count_3d, kind=wp )
2138                    ENDDO
2139                    DO  l = 0, 3
2140                       DO  m = 1, surf_usm_v(l)%ns
2141                          surf_usm_v(l)%t_surf_10cm_av(m) =                         &
2142                                          surf_usm_v(l)%t_surf_10cm_av(m) /         &
2143                                          REAL( average_count_3d, kind=wp )
2144                       ENDDO
2145                    ENDDO
2146                   
2147                CASE ( 'usm_t_wall' )
2148!--                 wall temperature for  iwl layer of walls and land
2149                    DO  m = 1, surf_usm_h%ns
2150                       surf_usm_h%t_wall_av(iwl,m) =                           &
2151                                          surf_usm_h%t_wall_av(iwl,m) /        &
2152                                          REAL( average_count_3d, kind=wp )
2153                    ENDDO
2154                    DO  l = 0, 3
2155                       DO  m = 1, surf_usm_v(l)%ns
2156                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
2157                                          surf_usm_v(l)%t_wall_av(iwl,m) /     &
2158                                          REAL( average_count_3d, kind=wp )
2159                       ENDDO
2160                    ENDDO
2161
2162                CASE ( 'usm_t_window' )
2163!--                 window temperature for  iwl layer of walls and land
2164                    DO  m = 1, surf_usm_h%ns
2165                       surf_usm_h%t_window_av(iwl,m) =                           &
2166                                          surf_usm_h%t_window_av(iwl,m) /        &
2167                                          REAL( average_count_3d, kind=wp )
2168                    ENDDO
2169                    DO  l = 0, 3
2170                       DO  m = 1, surf_usm_v(l)%ns
2171                          surf_usm_v(l)%t_window_av(iwl,m) =                     &
2172                                          surf_usm_v(l)%t_window_av(iwl,m) /     &
2173                                          REAL( average_count_3d, kind=wp )
2174                       ENDDO
2175                    ENDDO
2176
2177                CASE ( 'usm_t_green' )
2178!--                 green temperature for  iwl layer of walls and land
2179                    DO  m = 1, surf_usm_h%ns
2180                       surf_usm_h%t_green_av(iwl,m) =                           &
2181                                          surf_usm_h%t_green_av(iwl,m) /        &
2182                                          REAL( average_count_3d, kind=wp )
2183                    ENDDO
2184                    DO  l = 0, 3
2185                       DO  m = 1, surf_usm_v(l)%ns
2186                          surf_usm_v(l)%t_green_av(iwl,m) =                     &
2187                                          surf_usm_v(l)%t_green_av(iwl,m) /     &
2188                                          REAL( average_count_3d, kind=wp )
2189                       ENDDO
2190                    ENDDO
2191
2192
2193           END SELECT
2194
2195        ENDIF
2196
2197    END SUBROUTINE usm_average_3d_data
2198
2199
2200
2201!------------------------------------------------------------------------------!
2202! Description:
2203! ------------
2204!> Set internal Neumann boundary condition at outer soil grid points
2205!> for temperature and humidity.
2206!------------------------------------------------------------------------------!
2207 SUBROUTINE usm_boundary_condition
2208 
2209    IMPLICIT NONE
2210
2211    INTEGER(iwp) :: i      !< grid index x-direction
2212    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
2213    INTEGER(iwp) :: j      !< grid index y-direction
2214    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
2215    INTEGER(iwp) :: k      !< grid index z-direction
2216    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
2217    INTEGER(iwp) :: l      !< running index surface-orientation
2218    INTEGER(iwp) :: m      !< running index surface elements
2219
2220    koff = surf_usm_h%koff
2221    DO  m = 1, surf_usm_h%ns
2222       i = surf_usm_h%i(m)
2223       j = surf_usm_h%j(m)
2224       k = surf_usm_h%k(m)
2225       pt(k+koff,j,i) = pt(k,j,i)
2226    ENDDO
2227
2228    DO  l = 0, 3
2229       ioff = surf_usm_v(l)%ioff
2230       joff = surf_usm_v(l)%joff
2231       DO  m = 1, surf_usm_v(l)%ns
2232          i = surf_usm_v(l)%i(m)
2233          j = surf_usm_v(l)%j(m)
2234          k = surf_usm_v(l)%k(m)
2235          pt(k,j+joff,i+ioff) = pt(k,j,i)
2236       ENDDO
2237    ENDDO
2238
2239 END SUBROUTINE usm_boundary_condition
2240
2241
2242!------------------------------------------------------------------------------!
2243!
2244! Description:
2245! ------------
2246!> Subroutine checks variables and assigns units.
2247!> It is called out from subroutine check_parameters.
2248!------------------------------------------------------------------------------!
2249    SUBROUTINE usm_check_data_output( variable, unit )
2250       
2251        IMPLICIT NONE
2252 
2253        CHARACTER (len=*),INTENT(IN)    ::  variable !:
2254        CHARACTER (len=*),INTENT(OUT)   ::  unit     !:
2255       
2256        CHARACTER (len=varnamelength)   :: var
2257
2258        var = TRIM(variable)
2259        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.        &
2260             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.    &
2261             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR. &
2262             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR. &
2263             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.     &
2264             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.     &
2265             var(1:11) == 'usm_rad_hf_'  .OR.                                             &
2266             var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
2267             var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
2268             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' )  THEN
2269            unit = 'W/m2'
2270        ELSE IF ( var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall' .OR.         &
2271                  var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR.  &
2272                  var(1:16) == 'usm_t_surf_green'  .OR.                                   &
2273                  var(1:16) == 'usm_t_surf_whole' .OR. var(1:11) == 'usm_t_green' .OR.    &
2274                  var(1:15) == 'usm_t_surf_10cm')  THEN
2275            unit = 'K'
2276        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              & 
2277                  var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
2278                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis')  THEN
2279            unit = '1'
2280        ELSE
2281            unit = 'illegal'
2282        ENDIF
2283
2284    END SUBROUTINE usm_check_data_output
2285
2286
2287!------------------------------------------------------------------------------!
2288! Description:
2289! ------------
2290!> Check parameters routine for urban surface model
2291!------------------------------------------------------------------------------!
2292    SUBROUTINE usm_check_parameters
2293   
2294       USE control_parameters,                                                 &
2295           ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing,   &
2296                  lsf_surf, topography
2297
2298!
2299!--    Dirichlet boundary conditions are required as the surface fluxes are
2300!--    calculated from the temperature/humidity gradients in the urban surface
2301!--    model
2302       IF ( bc_pt_b == 'neumann'   .OR.   bc_q_b == 'neumann' )  THEN
2303          message_string = 'urban surface model requires setting of '//        &
2304                           'bc_pt_b = "dirichlet" and '//                      &
2305                           'bc_q_b  = "dirichlet"'
2306          CALL message( 'check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
2307       ENDIF
2308
2309       IF ( .NOT.  constant_flux_layer )  THEN
2310          message_string = 'urban surface model requires '//                   &
2311                           'constant_flux_layer = .T.'
2312          CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
2313       ENDIF
2314!       
2315!--    Surface forcing has to be disabled for LSF in case of enabled
2316!--    urban surface module
2317       IF ( large_scale_forcing )  THEN
2318          lsf_surf = .FALSE.
2319       ENDIF
2320!
2321!--    Topography
2322       IF ( topography == 'flat' )  THEN
2323          message_string = 'topography /= "flat" is required '//               &
2324                           'when using the urban surface model'
2325          CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2326       ENDIF
2327
2328    END SUBROUTINE usm_check_parameters
2329
2330
2331!------------------------------------------------------------------------------!
2332!
2333! Description:
2334! ------------
2335!> Output of the 3D-arrays in netCDF and/or AVS format
2336!> for variables of urban_surface model.
2337!> It resorts the urban surface module output quantities from surf style
2338!> indexing into temporary 3D array with indices (i,j,k).
2339!> It is called from subroutine data_output_3d.
2340!------------------------------------------------------------------------------!
2341    SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2342       
2343        IMPLICIT NONE
2344
2345        INTEGER(iwp), INTENT(IN)       ::  av        !<
2346        CHARACTER (len=*), INTENT(IN)  ::  variable  !<
2347        INTEGER(iwp), INTENT(IN)       ::  nzb_do    !< lower limit of the data output (usually 0)
2348        INTEGER(iwp), INTENT(IN)       ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
2349        LOGICAL, INTENT(OUT)           ::  found     !<
2350        REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< sp - it has to correspond to module data_output_3d
2351        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)     ::  temp_pf    !< temp array for urban surface output procedure
2352       
2353        CHARACTER (len=varnamelength)                          :: var, surfid
2354        INTEGER(iwp), PARAMETER                                :: nd = 5
2355        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER         :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
2356        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
2357        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
2358        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
2359        INTEGER(iwp)                                           :: ids,isurf,isvf,isurfs,isurflt
2360        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
2361        INTEGER(iwp)                                           ::  k_topo    !< topography top index
2362
2363        dirstart = (/ startland, startwall, startwall, startwall, startwall /)
2364        dirend = (/ endland, endwall, endwall, endwall, endwall /)
2365
2366        found = .TRUE.
2367        temp_pf = -1._wp
2368       
2369        ids = -1
2370        var = TRIM(variable)
2371        DO i = 0, nd-1
2372            k = len(TRIM(var))
2373            j = len(TRIM(dirname(i)))
2374            IF ( var(k-j+1:k) == dirname(i) )  THEN
2375                ids = i
2376                var = var(:k-j)
2377                EXIT
2378            ENDIF
2379        ENDDO
2380        IF ( ids == -1 )  THEN
2381            var = TRIM(variable)
2382        ENDIF
2383        IF ( var(1:11) == 'usm_t_wall_'  .AND.  len(TRIM(var)) >= 12 )  THEN
2384!--         wall layers
2385            READ(var(12:12), '(I1)', iostat=istat ) iwl
2386            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
2387                var = var(1:10)
2388            ENDIF
2389        ENDIF
2390        IF ( var(1:13) == 'usm_t_window_'  .AND.  len(TRIM(var)) >= 14 )  THEN
2391!--         window layers
2392            READ(var(14:14), '(I1)', iostat=istat ) iwl
2393            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
2394                var = var(1:12)
2395            ENDIF
2396        ENDIF
2397        IF ( var(1:12) == 'usm_t_green_'  .AND.  len(TRIM(var)) >= 13 )  THEN
2398!--         green layers
2399            READ(var(13:13), '(I1)', iostat=istat ) iwl
2400            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
2401                var = var(1:11)
2402            ENDIF
2403        ENDIF
2404        IF ( (var(1:8) == 'usm_svf_'  .OR.  var(1:8) == 'usm_dif_')  .AND.  len(TRIM(var)) >= 13 )  THEN
2405!--         svf values to particular surface
2406            surfid = var(9:)
2407            i = index(surfid,'_')
2408            j = index(surfid(i+1:),'_')
2409            READ(surfid(1:i-1),*, iostat=istat ) is
2410            IF ( istat == 0 )  THEN
2411                READ(surfid(i+1:i+j-1),*, iostat=istat ) js
2412            ENDIF
2413            IF ( istat == 0 )  THEN
2414                READ(surfid(i+j+1:),*, iostat=istat ) ks
2415            ENDIF
2416            IF ( istat == 0 )  THEN
2417                var = var(1:7)
2418            ENDIF
2419        ENDIF
2420       
2421        SELECT CASE ( TRIM(var) )
2422
2423          CASE ( 'usm_surfz' )
2424!--           array of lw radiation falling to local surface after i-th reflection
2425              DO  m = 1, surf_usm_h%ns
2426                 i = surf_usm_h%i(m)
2427                 j = surf_usm_h%j(m)
2428                 k = surf_usm_h%k(m)
2429                 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
2430              ENDDO
2431              DO  l = 0, 3
2432                 DO  m = 1, surf_usm_v(l)%ns
2433                    i = surf_usm_v(l)%i(m)
2434                    j = surf_usm_v(l)%j(m)
2435                    k = surf_usm_v(l)%k(m)
2436                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
2437                 ENDDO
2438              ENDDO
2439
2440          CASE ( 'usm_surfcat' )
2441!--           surface category
2442              DO  m = 1, surf_usm_h%ns
2443                 i = surf_usm_h%i(m)
2444                 j = surf_usm_h%j(m)
2445                 k = surf_usm_h%k(m)
2446                 temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2447              ENDDO
2448              DO  l = 0, 3
2449                 DO  m = 1, surf_usm_v(l)%ns
2450                    i = surf_usm_v(l)%i(m)
2451                    j = surf_usm_v(l)%j(m)
2452                    k = surf_usm_v(l)%k(m)
2453                    temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
2454                 ENDDO
2455              ENDDO
2456             
2457          CASE ( 'usm_surfalb' )
2458!--           surface albedo, weighted average
2459              DO  m = 1, surf_usm_h%ns
2460                 i = surf_usm_h%i(m)
2461                 j = surf_usm_h%j(m)
2462                 k = surf_usm_h%k(m)
2463                 temp_pf(k,j,i) = surf_usm_h%frac(0,m)     *                   &
2464                                  surf_usm_h%albedo(0,m) +                     &
2465                                  surf_usm_h%frac(1,m)     *                   &
2466                                  surf_usm_h%albedo(1,m) +                     &
2467                                  surf_usm_h%frac(2,m)     *                   &
2468                                  surf_usm_h%albedo(2,m)
2469              ENDDO
2470              DO  l = 0, 3
2471                 DO  m = 1, surf_usm_v(l)%ns
2472                    i = surf_usm_v(l)%i(m)
2473                    j = surf_usm_v(l)%j(m)
2474                    k = surf_usm_v(l)%k(m)
2475                    temp_pf(k,j,i) = surf_usm_v(l)%frac(0,m)     *             &
2476                                     surf_usm_v(l)%albedo(0,m) +               &
2477                                     surf_usm_v(l)%frac(1,m)     *             &
2478                                     surf_usm_v(l)%albedo(1,m) +               &
2479                                     surf_usm_v(l)%frac(2,m)     *             &
2480                                     surf_usm_v(l)%albedo(2,m)
2481                 ENDDO
2482              ENDDO
2483             
2484          CASE ( 'usm_surfemis' )
2485!--           surface emissivity, weighted average
2486              DO  m = 1, surf_usm_h%ns
2487                 i = surf_usm_h%i(m)
2488                 j = surf_usm_h%j(m)
2489                 k = surf_usm_h%k(m)
2490                 temp_pf(k,j,i) =  surf_usm_h%frac(0,m)     *                  &
2491                                   surf_usm_h%emissivity(0,m) +                &
2492                                   surf_usm_h%frac(1,m)     *                  &
2493                                   surf_usm_h%emissivity(1,m) +                &
2494                                   surf_usm_h%frac(2,m)     *                  &
2495                                   surf_usm_h%emissivity(2,m)
2496              ENDDO
2497              DO  l = 0, 3
2498                 DO  m = 1, surf_usm_v(l)%ns
2499                    i = surf_usm_v(l)%i(m)
2500                    j = surf_usm_v(l)%j(m)
2501                    k = surf_usm_v(l)%k(m)
2502                    temp_pf(k,j,i) = surf_usm_v(l)%frac(0,m)       *           &
2503                                     surf_usm_v(l)%emissivity(0,m) +           &
2504                                     surf_usm_v(l)%frac(1,m) *                 &
2505                                     surf_usm_v(l)%emissivity(1,m) +           &
2506                                     surf_usm_v(l)%frac(2,m)     *             &
2507                                     surf_usm_v(l)%emissivity(2,m)
2508                 ENDDO
2509              ENDDO
2510
2511          CASE ( 'usm_surfwintrans' )
2512!--           transmissivity window tiles
2513              DO  m = 1, surf_usm_h%ns
2514                 i = surf_usm_h%i(m)
2515                 j = surf_usm_h%j(m)
2516                 k = surf_usm_h%k(m)
2517                 temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
2518              ENDDO
2519              DO  l = 0, 3
2520                 DO  m = 1, surf_usm_v(l)%ns
2521                    i = surf_usm_v(l)%i(m)
2522                    j = surf_usm_v(l)%j(m)
2523                    k = surf_usm_v(l)%k(m)
2524                    temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
2525                 ENDDO
2526
2527              ENDDO
2528
2529!
2530!-- Not adjusted so far             
2531          CASE ( 'usm_svf', 'usm_dif' )
2532!--           shape view factors or iradiance factors to selected surface
2533              IF ( TRIM(var)=='usm_svf' )  THEN
2534                  k = 1
2535              ELSE
2536                  k = 2
2537              ENDIF
2538              DO isvf = 1, nsvfl
2539                  isurflt = svfsurf(1, isvf)
2540                  isurfs = svfsurf(2, isvf)
2541                             
2542                  IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND.       &
2543                       surf(iz,isurfs) == ks  .AND.  surf(id,isurfs) == ids )  THEN
2544  !--                 correct source surface
2545                      temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
2546                  ENDIF
2547              ENDDO
2548
2549          CASE ( 'usm_rad_net' )
2550!--           array of complete radiation balance
2551              IF ( av == 0 )  THEN
2552                 DO  m = 1, surf_usm_h%ns
2553                    i = surf_usm_h%i(m)
2554                    j = surf_usm_h%j(m)
2555                    k = surf_usm_h%k(m)
2556                    temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
2557                 ENDDO
2558                 DO  l = 0, 3
2559                    DO  m = 1, surf_usm_v(l)%ns
2560                       i = surf_usm_v(l)%i(m)
2561                       j = surf_usm_v(l)%j(m)
2562                       k = surf_usm_v(l)%k(m)
2563                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
2564                    ENDDO
2565                 ENDDO
2566              ELSE
2567                 DO  m = 1, surf_usm_h%ns
2568                    i = surf_usm_h%i(m)
2569                    j = surf_usm_h%j(m)
2570                    k = surf_usm_h%k(m)
2571                    temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
2572                 ENDDO
2573                 DO  l = 0, 3
2574                    DO  m = 1, surf_usm_v(l)%ns
2575                       i = surf_usm_v(l)%i(m)
2576                       j = surf_usm_v(l)%j(m)
2577                       k = surf_usm_v(l)%k(m)
2578                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
2579                    ENDDO
2580                 ENDDO
2581              ENDIF
2582
2583          CASE ( 'usm_rad_insw' )
2584!--           array of sw radiation falling to surface after i-th reflection
2585              DO isurf = dirstart(ids), dirend(ids)
2586                 IF ( surfl(id,isurf) == ids )  THEN
2587                   IF ( av == 0 )  THEN
2588                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
2589                   ELSE
2590                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
2591                   ENDIF
2592                 ENDIF
2593              ENDDO
2594
2595          CASE ( 'usm_rad_inlw' )
2596!--           array of lw radiation falling to surface after i-th reflection
2597              DO isurf = dirstart(ids), dirend(ids)
2598                 IF ( surfl(id,isurf) == ids )  THEN
2599                   IF ( av == 0 )  THEN
2600                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
2601                   ELSE
2602                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
2603                   ENDIF
2604                 ENDIF
2605              ENDDO
2606
2607          CASE ( 'usm_rad_inswdir' )
2608!--           array of direct sw radiation falling to surface from sun
2609              DO isurf = dirstart(ids), dirend(ids)
2610                 IF ( surfl(id,isurf) == ids )  THEN
2611                   IF ( av == 0 )  THEN
2612                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
2613                   ELSE
2614                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
2615                   ENDIF
2616                 ENDIF
2617              ENDDO
2618
2619          CASE ( 'usm_rad_inswdif' )
2620!--           array of difusion sw radiation falling to surface from sky and borders of the domain
2621              DO isurf = dirstart(ids), dirend(ids)
2622                 IF ( surfl(id,isurf) == ids )  THEN
2623                   IF ( av == 0 )  THEN
2624                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
2625                   ELSE
2626                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
2627                   ENDIF
2628                 ENDIF
2629              ENDDO
2630
2631          CASE ( 'usm_rad_inswref' )
2632!--           array of sw radiation falling to surface from reflections
2633              DO isurf = dirstart(ids), dirend(ids)
2634                 IF ( surfl(id,isurf) == ids )  THEN
2635                   IF ( av == 0 )  THEN
2636                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
2637                       surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
2638                   ELSE
2639                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
2640                   ENDIF
2641                 ENDIF
2642              ENDDO
2643
2644          CASE ( 'usm_rad_inlwref' )
2645!--           array of lw radiation falling to surface from reflections
2646              DO isurf = dirstart(ids), dirend(ids)
2647                 IF ( surfl(id,isurf) == ids )  THEN
2648                   IF ( av == 0 )  THEN
2649                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
2650                   ELSE
2651                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
2652                   ENDIF
2653                 ENDIF
2654              ENDDO
2655
2656          CASE ( 'usm_rad_outsw' )
2657!--           array of sw radiation emitted from surface after i-th reflection
2658              DO isurf = dirstart(ids), dirend(ids)
2659                 IF ( surfl(id,isurf) == ids )  THEN
2660                   IF ( av == 0 )  THEN
2661                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
2662                   ELSE
2663                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
2664                   ENDIF
2665                 ENDIF
2666              ENDDO
2667
2668          CASE ( 'usm_rad_outlw' )
2669!--           array of lw radiation emitted from surface after i-th reflection
2670              DO isurf = dirstart(ids), dirend(ids)
2671                 IF ( surfl(id,isurf) == ids )  THEN
2672                   IF ( av == 0 )  THEN
2673                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
2674                   ELSE
2675                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
2676                   ENDIF
2677                 ENDIF
2678              ENDDO
2679
2680          CASE ( 'usm_rad_ressw' )
2681!--           average of array of residua of sw radiation absorbed in surface after last reflection
2682              DO isurf = dirstart(ids), dirend(ids)
2683                 IF ( surfl(id,isurf) == ids )  THEN
2684                   IF ( av == 0 )  THEN
2685                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
2686                   ELSE
2687                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
2688                   ENDIF
2689                 ENDIF
2690              ENDDO
2691
2692          CASE ( 'usm_rad_reslw' )
2693!--           average of array of residua of lw radiation absorbed in surface after last reflection
2694              DO isurf = dirstart(ids), dirend(ids)
2695                 IF ( surfl(id,isurf) == ids )  THEN
2696                   IF ( av == 0 )  THEN
2697                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
2698                   ELSE
2699                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
2700                   ENDIF
2701                 ENDIF
2702              ENDDO
2703 
2704          CASE ( 'usm_rad_hf' )
2705!--           array of heat flux from radiation for surfaces after all reflections
2706              IF ( av == 0 )  THEN
2707                 DO  m = 1, surf_usm_h%ns
2708                    i = surf_usm_h%i(m)
2709                    j = surf_usm_h%j(m)
2710                    k = surf_usm_h%k(m)
2711                    temp_pf(k,j,i) = surf_usm_h%surfhf(m)
2712                 ENDDO
2713                 DO  l = 0, 3
2714                    DO  m = 1, surf_usm_v(l)%ns
2715                       i = surf_usm_v(l)%i(m)
2716                       j = surf_usm_v(l)%j(m)
2717                       k = surf_usm_v(l)%k(m)
2718                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
2719                    ENDDO
2720                 ENDDO
2721              ELSE
2722                 DO  m = 1, surf_usm_h%ns
2723                    i = surf_usm_h%i(m)
2724                    j = surf_usm_h%j(m)
2725                    k = surf_usm_h%k(m)
2726                    temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
2727                 ENDDO
2728                 DO  l = 0, 3
2729                    DO  m = 1, surf_usm_v(l)%ns
2730                       i = surf_usm_v(l)%i(m)
2731                       j = surf_usm_v(l)%j(m)
2732                       k = surf_usm_v(l)%k(m)
2733                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
2734                    ENDDO
2735                 ENDDO
2736              ENDIF
2737 
2738          CASE ( 'usm_wshf' )
2739!--           array of sensible heat flux from surfaces
2740              IF ( av == 0 )  THEN
2741                 DO  m = 1, surf_usm_h%ns
2742                    i = surf_usm_h%i(m)
2743                    j = surf_usm_h%j(m)
2744                    k = surf_usm_h%k(m)
2745                    temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2746                 ENDDO
2747                 DO  l = 0, 3
2748                    DO  m = 1, surf_usm_v(l)%ns
2749                       i = surf_usm_v(l)%i(m)
2750                       j = surf_usm_v(l)%j(m)
2751                       k = surf_usm_v(l)%k(m)
2752                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2753                    ENDDO
2754                 ENDDO
2755              ELSE
2756                 DO  m = 1, surf_usm_h%ns
2757                    i = surf_usm_h%i(m)
2758                    j = surf_usm_h%j(m)
2759                    k = surf_usm_h%k(m)
2760                    temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2761                 ENDDO
2762                 DO  l = 0, 3
2763                    DO  m = 1, surf_usm_v(l)%ns
2764                       i = surf_usm_v(l)%i(m)
2765                       j = surf_usm_v(l)%j(m)
2766                       k = surf_usm_v(l)%k(m)
2767                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2768                    ENDDO
2769                 ENDDO
2770              ENDIF
2771
2772
2773          CASE ( 'usm_wghf' )
2774!--           array of heat flux from ground (land, wall, roof)
2775              IF ( av == 0 )  THEN
2776                 DO  m = 1, surf_usm_h%ns
2777                    i = surf_usm_h%i(m)
2778                    j = surf_usm_h%j(m)
2779                    k = surf_usm_h%k(m)
2780                    temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2781                 ENDDO
2782                 DO  l = 0, 3
2783                    DO  m = 1, surf_usm_v(l)%ns
2784                       i = surf_usm_v(l)%i(m)
2785                       j = surf_usm_v(l)%j(m)
2786                       k = surf_usm_v(l)%k(m)
2787                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2788                    ENDDO
2789                 ENDDO
2790              ELSE
2791                 DO  m = 1, surf_usm_h%ns
2792                    i = surf_usm_h%i(m)
2793                    j = surf_usm_h%j(m)
2794                    k = surf_usm_h%k(m)
2795                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2796                 ENDDO
2797                 DO  l = 0, 3
2798                    DO  m = 1, surf_usm_v(l)%ns
2799                       i = surf_usm_v(l)%i(m)
2800                       j = surf_usm_v(l)%j(m)
2801                       k = surf_usm_v(l)%k(m)
2802                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2803                    ENDDO
2804                 ENDDO
2805              ENDIF
2806
2807          CASE ( 'usm_wghf_window' )
2808!--           array of heat flux from window ground (land, wall, roof)
2809
2810              IF ( av == 0 )  THEN
2811                 DO  m = 1, surf_usm_h%ns
2812                    i = surf_usm_h%i(m)
2813                    j = surf_usm_h%j(m)
2814                    k = surf_usm_h%k(m)
2815                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
2816                 ENDDO
2817                 DO  l = 0, 3
2818                    DO  m = 1, surf_usm_v(l)%ns
2819                       i = surf_usm_v(l)%i(m)
2820                       j = surf_usm_v(l)%j(m)
2821                       k = surf_usm_v(l)%k(m)
2822                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m)
2823                    ENDDO
2824                 ENDDO
2825              ELSE
2826                 DO  m = 1, surf_usm_h%ns
2827                    i = surf_usm_h%i(m)
2828                    j = surf_usm_h%j(m)
2829                    k = surf_usm_h%k(m)
2830                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
2831                 ENDDO
2832                 DO  l = 0, 3
2833                    DO  m = 1, surf_usm_v(l)%ns
2834                       i = surf_usm_v(l)%i(m)
2835                       j = surf_usm_v(l)%j(m)
2836                       k = surf_usm_v(l)%k(m)
2837                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
2838                    ENDDO
2839                 ENDDO
2840              ENDIF
2841
2842          CASE ( 'usm_wghf_green' )
2843!--           array of heat flux from green ground (land, wall, roof)
2844
2845              IF ( av == 0 )  THEN
2846                 DO  m = 1, surf_usm_h%ns
2847                    i = surf_usm_h%i(m)
2848                    j = surf_usm_h%j(m)
2849                    k = surf_usm_h%k(m)
2850                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
2851                 ENDDO
2852                 DO  l = 0, 3
2853                    DO  m = 1, surf_usm_v(l)%ns
2854                       i = surf_usm_v(l)%i(m)
2855                       j = surf_usm_v(l)%j(m)
2856                       k = surf_usm_v(l)%k(m)
2857                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
2858                    ENDDO
2859                 ENDDO
2860              ELSE
2861                 DO  m = 1, surf_usm_h%ns
2862                    i = surf_usm_h%i(m)
2863                    j = surf_usm_h%j(m)
2864                    k = surf_usm_h%k(m)
2865                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
2866                 ENDDO
2867                 DO  l = 0, 3
2868                    DO  m = 1, surf_usm_v(l)%ns
2869                       i = surf_usm_v(l)%i(m)
2870                       j = surf_usm_v(l)%j(m)
2871                       k = surf_usm_v(l)%k(m)
2872                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m)
2873                    ENDDO
2874                 ENDDO
2875              ENDIF
2876
2877          CASE ( 'usm_iwghf' )
2878!--           array of heat flux from indoor ground (land, wall, roof)
2879              IF ( av == 0 )  THEN
2880                 DO  m = 1, surf_usm_h%ns
2881                    i = surf_usm_h%i(m)
2882                    j = surf_usm_h%j(m)
2883                    k = surf_usm_h%k(m)
2884                    temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
2885                 ENDDO
2886                 DO  l = 0, 3
2887                    DO  m = 1, surf_usm_v(l)%ns
2888                       i = surf_usm_v(l)%i(m)
2889                       j = surf_usm_v(l)%j(m)
2890                       k = surf_usm_v(l)%k(m)
2891                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m)
2892                    ENDDO
2893                 ENDDO
2894              ELSE
2895                 DO  m = 1, surf_usm_h%ns
2896                    i = surf_usm_h%i(m)
2897                    j = surf_usm_h%j(m)
2898                    k = surf_usm_h%k(m)
2899                    temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
2900                 ENDDO
2901                 DO  l = 0, 3
2902                    DO  m = 1, surf_usm_v(l)%ns
2903                       i = surf_usm_v(l)%i(m)
2904                       j = surf_usm_v(l)%j(m)
2905                       k = surf_usm_v(l)%k(m)
2906                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m)
2907                    ENDDO
2908                 ENDDO
2909              ENDIF
2910
2911          CASE ( 'usm_iwghf_window' )
2912!--           array of heat flux from indoor window ground (land, wall, roof)
2913
2914              IF ( av == 0 )  THEN
2915                 DO  m = 1, surf_usm_h%ns
2916                    i = surf_usm_h%i(m)
2917                    j = surf_usm_h%j(m)
2918                    k = surf_usm_h%k(m)
2919                    temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
2920                 ENDDO
2921                 DO  l = 0, 3
2922                    DO  m = 1, surf_usm_v(l)%ns
2923                       i = surf_usm_v(l)%i(m)
2924                       j = surf_usm_v(l)%j(m)
2925                       k = surf_usm_v(l)%k(m)
2926                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m)
2927                    ENDDO
2928                 ENDDO
2929              ELSE
2930                 DO  m = 1, surf_usm_h%ns
2931                    i = surf_usm_h%i(m)
2932                    j = surf_usm_h%j(m)
2933                    k = surf_usm_h%k(m)
2934                    temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
2935                 ENDDO
2936                 DO  l = 0, 3
2937                    DO  m = 1, surf_usm_v(l)%ns
2938                       i = surf_usm_v(l)%i(m)
2939                       j = surf_usm_v(l)%j(m)
2940                       k = surf_usm_v(l)%k(m)
2941                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m)
2942                    ENDDO
2943                 ENDDO
2944              ENDIF
2945             
2946          CASE ( 'usm_t_surf' )
2947!--           surface temperature for surfaces
2948              IF ( av == 0 )  THEN
2949                 DO  m = 1, surf_usm_h%ns
2950                    i = surf_usm_h%i(m)
2951                    j = surf_usm_h%j(m)
2952                    k = surf_usm_h%k(m)
2953                    temp_pf(k,j,i) = t_surf_h(m)
2954                 ENDDO
2955                 DO  l = 0, 3
2956                    DO  m = 1, surf_usm_v(l)%ns
2957                       i = surf_usm_v(l)%i(m)
2958                       j = surf_usm_v(l)%j(m)
2959                       k = surf_usm_v(l)%k(m)
2960                       temp_pf(k,j,i) = t_surf_v(l)%t(m)
2961                    ENDDO
2962                 ENDDO
2963              ELSE
2964                 DO  m = 1, surf_usm_h%ns
2965                    i = surf_usm_h%i(m)
2966                    j = surf_usm_h%j(m)
2967                    k = surf_usm_h%k(m)
2968                    temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
2969                 ENDDO
2970                 DO  l = 0, 3
2971                    DO  m = 1, surf_usm_v(l)%ns
2972                       i = surf_usm_v(l)%i(m)
2973                       j = surf_usm_v(l)%j(m)
2974                       k = surf_usm_v(l)%k(m)
2975                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
2976                    ENDDO
2977                 ENDDO
2978              ENDIF
2979             
2980          CASE ( 'usm_t_surf_window' )
2981!--           surface temperature for window surfaces
2982
2983              IF ( av == 0 )  THEN
2984                 DO  m = 1, surf_usm_h%ns
2985                    i = surf_usm_h%i(m)
2986                    j = surf_usm_h%j(m)
2987                    k = surf_usm_h%k(m)
2988                    temp_pf(k,j,i) = t_surf_window_h(m)
2989                 ENDDO
2990                 DO  l = 0, 3
2991                    DO  m = 1, surf_usm_v(l)%ns
2992                       i = surf_usm_v(l)%i(m)
2993                       j = surf_usm_v(l)%j(m)
2994                       k = surf_usm_v(l)%k(m)
2995                       temp_pf(k,j,i) = t_surf_window_v(l)%t(m)
2996                    ENDDO
2997                 ENDDO
2998
2999              ELSE
3000                 DO  m = 1, surf_usm_h%ns
3001                    i = surf_usm_h%i(m)
3002                    j = surf_usm_h%j(m)
3003                    k = surf_usm_h%k(m)
3004                    temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
3005                 ENDDO
3006                 DO  l = 0, 3
3007                    DO  m = 1, surf_usm_v(l)%ns
3008                       i = surf_usm_v(l)%i(m)
3009                       j = surf_usm_v(l)%j(m)
3010                       k = surf_usm_v(l)%k(m)
3011                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_window_av(m)
3012                    ENDDO
3013
3014                 ENDDO
3015
3016              ENDIF
3017
3018          CASE ( 'usm_t_surf_green' )
3019!--           surface temperature for green surfaces
3020
3021              IF ( av == 0 )  THEN
3022                 DO  m = 1, surf_usm_h%ns
3023                    i = surf_usm_h%i(m)
3024                    j = surf_usm_h%j(m)
3025                    k = surf_usm_h%k(m)
3026                    temp_pf(k,j,i) = t_surf_green_h(m)
3027                 ENDDO
3028                 DO  l = 0, 3
3029                    DO  m = 1, surf_usm_v(l)%ns
3030                       i = surf_usm_v(l)%i(m)
3031                       j = surf_usm_v(l)%j(m)
3032                       k = surf_usm_v(l)%k(m)
3033                       temp_pf(k,j,i) = t_surf_green_v(l)%t(m)
3034                    ENDDO
3035                 ENDDO
3036
3037              ELSE
3038                 DO  m = 1, surf_usm_h%ns
3039                    i = surf_usm_h%i(m)
3040                    j = surf_usm_h%j(m)
3041                    k = surf_usm_h%k(m)
3042                    temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
3043                 ENDDO
3044                 DO  l = 0, 3
3045                    DO  m = 1, surf_usm_v(l)%ns
3046                       i = surf_usm_v(l)%i(m)
3047                       j = surf_usm_v(l)%j(m)
3048                       k = surf_usm_v(l)%k(m)
3049                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_green_av(m)
3050                    ENDDO
3051
3052                 ENDDO
3053
3054              ENDIF
3055
3056          CASE ( 'usm_t_surf_whole' )
3057!--           surface temperature for whole surfaces
3058
3059              IF ( av == 0 )  THEN
3060                 DO  m = 1, surf_usm_h%ns
3061                    i = surf_usm_h%i(m)
3062                    j = surf_usm_h%j(m)
3063                    k = surf_usm_h%k(m)
3064                    temp_pf(k,j,i) = t_surf_whole_h(m)
3065                 ENDDO
3066                 DO  l = 0, 3
3067                    DO  m = 1, surf_usm_v(l)%ns
3068                       i = surf_usm_v(l)%i(m)
3069                       j = surf_usm_v(l)%j(m)
3070                       k = surf_usm_v(l)%k(m)
3071                       temp_pf(k,j,i) = t_surf_whole_v(l)%t(m)
3072                    ENDDO
3073                 ENDDO
3074
3075              ELSE
3076                 DO  m = 1, surf_usm_h%ns
3077                    i = surf_usm_h%i(m)
3078                    j = surf_usm_h%j(m)
3079                    k = surf_usm_h%k(m)
3080                    temp_pf(k,j,i) = surf_usm_h%t_surf_whole_av(m)
3081                 ENDDO
3082                 DO  l = 0, 3
3083                    DO  m = 1, surf_usm_v(l)%ns
3084                       i = surf_usm_v(l)%i(m)
3085                       j = surf_usm_v(l)%j(m)
3086                       k = surf_usm_v(l)%k(m)
3087                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_whole_av(m)
3088                    ENDDO
3089
3090                 ENDDO
3091
3092              ENDIF
3093             
3094          CASE ( 'usm_t_surf_10cm' )
3095!--           near surface temperature for whole surfaces
3096
3097              IF ( av == 0 )  THEN
3098                 DO  m = 1, surf_usm_h%ns
3099                    i = surf_usm_h%i(m)
3100                    j = surf_usm_h%j(m)
3101                    k = surf_usm_h%k(m)
3102                    temp_pf(k,j,i) = t_surf_10cm_h(m)
3103                 ENDDO
3104                 DO  l = 0, 3
3105                    DO  m = 1, surf_usm_v(l)%ns
3106                       i = surf_usm_v(l)%i(m)
3107                       j = surf_usm_v(l)%j(m)
3108                       k = surf_usm_v(l)%k(m)
3109                       temp_pf(k,j,i) = t_surf_10cm_v(l)%t(m)
3110                    ENDDO
3111                 ENDDO
3112
3113              ELSE
3114                 DO  m = 1, surf_usm_h%ns
3115                    i = surf_usm_h%i(m)
3116                    j = surf_usm_h%j(m)
3117                    k = surf_usm_h%k(m)
3118                    temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m)
3119                 ENDDO
3120                 DO  l = 0, 3
3121                    DO  m = 1, surf_usm_v(l)%ns
3122                       i = surf_usm_v(l)%i(m)
3123                       j = surf_usm_v(l)%j(m)
3124                       k = surf_usm_v(l)%k(m)
3125                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_10cm_av(m)
3126                    ENDDO
3127
3128                 ENDDO
3129
3130              ENDIF
3131
3132             
3133          CASE ( 'usm_t_wall' )
3134!--           wall temperature for  iwl layer of walls and land
3135              IF ( av == 0 )  THEN
3136                 DO  m = 1, surf_usm_h%ns
3137                    i = surf_usm_h%i(m)
3138                    j = surf_usm_h%j(m)
3139                    k = surf_usm_h%k(m)
3140                    temp_pf(k,j,i) = t_wall_h(iwl,m)
3141                 ENDDO
3142                 DO  l = 0, 3
3143                    DO  m = 1, surf_usm_v(l)%ns
3144                       i = surf_usm_v(l)%i(m)
3145                       j = surf_usm_v(l)%j(m)
3146                       k = surf_usm_v(l)%k(m)
3147                       temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
3148                    ENDDO
3149                 ENDDO
3150              ELSE
3151                 DO  m = 1, surf_usm_h%ns
3152                    i = surf_usm_h%i(m)
3153                    j = surf_usm_h%j(m)
3154                    k = surf_usm_h%k(m)
3155                    temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
3156                 ENDDO
3157                 DO  l = 0, 3
3158                    DO  m = 1, surf_usm_v(l)%ns
3159                       i = surf_usm_v(l)%i(m)
3160                       j = surf_usm_v(l)%j(m)
3161                       k = surf_usm_v(l)%k(m)
3162                       temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
3163                    ENDDO
3164                 ENDDO
3165              ENDIF
3166             
3167          CASE ( 'usm_t_window' )
3168!--           window temperature for iwl layer of walls and land
3169              IF ( av == 0 )  THEN
3170                 DO  m = 1, surf_usm_h%ns
3171                    i = surf_usm_h%i(m)
3172                    j = surf_usm_h%j(m)
3173                    k = surf_usm_h%k(m)
3174                    temp_pf(k,j,i) = t_window_h(iwl,m)
3175                 ENDDO
3176                 DO  l = 0, 3
3177                    DO  m = 1, surf_usm_v(l)%ns
3178                       i = surf_usm_v(l)%i(m)
3179                       j = surf_usm_v(l)%j(m)
3180                       k = surf_usm_v(l)%k(m)
3181                       temp_pf(k,j,i) = t_window_v(l)%t(iwl,m)
3182                    ENDDO
3183                 ENDDO
3184              ELSE
3185                 DO  m = 1, surf_usm_h%ns
3186                    i = surf_usm_h%i(m)
3187                    j = surf_usm_h%j(m)
3188                    k = surf_usm_h%k(m)
3189                    temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
3190                 ENDDO
3191                 DO  l = 0, 3
3192                    DO  m = 1, surf_usm_v(l)%ns
3193                       i = surf_usm_v(l)%i(m)
3194                       j = surf_usm_v(l)%j(m)
3195                       k = surf_usm_v(l)%k(m)
3196                       temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m)
3197                    ENDDO
3198                 ENDDO
3199              ENDIF
3200
3201          CASE ( 'usm_t_green' )
3202!--           green temperature for  iwl layer of walls and land
3203              IF ( av == 0 )  THEN
3204                 DO  m = 1, surf_usm_h%ns
3205                    i = surf_usm_h%i(m)
3206                    j = surf_usm_h%j(m)
3207                    k = surf_usm_h%k(m)
3208                    temp_pf(k,j,i) = t_green_h(iwl,m)
3209                 ENDDO
3210                 DO  l = 0, 3
3211                    DO  m = 1, surf_usm_v(l)%ns
3212                       i = surf_usm_v(l)%i(m)
3213                       j = surf_usm_v(l)%j(m)
3214                       k = surf_usm_v(l)%k(m)
3215                       temp_pf(k,j,i) = t_green_v(l)%t(iwl,m)
3216                    ENDDO
3217                 ENDDO
3218              ELSE
3219                 DO  m = 1, surf_usm_h%ns
3220                    i = surf_usm_h%i(m)
3221                    j = surf_usm_h%j(m)
3222                    k = surf_usm_h%k(m)
3223                    temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
3224                 ENDDO
3225                 DO  l = 0, 3
3226                    DO  m = 1, surf_usm_v(l)%ns
3227                       i = surf_usm_v(l)%i(m)
3228                       j = surf_usm_v(l)%j(m)
3229                       k = surf_usm_v(l)%k(m)
3230                       temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m)
3231                    ENDDO
3232                 ENDDO
3233              ENDIF
3234
3235             
3236          CASE DEFAULT
3237              found = .FALSE.
3238             
3239        END SELECT
3240
3241!
3242!--     Rearrange dimensions for NetCDF output
3243        DO  j = nys, nyn
3244            DO  i = nxl, nxr
3245                DO  k = nzb_do, nzt_do
3246                    local_pf(i,j,k) = temp_pf(k,j,i)
3247                ENDDO
3248            ENDDO
3249        ENDDO
3250       
3251    END SUBROUTINE usm_data_output_3d
3252   
3253
3254!------------------------------------------------------------------------------!
3255!
3256! Description:
3257! ------------
3258!> Soubroutine defines appropriate grid for netcdf variables.
3259!> It is called out from subroutine netcdf.
3260!------------------------------------------------------------------------------!
3261    SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
3262   
3263        IMPLICIT NONE
3264
3265        CHARACTER (len=*), INTENT(IN)  ::  variable    !<
3266        LOGICAL, INTENT(OUT)           ::  found       !<
3267        CHARACTER (len=*), INTENT(OUT) ::  grid_x      !<
3268        CHARACTER (len=*), INTENT(OUT) ::  grid_y      !<
3269        CHARACTER (len=*), INTENT(OUT) ::  grid_z      !<
3270
3271        CHARACTER (len=varnamelength)  :: var
3272
3273        var = TRIM(variable)
3274        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.          &
3275             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.      &
3276             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR.   &
3277             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR.   &
3278             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.       &
3279             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.       &
3280             var(1:11) == 'usm_rad_hf_'  .OR.                                               &
3281             var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
3282             var(1:16) == 'usm_wghf_window_'  .OR. var(1:15) == 'usm_wghf_green_' .OR.      &
3283             var(1:10) == 'usm_iwghf_'  .OR. var(1:17) == 'usm_iwghf_window_' .OR.          &
3284             var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall'  .OR.               &
3285             var(1:17) == 'usm_t_surf_window'  .OR.  var(1:12) == 'usm_t_window'  .OR.      &
3286             var(1:16) == 'usm_t_surf_green'  .OR.  var(1:16) == 'usm_t_surf_whole' .OR.    &
3287             var(1:15) == 'usm_t_surf_10cm' .OR.                                            &
3288             var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                     & 
3289             var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
3290             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
3291             var(1:16) == 'usm_surfwintrans' )  THEN
3292
3293            found = .TRUE.
3294            grid_x = 'x'
3295            grid_y = 'y'
3296            grid_z = 'zu'
3297        ELSE
3298            found  = .FALSE.
3299            grid_x = 'none'
3300            grid_y = 'none'
3301            grid_z = 'none'
3302        ENDIF
3303
3304    END SUBROUTINE usm_define_netcdf_grid
3305   
3306
3307!------------------------------------------------------------------------------!
3308! Description:
3309! ------------
3310!> Initialization of the wall surface model
3311!------------------------------------------------------------------------------!
3312    SUBROUTINE usm_init_material_model
3313
3314        IMPLICIT NONE
3315
3316        INTEGER(iwp) ::  k, l, m            !< running indices
3317       
3318        CALL location_message( '    initialization of wall surface model', .TRUE. )
3319       
3320!--     Calculate wall grid spacings.
3321!--     Temperature is defined at the center of the wall layers,
3322!--     whereas gradients/fluxes are defined at the edges (_stag)     
3323!--     apply for all particular surface grids. First for horizontal surfaces
3324        DO  m = 1, surf_usm_h%ns
3325
3326           surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
3327           DO k = nzb_wall+1, nzt_wall
3328               surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) -                  &
3329                                         surf_usm_h%zw(k-1,m)
3330           ENDDO
3331           surf_usm_h%dz_window(nzb_wall,m) = surf_usm_h%zw_window(nzb_wall,m)
3332           DO k = nzb_wall+1, nzt_wall
3333               surf_usm_h%dz_window(k,m) = surf_usm_h%zw_window(k,m) -         &
3334                                         surf_usm_h%zw_window(k-1,m)
3335           ENDDO
3336           surf_usm_h%dz_green(nzb_wall,m) = surf_usm_h%zw_green(nzb_wall,m)
3337           DO k = nzb_wall+1, nzt_wall
3338               surf_usm_h%dz_green(k,m) = surf_usm_h%zw_green(k,m) -           &
3339                                         surf_usm_h%zw_green(k-1,m)
3340           ENDDO
3341           
3342           surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
3343
3344           DO k = nzb_wall, nzt_wall-1
3345               surf_usm_h%dz_wall_stag(k,m) = 0.5 * (                          &
3346                           surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
3347           ENDDO
3348           surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
3349           
3350           surf_usm_h%dz_window(nzt_wall+1,m) = surf_usm_h%dz_window(nzt_wall,m)
3351
3352           DO k = nzb_wall, nzt_wall-1
3353               surf_usm_h%dz_window_stag(k,m) = 0.5 * (                        &
3354                           surf_usm_h%dz_window(k+1,m) + surf_usm_h%dz_window(k,m) )
3355           ENDDO
3356           surf_usm_h%dz_window_stag(nzt_wall,m) = surf_usm_h%dz_window(nzt_wall,m)
3357
3358           surf_usm_h%dz_green(nzt_wall+1,m) = surf_usm_h%dz_green(nzt_wall,m)
3359
3360           DO k = nzb_wall, nzt_wall-1
3361               surf_usm_h%dz_green_stag(k,m) = 0.5 * (                         &
3362                           surf_usm_h%dz_green(k+1,m) + surf_usm_h%dz_green(k,m) )
3363           ENDDO
3364           surf_usm_h%dz_green_stag(nzt_wall,m) = surf_usm_h%dz_green(nzt_wall,m)
3365        ENDDO
3366        surf_usm_h%ddz_wall        = 1.0_wp / surf_usm_h%dz_wall
3367        surf_usm_h%ddz_wall_stag   = 1.0_wp / surf_usm_h%dz_wall_stag
3368        surf_usm_h%ddz_window      = 1.0_wp / surf_usm_h%dz_window
3369        surf_usm_h%ddz_window_stag = 1.0_wp / surf_usm_h%dz_window_stag
3370        surf_usm_h%ddz_green       = 1.0_wp / surf_usm_h%dz_green
3371        surf_usm_h%ddz_green_stag  = 1.0_wp / surf_usm_h%dz_green_stag
3372!       
3373!--     For vertical surfaces
3374        DO  l = 0, 3
3375           DO  m = 1, surf_usm_v(l)%ns
3376              surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
3377              DO k = nzb_wall+1, nzt_wall
3378                  surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) -         &
3379                                               surf_usm_v(l)%zw(k-1,m)
3380              ENDDO
3381              surf_usm_v(l)%dz_window(nzb_wall,m) = surf_usm_v(l)%zw_window(nzb_wall,m)
3382              DO k = nzb_wall+1, nzt_wall
3383                  surf_usm_v(l)%dz_window(k,m) = surf_usm_v(l)%zw_window(k,m) - &
3384                                               surf_usm_v(l)%zw_window(k-1,m)
3385              ENDDO
3386              surf_usm_v(l)%dz_green(nzb_wall,m) = surf_usm_v(l)%zw_green(nzb_wall,m)
3387              DO k = nzb_wall+1, nzt_wall
3388                  surf_usm_v(l)%dz_green(k,m) = surf_usm_v(l)%zw_green(k,m) - &
3389                                               surf_usm_v(l)%zw_green(k-1,m)
3390              ENDDO
3391           
3392              surf_usm_v(l)%dz_wall(nzt_wall+1,m) =                            &
3393                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
3394
3395              DO k = nzb_wall, nzt_wall-1
3396                  surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * (                    &
3397                                                surf_usm_v(l)%dz_wall(k+1,m) + &
3398                                                surf_usm_v(l)%dz_wall(k,m) )
3399              ENDDO
3400              surf_usm_v(l)%dz_wall_stag(nzt_wall,m) =                         &
3401                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
3402              surf_usm_v(l)%dz_window(nzt_wall+1,m) =                            &
3403                                              surf_usm_v(l)%dz_window(nzt_wall,m)
3404
3405              DO k = nzb_wall, nzt_wall-1
3406                  surf_usm_v(l)%dz_window_stag(k,m) = 0.5 * (                    &
3407                                                surf_usm_v(l)%dz_window(k+1,m) + &
3408                                                surf_usm_v(l)%dz_window(k,m) )
3409              ENDDO
3410              surf_usm_v(l)%dz_window_stag(nzt_wall,m) =                         &
3411                                              surf_usm_v(l)%dz_window(nzt_wall,m)
3412              surf_usm_v(l)%dz_green(nzt_wall+1,m) =                            &
3413                                              surf_usm_v(l)%dz_green(nzt_wall,m)
3414
3415              DO k = nzb_wall, nzt_wall-1
3416                  surf_usm_v(l)%dz_green_stag(k,m) = 0.5 * (                    &
3417                                                surf_usm_v(l)%dz_green(k+1,m) + &
3418                                                surf_usm_v(l)%dz_green(k,m) )
3419              ENDDO
3420              surf_usm_v(l)%dz_green_stag(nzt_wall,m) =                         &
3421                                              surf_usm_v(l)%dz_green(nzt_wall,m)
3422           ENDDO
3423           surf_usm_v(l)%ddz_wall        = 1.0_wp / surf_usm_v(l)%dz_wall
3424           surf_usm_v(l)%ddz_wall_stag   = 1.0_wp / surf_usm_v(l)%dz_wall_stag
3425           surf_usm_v(l)%ddz_window      = 1.0_wp / surf_usm_v(l)%dz_window
3426           surf_usm_v(l)%ddz_window_stag = 1.0_wp / surf_usm_v(l)%dz_window_stag
3427           surf_usm_v(l)%ddz_green       = 1.0_wp / surf_usm_v(l)%dz_green
3428           surf_usm_v(l)%ddz_green_stag  = 1.0_wp / surf_usm_v(l)%dz_green_stag
3429        ENDDO     
3430
3431       
3432        CALL location_message( '    wall structures filed out', .TRUE. )
3433
3434        CALL location_message( '    initialization of wall surface model finished', .TRUE. )
3435
3436    END SUBROUTINE usm_init_material_model
3437
3438 
3439!------------------------------------------------------------------------------!
3440! Description:
3441! ------------
3442!> Initialization of the urban surface model
3443!------------------------------------------------------------------------------!
3444    SUBROUTINE usm_init_urban_surface
3445
3446        USE arrays_3d,                                                         &
3447            ONLY:  zw
3448
3449        USE netcdf_data_input_mod,                                             &
3450            ONLY:  building_pars_f, building_type_f, terrain_height_f
3451   
3452        IMPLICIT NONE
3453
3454        INTEGER(iwp) ::  i                   !< loop index x-dirction
3455        INTEGER(iwp) ::  ind_emis_wall       !< index in input list for wall emissivity
3456        INTEGER(iwp) ::  ind_emis_green      !< index in input list for green emissivity
3457        INTEGER(iwp) ::  ind_emis_win        !< index in input list for window emissivity
3458        INTEGER(iwp) ::  ind_green_frac_w    !< index in input list for green fraction on wall
3459        INTEGER(iwp) ::  ind_green_frac_r    !< index in input list for green fraction on roof
3460        INTEGER(iwp) ::  ind_hc1             !< index in input list for heat capacity at first wall layer
3461        INTEGER(iwp) ::  ind_hc2             !< index in input list for heat capacity at second wall layer
3462        INTEGER(iwp) ::  ind_hc3             !< index in input list for heat capacity at third wall layer
3463        INTEGER(iwp) ::  ind_lai_r           !< index in input list for LAI on roof
3464        INTEGER(iwp) ::  ind_lai_w           !< index in input list for LAI on wall
3465        INTEGER(iwp) ::  ind_tc1             !< index in input list for thermal conductivity at first wall layer
3466        INTEGER(iwp) ::  ind_tc2             !< index in input list for thermal conductivity at second wall layer
3467        INTEGER(iwp) ::  ind_tc3             !< index in input list for thermal conductivity at third wall layer
3468        INTEGER(iwp) ::  ind_trans           !< index in input list for window transmissivity
3469        INTEGER(iwp) ::  ind_wall_frac       !< index in input list for wall fraction
3470        INTEGER(iwp) ::  ind_win_frac        !< index in input list for window fraction
3471        INTEGER(iwp) ::  ind_z0              !< index in input list for z0
3472        INTEGER(iwp) ::  ind_z0qh            !< index in input list for z0h / z0q
3473        INTEGER(iwp) ::  j                   !< loop index y-dirction
3474        INTEGER(iwp) ::  k                   !< loop index z-dirction
3475        INTEGER(iwp) ::  l                   !< loop index surface orientation
3476        INTEGER(iwp) ::  m                   !< loop index surface element
3477        INTEGER(iwp) ::  st                  !< dummy 
3478
3479        REAL(wp)     ::  c, d, tin, twin, exn
3480        REAL(wp)     ::  ground_floor_level_l         !< local height of ground floor level
3481        REAL(wp)     ::  z_agl                        !< height above ground
3482
3483!
3484!-- NOPOINTER version not implemented yet
3485#if defined( __nopointer )
3486    message_string = 'The urban surface module only runs with POINTER version'
3487    CALL message( 'urban_surface_mod', 'PA0452', 1, 2, 0, 6, 0 )
3488#endif
3489
3490        CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
3491!--     surface forcing have to be disabled for LSF
3492!--     in case of enabled urban surface module
3493        IF ( large_scale_forcing )  THEN
3494            lsf_surf = .FALSE.
3495        ENDIF
3496
3497!
3498!--     Flag surface elements belonging to the ground floor level. Therefore,
3499!--     use terrain height array from file, if available. This flag is later used
3500!--     to control initialization of surface attributes.
3501        surf_usm_h%ground_level = .FALSE. 
3502        DO  m = 1, surf_usm_h%ns
3503           i = surf_usm_h%i(m)
3504           j = surf_usm_h%j(m)
3505           k = surf_usm_h%k(m)
3506!
3507!--        Get local ground level. If no ground level is given in input file,
3508!--        use default value.
3509           ground_floor_level_l = ground_floor_level
3510           IF ( building_pars_f%from_file )  THEN
3511              IF ( building_pars_f%pars_xy(ind_gflh,j,i) /=                    &
3512                   building_pars_f%fill )  &
3513                 ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)         
3514           ENDIF
3515!
3516!--        Determine height of surface element above ground level
3517           IF (  terrain_height_f%from_file )  THEN
3518              z_agl = zw(k) - terrain_height_f%var(j,i)
3519           ELSE
3520              z_agl = zw(k)
3521           ENDIF
3522!
3523!--        Set flag for ground level
3524           IF ( z_agl <= ground_floor_level_l )                                &
3525              surf_usm_h%ground_level(m) = .TRUE.
3526        ENDDO
3527
3528        DO  l = 0, 3
3529           surf_usm_v(l)%ground_level = .FALSE.
3530           DO  m = 1, surf_usm_v(l)%ns
3531              i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
3532              j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
3533              k = surf_usm_v(l)%k(m)
3534!
3535!--           Get local ground level. If no ground level is given in input file,
3536!--           use default value.
3537              ground_floor_level_l = ground_floor_level
3538              IF ( building_pars_f%from_file )  THEN
3539                 IF ( building_pars_f%pars_xy(ind_gflh,j,i) /=                 &
3540                      building_pars_f%fill ) &
3541                    ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)
3542              ENDIF
3543!
3544!--           Determine height of surface element above ground level. Please
3545!--           note, height of surface element is determined with respect to
3546!--           its height of the adjoing atmospheric grid point.
3547              IF (  terrain_height_f%from_file )  THEN
3548                 z_agl = zw(k) - terrain_height_f%var(j-surf_usm_v(l)%joff,    &
3549                                                      i-surf_usm_v(l)%ioff)
3550              ELSE
3551                 z_agl = zw(k)
3552              ENDIF
3553!
3554!--           Set flag for ground level
3555              IF ( z_agl <= ground_floor_level_l )                                &
3556                 surf_usm_v(l)%ground_level(m) = .TRUE.
3557
3558           ENDDO
3559        ENDDO
3560!
3561!--     Initialize urban-type surface attribute. According to initialization in
3562!--     land-surface model, follow a 3-level approach.
3563!--     Level 1 - initialization via default attributes
3564        DO  m = 1, surf_usm_h%ns
3565!
3566!--        Now, all horizontal surfaces are roof surfaces (?)
3567           surf_usm_h%isroof_surf(m)   = .TRUE.
3568           surf_usm_h%surface_types(m) = roof_category         !< default category for root surface
3569!
3570!--        In order to distinguish between ground floor level and
3571!--        above-ground-floor level surfaces, set input indices.
3572           ind_wall_frac    = MERGE( ind_wall_frac_gfl,    ind_wall_frac_agfl,    &
3573                                     surf_usm_h%ground_level(m) )
3574           ind_win_frac     = MERGE( ind_win_frac_gfl,     ind_win_frac_agfl,     &
3575                                     surf_usm_h%ground_level(m) )
3576           ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3577                                     surf_usm_h%ground_level(m) )
3578           ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3579                                     surf_usm_h%ground_level(m) )
3580           ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
3581                                     surf_usm_h%ground_level(m) )
3582           ind_lai_w        = MERGE( ind_lai_w_gfl,        ind_lai_w_agfl,        &
3583                                     surf_usm_h%ground_level(m) )
3584           ind_hc1          = MERGE( ind_hc1_gfl,          ind_hc1_agfl,          &
3585                                     surf_usm_h%ground_level(m) )
3586           ind_hc2          = MERGE( ind_hc2_gfl,          ind_hc2_agfl,          &
3587                                     surf_usm_h%ground_level(m) )
3588           ind_hc3          = MERGE( ind_hc3_gfl,          ind_hc3_agfl,          &
3589                                     surf_usm_h%ground_level(m) )
3590           ind_tc1          = MERGE( ind_tc1_gfl,          ind_tc1_agfl,          &
3591                                     surf_usm_h%ground_level(m) )
3592           ind_tc2          = MERGE( ind_tc2_gfl,          ind_tc2_agfl,          &
3593                                     surf_usm_h%ground_level(m) )
3594           ind_tc3          = MERGE( ind_tc3_gfl,          ind_tc3_agfl,          &
3595                                     surf_usm_h%ground_level(m) )
3596           ind_emis_wall    = MERGE( ind_emis_wall_gfl,    ind_emis_wall_agfl,    &
3597                                     surf_usm_h%ground_level(m) )
3598           ind_emis_green   = MERGE( ind_emis_green_gfl,   ind_emis_green_agfl,   &
3599                                     surf_usm_h%ground_level(m) )
3600           ind_emis_win     = MERGE( ind_emis_win_gfl,     ind_emis_win_agfl,     &
3601                                     surf_usm_h%ground_level(m) )
3602           ind_trans        = MERGE( ind_trans_gfl,        ind_trans_agfl,        &
3603                                     surf_usm_h%ground_level(m) )
3604           ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
3605                                     surf_usm_h%ground_level(m) )
3606           ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
3607                                     surf_usm_h%ground_level(m) )
3608!
3609!--        Initialize relatvie wall- (0), green- (1) and window (2) fractions
3610           surf_usm_h%frac(0,m)  = building_pars(ind_wall_frac,building_type)   
3611           surf_usm_h%frac(1,m)  = building_pars(ind_green_frac_r,building_type) 
3612           surf_usm_h%frac(2,m)  = building_pars(ind_win_frac,building_type) 
3613           surf_usm_h%lai(m)     = building_pars(ind_green_frac_r,building_type) 
3614
3615           surf_usm_h%rho_c_wall(nzb_wall,m)   = building_pars(ind_hc1,building_type) 
3616           surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3617           surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3618           surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type)   
3619           surf_usm_h%lambda_h(nzb_wall,m)   = building_pars(ind_tc1,building_type) 
3620           surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 
3621           surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3622           surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type)   
3623           surf_usm_h%rho_c_green(nzb_wall,m)   = building_pars(ind_hc1,building_type) 
3624           surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3625           surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3626           surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,building_type)   
3627           surf_usm_h%lambda_h_green(nzb_wall,m)   = building_pars(ind_tc1,building_type) 
3628           surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 
3629           surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3630           surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3631           surf_usm_h%rho_c_window(nzb_wall,m)   = building_pars(ind_hc1,building_type) 
3632           surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3633           surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3634           surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,building_type)   
3635           surf_usm_h%lambda_h_window(nzb_wall,m)   = building_pars(ind_tc1,building_type) 
3636           surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 
3637           surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3638           surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,building_type)   
3639
3640           surf_usm_h%target_temp_summer(m)  = building_pars(12,building_type)   
3641           surf_usm_h%target_temp_winter(m)  = building_pars(13,building_type)   
3642!
3643!--        emissivity of wall-, green- and window fraction
3644           surf_usm_h%emissivity(0,m)        = building_pars(ind_emis_wall,building_type)
3645           surf_usm_h%emissivity(1,m)        = building_pars(ind_emis_green,building_type)
3646           surf_usm_h%emissivity(2,m)        = building_pars(ind_emis_win,building_type)
3647
3648           surf_usm_h%transmissivity(m)      = building_pars(ind_trans,building_type)
3649
3650           surf_usm_h%z0(m)                  = building_pars(ind_z0,building_type)
3651           surf_usm_h%z0h(m)                 = building_pars(ind_z0qh,building_type)
3652           surf_usm_h%z0q(m)                 = building_pars(ind_z0qh,building_type)
3653!
3654!--        albedo type for wall fraction, green fraction, window fraction
3655           surf_usm_h%albedo_type(0,m)       = INT( building_pars(ind_alb_wall,building_type)  )
3656           surf_usm_h%albedo_type(1,m)       = INT( building_pars(ind_alb_green,building_type) )
3657           surf_usm_h%albedo_type(2,m)       = INT( building_pars(ind_alb_win,building_type)   )
3658
3659           surf_usm_h%zw(nzb_wall,m)         = building_pars(ind_thick_1,building_type)
3660           surf_usm_h%zw(nzb_wall+1,m)       = building_pars(ind_thick_2,building_type)
3661           surf_usm_h%zw(nzb_wall+2,m)       = building_pars(ind_thick_3,building_type)
3662           surf_usm_h%zw(nzb_wall+3,m)       = building_pars(ind_thick_4,building_type)
3663           
3664           surf_usm_h%zw_green(nzb_wall,m)         = building_pars(ind_thick_1,building_type)
3665           surf_usm_h%zw_green(nzb_wall+1,m)       = building_pars(ind_thick_2,building_type)
3666           surf_usm_h%zw_green(nzb_wall+2,m)       = building_pars(ind_thick_3,building_type)
3667           surf_usm_h%zw_green(nzb_wall+3,m)       = building_pars(ind_thick_4,building_type)
3668           
3669           surf_usm_h%zw_window(nzb_wall,m)         = building_pars(ind_thick_1,building_type)
3670           surf_usm_h%zw_window(nzb_wall+1,m)       = building_pars(ind_thick_2,building_type)
3671           surf_usm_h%zw_window(nzb_wall+2,m)       = building_pars(ind_thick_3,building_type)
3672           surf_usm_h%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4,building_type)
3673
3674           surf_usm_h%c_surface(m)           = building_pars(45,building_type) 
3675           surf_usm_h%lambda_surf(m)         = building_pars(46,building_type) 
3676           surf_usm_h%c_surface_green(m)     = building_pars(45,building_type) 
3677           surf_usm_h%lambda_surf_green(m)   = building_pars(46,building_type) 
3678           surf_usm_h%c_surface_window(m)    = building_pars(45,building_type) 
3679           surf_usm_h%lambda_surf_window(m)  = building_pars(46,building_type) 
3680
3681        ENDDO
3682
3683        DO  l = 0, 3
3684           DO  m = 1, surf_usm_v(l)%ns
3685
3686              surf_usm_v(l)%surface_types(m) = wall_category         !< default category for root surface
3687!
3688!--           In order to distinguish between ground floor level and
3689!--           above-ground-floor level surfaces, set input indices.
3690              ind_wall_frac    = MERGE( ind_wall_frac_gfl,    ind_wall_frac_agfl,    &
3691                                        surf_usm_v(l)%ground_level(m) )
3692              ind_win_frac     = MERGE( ind_win_frac_gfl,     ind_win_frac_agfl,     &
3693                                        surf_usm_v(l)%ground_level(m) )
3694              ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3695                                        surf_usm_v(l)%ground_level(m) )
3696              ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3697                                        surf_usm_v(l)%ground_level(m) )
3698              ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
3699                                        surf_usm_v(l)%ground_level(m) )
3700              ind_lai_w        = MERGE( ind_lai_w_gfl,        ind_lai_w_agfl,        &
3701                                        surf_usm_v(l)%ground_level(m) )
3702              ind_hc1          = MERGE( ind_hc1_gfl,          ind_hc1_agfl,          &
3703                                        surf_usm_v(l)%ground_level(m) )
3704              ind_hc2          = MERGE( ind_hc2_gfl,          ind_hc2_agfl,          &
3705                                        surf_usm_v(l)%ground_level(m) )
3706              ind_hc3          = MERGE( ind_hc3_gfl,          ind_hc3_agfl,          &
3707                                        surf_usm_v(l)%ground_level(m) )
3708              ind_tc1          = MERGE( ind_tc1_gfl,          ind_tc1_agfl,          &
3709                                        surf_usm_v(l)%ground_level(m) )
3710              ind_tc2          = MERGE( ind_tc2_gfl,          ind_tc2_agfl,          &
3711                                        surf_usm_v(l)%ground_level(m) )
3712              ind_tc3          = MERGE( ind_tc3_gfl,          ind_tc3_agfl,          &
3713                                        surf_usm_v(l)%ground_level(m) )
3714              ind_emis_wall    = MERGE( ind_emis_wall_gfl,    ind_emis_wall_agfl,    &
3715                                        surf_usm_v(l)%ground_level(m) )
3716              ind_emis_green   = MERGE( ind_emis_green_gfl,   ind_emis_green_agfl,   &
3717                                        surf_usm_v(l)%ground_level(m) )
3718              ind_emis_win     = MERGE( ind_emis_win_gfl,     ind_emis_win_agfl,     &
3719                                        surf_usm_v(l)%ground_level(m) )
3720              ind_trans        = MERGE( ind_trans_gfl,       ind_trans_agfl,         &
3721                                        surf_usm_v(l)%ground_level(m) )
3722              ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
3723                                        surf_usm_v(l)%ground_level(m) )
3724              ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
3725                                        surf_usm_v(l)%ground_level(m) )
3726
3727!
3728!--           Initialize relatvie wall- (0), green- (1) and window (2) fractions
3729              surf_usm_v(l)%frac(0,m)              = building_pars(ind_wall_frac,building_type)   
3730              surf_usm_v(l)%frac(1,m)              = building_pars(ind_green_frac_w,building_type) 
3731              surf_usm_v(l)%frac(2,m)              = building_pars(ind_win_frac,building_type) 
3732              surf_usm_v(l)%lai(m)                 = building_pars(ind_lai_w,building_type) 
3733
3734              surf_usm_v(l)%rho_c_wall(nzb_wall,m)   = building_pars(ind_hc1,building_type) 
3735              surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3736              surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3737              surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type)   
3738             
3739              surf_usm_v(l)%rho_c_green(nzb_wall,m)   = building_pars(ind_hc1,building_type) 
3740              surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3741              surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3742              surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,building_type)   
3743             
3744              surf_usm_v(l)%rho_c_window(nzb_wall,m)   = building_pars(ind_hc1,building_type) 
3745              surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3746              surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3747              surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,building_type)   
3748
3749              surf_usm_v(l)%lambda_h(nzb_wall,m)   = building_pars(ind_tc1,building_type) 
3750              surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 
3751              surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3752              surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type)   
3753             
3754              surf_usm_v(l)%lambda_h_green(nzb_wall,m)   = building_pars(ind_tc1,building_type) 
3755              surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 
3756              surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3757              surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,building_type)   
3758
3759              surf_usm_v(l)%lambda_h_window(nzb_wall,m)   = building_pars(ind_tc1,building_type) 
3760              surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,building_type) 
3761              surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3762              surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,building_type)   
3763
3764              surf_usm_v(l)%target_temp_summer(m)  = building_pars(12,building_type)   
3765              surf_usm_v(l)%target_temp_winter(m)  = building_pars(13,building_type)   
3766!
3767!--           emissivity of wall-, green- and window fraction
3768              surf_usm_v(l)%emissivity(0,m)        = building_pars(ind_emis_wall,building_type)
3769              surf_usm_v(l)%emissivity(1,m)        = building_pars(ind_emis_green,building_type)
3770              surf_usm_v(l)%emissivity(2,m)        = building_pars(ind_emis_win,building_type)
3771
3772              surf_usm_v(l)%transmissivity(m)      = building_pars(ind_trans,building_type)
3773
3774              surf_usm_v(l)%z0(m)                  = building_pars(ind_z0,building_type)
3775              surf_usm_v(l)%z0h(m)                 = building_pars(ind_z0qh,building_type)
3776              surf_usm_v(l)%z0q(m)                 = building_pars(ind_z0qh,building_type)
3777
3778              surf_usm_v(l)%albedo_type(0,m)       = INT( building_pars(ind_alb_wall,building_type) )
3779              surf_usm_v(l)%albedo_type(1,m)       = INT( building_pars(ind_alb_green,building_type) )
3780              surf_usm_v(l)%albedo_type(2,m)       = INT( building_pars(ind_alb_win,building_type) )
3781
3782              surf_usm_v(l)%zw(nzb_wall,m)         = building_pars(ind_thick_1,building_type)
3783              surf_usm_v(l)%zw(nzb_wall+1,m)       = building_pars(ind_thick_2,building_type)
3784              surf_usm_v(l)%zw(nzb_wall+2,m)       = building_pars(ind_thick_3,building_type)
3785              surf_usm_v(l)%zw(nzb_wall+3,m)       = building_pars(ind_thick_4,building_type)
3786             
3787              surf_usm_v(l)%zw_green(nzb_wall,m)         = building_pars(ind_thick_1,building_type)
3788              surf_usm_v(l)%zw_green(nzb_wall+1,m)       = building_pars(ind_thick_2,building_type)
3789              surf_usm_v(l)%zw_green(nzb_wall+2,m)       = building_pars(ind_thick_3,building_type)
3790              surf_usm_v(l)%zw_green(nzb_wall+3,m)       = building_pars(ind_thick_4,building_type)
3791
3792              surf_usm_v(l)%zw_window(nzb_wall,m)         = building_pars(ind_thick_1,building_type)
3793              surf_usm_v(l)%zw_window(nzb_wall+1,m)       = building_pars(ind_thick_2,building_type)
3794              surf_usm_v(l)%zw_window(nzb_wall+2,m)       = building_pars(ind_thick_3,building_type)
3795              surf_usm_v(l)%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4,building_type)
3796
3797              surf_usm_v(l)%c_surface(m)           = building_pars(45,building_type) 
3798              surf_usm_v(l)%lambda_surf(m)         = building_pars(46,building_type)
3799              surf_usm_v(l)%c_surface_green(m)     = building_pars(45,building_type) 
3800              surf_usm_v(l)%lambda_surf_green(m)   = building_pars(46,building_type)
3801              surf_usm_v(l)%c_surface_window(m)    = building_pars(45,building_type) 
3802              surf_usm_v(l)%lambda_surf_window(m)  = building_pars(46,building_type)
3803
3804           ENDDO
3805        ENDDO
3806!
3807!--     Level 2 - initialization via building type read from file
3808        IF ( building_type_f%from_file )  THEN
3809           DO  m = 1, surf_usm_h%ns
3810              i = surf_usm_h%i(m)
3811              j = surf_usm_h%j(m)
3812!
3813!--           For the moment, limit building type to 6 (to overcome errors in input file).
3814              st = building_type_f%var(j,i)
3815              IF ( st /= building_type_f%fill )  THEN
3816
3817!
3818!--              In order to distinguish between ground floor level and
3819!--              above-ground-floor level surfaces, set input indices.
3820                 ind_wall_frac    = MERGE( ind_wall_frac_gfl,    ind_wall_frac_agfl,    &
3821                                           surf_usm_h%ground_level(m) )
3822                 ind_win_frac     = MERGE( ind_win_frac_gfl,     ind_win_frac_agfl,     &
3823                                           surf_usm_h%ground_level(m) )
3824                 ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3825                                           surf_usm_h%ground_level(m) )
3826                 ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3827                                           surf_usm_h%ground_level(m) )
3828                 ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
3829                                           surf_usm_h%ground_level(m) )
3830                 ind_lai_w        = MERGE( ind_lai_w_gfl,        ind_lai_w_agfl,        &
3831                                           surf_usm_h%ground_level(m) )
3832                 ind_hc1          = MERGE( ind_hc1_gfl,          ind_hc1_agfl,          &
3833                                           surf_usm_h%ground_level(m) )
3834                 ind_hc2          = MERGE( ind_hc2_gfl,          ind_hc2_agfl,          &
3835                                           surf_usm_h%ground_level(m) )
3836                 ind_hc3          = MERGE( ind_hc3_gfl,          ind_hc3_agfl,          &
3837                                           surf_usm_h%ground_level(m) )
3838                 ind_tc1          = MERGE( ind_tc1_gfl,          ind_tc1_agfl,          &
3839                                           surf_usm_h%ground_level(m) )
3840                 ind_tc2          = MERGE( ind_tc2_gfl,          ind_tc2_agfl,          &
3841                                           surf_usm_h%ground_level(m) )
3842                 ind_tc3          = MERGE( ind_tc3_gfl,          ind_tc3_agfl,          &
3843                                           surf_usm_h%ground_level(m) )
3844                 ind_emis_wall    = MERGE( ind_emis_wall_gfl,    ind_emis_wall_agfl,    &
3845                                           surf_usm_h%ground_level(m) )
3846                 ind_emis_green   = MERGE( ind_emis_green_gfl,   ind_emis_green_agfl,   &
3847                                           surf_usm_h%ground_level(m) )
3848                 ind_emis_win     = MERGE( ind_emis_win_gfl,     ind_emis_win_agfl,     &
3849                                           surf_usm_h%ground_level(m) )
3850                 ind_trans        = MERGE( ind_trans_gfl,        ind_trans_agfl,        &
3851                                           surf_usm_h%ground_level(m) )
3852                 ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
3853                                           surf_usm_h%ground_level(m) )
3854                 ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
3855                                           surf_usm_h%ground_level(m) )
3856
3857!
3858!--              Initialize relatvie wall- (0), green- (1) and window (2) fractions
3859                 surf_usm_h%frac(0,m)   = building_pars(ind_wall_frac,st)   
3860                 surf_usm_h%frac(1,m)   = building_pars(ind_green_frac_r,st) 
3861                 surf_usm_h%frac(2,m)   = building_pars(ind_win_frac,st) 
3862                 surf_usm_h%lai(m)      = building_pars(ind_green_frac_r,st) 
3863
3864                 surf_usm_h%rho_c_wall(nzb_wall,m)   = building_pars(ind_hc1,st) 
3865                 surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st)
3866                 surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st)
3867                 surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st)   
3868                 surf_usm_h%lambda_h(nzb_wall,m)   = building_pars(ind_tc1,st) 
3869                 surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st) 
3870                 surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st)
3871                 surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st)   
3872                 
3873                 surf_usm_h%rho_c_green(nzb_wall,m)   = building_pars(ind_hc1,st) 
3874                 surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,st)
3875                 surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,st)
3876                 surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,st)   
3877                 surf_usm_h%lambda_h_green(nzb_wall,m)   = building_pars(ind_tc1,st) 
3878                 surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,st) 
3879                 surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,st)
3880                 surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,st)   
3881               
3882                 surf_usm_h%rho_c_window(nzb_wall,m)   = building_pars(ind_hc1,st) 
3883                 surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,st)
3884                 surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,st)
3885                 surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,st)   
3886                 surf_usm_h%lambda_h_window(nzb_wall,m)   = building_pars(ind_tc1,st) 
3887                 surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,st) 
3888                 surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,st)
3889                 surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,st)   
3890
3891                 surf_usm_h%target_temp_summer(m)  = building_pars(12,st)   
3892                 surf_usm_h%target_temp_winter(m)  = building_pars(13,st)   
3893!
3894!--              emissivity of wall-, green- and window fraction
3895                 surf_usm_h%emissivity(0,m)        = building_pars(ind_emis_wall,st)
3896                 surf_usm_h%emissivity(1,m)        = building_pars(ind_emis_green,st)
3897                 surf_usm_h%emissivity(2,m)        = building_pars(ind_emis_win,st)
3898
3899                 surf_usm_h%transmissivity(m)      = building_pars(ind_trans,st)
3900
3901                 surf_usm_h%z0(m)                  = building_pars(ind_z0,st)
3902                 surf_usm_h%z0h(m)                 = building_pars(ind_z0qh,st)
3903                 surf_usm_h%z0q(m)                 = building_pars(ind_z0qh,st)
3904!
3905!--              albedo type for wall fraction, green fraction, window fraction
3906                 surf_usm_h%albedo_type(0,m)       = INT( building_pars(ind_alb_wall,st) )
3907                 surf_usm_h%albedo_type(1,m)       = INT( building_pars(ind_alb_green,st) )
3908                 surf_usm_h%albedo_type(2,m)       = INT( building_pars(ind_alb_win,st) )
3909
3910                 surf_usm_h%zw(nzb_wall,m)         = building_pars(ind_thick_1,st)
3911                 surf_usm_h%zw(nzb_wall+1,m)       = building_pars(ind_thick_2,st)
3912                 surf_usm_h%zw(nzb_wall+2,m)       = building_pars(ind_thick_3,st)
3913                 surf_usm_h%zw(nzb_wall+3,m)       = building_pars(ind_thick_4,st)
3914                 
3915                 surf_usm_h%zw_green(nzb_wall,m)         = building_pars(ind_thick_1,st)
3916                 surf_usm_h%zw_green(nzb_wall+1,m)       = building_pars(ind_thick_2,st)
3917                 surf_usm_h%zw_green(nzb_wall+2,m)       = building_pars(ind_thick_3,st)
3918                 surf_usm_h%zw_green(nzb_wall+3,m)       = building_pars(ind_thick_4,st)
3919
3920                 surf_usm_h%zw_window(nzb_wall,m)         = building_pars(ind_thick_1,st)
3921                 surf_usm_h%zw_window(nzb_wall+1,m)       = building_pars(ind_thick_2,st)
3922                 surf_usm_h%zw_window(nzb_wall+2,m)       = building_pars(ind_thick_3,st)
3923                 surf_usm_h%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4,st)
3924
3925                 surf_usm_h%c_surface(m)           = building_pars(45,st) 
3926                 surf_usm_h%lambda_surf(m)         = building_pars(46,st)
3927                 surf_usm_h%c_surface_green(m)     = building_pars(45,st) 
3928                 surf_usm_h%lambda_surf_green(m)   = building_pars(46,st)
3929                 surf_usm_h%c_surface_window(m)    = building_pars(45,st) 
3930                 surf_usm_h%lambda_surf_window(m)  = building_pars(46,st)
3931
3932              ENDIF
3933           ENDDO
3934
3935           DO  l = 0, 3
3936              DO  m = 1, surf_usm_v(l)%ns
3937                 i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
3938                 j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
3939!
3940!--              For the moment, limit building type to 6 (to overcome errors in input file).
3941
3942                 st = building_type_f%var(j,i)
3943                 IF ( st /= building_type_f%fill )  THEN
3944
3945!
3946!--                 In order to distinguish between ground floor level and
3947!--                 above-ground-floor level surfaces, set input indices.
3948                    ind_wall_frac    = MERGE( ind_wall_frac_gfl,    ind_wall_frac_agfl,    &
3949                                              surf_usm_v(l)%ground_level(m) )
3950                    ind_win_frac     = MERGE( ind_win_frac_gfl,     ind_win_frac_agfl,     &
3951                                              surf_usm_v(l)%ground_level(m) )
3952                    ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3953                                              surf_usm_v(l)%ground_level(m) )
3954                    ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3955                                              surf_usm_v(l)%ground_level(m) )
3956                    ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
3957                                              surf_usm_v(l)%ground_level(m) )
3958                    ind_lai_w        = MERGE( ind_lai_w_gfl,        ind_lai_w_agfl,        &
3959                                              surf_usm_v(l)%ground_level(m) )
3960                    ind_hc1          = MERGE( ind_hc1_gfl,          ind_hc1_agfl,          &
3961                                              surf_usm_v(l)%ground_level(m) )
3962                    ind_hc2          = MERGE( ind_hc2_gfl,          ind_hc2_agfl,          &
3963                                              surf_usm_v(l)%ground_level(m) )
3964                    ind_hc3          = MERGE( ind_hc3_gfl,          ind_hc3_agfl,          &
3965                                              surf_usm_v(l)%ground_level(m) )
3966                    ind_tc1          = MERGE( ind_tc1_gfl,          ind_tc1_agfl,          &
3967                                              surf_usm_v(l)%ground_level(m) )
3968                    ind_tc2          = MERGE( ind_tc2_gfl,          ind_tc2_agfl,          &
3969                                              surf_usm_v(l)%ground_level(m) )
3970                    ind_tc3          = MERGE( ind_tc3_gfl,          ind_tc3_agfl,          &
3971                                              surf_usm_v(l)%ground_level(m) )
3972                    ind_emis_wall    = MERGE( ind_emis_wall_gfl,    ind_emis_wall_agfl,    &
3973                                              surf_usm_v(l)%ground_level(m) )
3974                    ind_emis_green   = MERGE( ind_emis_green_gfl,   ind_emis_green_agfl,   &
3975                                              surf_usm_v(l)%ground_level(m) )
3976                    ind_emis_win     = MERGE( ind_emis_win_gfl,     ind_emis_win_agfl,     &
3977                                              surf_usm_v(l)%ground_level(m) )
3978                    ind_trans        = MERGE( ind_trans_gfl,       ind_trans_agfl,         &
3979                                        surf_usm_v(l)%ground_level(m) )
3980                    ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
3981                                              surf_usm_v(l)%ground_level(m) )
3982                    ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
3983                                              surf_usm_v(l)%ground_level(m) )
3984
3985!
3986!--                 Initialize relatvie wall- (0), green- (1) and window (2) fractions
3987                    surf_usm_v(l)%frac(0,m) = building_pars(ind_wall_frac,st)   
3988                    surf_usm_v(l)%frac(1,m) = building_pars(ind_green_frac_w,st) 
3989                    surf_usm_v(l)%frac(2,m) = building_pars(ind_win_frac,st)   
3990                    surf_usm_v(l)%lai(m)    = building_pars(ind_lai_w,st) 
3991
3992                    surf_usm_v(l)%rho_c_wall(nzb_wall,m)   = building_pars(ind_hc1,st) 
3993                    surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st)
3994                    surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st)
3995                    surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st)
3996                   
3997                    surf_usm_v(l)%rho_c_green(nzb_wall,m)   = building_pars(ind_hc1,st) 
3998                    surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,st)
3999                    surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,st)
4000                    surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,st)
4001                   
4002                    surf_usm_v(l)%rho_c_window(nzb_wall,m)   = building_pars(ind_hc1,st) 
4003                    surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,st)
4004                    surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,st)
4005                    surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,st)
4006
4007                    surf_usm_v(l)%lambda_h(nzb_wall,m)   = building_pars(ind_tc1,st) 
4008                    surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st) 
4009                    surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st)
4010                    surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st) 
4011                   
4012                    surf_usm_v(l)%lambda_h_green(nzb_wall,m)   = building_pars(ind_tc1,st) 
4013                    surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,st) 
4014                    surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,st)
4015                    surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,st) 
4016                   
4017                    surf_usm_v(l)%lambda_h_window(nzb_wall,m)   = building_pars(ind_tc1,st) 
4018                    surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,st) 
4019                    surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,st)
4020                    surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,st) 
4021
4022                    surf_usm_v(l)%target_temp_summer(m)  = building_pars(12,st)   
4023                    surf_usm_v(l)%target_temp_winter(m)  = building_pars(13,st)   
4024!
4025!--                 emissivity of wall-, green- and window fraction
4026                    surf_usm_v(l)%emissivity(0,m)        = building_pars(ind_emis_wall,st)
4027                    surf_usm_v(l)%emissivity(1,m)        = building_pars(ind_emis_green,st)
4028                    surf_usm_v(l)%emissivity(2,m)        = building_pars(ind_emis_win,st)
4029
4030                    surf_usm_v(l)%transmissivity(m)      = building_pars(ind_trans,st)
4031
4032                    surf_usm_v(l)%z0(m)                  = building_pars(ind_z0,st)
4033                    surf_usm_v(l)%z0h(m)                 = building_pars(ind_z0qh,st)
4034                    surf_usm_v(l)%z0q(m)                 = building_pars(ind_z0qh,st)
4035
4036                    surf_usm_v(l)%albedo_type(0,m)       = INT( building_pars(ind_alb_wall,st) )
4037                    surf_usm_v(l)%albedo_type(1,m)       = INT( building_pars(ind_alb_green,st) )
4038                    surf_usm_v(l)%albedo_type(2,m)       = INT( building_pars(ind_alb_win,st) )
4039
4040                    surf_usm_v(l)%zw(nzb_wall,m)         = building_pars(ind_thick_1,st)
4041                    surf_usm_v(l)%zw(nzb_wall+1,m)       = building_pars(ind_thick_2,st)
4042                    surf_usm_v(l)%zw(nzb_wall+2,m)       = building_pars(ind_thick_3,st)
4043                    surf_usm_v(l)%zw(nzb_wall+3,m)       = building_pars(ind_thick_4,st)
4044                   
4045                    surf_usm_v(l)%zw_green(nzb_wall,m)         = building_pars(ind_thick_1,st)
4046                    surf_usm_v(l)%zw_green(nzb_wall+1,m)       = building_pars(ind_thick_2,st)
4047                    surf_usm_v(l)%zw_green(nzb_wall+2,m)       = building_pars(ind_thick_3,st)
4048                    surf_usm_v(l)%zw_green(nzb_wall+3,m)       = building_pars(ind_thick_4,st)
4049                   
4050                    surf_usm_v(l)%zw_window(nzb_wall,m)         = building_pars(ind_thick_1,st)
4051                    surf_usm_v(l)%zw_window(nzb_wall+1,m)       = building_pars(ind_thick_2,st)
4052                    surf_usm_v(l)%zw_window(nzb_wall+2,m)       = building_pars(ind_thick_3,st)
4053                    surf_usm_v(l)%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4,st)
4054
4055                    surf_usm_v(l)%c_surface(m)           = building_pars(45,st) 
4056                    surf_usm_v(l)%lambda_surf(m)         = building_pars(46,st) 
4057                    surf_usm_v(l)%c_surface_green(m)     = building_pars(45,st) 
4058                    surf_usm_v(l)%lambda_surf_green(m)   = building_pars(46,st) 
4059                    surf_usm_v(l)%c_surface_window(m)    = building_pars(45,st) 
4060                    surf_usm_v(l)%lambda_surf_window(m)  = building_pars(46,st) 
4061
4062
4063                 ENDIF
4064              ENDDO
4065           ENDDO
4066        ENDIF
4067
4068!
4069!--     Level 3 - initialization via building_pars read from file
4070        IF ( building_pars_f%from_file )  THEN
4071           DO  m = 1, surf_usm_h%ns
4072              i = surf_usm_h%i(m)
4073              j = surf_usm_h%j(m)
4074
4075!
4076!--           In order to distinguish between ground floor level and
4077!--           above-ground-floor level surfaces, set input indices.
4078              ind_wall_frac    = MERGE( ind_wall_frac_gfl,    ind_wall_frac_agfl,    &
4079                                        surf_usm_h%ground_level(m) )
4080              ind_win_frac     = MERGE( ind_win_frac_gfl,     ind_win_frac_agfl,     &
4081                                        surf_usm_h%ground_level(m) )
4082              ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
4083                                        surf_usm_h%ground_level(m) )
4084              ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
4085                                        surf_usm_h%ground_level(m) )
4086              ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
4087                                        surf_usm_h%ground_level(m) )
4088              ind_lai_w        = MERGE( ind_lai_w_gfl,        ind_lai_w_agfl,        &
4089                                        surf_usm_h%ground_level(m) )
4090              ind_hc1          = MERGE( ind_hc1_gfl,          ind_hc1_agfl,          &
4091                                        surf_usm_h%ground_level(m) )
4092              ind_hc2          = MERGE( ind_hc2_gfl,          ind_hc2_agfl,          &
4093                                        surf_usm_h%ground_level(m) )
4094              ind_hc3          = MERGE( ind_hc3_gfl,          ind_hc3_agfl,          &
4095                                        surf_usm_h%ground_level(m) )
4096              ind_tc1          = MERGE( ind_tc1_gfl,          ind_tc1_agfl,          &
4097                                        surf_usm_h%ground_level(m) )
4098              ind_tc2          = MERGE( ind_tc2_gfl,          ind_tc2_agfl,          &
4099                                        surf_usm_h%ground_level(m) )
4100              ind_tc3          = MERGE( ind_tc3_gfl,          ind_tc3_agfl,          &
4101                                        surf_usm_h%ground_level(m) )
4102              ind_emis_wall    = MERGE( ind_emis_wall_gfl,    ind_emis_wall_agfl,    &
4103                                        surf_usm_h%ground_level(m) )
4104              ind_emis_green   = MERGE( ind_emis_green_gfl,   ind_emis_green_agfl,   &
4105                                        surf_usm_h%ground_level(m) )
4106              ind_emis_win     = MERGE( ind_emis_win_gfl,     ind_emis_win_agfl,     &
4107                                        surf_usm_h%ground_level(m) )
4108              ind_trans        = MERGE( ind_trans_gfl,        ind_trans_agfl,        &
4109                                        surf_usm_h%ground_level(m) )
4110              ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
4111                                        surf_usm_h%ground_level(m) )
4112              ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
4113                                        surf_usm_h%ground_level(m) )
4114
4115!
4116!--           Initialize relatvie wall- (0), green- (1) and window (2) fractions
4117              IF ( building_pars_f%pars_xy(ind_wall_frac,j,i) /= building_pars_f%fill )    &
4118                 surf_usm_h%frac(0,m)  = building_pars_f%pars_xy(ind_wall_frac,j,i)   
4119              IF ( building_pars_f%pars_xy(ind_green_frac_r,j,i) /= building_pars_f%fill ) & 
4120                 surf_usm_h%frac(1,m)  = building_pars_f%pars_xy(ind_green_frac_r,j,i) 
4121              IF ( building_pars_f%pars_xy(ind_win_frac,j,i) /= building_pars_f%fill )     & 
4122                 surf_usm_h%frac(2,m)  = building_pars_f%pars_xy(ind_win_frac,j,i)
4123
4124 
4125              IF ( building_pars_f%pars_xy(ind_lai_r,j,i) /= building_pars_f%fill )        &
4126                 surf_usm_h%lai(m)             = building_pars_f%pars_xy(ind_lai_r,j,i)
4127
4128              IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill )  THEN
4129                 surf_usm_h%rho_c_wall(nzb_wall,m)   = building_pars_f%pars_xy(ind_hc1,j,i) 
4130                 surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4131              ENDIF
4132              IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill )    &
4133                 surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i)
4134              IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill )    & 
4135                 surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i)
4136              IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill )  THEN
4137                 surf_usm_h%rho_c_green(nzb_wall,m)   = building_pars_f%pars_xy(ind_hc1,j,i) 
4138                 surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4139              ENDIF
4140              IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill )    &
4141                 surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i)
4142              IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill )    & 
4143                 surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i)
4144              IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill )  THEN
4145                 surf_usm_h%rho_c_window(nzb_wall,m)   = building_pars_f%pars_xy(ind_hc1,j,i) 
4146                 surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4147              ENDIF
4148              IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill )    &
4149                 surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i)
4150              IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill )    & 
4151                 surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i)
4152
4153              IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill )  THEN
4154                 surf_usm_h%lambda_h(nzb_wall,m)   = building_pars_f%pars_xy(ind_tc1,j,i)         
4155                 surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i)       
4156              ENDIF
4157              IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill )    &
4158                 surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i)
4159              IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill )    & 
4160                 surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i)   
4161              IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill )  THEN
4162                 surf_usm_h%lambda_h_green(nzb_wall,m)   = building_pars_f%pars_xy(ind_tc1,j,i)         
4163                 surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i)       
4164              ENDIF
4165              IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill )    &
4166                 surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i)
4167              IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill )    & 
4168                 surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i)   
4169              IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill )  THEN
4170                 surf_usm_h%lambda_h_window(nzb_wall,m)   = building_pars_f%pars_xy(ind_tc1,j,i)         
4171                 surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i)       
4172              ENDIF
4173              IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill )    &
4174                 surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i)
4175              IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill )    & 
4176                 surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i)   
4177
4178              IF ( building_pars_f%pars_xy(12,j,i) /= building_pars_f%fill )         & 
4179                 surf_usm_h%target_temp_summer(m)  = building_pars_f%pars_xy(12,j,i)   
4180              IF ( building_pars_f%pars_xy(13,j,i) /= building_pars_f%fill )         & 
4181                 surf_usm_h%target_temp_winter(m)  = building_pars_f%pars_xy(13,j,i)   
4182
4183              IF ( building_pars_f%pars_xy(ind_emis_wall,j,i) /= building_pars_f%fill ) & 
4184                 surf_usm_h%emissivity(0,m)        = building_pars_f%pars_xy(ind_emis_wall,j,i)
4185              IF ( building_pars_f%pars_xy(ind_emis_green,j,i) /= building_pars_f%fill )& 
4186                 surf_usm_h%emissivity(1,m)        = building_pars_f%pars_xy(ind_emis_green,j,i)
4187              IF ( building_pars_f%pars_xy(ind_emis_win,j,i) /= building_pars_f%fill )  & 
4188                 surf_usm_h%emissivity(2,m)        = building_pars_f%pars_xy(ind_emis_win,j,i)
4189
4190              IF ( building_pars_f%pars_xy(ind_trans,j,i) /= building_pars_f%fill )    & 
4191                 surf_usm_h%transmissivity(m)      = building_pars_f%pars_xy(ind_trans,j,i)
4192
4193              IF ( building_pars_f%pars_xy(ind_z0,j,i) /= building_pars_f%fill )    & 
4194                 surf_usm_h%z0(m)                  = building_pars_f%pars_xy(ind_z0,j,i)
4195              IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill )  & 
4196                 surf_usm_h%z0h(m)                 = building_pars_f%pars_xy(ind_z0qh,j,i)
4197              IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill )  & 
4198                 surf_usm_h%z0q(m)                 = building_pars_f%pars_xy(ind_z0qh,j,i)
4199
4200              IF ( building_pars_f%pars_xy(ind_alb_wall,j,i) /= building_pars_f%fill )    & 
4201                 surf_usm_h%albedo_type(0,m)       = building_pars_f%pars_xy(ind_alb_wall,j,i)
4202              IF ( building_pars_f%pars_xy(ind_alb_green,j,i) /= building_pars_f%fill )    & 
4203                 surf_usm_h%albedo_type(1,m)       = building_pars_f%pars_xy(ind_alb_green,j,i)
4204              IF ( building_pars_f%pars_xy(ind_alb_win,j,i) /= building_pars_f%fill )    & 
4205                 surf_usm_h%albedo_type(2,m)       = building_pars_f%pars_xy(ind_alb_win,j,i)
4206
4207              IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill )    & 
4208                 surf_usm_h%zw(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i)
4209              IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill )    & 
4210                 surf_usm_h%zw(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i)
4211              IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill )    & 
4212                 surf_usm_h%zw(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i)
4213              IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill )    & 
4214                 surf_usm_h%zw(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i)
4215              IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill )    & 
4216                 surf_usm_h%zw_green(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i)
4217              IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill )    & 
4218                 surf_usm_h%zw_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i)
4219              IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill )    & 
4220                 surf_usm_h%zw_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i)
4221              IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill )    & 
4222                 surf_usm_h%zw_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i)
4223              IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill )    & 
4224                 surf_usm_h%zw_window(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i)
4225              IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /=