Home

Context Navigation

source: palm/trunk/SOURCE/surface_mod.f90 @ 2965

Last change on this file since 2965 was 2963, checked in by suehring, 7 years ago
Minor revision of static input file checks, bugfix in initialization of surface-fractions in LSM; minor bugfix in initialization of albedo at window-surfaces; for clearer access of albedo and emissivity introduce index for vegetation/wall, pavement/green-wall and water/window surfaces
Property svn:keywords set to `Id`
File size: 210.0 KB

Line
1	!> @file 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 1997-2018 Leibniz Universitaet Hannover
18	!
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: surface_mod.f90 2963 2018-04-12 14:47:44Z scharf $
28	! Introduce index for vegetation/wall, pavement/green-wall and water/window
29	! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
30	!
31	! 2942 2018-04-03 13:51:09Z suehring
32	! Bugfix in assigning surface element data after restart
33	!
34	! 2940 2018-04-03 11:22:42Z suehring
35	! Bugfix in reading restart data of vertical surface elements
36	!
37	! 2920 2018-03-22 11:22:01Z kanani
38	! Correct comment for surface directions
39	!
40	! 2894 2018-03-15 09:17:58Z Giersch
41	! Calculations of the index range of the subdomain on file which overlaps with
42	! the current subdomain are already done in read_restart_data_mod,
43	! surface_read/write_restart_data was renamed to surface_r/wrd_local, variable
44	! named found has been introduced for checking if restart data was found,
45	! reading of restart strings has been moved completely to
46	! read_restart_data_mod, surface_rrd_local is already inside the overlap loop
47	! programmed in read_restart_data_mod, SAVE attribute added where necessary,
48	! deallocation and allocation of some arrays have been changed to take care of
49	! different restart files that can be opened (index i), the marker *** end
50	! surf *** is not necessary anymore, strings and their respective lengths are
51	! written out and read now in case of restart runs to get rid of prescribed
52	! character lengths (Giersch)
53	!
54	! 2813 2018-02-16 16:28:14Z suehring
55	! Some more bugfixes concerning restart runs
56	!
57	! 2812 2018-02-16 13:40:25Z hellstea
58	! Entries 'u_out', 'v_out' and 'w_out' removed from the functions
59	! get_topography_top_index and get_topography_top_index_ji
60	!
61	! 2805 2018-02-14 17:00:09Z suehring
62	! Bugfix in re-mapping surface elements in case of restart runs
63	!
64	! 2766 2018-01-22 17:17:47Z kanani
65	! Removed preprocessor directive __chem
66	!
67	! 2759 2018-01-17 16:24:59Z suehring
68	! Bugfix, consider density in vertical fluxes of passive scalar as well as
69	! chemical species.
70	!
71	! 2753 2018-01-16 14:16:49Z suehring
72	! +r_a_green, r_a_window
73	!
74	! 2718 2018-01-02 08:49:38Z maronga
75	! Changes from last commit documented
76	!
77	! 2706 2017-12-18 18:33:49Z suehring
78	! In case of restarts read and write pt_surface
79	!
80	! 2698 2017-12-14 18:46:24Z suehring
81	!
82	! 2696 2017-12-14 17:12:51Z kanani
83	! - Change in file header (GPL part)
84	! - Missing code for css added to surf_*, handling of surface_csflux updated (FK)
85	! - Bugfixes in reading/writing restart data in case several surface types are
86	! present at the same time (MS)
87	! - Implementation of chemistry module (FK)
88	! - Allocation of pt1 and qv1 now done for all surface types (MS)
89	! - Revised classification of surface types
90	! - Introduce offset values to simplify index determination of surface elements
91	! - Dimensions of function get_topo_top_index (MS)
92	! - added variables for USM
93	! - adapted to changes in USM (RvT)
94	!
95	! 2688 2017-12-12 17:27:04Z Giersch
96	! Allocation and initialization of the latent heat flux (qsws) at the top of
97	! the ocean domain in case of coupled runs. In addtion, a double if-query has
98	! been removed.
99	!
100	! 2638 2017-11-23 12:44:23Z raasch
101	! bugfix for constant top momentumflux
102	!
103	! 2575 2017-10-24 09:57:58Z maronga
104	! Pavement parameterization revised
105	!
106	! 2547 2017-10-16 12:41:56Z schwenkel
107	! extended by cloud_droplets option
108	!
109	! 2508 2017-10-02 08:57:09Z suehring
110	! Minor formatting adjustment
111	!
112	! 2478 2017-09-18 13:37:24Z suehring
113	! Bugfixes in initializing model top
114	!
115	! 2378 2017-08-31 13:57:27Z suehring
116	! Bugfix in write restart data
117	!
118	! 2339 2017-08-07 13:55:26Z gronemeier
119	! corrected timestamp in header
120	!
121	! 2338 2017-08-07 12:15:38Z gronemeier
122	! Modularize 1D model
123	!
124	! 2318 2017-07-20 17:27:44Z suehring
125	! New function to obtain topography top index.
126	!
127	! 2317 2017-07-20 17:27:19Z suehring
128	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
129	! includes two more prognostic equations for cloud drop concentration (nc)
130	! and cloud water content (qc).
131	!
132	! 2270 2017-06-09 12:18:47Z maronga
133	! Parameters removed/added due to changes in the LSM
134	!
135	! 2269 2017-06-09 11:57:32Z suehring
136	! Formatting and description adjustments
137	!
138	! 2256 2017-06-07 13:58:08Z suehring
139	! Enable heating at downward-facing surfaces
140	!
141	! 2233 2017-05-30 18:08:54Z suehring
142	! Initial revision
143	!
144	!
145	! Description:
146	! ------------
147	!> Surface module defines derived data structures to treat surface-
148	!> bounded grid cells. Three different types of surfaces are defined:
149	!> default surfaces, natural surfaces, and urban surfaces. The module
150	!> encompasses the allocation and initialization of surface arrays, and handles
151	!> reading and writing restart data.
152	!> In addition, a further derived data structure is defined, in order to set
153	!> boundary conditions at surfaces.
154	!> @todo For the moment, downward-facing surfaces are only classified as
155	!> default type
156	!> @todo Clean up urban-surface variables (some of them are not used any more)
157	!> @todo Revise chemistry surface flux part (reduce loops?!)
158	!------------------------------------------------------------------------------!
159	MODULE surface_mod
160
161	USE arrays_3d, &
162	ONLY: heatflux_input_conversion, momentumflux_input_conversion, &
163	rho_air, rho_air_zw, zu, zw, waterflux_input_conversion
164
165	USE chem_modules
166
167	USE control_parameters
168
169	USE indices, &
170	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, wall_flags_0
171
172	USE grid_variables, &
173	ONLY: dx, dy
174
175	USE kinds
176
177	USE model_1d_mod, &
178	ONLY: rif1d, us1d, usws1d, vsws1d
179
180
181	IMPLICIT NONE
182
183	!
184	!-- Data type used to identify grid-points where horizontal boundary conditions
185	!-- are applied
186	TYPE bc_type
187
188	INTEGER(iwp) :: ns !< number of surface elements on the PE
189
190	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid
191	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid
192	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid
193
194	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
195	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< end index within surface data type for given (j,i)
196
197	END TYPE bc_type
198	!
199	!-- Data type used to identify and treat surface-bounded grid points
200	TYPE surf_type
201
202	INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
203	INTEGER(iwp) :: joff !< offset value in y-direction, used to determine index of surface element
204	INTEGER(iwp) :: koff !< offset value in z-direction, used to determine index of surface element
205	INTEGER(iwp) :: ns !< number of surface elements on the PE
206
207	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid
208	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid
209	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid
210
211	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: facing !< Bit indicating surface orientation
212
213	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i)
214	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< End index within surface data type for given (j,i)
215
216	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_mo !< surface-layer height
217	REAL(wp), DIMENSION(:), ALLOCATABLE :: uvw_abs !< absolute surface-parallel velocity
218	REAL(wp), DIMENSION(:), ALLOCATABLE :: us !< friction velocity
219	REAL(wp), DIMENSION(:), ALLOCATABLE :: ts !< scaling parameter temerature
220	REAL(wp), DIMENSION(:), ALLOCATABLE :: qs !< scaling parameter humidity
221	REAL(wp), DIMENSION(:), ALLOCATABLE :: ss !< scaling parameter passive scalar
222	REAL(wp), DIMENSION(:), ALLOCATABLE :: qcs !< scaling parameter qc
223	REAL(wp), DIMENSION(:), ALLOCATABLE :: ncs !< scaling parameter nc
224	REAL(wp), DIMENSION(:), ALLOCATABLE :: qrs !< scaling parameter qr
225	REAL(wp), DIMENSION(:), ALLOCATABLE :: nrs !< scaling parameter nr
226
227	REAL(wp), DIMENSION(:), ALLOCATABLE :: ol !< Obukhov length
228	REAL(wp), DIMENSION(:), ALLOCATABLE :: rib !< Richardson bulk number
229
230	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0 !< roughness length for momentum
231	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0h !< roughness length for heat
232	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0q !< roughness length for humidity
233
234	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt1 !< Potential temperature at first grid level
235	REAL(wp), DIMENSION(:), ALLOCATABLE :: qv1 !< Specific humidity at first grid level
236	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: css !< scaling parameter chemical species
237	!
238	!-- Define arrays for surface fluxes
239	REAL(wp), DIMENSION(:), ALLOCATABLE :: usws !< vertical momentum flux for u-component at horizontal surfaces
240	REAL(wp), DIMENSION(:), ALLOCATABLE :: vsws !< vertical momentum flux for v-component at horizontal surfaces
241
242	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf !< surface flux sensible heat
243	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws !< surface flux latent heat
244	REAL(wp), DIMENSION(:), ALLOCATABLE :: ssws !< surface flux passive scalar
245	REAL(wp), DIMENSION(:), ALLOCATABLE :: qcsws !< surface flux qc
246	REAL(wp), DIMENSION(:), ALLOCATABLE :: ncsws !< surface flux nc
247	REAL(wp), DIMENSION(:), ALLOCATABLE :: qrsws !< surface flux qr
248	REAL(wp), DIMENSION(:), ALLOCATABLE :: nrsws !< surface flux nr
249	REAL(wp), DIMENSION(:), ALLOCATABLE :: sasws !< surface flux salinity
250
251	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: cssws !< surface flux chemical species
252	!
253	!-- Required for horizontal walls in production_e
254	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_0 !< virtual velocity component (see production_e_init for further explanation)
255	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_0 !< virtual velocity component (see production_e_init for further explanation)
256
257	REAL(wp), DIMENSION(:), ALLOCATABLE :: mom_flux_uv !< momentum flux usvs and vsus at vertical surfaces (used in diffusion_u and diffusion_v)
258	REAL(wp), DIMENSION(:), ALLOCATABLE :: mom_flux_w !< momentum flux wsus and wsvs at vertical surfaces (used in diffusion_w)
259	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mom_flux_tke !< momentum flux usvs, vsus, wsus, wsvs at vertical surfaces at grid center (used in production_e)
260	!
261	!-- Variables required for LSM as well as for USM
262	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzt_pavement !< top index for pavement in soil
263	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: albedo_type !< albedo type, for each fraction (wall,green,window or vegetation,pavement water)
264
265	LOGICAL, DIMENSION(:), ALLOCATABLE :: building_surface !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
266	LOGICAL, DIMENSION(:), ALLOCATABLE :: building_covered !< flag indicating that buildings are on top of orography, only used for vertical surfaces in LSM
267	LOGICAL, DIMENSION(:), ALLOCATABLE :: pavement_surface !< flag parameter for pavements
268	LOGICAL, DIMENSION(:), ALLOCATABLE :: water_surface !< flag parameter for water surfaces
269	LOGICAL, DIMENSION(:), ALLOCATABLE :: vegetation_surface !< flag parameter for natural land surfaces
270
271	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: albedo !< broadband albedo for each surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
272	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: emissivity !< emissivity of the surface, for each fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
273	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: frac !< relative surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
274
275	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif !< albedo for longwave diffusive radiation, solar angle of 60Â°
276	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldir !< albedo for longwave direct radiation, solar angle of 60Â°
277	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: asdif !< albedo for shortwave diffusive radiation, solar angle of 60Â°
278	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: asdir !< albedo for shortwave direct radiation, solar angle of 60Â°
279	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_aldif !< albedo for longwave diffusive radiation, solar angle of 60Â°
280	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_aldir !< albedo for longwave direct radiation, solar angle of 60Â°
281	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_asdif !< albedo for shortwave diffusive radiation, solar angle of 60Â°
282	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_asdir !< albedo for shortwave direct radiation, solar angle of 60Â°
283
284	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_surface !< skin-surface temperature
285	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net !< net radiation
286	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< net radiation, used in USM
287	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h !< heat conductivity of soil/ wall (W/m/K)
288	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h_green !< heat conductivity of green soil (W/m/K)
289	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h_window !< heat conductivity of windows (W/m/K)
290	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h_def !< default heat conductivity of soil (W/m/K)
291
292	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_lw_in !< incoming longwave radiation
293	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_lw_out !< emitted longwave radiation
294	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_sw_in !< incoming shortwave radiation
295	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_sw_out !< emitted shortwave radiation
296
297
298
299	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_liq !< liquid water coverage (of vegetated area)
300	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_veg !< vegetation coverage
301	REAL(wp), DIMENSION(:), ALLOCATABLE :: f_sw_in !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
302	REAL(wp), DIMENSION(:), ALLOCATABLE :: ghf !< ground heat flux
303	REAL(wp), DIMENSION(:), ALLOCATABLE :: g_d !< coefficient for dependence of r_canopy on water vapour pressure deficit
304	REAL(wp), DIMENSION(:), ALLOCATABLE :: lai !< leaf area index
305	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surface_u !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
306	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surface_s !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
307	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_liq !< surface flux of latent heat (liquid water portion)
308	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_soil !< surface flux of latent heat (soil portion)
309	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_veg !< surface flux of latent heat (vegetation portion)
310	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_a !< aerodynamic resistance
311	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_a_green !< aerodynamic resistance at green fraction
312	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_a_window !< aerodynamic resistance at window fraction
313	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_canopy !< canopy resistance
314	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_soil !< soil resistance
315	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_soil_min !< minimum soil resistance
316	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_s !< total surface resistance (combination of r_soil and r_canopy)
317	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_canopy_min !< minimum canopy (stomatal) resistance
318
319	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: alpha_vg !< coef. of Van Genuchten
320	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_w !< hydraulic diffusivity of soil (?)
321	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gamma_w !< hydraulic conductivity of soil (W/m/K)
322	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gamma_w_sat !< hydraulic conductivity at saturation
323	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: l_vg !< coef. of Van Genuchten
324	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_fc !< soil moisture at field capacity (m3/m3)
325	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_res !< residual soil moisture
326	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_sat !< saturation soil moisture (m3/m3)
327	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_wilt !< soil moisture at permanent wilting point (m3/m3)
328	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_vg !< coef. Van Genuchten
329	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_total_def !< default volumetric heat capacity of the (soil) layer (J/m3/K)
330	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_total !< volumetric heat capacity of the actual soil matrix (J/m3/K)
331	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: root_fr !< root fraction within the soil layers
332	!
333	!-- Urban surface variables
334	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_types !< array of types of wall parameters
335
336	LOGICAL, DIMENSION(:), ALLOCATABLE :: isroof_surf !< flag indicating roof surfaces
337	LOGICAL, DIMENSION(:), ALLOCATABLE :: ground_level !< flag indicating ground floor level surfaces
338
339	REAL(wp), DIMENSION(:), ALLOCATABLE :: target_temp_summer !< indoor target temperature summer
340	REAL(wp), DIMENSION(:), ALLOCATABLE :: target_temp_winter !< indoor target temperature summer
341
342	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface !< heat capacity of the wall surface skin (J/m2/K)
343	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface_green !< heat capacity of the green surface skin (J/m2/K)
344	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface_window !< heat capacity of the window surface skin (J/m2/K)
345	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf !< heat conductivity between air and surface (W/m2/K)
346	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf_green !< heat conductivity between air and green surface (W/m2/K)
347	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf_window !< heat conductivity between air and window surface (W/m2/K)
348	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_wall !< thickness of the wall, roof and soil layers
349	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_green !< thickness of the green wall, roof and soil layers
350	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_window !< thickness of the window wall, roof and soil layers
351	REAL(wp), DIMENSION(:), ALLOCATABLE :: transmissivity !< transmissivity of windows
352
353	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< reflected shortwave radiation for local surface in i-th reflection
354	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< reflected + emitted longwave radiation for local surface in i-th reflection
355	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf !< total radiation flux incoming to minus outgoing from local surface
356
357	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K)
358	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_window_m !< window surface temperature tendency (K)
359	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_green_m !< green surface temperature tendency (K)
360	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf !< kinematic wall heat flux of sensible heat (actually no longer needed)
361	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb !< wall heat flux of sensible heat in wall normal direction
362
363	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb !< wall ground heat flux
364	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_window !< window ground heat flux
365	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_green !< green ground heat flux
366	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb !< indoor wall ground heat flux
367	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb_window !< indoor window ground heat flux
368
369	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_lw_out_change_0
370
371	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< shortwave radiation falling to local surface including radiation from reflections
372	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< total shortwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
373	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< longwave radiation falling to local surface including radiation from reflections
374	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< total longwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
375
376
377	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_wall !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)
378	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall !< wall grid spacing (center-center)
379	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall !< 1/dz_wall
380	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall_stag !< wall grid spacing (edge-edge)
381	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall_stag !< 1/dz_wall_stag
382	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array
383	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw !< wall layer depths (m)
384	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_window !< volumetric heat capacity of the window material ( J m-3 K-1 ) (= 2.19E6)
385	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_window !< window grid spacing (center-center)
386	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_window !< 1/dz_window
387	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_window_stag !< window grid spacing (edge-edge)
388	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_window_stag !< 1/dz_window_stag
389	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_window_m !< t_window prognostic array
390	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw_window !< window layer depths (m)
391	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_green !< volumetric heat capacity of the green material ( J m-3 K-1 ) (= 2.19E6)
392	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_green !< green grid spacing (center-center)
393	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_green !< 1/dz_green
394	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_green_stag !< green grid spacing (edge-edge)
395	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_green_stag !< 1/dz_green_stag
396	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_green_m !< t_green prognostic array
397	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw_green !< green layer depths (m)
398
399
400	!-- arrays for time averages
401	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l
402	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
403	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
404	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
405	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
406	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
407	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
408	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
409	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
410	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
411	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
412	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
413	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface
414	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_av !< average of wghf_eb
415	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_window_av !< average of wghf_eb window
416	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_green_av !< average of wghf_eb window
417	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb_av !< indoor average of wghf_eb
418	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb_window_av !< indoor average of wghf_eb window
419	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb_av !< average of wshf_eb
420	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_av !< average of wall surface temperature (K)
421	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_window_av !< average of window surface temperature (K)
422	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_green_av !< average of green wall surface temperature (K)
423	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_whole_av !< average of whole wall surface temperature (K)
424	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_10cm_av !< average of the near surface temperature (K)
425
426	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_wall_av !< Average of t_wall
427	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_window_av !< Average of t_window
428	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_green_av !< Average of t_green
429
430	END TYPE surf_type
431
432	TYPE (bc_type), DIMENSION(0:1) :: bc_h !< boundary condition data type, horizontal upward- and downward facing surfaces
433
434	TYPE (surf_type), DIMENSION(0:2), TARGET :: surf_def_h !< horizontal default surfaces (Up, Down, and Top)
435	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_def_v !< vertical default surfaces (North, South, East, West)
436	TYPE (surf_type) , TARGET :: surf_lsm_h !< horizontal natural land surfaces, so far only upward-facing
437	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_lsm_v !< vertical land surfaces (North, South, East, West)
438	TYPE (surf_type) , TARGET :: surf_usm_h !< horizontal urban surfaces, so far only upward-facing
439	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_usm_v !< vertical urban surfaces (North, South, East, West)
440
441	INTEGER(iwp), PARAMETER :: ind_veg_wall = 0 !< index for vegetation / wall-surface fraction, used for access of albedo, emissivity, etc., for each surface type
442	INTEGER(iwp), PARAMETER :: ind_pav_green = 1 !< index for pavement / green-wall surface fraction, used for access of albedo, emissivity, etc., for each surface type
443	INTEGER(iwp), PARAMETER :: ind_wat_win = 2 !< index for water / window-surface fraction, used for access of albedo, emissivity, etc., for each surface type
444
445	INTEGER(iwp) :: ns_h_on_file(0:2) !< total number of horizontal surfaces with the same facing, required for writing restart data
446	INTEGER(iwp) :: ns_v_on_file(0:3) !< total number of vertical surfaces with the same facing, required for writing restart data
447
448
449	SAVE
450
451	PRIVATE
452
453	INTERFACE init_bc
454	MODULE PROCEDURE init_bc
455	END INTERFACE init_bc
456
457	INTERFACE init_surfaces
458	MODULE PROCEDURE init_surfaces
459	END INTERFACE init_surfaces
460
461	INTERFACE init_surface_arrays
462	MODULE PROCEDURE init_surface_arrays
463	END INTERFACE init_surface_arrays
464
465	INTERFACE surface_rrd_local
466	MODULE PROCEDURE surface_rrd_local
467	END INTERFACE surface_rrd_local
468
469	INTERFACE surface_wrd_local
470	MODULE PROCEDURE surface_wrd_local
471	END INTERFACE surface_wrd_local
472
473	INTERFACE surface_last_actions
474	MODULE PROCEDURE surface_last_actions
475	END INTERFACE surface_last_actions
476
477	INTERFACE surface_restore_elements
478	MODULE PROCEDURE surface_restore_elements_1d
479	MODULE PROCEDURE surface_restore_elements_2d
480	END INTERFACE surface_restore_elements
481
482	!
483	!-- Public variables
484	PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, &
485	ns_v_on_file, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, &
486	surf_usm_h, surf_usm_v, surf_type
487	!
488	!-- Public subroutines and functions
489	PUBLIC get_topography_top_index, get_topography_top_index_ji, init_bc, &
490	init_surfaces, &
491	init_surface_arrays, surface_rrd_local, &
492	surface_restore_elements, surface_wrd_local, &
493	surface_last_actions
494
495
496	CONTAINS
497
498	!------------------------------------------------------------------------------!
499	! Description:
500	! ------------
501	!> Initialize data type for setting boundary conditions at horizontal surfaces.
502	!------------------------------------------------------------------------------!
503	SUBROUTINE init_bc
504
505	IMPLICIT NONE
506
507	INTEGER(iwp) :: i !<
508	INTEGER(iwp) :: j !<
509	INTEGER(iwp) :: k !<
510
511	INTEGER(iwp), DIMENSION(0:1) :: num_h !<
512	INTEGER(iwp), DIMENSION(0:1) :: num_h_kji !<
513	INTEGER(iwp), DIMENSION(0:1) :: start_index_h !<
514
515	!
516	!-- First of all, count the number of upward- and downward-facing surfaces
517	num_h = 0
518	DO i = nxlg, nxrg
519	DO j = nysg, nyng
520	DO k = nzb+1, nzt
521	!
522	!-- Check if current gridpoint belongs to the atmosphere
523	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
524	!
525	!-- Upward-facing
526	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) &
527	num_h(0) = num_h(0) + 1
528	!
529	!-- Downward-facing
530	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
531	num_h(1) = num_h(1) + 1
532	ENDIF
533	ENDDO
534	ENDDO
535	ENDDO
536	!
537	!-- Save the number of surface elements
538	bc_h(0)%ns = num_h(0)
539	bc_h(1)%ns = num_h(1)
540	!
541	!-- ALLOCATE data type variables
542	!-- Upward facing
543	ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
544	ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
545	ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
546	ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
547	ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg) )
548	bc_h(0)%start_index = 1
549	bc_h(0)%end_index = 0
550	!
551	!-- Downward facing
552	ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
553	ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
554	ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
555	ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
556	ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg) )
557	bc_h(1)%start_index = 1
558	bc_h(1)%end_index = 0
559	!
560	!-- Store the respective indices on data type
561	num_h(0:1) = 1
562	start_index_h(0:1) = 1
563	DO i = nxlg, nxrg
564	DO j = nysg, nyng
565
566	num_h_kji(0:1) = 0
567	DO k = nzb+1, nzt
568	!
569	!-- Check if current gridpoint belongs to the atmosphere
570	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
571	!
572	!-- Upward-facing
573	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
574	bc_h(0)%i(num_h(0)) = i
575	bc_h(0)%j(num_h(0)) = j
576	bc_h(0)%k(num_h(0)) = k
577	num_h_kji(0) = num_h_kji(0) + 1
578	num_h(0) = num_h(0) + 1
579	ENDIF
580	!
581	!-- Downward-facing
582	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
583	bc_h(1)%i(num_h(1)) = i
584	bc_h(1)%j(num_h(1)) = j
585	bc_h(1)%k(num_h(1)) = k
586	num_h_kji(1) = num_h_kji(1) + 1
587	num_h(1) = num_h(1) + 1
588	ENDIF
589	ENDIF
590	ENDDO
591	bc_h(0)%start_index(j,i) = start_index_h(0)
592	bc_h(0)%end_index(j,i) = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
593	start_index_h(0) = bc_h(0)%end_index(j,i) + 1
594
595	bc_h(1)%start_index(j,i) = start_index_h(1)
596	bc_h(1)%end_index(j,i) = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
597	start_index_h(1) = bc_h(1)%end_index(j,i) + 1
598	ENDDO
599	ENDDO
600
601
602	END SUBROUTINE init_bc
603
604
605	!------------------------------------------------------------------------------!
606	! Description:
607	! ------------
608	!> Initialize horizontal and vertical surfaces. Counts the number of default-,
609	!> natural and urban surfaces and allocates memory, respectively.
610	!------------------------------------------------------------------------------!
611	SUBROUTINE init_surface_arrays
612
613	IMPLICIT NONE
614
615	INTEGER(iwp) :: i !< running index x-direction
616	INTEGER(iwp) :: j !< running index y-direction
617	INTEGER(iwp) :: k !< running index z-direction
618	INTEGER(iwp) :: l !< index variable for surface facing
619	INTEGER(iwp) :: num_lsm_h !< number of horizontally-aligned natural surfaces
620	INTEGER(iwp) :: num_usm_h !< number of horizontally-aligned urban surfaces
621
622	INTEGER(iwp), DIMENSION(0:2) :: num_def_h !< number of horizontally-aligned default surfaces
623	INTEGER(iwp), DIMENSION(0:3) :: num_def_v !< number of vertically-aligned default surfaces
624	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< number of vertically-aligned natural surfaces
625	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< number of vertically-aligned urban surfaces
626
627	LOGICAL :: building !< flag indicating building grid point
628	LOGICAL :: terrain !< flag indicating natural terrain grid point
629
630	num_def_h = 0
631	num_def_v = 0
632	num_lsm_h = 0
633	num_lsm_v = 0
634	num_usm_h = 0
635	num_usm_v = 0
636	!
637	!-- Surfaces are classified according to the input data read from static
638	!-- input file. If no input file is present, all surfaces are classified
639	!-- either as natural, urban, or default, depending on the setting of
640	!-- land_surface and urban_surface. To control this, use the control
641	!-- flag topo_no_distinct
642	!
643	!-- Count number of horizontal surfaces on local domain
644	DO i = nxl, nxr
645	DO j = nys, nyn
646	DO k = nzb+1, nzt
647	!
648	!-- Check if current gridpoint belongs to the atmosphere
649	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
650	!
651	!-- Check if grid point adjoins to any upward-facing horizontal
652	!-- surface, e.g. the Earth surface, plane roofs, or ceilings.
653
654	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
655	!
656	!-- Determine flags indicating terrain or building.
657	terrain = BTEST( wall_flags_0(k-1,j,i), 5 ) .OR. &
658	topo_no_distinct
659	building = BTEST( wall_flags_0(k-1,j,i), 6 ) .OR. &
660	topo_no_distinct
661	!
662	!-- Land-surface type
663	IF ( land_surface .AND. terrain ) THEN
664	num_lsm_h = num_lsm_h + 1
665	!
666	!-- Urban surface tpye
667	ELSEIF ( urban_surface .AND. building ) THEN
668	num_usm_h = num_usm_h + 1
669	!
670	!-- Default-surface type
671	ELSEIF ( .NOT. land_surface .AND. &
672	.NOT. urban_surface ) THEN
673
674	num_def_h(0) = num_def_h(0) + 1
675	!
676	!-- Unclassifified surface-grid point. Give error message.
677	ELSE
678	WRITE( message_string, * ) &
679	'Unclassified upward-facing ' // &
680	'surface element at '// &
681	'grid point (k,j,i) = ', k, j, i
682	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
683	ENDIF
684
685	ENDIF
686	!
687	!-- Check for top-fluxes
688	IF ( k == nzt .AND. use_top_fluxes ) THEN
689	num_def_h(2) = num_def_h(2) + 1
690	!
691	!-- Check for any other downward-facing surface. So far only for
692	!-- default surface type.
693	ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
694	num_def_h(1) = num_def_h(1) + 1
695	ENDIF
696
697	ENDIF
698	ENDDO
699	ENDDO
700	ENDDO
701	!
702	!-- Count number of vertical surfaces on local domain
703	DO i = nxl, nxr
704	DO j = nys, nyn
705	DO k = nzb+1, nzt
706	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
707	!
708	!-- Northward-facing
709	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) ) THEN
710	!
711	!-- Determine flags indicating terrain or building
712
713	terrain = BTEST( wall_flags_0(k,j-1,i), 5 ) .OR. &
714	topo_no_distinct
715	building = BTEST( wall_flags_0(k,j-1,i), 6 ) .OR. &
716	topo_no_distinct
717	IF ( land_surface .AND. terrain ) THEN
718	num_lsm_v(0) = num_lsm_v(0) + 1
719	ELSEIF ( urban_surface .AND. building ) THEN
720	num_usm_v(0) = num_usm_v(0) + 1
721	!
722	!-- Default-surface type
723	ELSEIF ( .NOT. land_surface .AND. &
724	.NOT. urban_surface ) THEN
725	num_def_v(0) = num_def_v(0) + 1
726	!
727	!-- Unclassifified surface-grid point. Give error message.
728	ELSE
729	WRITE( message_string, * ) &
730	'Unclassified northward-facing ' // &
731	'surface element at '// &
732	'grid point (k,j,i) = ', k, j, i
733	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
734
735	ENDIF
736	ENDIF
737	!
738	!-- Southward-facing
739	IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) THEN
740	!
741	!-- Determine flags indicating terrain or building
742	terrain = BTEST( wall_flags_0(k,j+1,i), 5 ) .OR. &
743	topo_no_distinct
744	building = BTEST( wall_flags_0(k,j+1,i), 6 ) .OR. &
745	topo_no_distinct
746	IF ( land_surface .AND. terrain ) THEN
747	num_lsm_v(1) = num_lsm_v(1) + 1
748	ELSEIF ( urban_surface .AND. building ) THEN
749	num_usm_v(1) = num_usm_v(1) + 1
750	!
751	!-- Default-surface type
752	ELSEIF ( .NOT. land_surface .AND. &
753	.NOT. urban_surface ) THEN
754	num_def_v(1) = num_def_v(1) + 1
755	!
756	!-- Unclassifified surface-grid point. Give error message.
757	ELSE
758	WRITE( message_string, * ) &
759	'Unclassified southward-facing ' // &
760	'surface element at '// &
761	'grid point (k,j,i) = ', k, j, i
762	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
763
764	ENDIF
765	ENDIF
766	!
767	!-- Eastward-facing
768	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) ) THEN
769	!
770	!-- Determine flags indicating terrain or building
771	terrain = BTEST( wall_flags_0(k,j,i-1), 5 ) .OR. &
772	topo_no_distinct
773	building = BTEST( wall_flags_0(k,j,i-1), 6 ) .OR. &
774	topo_no_distinct
775	IF ( land_surface .AND. terrain ) THEN
776	num_lsm_v(2) = num_lsm_v(2) + 1
777	ELSEIF ( urban_surface .AND. building ) THEN
778	num_usm_v(2) = num_usm_v(2) + 1
779	!
780	!-- Default-surface type
781	ELSEIF ( .NOT. land_surface .AND. &
782	.NOT. urban_surface ) THEN
783	num_def_v(2) = num_def_v(2) + 1
784	!
785	!-- Unclassifified surface-grid point. Give error message.
786	ELSE
787	WRITE( message_string, * ) &
788	'Unclassified eastward-facing ' // &
789	'surface element at '// &
790	'grid point (k,j,i) = ', k, j, i
791	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
792
793	ENDIF
794	ENDIF
795	!
796	!-- Westward-facing
797	IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) THEN
798	!
799	!-- Determine flags indicating terrain or building
800	terrain = BTEST( wall_flags_0(k,j,i+1), 5 ) .OR. &
801	topo_no_distinct
802	building = BTEST( wall_flags_0(k,j,i+1), 6 ) .OR. &
803	topo_no_distinct
804	IF ( land_surface .AND. terrain ) THEN
805	num_lsm_v(3) = num_lsm_v(3) + 1
806	ELSEIF ( urban_surface .AND. building ) THEN
807	num_usm_v(3) = num_usm_v(3) + 1
808	!
809	!-- Default-surface type
810	ELSEIF ( .NOT. land_surface .AND. &
811	.NOT. urban_surface ) THEN
812	num_def_v(3) = num_def_v(3) + 1
813	!
814	!-- Unclassifified surface-grid point. Give error message.
815	ELSE
816	WRITE( message_string, * ) &
817	'Unclassified westward-facing ' // &
818	'surface element at '// &
819	'grid point (k,j,i) = ', k, j, i
820	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
821
822	ENDIF
823	ENDIF
824	ENDIF
825	ENDDO
826	ENDDO
827	ENDDO
828
829	!
830	!-- Store number of surfaces per core.
831	!-- Horizontal surface, default type, upward facing
832	surf_def_h(0)%ns = num_def_h(0)
833	!
834	!-- Horizontal surface, default type, downward facing
835	surf_def_h(1)%ns = num_def_h(1)
836	!
837	!-- Horizontal surface, default type, top downward facing
838	surf_def_h(2)%ns = num_def_h(2)
839	!
840	!-- Horizontal surface, natural type, so far only upward-facing
841	surf_lsm_h%ns = num_lsm_h
842	!
843	!-- Horizontal surface, urban type, so far only upward-facing
844	surf_usm_h%ns = num_usm_h
845	!
846	!-- Vertical surface, default type, northward facing
847	surf_def_v(0)%ns = num_def_v(0)
848	!
849	!-- Vertical surface, default type, southward facing
850	surf_def_v(1)%ns = num_def_v(1)
851	!
852	!-- Vertical surface, default type, eastward facing
853	surf_def_v(2)%ns = num_def_v(2)
854	!
855	!-- Vertical surface, default type, westward facing
856	surf_def_v(3)%ns = num_def_v(3)
857	!
858	!-- Vertical surface, natural type, northward facing
859	surf_lsm_v(0)%ns = num_lsm_v(0)
860	!
861	!-- Vertical surface, natural type, southward facing
862	surf_lsm_v(1)%ns = num_lsm_v(1)
863	!
864	!-- Vertical surface, natural type, eastward facing
865	surf_lsm_v(2)%ns = num_lsm_v(2)
866	!
867	!-- Vertical surface, natural type, westward facing
868	surf_lsm_v(3)%ns = num_lsm_v(3)
869	!
870	!-- Vertical surface, urban type, northward facing
871	surf_usm_v(0)%ns = num_usm_v(0)
872	!
873	!-- Vertical surface, urban type, southward facing
874	surf_usm_v(1)%ns = num_usm_v(1)
875	!
876	!-- Vertical surface, urban type, eastward facing
877	surf_usm_v(2)%ns = num_usm_v(2)
878	!
879	!-- Vertical surface, urban type, westward facing
880	surf_usm_v(3)%ns = num_usm_v(3)
881	!
882	!-- Allocate required attributes for horizontal surfaces - default type.
883	!-- Upward-facing (l=0) and downward-facing (l=1).
884	DO l = 0, 1
885	CALL allocate_surface_attributes_h ( surf_def_h(l), nys, nyn, nxl, nxr )
886	ENDDO
887	!
888	!-- Allocate required attributes for model top
889	CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
890	!
891	!-- Allocate required attributes for horizontal surfaces - natural type.
892	CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
893	!
894	!-- Allocate required attributes for horizontal surfaces - urban type.
895	CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
896
897	!
898	!-- Allocate required attributes for vertical surfaces.
899	!-- Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2)
900	!-- and westward-facing (l=3).
901	!-- Default type.
902	DO l = 0, 3
903	CALL allocate_surface_attributes_v ( surf_def_v(l), &
904	nys, nyn, nxl, nxr )
905	ENDDO
906	!
907	!-- Natural type
908	DO l = 0, 3
909	CALL allocate_surface_attributes_v ( surf_lsm_v(l), &
910	nys, nyn, nxl, nxr )
911	ENDDO
912	!
913	!-- Urban type
914	DO l = 0, 3
915	CALL allocate_surface_attributes_v ( surf_usm_v(l), &
916	nys, nyn, nxl, nxr )
917	ENDDO
918
919
920	END SUBROUTINE init_surface_arrays
921
922
923	!------------------------------------------------------------------------------!
924	! Description:
925	! ------------
926	!> Deallocating memory for upward and downward-facing horizontal surface types,
927	!> except for top fluxes.
928	!------------------------------------------------------------------------------!
929	SUBROUTINE deallocate_surface_attributes_h( surfaces )
930
931	IMPLICIT NONE
932
933
934	TYPE(surf_type) :: surfaces !< respective surface type
935
936
937	DEALLOCATE ( surfaces%start_index )
938	DEALLOCATE ( surfaces%end_index )
939	!
940	!-- Indices to locate surface element
941	DEALLOCATE ( surfaces%i )
942	DEALLOCATE ( surfaces%j )
943	DEALLOCATE ( surfaces%k )
944	!
945	!-- Surface-layer height
946	DEALLOCATE ( surfaces%z_mo )
947	!
948	!-- Surface orientation
949	DEALLOCATE ( surfaces%facing )
950	!
951	!-- Surface-parallel wind velocity
952	DEALLOCATE ( surfaces%uvw_abs )
953	!
954	!-- Roughness
955	DEALLOCATE ( surfaces%z0 )
956	DEALLOCATE ( surfaces%z0h )
957	DEALLOCATE ( surfaces%z0q )
958	!
959	!-- Friction velocity
960	DEALLOCATE ( surfaces%us )
961	!
962	!-- Stability parameter
963	DEALLOCATE ( surfaces%ol )
964	!
965	!-- Bulk Richardson number
966	DEALLOCATE ( surfaces%rib )
967	!
968	!-- Vertical momentum fluxes of u and v
969	DEALLOCATE ( surfaces%usws )
970	DEALLOCATE ( surfaces%vsws )
971	!
972	!-- Required in production_e
973	IF ( .NOT. constant_diffusion ) THEN
974	DEALLOCATE ( surfaces%u_0 )
975	DEALLOCATE ( surfaces%v_0 )
976	ENDIF
977	!
978	!-- Characteristic temperature and surface flux of sensible heat
979	DEALLOCATE ( surfaces%ts )
980	DEALLOCATE ( surfaces%shf )
981	!
982	!-- surface temperature
983	DEALLOCATE ( surfaces%pt_surface )
984	!
985	!-- Characteristic humidity and surface flux of latent heat
986	IF ( humidity ) THEN
987	DEALLOCATE ( surfaces%qs )
988	DEALLOCATE ( surfaces%qsws )
989	ENDIF
990	!
991	!-- Characteristic scalar and surface flux of scalar
992	IF ( passive_scalar ) THEN
993	DEALLOCATE ( surfaces%ss )
994	DEALLOCATE ( surfaces%ssws )
995	ENDIF
996	!
997	!-- Scaling parameter (cs*) and surface flux of chemical species
998	IF ( air_chemistry ) THEN
999	DEALLOCATE ( surfaces%css )
1000	DEALLOCATE ( surfaces%cssws )
1001	ENDIF
1002	!
1003	!-- Arrays for storing potential temperature and
1004	!-- specific humidity at first grid level
1005	DEALLOCATE ( surfaces%pt1 )
1006	DEALLOCATE ( surfaces%qv1 )
1007	!
1008	!--
1009	IF ( cloud_physics .AND. microphysics_morrison) THEN
1010	DEALLOCATE ( surfaces%qcs )
1011	DEALLOCATE ( surfaces%ncs )
1012	DEALLOCATE ( surfaces%qcsws )
1013	DEALLOCATE ( surfaces%ncsws )
1014	ENDIF
1015	!
1016	!--
1017	IF ( cloud_physics .AND. microphysics_seifert) THEN
1018	DEALLOCATE ( surfaces%qrs )
1019	DEALLOCATE ( surfaces%nrs )
1020	DEALLOCATE ( surfaces%qrsws )
1021	DEALLOCATE ( surfaces%nrsws )
1022	ENDIF
1023	!
1024	!-- Salinity surface flux
1025	IF ( ocean ) DEALLOCATE ( surfaces%sasws )
1026
1027	END SUBROUTINE deallocate_surface_attributes_h
1028
1029
1030	!------------------------------------------------------------------------------!
1031	! Description:
1032	! ------------
1033	!> Allocating memory for upward and downward-facing horizontal surface types,
1034	!> except for top fluxes.
1035	!------------------------------------------------------------------------------!
1036	SUBROUTINE allocate_surface_attributes_h( surfaces, &
1037	nys_l, nyn_l, nxl_l, nxr_l )
1038
1039	IMPLICIT NONE
1040
1041	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1042	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1043	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1044	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1045
1046	TYPE(surf_type) :: surfaces !< respective surface type
1047
1048	!
1049	!-- Allocate arrays for start and end index of horizontal surface type
1050	!-- for each (j,i)-grid point. This is required e.g. in diffion_x, which is
1051	!-- called for each (j,i). In order to find the location where the
1052	!-- respective flux is store within the surface-type, start- and end-
1053	!-- index are stored for each (j,i). For example, each (j,i) can have
1054	!-- several entries where fluxes for horizontal surfaces might be stored,
1055	!-- e.g. for overhanging structures where several upward-facing surfaces
1056	!-- might exist for given (j,i).
1057	!-- If no surface of respective type exist at current (j,i), set indicies
1058	!-- such that loop in diffusion routines will not be entered.
1059	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1060	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1061	surfaces%start_index = 0
1062	surfaces%end_index = -1
1063	!
1064	!-- Indices to locate surface element
1065	ALLOCATE ( surfaces%i(1:surfaces%ns) )
1066	ALLOCATE ( surfaces%j(1:surfaces%ns) )
1067	ALLOCATE ( surfaces%k(1:surfaces%ns) )
1068	!
1069	!-- Surface-layer height
1070	ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
1071	!
1072	!-- Surface orientation
1073	ALLOCATE ( surfaces%facing(1:surfaces%ns) )
1074	!
1075	!-- Surface-parallel wind velocity
1076	ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
1077	!
1078	!-- Roughness
1079	ALLOCATE ( surfaces%z0(1:surfaces%ns) )
1080	ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
1081	ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
1082	!
1083	!-- Friction velocity
1084	ALLOCATE ( surfaces%us(1:surfaces%ns) )
1085	!
1086	!-- Stability parameter
1087	ALLOCATE ( surfaces%ol(1:surfaces%ns) )
1088	!
1089	!-- Bulk Richardson number
1090	ALLOCATE ( surfaces%rib(1:surfaces%ns) )
1091	!
1092	!-- Vertical momentum fluxes of u and v
1093	ALLOCATE ( surfaces%usws(1:surfaces%ns) )
1094	ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
1095	!
1096	!-- Required in production_e
1097	IF ( .NOT. constant_diffusion ) THEN
1098	ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
1099	ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
1100	ENDIF
1101	!
1102	!-- Characteristic temperature and surface flux of sensible heat
1103	ALLOCATE ( surfaces%ts(1:surfaces%ns) )
1104	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
1105	!
1106	!-- surface temperature
1107	ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
1108	!
1109	!-- Characteristic humidity and surface flux of latent heat
1110	IF ( humidity ) THEN
1111	ALLOCATE ( surfaces%qs(1:surfaces%ns) )
1112	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1113	ENDIF
1114	!
1115	!-- Characteristic scalar and surface flux of scalar
1116	IF ( passive_scalar ) THEN
1117	ALLOCATE ( surfaces%ss(1:surfaces%ns) )
1118	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1119	ENDIF
1120	!
1121	!-- Scaling parameter (cs*) and surface flux of chemical species
1122	IF ( air_chemistry ) THEN
1123	ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
1124	ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1125	ENDIF
1126	!
1127	!-- Arrays for storing potential temperature and
1128	!-- specific humidity at first grid level
1129	ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
1130	ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
1131	!
1132	!--
1133	IF ( cloud_physics .AND. microphysics_morrison) THEN
1134	ALLOCATE ( surfaces%qcs(1:surfaces%ns) )
1135	ALLOCATE ( surfaces%ncs(1:surfaces%ns) )
1136	ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1137	ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1138	ENDIF
1139	!
1140	!--
1141	IF ( cloud_physics .AND. microphysics_seifert) THEN
1142	ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
1143	ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
1144	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1145	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1146	ENDIF
1147	!
1148	!-- Salinity surface flux
1149	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
1150
1151	END SUBROUTINE allocate_surface_attributes_h
1152
1153
1154	!------------------------------------------------------------------------------!
1155	! Description:
1156	! ------------
1157	!> Deallocating memory for model-top fluxes
1158	!------------------------------------------------------------------------------!
1159	SUBROUTINE deallocate_surface_attributes_h_top( surfaces )
1160
1161	IMPLICIT NONE
1162
1163
1164	TYPE(surf_type) :: surfaces !< respective surface type
1165
1166	DEALLOCATE ( surfaces%start_index )
1167	DEALLOCATE ( surfaces%end_index )
1168	!
1169	!-- Indices to locate surface (model-top) element
1170	DEALLOCATE ( surfaces%i )
1171	DEALLOCATE ( surfaces%j )
1172	DEALLOCATE ( surfaces%k )
1173
1174	IF ( .NOT. constant_diffusion ) THEN
1175	DEALLOCATE ( surfaces%u_0 )
1176	DEALLOCATE ( surfaces%v_0 )
1177	ENDIF
1178	!
1179	!-- Vertical momentum fluxes of u and v
1180	DEALLOCATE ( surfaces%usws )
1181	DEALLOCATE ( surfaces%vsws )
1182	!
1183	!-- Sensible heat flux
1184	DEALLOCATE ( surfaces%shf )
1185	!
1186	!-- Latent heat flux
1187	IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere') THEN
1188	DEALLOCATE ( surfaces%qsws )
1189	ENDIF
1190	!
1191	!-- Scalar flux
1192	IF ( passive_scalar ) THEN
1193	DEALLOCATE ( surfaces%ssws )
1194	ENDIF
1195	!
1196	!-- Chemical species flux
1197	IF ( air_chemistry ) THEN
1198	DEALLOCATE ( surfaces%cssws )
1199	ENDIF
1200	!
1201	!--
1202	IF ( cloud_physics .AND. microphysics_morrison) THEN
1203	DEALLOCATE ( surfaces%qcsws )
1204	DEALLOCATE ( surfaces%ncsws )
1205	ENDIF
1206	!
1207	!--
1208	IF ( cloud_physics .AND. microphysics_seifert) THEN
1209	DEALLOCATE ( surfaces%qrsws )
1210	DEALLOCATE ( surfaces%nrsws )
1211	ENDIF
1212	!
1213	!-- Salinity flux
1214	IF ( ocean ) DEALLOCATE ( surfaces%sasws )
1215
1216	END SUBROUTINE deallocate_surface_attributes_h_top
1217
1218
1219	!------------------------------------------------------------------------------!
1220	! Description:
1221	! ------------
1222	!> Allocating memory for model-top fluxes
1223	!------------------------------------------------------------------------------!
1224	SUBROUTINE allocate_surface_attributes_h_top( surfaces, &
1225	nys_l, nyn_l, nxl_l, nxr_l )
1226
1227	IMPLICIT NONE
1228
1229	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1230	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1231	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1232	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1233
1234	TYPE(surf_type) :: surfaces !< respective surface type
1235
1236	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1237	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1238	surfaces%start_index = 0
1239	surfaces%end_index = -1
1240	!
1241	!-- Indices to locate surface (model-top) element
1242	ALLOCATE ( surfaces%i(1:surfaces%ns) )
1243	ALLOCATE ( surfaces%j(1:surfaces%ns) )
1244	ALLOCATE ( surfaces%k(1:surfaces%ns) )
1245
1246	IF ( .NOT. constant_diffusion ) THEN
1247	ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
1248	ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
1249	ENDIF
1250	!
1251	!-- Vertical momentum fluxes of u and v
1252	ALLOCATE ( surfaces%usws(1:surfaces%ns) )
1253	ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
1254	!
1255	!-- Sensible heat flux
1256	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
1257	!
1258	!-- Latent heat flux
1259	IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere') THEN
1260	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1261	ENDIF
1262	!
1263	!-- Scalar flux
1264	IF ( passive_scalar ) THEN
1265	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1266	ENDIF
1267	!
1268	!-- Chemical species flux
1269	IF ( air_chemistry ) THEN
1270	ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1271	ENDIF
1272	!
1273	!--
1274	IF ( cloud_physics .AND. microphysics_morrison) THEN
1275	ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1276	ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1277	ENDIF
1278	!
1279	!--
1280	IF ( cloud_physics .AND. microphysics_seifert) THEN
1281	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1282	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1283	ENDIF
1284	!
1285	!-- Salinity flux
1286	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
1287
1288	END SUBROUTINE allocate_surface_attributes_h_top
1289
1290
1291	!------------------------------------------------------------------------------!
1292	! Description:
1293	! ------------
1294	!> Deallocating memory for vertical surface types.
1295	!------------------------------------------------------------------------------!
1296	SUBROUTINE deallocate_surface_attributes_v( surfaces )
1297
1298	IMPLICIT NONE
1299
1300
1301	TYPE(surf_type) :: surfaces !< respective surface type
1302
1303	!
1304	!-- Allocate arrays for start and end index of vertical surface type
1305	!-- for each (j,i)-grid point. This is required in diffion_x, which is
1306	!-- called for each (j,i). In order to find the location where the
1307	!-- respective flux is store within the surface-type, start- and end-
1308	!-- index are stored for each (j,i). For example, each (j,i) can have
1309	!-- several entries where fluxes for vertical surfaces might be stored.
1310	!-- In the flat case, where no vertical walls exit, set indicies such
1311	!-- that loop in diffusion routines will not be entered.
1312	DEALLOCATE ( surfaces%start_index )
1313	DEALLOCATE ( surfaces%end_index )
1314	!
1315	!-- Indices to locate surface element.
1316	DEALLOCATE ( surfaces%i )
1317	DEALLOCATE ( surfaces%j )
1318	DEALLOCATE ( surfaces%k )
1319	!
1320	!-- Surface-layer height
1321	DEALLOCATE ( surfaces%z_mo )
1322	!
1323	!-- Surface orientation
1324	DEALLOCATE ( surfaces%facing )
1325	!
1326	!-- Surface parallel wind velocity
1327	DEALLOCATE ( surfaces%uvw_abs )
1328	!
1329	!-- Roughness
1330	DEALLOCATE ( surfaces%z0 )
1331	DEALLOCATE ( surfaces%z0h )
1332	DEALLOCATE ( surfaces%z0q )
1333
1334	!
1335	!-- Friction velocity
1336	DEALLOCATE ( surfaces%us )
1337	!
1338	!-- Allocate Obukhov length and bulk Richardson number. Actually, at
1339	!-- vertical surfaces these are only required for natural surfaces.
1340	!-- for natural land surfaces
1341	DEALLOCATE( surfaces%ol )
1342	DEALLOCATE( surfaces%rib )
1343	!
1344	!-- Allocate arrays for surface momentum fluxes for u and v. For u at north-
1345	!-- and south-facing surfaces, for v at east- and west-facing surfaces.
1346	DEALLOCATE ( surfaces%mom_flux_uv )
1347	!
1348	!-- Allocate array for surface momentum flux for w - wsus and wsvs
1349	DEALLOCATE ( surfaces%mom_flux_w )
1350	!
1351	!-- Allocate array for surface momentum flux for subgrid-scale tke wsus and
1352	!-- wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
1353	!-- on surface.
1354	DEALLOCATE ( surfaces%mom_flux_tke )
1355	!
1356	!-- Characteristic temperature and surface flux of sensible heat
1357	DEALLOCATE ( surfaces%ts )
1358	DEALLOCATE ( surfaces%shf )
1359	!
1360	!-- surface temperature
1361	DEALLOCATE ( surfaces%pt_surface )
1362	!
1363	!-- Characteristic humidity and surface flux of latent heat
1364	IF ( humidity ) THEN
1365	DEALLOCATE ( surfaces%qs )
1366	DEALLOCATE ( surfaces%qsws )
1367	ENDIF
1368	!
1369	!-- Characteristic scalar and surface flux of scalar
1370	IF ( passive_scalar ) THEN
1371	DEALLOCATE ( surfaces%ss )
1372	DEALLOCATE ( surfaces%ssws )
1373	ENDIF
1374	!
1375	!-- Scaling parameter (cs*) and surface flux of chemical species
1376	IF ( air_chemistry ) THEN
1377	DEALLOCATE ( surfaces%css )
1378	DEALLOCATE ( surfaces%cssws )
1379	ENDIF
1380	!
1381	!-- Arrays for storing potential temperature and
1382	!-- specific humidity at first grid level
1383	DEALLOCATE ( surfaces%pt1 )
1384	DEALLOCATE ( surfaces%qv1 )
1385
1386	IF ( cloud_physics .AND. microphysics_seifert) THEN
1387	DEALLOCATE ( surfaces%qcs )
1388	DEALLOCATE ( surfaces%ncs )
1389	DEALLOCATE ( surfaces%qcsws )
1390	DEALLOCATE ( surfaces%ncsws )
1391	ENDIF
1392
1393	IF ( cloud_physics .AND. microphysics_seifert) THEN
1394	DEALLOCATE ( surfaces%qrs )
1395	DEALLOCATE ( surfaces%nrs )
1396	DEALLOCATE ( surfaces%qrsws )
1397	DEALLOCATE ( surfaces%nrsws )
1398	ENDIF
1399	!
1400	!-- Salinity surface flux
1401	IF ( ocean ) DEALLOCATE ( surfaces%sasws )
1402
1403	END SUBROUTINE deallocate_surface_attributes_v
1404
1405
1406	!------------------------------------------------------------------------------!
1407	! Description:
1408	! ------------
1409	!> Allocating memory for vertical surface types.
1410	!------------------------------------------------------------------------------!
1411	SUBROUTINE allocate_surface_attributes_v( surfaces, &
1412	nys_l, nyn_l, nxl_l, nxr_l )
1413
1414	IMPLICIT NONE
1415
1416	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1417	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1418	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1419	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1420
1421	TYPE(surf_type) :: surfaces !< respective surface type
1422
1423	!
1424	!-- Allocate arrays for start and end index of vertical surface type
1425	!-- for each (j,i)-grid point. This is required in diffion_x, which is
1426	!-- called for each (j,i). In order to find the location where the
1427	!-- respective flux is store within the surface-type, start- and end-
1428	!-- index are stored for each (j,i). For example, each (j,i) can have
1429	!-- several entries where fluxes for vertical surfaces might be stored.
1430	!-- In the flat case, where no vertical walls exit, set indicies such
1431	!-- that loop in diffusion routines will not be entered.
1432	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1433	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1434	surfaces%start_index = 0
1435	surfaces%end_index = -1
1436	!
1437	!-- Indices to locate surface element.
1438	ALLOCATE ( surfaces%i(1:surfaces%ns) )
1439	ALLOCATE ( surfaces%j(1:surfaces%ns) )
1440	ALLOCATE ( surfaces%k(1:surfaces%ns) )
1441	!
1442	!-- Surface-layer height
1443	ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
1444	!
1445	!-- Surface orientation
1446	ALLOCATE ( surfaces%facing(1:surfaces%ns) )
1447	!
1448	!-- Surface parallel wind velocity
1449	ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
1450	!
1451	!-- Roughness
1452	ALLOCATE ( surfaces%z0(1:surfaces%ns) )
1453	ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
1454	ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
1455
1456	!
1457	!-- Friction velocity
1458	ALLOCATE ( surfaces%us(1:surfaces%ns) )
1459	!
1460	!-- Allocate Obukhov length and bulk Richardson number. Actually, at
1461	!-- vertical surfaces these are only required for natural surfaces.
1462	!-- for natural land surfaces
1463	ALLOCATE( surfaces%ol(1:surfaces%ns) )
1464	ALLOCATE( surfaces%rib(1:surfaces%ns) )
1465	!
1466	!-- Allocate arrays for surface momentum fluxes for u and v. For u at north-
1467	!-- and south-facing surfaces, for v at east- and west-facing surfaces.
1468	ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
1469	!
1470	!-- Allocate array for surface momentum flux for w - wsus and wsvs
1471	ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) )
1472	!
1473	!-- Allocate array for surface momentum flux for subgrid-scale tke wsus and
1474	!-- wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
1475	!-- on surface.
1476	ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )
1477	!
1478	!-- Characteristic temperature and surface flux of sensible heat
1479	ALLOCATE ( surfaces%ts(1:surfaces%ns) )
1480	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
1481	!
1482	!-- surface temperature
1483	ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
1484	!
1485	!-- Characteristic humidity and surface flux of latent heat
1486	IF ( humidity ) THEN
1487	ALLOCATE ( surfaces%qs(1:surfaces%ns) )
1488	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1489	ENDIF
1490	!
1491	!-- Characteristic scalar and surface flux of scalar
1492	IF ( passive_scalar ) THEN
1493	ALLOCATE ( surfaces%ss(1:surfaces%ns) )
1494	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1495	ENDIF
1496	!
1497	!-- Scaling parameter (cs*) and surface flux of chemical species
1498	IF ( air_chemistry ) THEN
1499	ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
1500	ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1501	ENDIF
1502	!
1503	!-- Arrays for storing potential temperature and
1504	!-- specific humidity at first grid level
1505	ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
1506	ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
1507
1508	IF ( cloud_physics .AND. microphysics_seifert) THEN
1509	ALLOCATE ( surfaces%qcs(1:surfaces%ns) )
1510	ALLOCATE ( surfaces%ncs(1:surfaces%ns) )
1511	ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1512	ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1513	ENDIF
1514
1515	IF ( cloud_physics .AND. microphysics_seifert) THEN
1516	ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
1517	ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
1518	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1519	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1520	ENDIF
1521	!
1522	!-- Salinity surface flux
1523	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
1524
1525	END SUBROUTINE allocate_surface_attributes_v
1526
1527
1528	!------------------------------------------------------------------------------!
1529	! Description:
1530	! ------------
1531	!> Initialize surface elements, i.e. set initial values for surface fluxes,
1532	!> friction velocity, calcuation of start/end indices, etc. .
1533	!> Please note, further initialization concerning
1534	!> special surface characteristics, e.g. soil- and vegatation type,
1535	!> building type, etc., is done in the land-surface and urban-surface module,
1536	!> respectively.
1537	!------------------------------------------------------------------------------!
1538	SUBROUTINE init_surfaces
1539
1540	IMPLICIT NONE
1541
1542	INTEGER(iwp) :: i !< running index x-direction
1543	INTEGER(iwp) :: j !< running index y-direction
1544	INTEGER(iwp) :: k !< running index z-direction
1545	INTEGER(iwp) :: l !< index variable used to distinguish surface facing
1546	INTEGER(iwp) :: m !< running index surface elements
1547
1548	INTEGER(iwp) :: start_index_lsm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal natural surfaces
1549	INTEGER(iwp) :: start_index_usm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal urban surfaces
1550
1551	INTEGER(iwp) :: num_lsm_h !< current number of horizontal surface element, natural type
1552	INTEGER(iwp) :: num_lsm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal natural surfaces
1553	INTEGER(iwp) :: num_usm_h !< current number of horizontal surface element, urban type
1554	INTEGER(iwp) :: num_usm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban surfaces
1555
1556	INTEGER(iwp), DIMENSION(0:2) :: num_def_h !< current number of horizontal surface element, default type
1557	INTEGER(iwp), DIMENSION(0:2) :: num_def_h_kji !< dummy to determing local end index in surface type for given (j,i), for horizonal default surfaces
1558	INTEGER(iwp), DIMENSION(0:2) :: start_index_def_h !< dummy to determing local start index in surface type for given (j,i), for horizontal default surfaces
1559
1560	INTEGER(iwp), DIMENSION(0:3) :: num_def_v !< current number of vertical surface element, default type
1561	INTEGER(iwp), DIMENSION(0:3) :: num_def_v_kji !< dummy to determing local end index in surface type for given (j,i), for vertical default surfaces
1562	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< current number of vertical surface element, natural type
1563	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v_kji !< dummy to determing local end index in surface type for given (j,i), for vertical natural surfaces
1564	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< current number of vertical surface element, urban type
1565	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v_kji !< dummy to determing local end index in surface type for given (j,i), for vertical urban surfaces
1566
1567	INTEGER(iwp), DIMENSION(0:3) :: start_index_def_v !< dummy to determing local start index in surface type for given (j,i), for vertical default surfaces
1568	INTEGER(iwp), DIMENSION(0:3) :: start_index_lsm_v !< dummy to determing local start index in surface type for given (j,i), for vertical natural surfaces
1569	INTEGER(iwp), DIMENSION(0:3) :: start_index_usm_v !< dummy to determing local start index in surface type for given (j,i), for vertical urban surfaces
1570
1571	LOGICAL :: building !< flag indicating building grid point
1572	LOGICAL :: terrain !< flag indicating natural terrain grid point
1573
1574	!
1575	!-- Set offset indices, i.e. index difference between surface element and
1576	!-- surface-bounded grid point.
1577	!-- Upward facing - no horizontal offsets
1578	surf_def_h(0:2)%ioff = 0
1579	surf_def_h(0:2)%joff = 0
1580
1581	surf_lsm_h%ioff = 0
1582	surf_lsm_h%joff = 0
1583
1584	surf_usm_h%ioff = 0
1585	surf_usm_h%joff = 0
1586	!
1587	!-- Upward facing vertical offsets
1588	surf_def_h(0)%koff = -1
1589	surf_lsm_h%koff = -1
1590	surf_usm_h%koff = -1
1591	!
1592	!-- Downward facing vertical offset
1593	surf_def_h(1:2)%koff = 1
1594	!
1595	!-- Vertical surfaces - no vertical offset
1596	surf_def_v(0:3)%koff = 0
1597	surf_lsm_v(0:3)%koff = 0
1598	surf_usm_v(0:3)%koff = 0
1599	!
1600	!-- North- and southward facing - no offset in x
1601	surf_def_v(0:1)%ioff = 0
1602	surf_lsm_v(0:1)%ioff = 0
1603	surf_usm_v(0:1)%ioff = 0
1604	!
1605	!-- Northward facing offset in y
1606	surf_def_v(0)%joff = -1
1607	surf_lsm_v(0)%joff = -1
1608	surf_usm_v(0)%joff = -1
1609	!
1610	!-- Southward facing offset in y
1611	surf_def_v(1)%joff = 1
1612	surf_lsm_v(1)%joff = 1
1613	surf_usm_v(1)%joff = 1
1614
1615	!
1616	!-- East- and westward facing - no offset in y
1617	surf_def_v(2:3)%joff = 0
1618	surf_lsm_v(2:3)%joff = 0
1619	surf_usm_v(2:3)%joff = 0
1620	!
1621	!-- Eastward facing offset in x
1622	surf_def_v(2)%ioff = -1
1623	surf_lsm_v(2)%ioff = -1
1624	surf_usm_v(2)%ioff = -1
1625	!
1626	!-- Westward facing offset in y
1627	surf_def_v(3)%ioff = 1
1628	surf_lsm_v(3)%ioff = 1
1629	surf_usm_v(3)%ioff = 1
1630
1631	!
1632	!-- Initialize surface attributes, store indicies, surfaces orientation, etc.,
1633	num_def_h(0:2) = 1
1634	num_def_v(0:3) = 1
1635
1636	num_lsm_h = 1
1637	num_lsm_v(0:3) = 1
1638
1639	num_usm_h = 1
1640	num_usm_v(0:3) = 1
1641
1642	start_index_def_h(0:2) = 1
1643	start_index_def_v(0:3) = 1
1644
1645	start_index_lsm_h = 1
1646	start_index_lsm_v(0:3) = 1
1647
1648	start_index_usm_h = 1
1649	start_index_usm_v(0:3) = 1
1650
1651	DO i = nxl, nxr
1652	DO j = nys, nyn
1653
1654	num_def_h_kji = 0
1655	num_def_v_kji = 0
1656	num_lsm_h_kji = 0
1657	num_lsm_v_kji = 0
1658	num_usm_h_kji = 0
1659	num_usm_v_kji = 0
1660
1661	DO k = nzb+1, nzt
1662	!
1663	!-- Check if current gridpoint belongs to the atmosphere
1664	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1665	!
1666	!-- Upward-facing surface. Distinguish between differet surface types.
1667	!-- To do, think about method to flag natural and non-natural
1668	!-- surfaces.
1669	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
1670	!
1671	!-- Determine flags indicating terrain or building
1672	terrain = BTEST( wall_flags_0(k-1,j,i), 5 ) .OR. &
1673	topo_no_distinct
1674	building = BTEST( wall_flags_0(k-1,j,i), 6 ) .OR. &
1675	topo_no_distinct
1676	!
1677	!-- Natural surface type
1678	IF ( land_surface .AND. terrain ) THEN
1679	CALL initialize_horizontal_surfaces( k, j, i, &
1680	surf_lsm_h, &
1681	num_lsm_h, &
1682	num_lsm_h_kji, &
1683	.TRUE., .FALSE. )
1684	!
1685	!-- Urban surface tpye
1686	ELSEIF ( urban_surface .AND. building ) THEN
1687	CALL initialize_horizontal_surfaces( k, j, i, &
1688	surf_usm_h, &
1689	num_usm_h, &
1690	num_usm_h_kji, &
1691	.TRUE., .FALSE. )
1692	!
1693	!-- Default surface type
1694	ELSE
1695	CALL initialize_horizontal_surfaces( k, j, i, &
1696	surf_def_h(0), &
1697	num_def_h(0), &
1698	num_def_h_kji(0),&
1699	.TRUE., .FALSE. )
1700	ENDIF
1701	ENDIF
1702	!
1703	!-- downward-facing surface, first, model top. Please note,
1704	!-- for the moment, downward-facing surfaces are always of
1705	!-- default type
1706	IF ( k == nzt .AND. use_top_fluxes ) THEN
1707	CALL initialize_top( k, j, i, surf_def_h(2), &
1708	num_def_h(2), num_def_h_kji(2) )
1709	!
1710	!-- Check for any other downward-facing surface. So far only for
1711	!-- default surface type.
1712	ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
1713	CALL initialize_horizontal_surfaces( k, j, i, &
1714	surf_def_h(1), &
1715	num_def_h(1), &
1716	num_def_h_kji(1), &
1717	.FALSE., .TRUE. )
1718	ENDIF
1719	!
1720	!-- Check for vertical walls and, if required, initialize it.
1721	! Start with northward-facing surface.
1722	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) ) THEN
1723	!
1724	!-- Determine flags indicating terrain or building
1725	terrain = BTEST( wall_flags_0(k,j-1,i), 5 ) .OR. &
1726	topo_no_distinct
1727	building = BTEST( wall_flags_0(k,j-1,i), 6 ) .OR. &
1728	topo_no_distinct
1729	IF ( urban_surface .AND. building ) THEN
1730	CALL initialize_vertical_surfaces( 0, k, j, i, &
1731	surf_usm_v(0), &
1732	num_usm_v(0), &
1733	num_usm_v_kji(0), &
1734	.FALSE., .FALSE., &
1735	.FALSE., .TRUE. )
1736	ELSEIF ( land_surface .AND. terrain ) THEN
1737	CALL initialize_vertical_surfaces( 0, k, j, i, &
1738	surf_lsm_v(0), &
1739	num_lsm_v(0), &
1740	num_lsm_v_kji(0), &
1741	.FALSE., .FALSE., &
1742	.FALSE., .TRUE. )
1743	ELSE
1744	CALL initialize_vertical_surfaces( 0, k, j, i, &
1745	surf_def_v(0), &
1746	num_def_v(0), &
1747	num_def_v_kji(0), &
1748	.FALSE., .FALSE., &
1749	.FALSE., .TRUE. )
1750	ENDIF
1751	ENDIF
1752	!
1753	!-- southward-facing surface
1754	IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) THEN
1755	!
1756	!-- Determine flags indicating terrain or building
1757	terrain = BTEST( wall_flags_0(k,j+1,i), 5 ) .OR. &
1758	topo_no_distinct
1759	building = BTEST( wall_flags_0(k,j+1,i), 6 ) .OR. &
1760	topo_no_distinct
1761	IF ( urban_surface .AND. building ) THEN
1762	CALL initialize_vertical_surfaces( 1, k, j, i, &
1763	surf_usm_v(1), &
1764	num_usm_v(1), &
1765	num_usm_v_kji(1), &
1766	.FALSE., .FALSE., &
1767	.TRUE., .FALSE. )
1768	ELSEIF ( land_surface .AND. terrain ) THEN
1769	CALL initialize_vertical_surfaces( 1, k, j, i, &
1770	surf_lsm_v(1), &
1771	num_lsm_v(1), &
1772	num_lsm_v_kji(1), &
1773	.FALSE., .FALSE., &
1774	.TRUE., .FALSE. )
1775	ELSE
1776	CALL initialize_vertical_surfaces( 1, k, j, i, &
1777	surf_def_v(1), &
1778	num_def_v(1), &
1779	num_def_v_kji(1), &
1780	.FALSE., .FALSE., &
1781	.TRUE., .FALSE. )
1782	ENDIF
1783	ENDIF
1784	!
1785	!-- eastward-facing surface
1786	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) ) THEN
1787	!
1788	!-- Determine flags indicating terrain or building
1789	terrain = BTEST( wall_flags_0(k,j,i-1), 5 ) .OR. &
1790	topo_no_distinct
1791	building = BTEST( wall_flags_0(k,j,i-1), 6 ) .OR. &
1792	topo_no_distinct
1793	IF ( urban_surface .AND. building ) THEN
1794	CALL initialize_vertical_surfaces( 2, k, j, i, &
1795	surf_usm_v(2), &
1796	num_usm_v(2), &
1797	num_usm_v_kji(2), &
1798	.TRUE., .FALSE., &
1799	.FALSE., .FALSE. )
1800	ELSEIF ( land_surface .AND. terrain ) THEN
1801	CALL initialize_vertical_surfaces( 2, k, j, i, &
1802	surf_lsm_v(2), &
1803	num_lsm_v(2), &
1804	num_lsm_v_kji(2), &
1805	.TRUE., .FALSE., &
1806	.FALSE., .FALSE. )
1807	ELSE
1808	CALL initialize_vertical_surfaces( 2, k, j, i, &
1809	surf_def_v(2), &
1810	num_def_v(2), &
1811	num_def_v_kji(2), &
1812	.TRUE., .FALSE., &
1813	.FALSE., .FALSE. )
1814	ENDIF
1815	ENDIF
1816	!
1817	!-- westward-facing surface
1818	IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) THEN
1819	!
1820	!-- Determine flags indicating terrain or building
1821	terrain = BTEST( wall_flags_0(k,j,i+1), 5 ) .OR. &
1822	topo_no_distinct
1823	building = BTEST( wall_flags_0(k,j,i+1), 6 ) .OR. &
1824	topo_no_distinct
1825	IF ( urban_surface .AND. building ) THEN
1826	CALL initialize_vertical_surfaces( 3, k, j, i, &
1827	surf_usm_v(3), &
1828	num_usm_v(3), &
1829	num_usm_v_kji(3), &
1830	.FALSE., .TRUE., &
1831	.FALSE., .FALSE. )
1832	ELSEIF ( land_surface .AND. terrain ) THEN
1833	CALL initialize_vertical_surfaces( 3, k, j, i, &
1834	surf_lsm_v(3), &
1835	num_lsm_v(3), &
1836	num_lsm_v_kji(3), &
1837	.FALSE., .TRUE., &
1838	.FALSE., .FALSE. )
1839	ELSE
1840	CALL initialize_vertical_surfaces( 3, k, j, i, &
1841	surf_def_v(3), &
1842	num_def_v(3), &
1843	num_def_v_kji(3), &
1844	.FALSE., .TRUE., &
1845	.FALSE., .FALSE. )
1846	ENDIF
1847	ENDIF
1848	ENDIF
1849
1850
1851	ENDDO
1852	!
1853	!-- Determine start- and end-index at grid point (j,i). Also, for
1854	!-- horizontal surfaces more than 1 horizontal surface element can
1855	!-- exist at grid point (j,i) if overhanging structures are present.
1856	!-- Upward-facing surfaces
1857	surf_def_h(0)%start_index(j,i) = start_index_def_h(0)
1858	surf_def_h(0)%end_index(j,i) = surf_def_h(0)%start_index(j,i) + &
1859	num_def_h_kji(0) - 1
1860	start_index_def_h(0) = surf_def_h(0)%end_index(j,i) + 1
1861	!
1862	!-- Downward-facing surfaces, except model top
1863	surf_def_h(1)%start_index(j,i) = start_index_def_h(1)
1864	surf_def_h(1)%end_index(j,i) = surf_def_h(1)%start_index(j,i) + &
1865	num_def_h_kji(1) - 1
1866	start_index_def_h(1) = surf_def_h(1)%end_index(j,i) + 1
1867	!
1868	!-- Downward-facing surfaces -- model top fluxes
1869	surf_def_h(2)%start_index(j,i) = start_index_def_h(2)
1870	surf_def_h(2)%end_index(j,i) = surf_def_h(2)%start_index(j,i) + &
1871	num_def_h_kji(2) - 1
1872	start_index_def_h(2) = surf_def_h(2)%end_index(j,i) + 1
1873	!
1874	!-- Horizontal natural land surfaces
1875	surf_lsm_h%start_index(j,i) = start_index_lsm_h
1876	surf_lsm_h%end_index(j,i) = surf_lsm_h%start_index(j,i) + &
1877	num_lsm_h_kji - 1
1878	start_index_lsm_h = surf_lsm_h%end_index(j,i) + 1
1879	!
1880	!-- Horizontal urban surfaces
1881	surf_usm_h%start_index(j,i) = start_index_usm_h
1882	surf_usm_h%end_index(j,i) = surf_usm_h%start_index(j,i) + &
1883	num_usm_h_kji - 1
1884	start_index_usm_h = surf_usm_h%end_index(j,i) + 1
1885
1886	!
1887	!-- Vertical surfaces - Default type
1888	surf_def_v(0)%start_index(j,i) = start_index_def_v(0)
1889	surf_def_v(1)%start_index(j,i) = start_index_def_v(1)
1890	surf_def_v(2)%start_index(j,i) = start_index_def_v(2)
1891	surf_def_v(3)%start_index(j,i) = start_index_def_v(3)
1892	surf_def_v(0)%end_index(j,i) = start_index_def_v(0) + &
1893	num_def_v_kji(0) - 1
1894	surf_def_v(1)%end_index(j,i) = start_index_def_v(1) + &
1895	num_def_v_kji(1) - 1
1896	surf_def_v(2)%end_index(j,i) = start_index_def_v(2) + &
1897	num_def_v_kji(2) - 1
1898	surf_def_v(3)%end_index(j,i) = start_index_def_v(3) + &
1899	num_def_v_kji(3) - 1
1900	start_index_def_v(0) = surf_def_v(0)%end_index(j,i) + 1
1901	start_index_def_v(1) = surf_def_v(1)%end_index(j,i) + 1
1902	start_index_def_v(2) = surf_def_v(2)%end_index(j,i) + 1
1903	start_index_def_v(3) = surf_def_v(3)%end_index(j,i) + 1
1904	!
1905	!-- Natural type
1906	surf_lsm_v(0)%start_index(j,i) = start_index_lsm_v(0)
1907	surf_lsm_v(1)%start_index(j,i) = start_index_lsm_v(1)
1908	surf_lsm_v(2)%start_index(j,i) = start_index_lsm_v(2)
1909	surf_lsm_v(3)%start_index(j,i) = start_index_lsm_v(3)
1910	surf_lsm_v(0)%end_index(j,i) = start_index_lsm_v(0) + &
1911	num_lsm_v_kji(0) - 1
1912	surf_lsm_v(1)%end_index(j,i) = start_index_lsm_v(1) + &
1913	num_lsm_v_kji(1) - 1
1914	surf_lsm_v(2)%end_index(j,i) = start_index_lsm_v(2) + &
1915	num_lsm_v_kji(2) - 1
1916	surf_lsm_v(3)%end_index(j,i) = start_index_lsm_v(3) + &
1917	num_lsm_v_kji(3) - 1
1918	start_index_lsm_v(0) = surf_lsm_v(0)%end_index(j,i) + 1
1919	start_index_lsm_v(1) = surf_lsm_v(1)%end_index(j,i) + 1
1920	start_index_lsm_v(2) = surf_lsm_v(2)%end_index(j,i) + 1
1921	start_index_lsm_v(3) = surf_lsm_v(3)%end_index(j,i) + 1
1922	!
1923	!-- Urban type
1924	surf_usm_v(0)%start_index(j,i) = start_index_usm_v(0)
1925	surf_usm_v(1)%start_index(j,i) = start_index_usm_v(1)
1926	surf_usm_v(2)%start_index(j,i) = start_index_usm_v(2)
1927	surf_usm_v(3)%start_index(j,i) = start_index_usm_v(3)
1928	surf_usm_v(0)%end_index(j,i) = start_index_usm_v(0) + &
1929	num_usm_v_kji(0) - 1
1930	surf_usm_v(1)%end_index(j,i) = start_index_usm_v(1) + &
1931	num_usm_v_kji(1) - 1
1932	surf_usm_v(2)%end_index(j,i) = start_index_usm_v(2) + &
1933	num_usm_v_kji(2) - 1
1934	surf_usm_v(3)%end_index(j,i) = start_index_usm_v(3) + &
1935	num_usm_v_kji(3) - 1
1936	start_index_usm_v(0) = surf_usm_v(0)%end_index(j,i) + 1
1937	start_index_usm_v(1) = surf_usm_v(1)%end_index(j,i) + 1
1938	start_index_usm_v(2) = surf_usm_v(2)%end_index(j,i) + 1
1939	start_index_usm_v(3) = surf_usm_v(3)%end_index(j,i) + 1
1940
1941
1942	ENDDO
1943	ENDDO
1944
1945	CONTAINS
1946
1947	!------------------------------------------------------------------------------!
1948	! Description:
1949	! ------------
1950	!> Initialize horizontal surface elements, upward- and downward-facing.
1951	!> Note, horizontal surface type alsw comprises model-top fluxes, which are,
1952	!> initialized in a different routine.
1953	!------------------------------------------------------------------------------!
1954	SUBROUTINE initialize_horizontal_surfaces( k, j, i, surf, num_h, &
1955	num_h_kji, upward_facing, &
1956	downward_facing )
1957
1958	IMPLICIT NONE
1959
1960	INTEGER(iwp) :: i !< running index x-direction
1961	INTEGER(iwp) :: j !< running index y-direction
1962	INTEGER(iwp) :: k !< running index z-direction
1963	INTEGER(iwp) :: num_h !< current number of surface element
1964	INTEGER(iwp) :: num_h_kji !< dummy increment
1965	INTEGER(iwp) :: lsp !< running index chemical species
1966	INTEGER(iwp) :: lsp_pr !< running index chemical species??
1967
1968	LOGICAL :: upward_facing !< flag indicating upward-facing surface
1969	LOGICAL :: downward_facing !< flag indicating downward-facing surface
1970
1971	TYPE( surf_type ) :: surf !< respective surface type
1972	!
1973	!-- Store indices of respective surface element
1974	surf%i(num_h) = i
1975	surf%j(num_h) = j
1976	surf%k(num_h) = k
1977	!
1978	!-- Surface orientation, bit 0 is set to 1 for upward-facing surfaces,
1979	!-- bit 1 is for downward-facing surfaces.
1980	IF ( upward_facing ) surf%facing(num_h) = IBSET( surf%facing(num_h), 0 )
1981	IF ( downward_facing ) surf%facing(num_h) = IBSET( surf%facing(num_h), 1 )
1982	!
1983	!-- Initialize surface-layer height
1984	IF ( upward_facing ) THEN
1985	surf%z_mo(num_h) = zu(k) - zw(k-1)
1986	ELSE
1987	surf%z_mo(num_h) = zw(k) - zu(k)
1988	ENDIF
1989
1990	surf%z0(num_h) = roughness_length
1991	surf%z0h(num_h) = z0h_factor * roughness_length
1992	surf%z0q(num_h) = z0h_factor * roughness_length
1993	!
1994	!-- Initialization in case of 1D pre-cursor run
1995	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )&
1996	THEN
1997	IF ( .NOT. constant_diffusion ) THEN
1998	IF ( constant_flux_layer ) THEN
1999	surf%ol(num_h) = surf%z_mo(num_h) / &
2000	( rif1d(nzb+1) + 1.0E-20_wp )
2001	surf%us(num_h) = us1d
2002	surf%usws(num_h) = usws1d
2003	surf%vsws(num_h) = vsws1d
2004	ELSE
2005	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
2006	surf%us(num_h) = 0.0_wp
2007	surf%usws(num_h) = 0.0_wp
2008	surf%vsws(num_h) = 0.0_wp
2009	ENDIF
2010	ELSE
2011	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
2012	surf%us(num_h) = 0.0_wp
2013	surf%usws(num_h) = 0.0_wp
2014	surf%vsws(num_h) = 0.0_wp
2015	ENDIF
2016	!
2017	!-- Initialization in case of constant profiles
2018	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0&
2019	.OR. INDEX(initializing_actions, 'inifor' ) /= 0 ) THEN
2020
2021	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
2022	!
2023	!-- Very small number is required for calculation of Obukhov length
2024	!-- at first timestep
2025	surf%us(num_h) = 1E-30_wp
2026	surf%usws(num_h) = 0.0_wp
2027	surf%vsws(num_h) = 0.0_wp
2028
2029	ENDIF
2030
2031	surf%rib(num_h) = 0.0_wp
2032	surf%uvw_abs(num_h) = 0.0_wp
2033
2034	IF ( .NOT. constant_diffusion ) THEN
2035	surf%u_0(num_h) = 0.0_wp
2036	surf%v_0(num_h) = 0.0_wp
2037	ENDIF
2038
2039	surf%ts(num_h) = 0.0_wp
2040
2041	IF ( humidity ) THEN
2042	surf%qs(num_h) = 0.0_wp
2043	IF ( cloud_physics .AND. microphysics_morrison) THEN
2044	surf%qcs(num_h) = 0.0_wp
2045	surf%ncs(num_h) = 0.0_wp
2046
2047	surf%qcsws(num_h) = 0.0_wp
2048	surf%ncsws(num_h) = 0.0_wp
2049
2050	ENDIF
2051	IF ( cloud_physics .AND. microphysics_seifert) THEN
2052	surf%qrs(num_h) = 0.0_wp
2053	surf%nrs(num_h) = 0.0_wp
2054
2055	surf%qrsws(num_h) = 0.0_wp
2056	surf%nrsws(num_h) = 0.0_wp
2057
2058	surf%pt1(num_h) = 0.0_wp
2059	surf%qv1(num_h) = 0.0_wp
2060
2061	ENDIF
2062	ENDIF
2063
2064	IF ( passive_scalar ) surf%ss(num_h) = 0.0_wp
2065
2066	DO lsp = 1, nvar
2067	IF ( air_chemistry ) surf%css(lsp,num_h) = 0.0_wp
2068	ENDDO
2069	!
2070	!-- Set initial value for surface temperature
2071	surf%pt_surface(num_h) = pt_surface
2072	!
2073	!-- Inititalize surface fluxes of sensible and latent heat, as well as
2074	!-- passive scalar
2075	IF ( use_surface_fluxes ) THEN
2076
2077	IF ( upward_facing ) THEN
2078	IF ( constant_heatflux ) THEN
2079	!
2080	!-- Initialize surface heatflux. However, skip this for now if
2081	!-- if random_heatflux is set. This case, shf is initialized later.
2082	IF ( .NOT. random_heatflux ) THEN
2083	surf%shf(num_h) = surface_heatflux * &
2084	heatflux_input_conversion(k-1)
2085	!
2086	!-- Check if surface heat flux might be replaced by
2087	!-- prescribed wall heatflux
2088	IF ( k-1 /= 0 ) THEN
2089	surf%shf(num_h) = wall_heatflux(0) * &
2090	heatflux_input_conversion(k-1)
2091	ENDIF
2092	!
2093	!-- Initialize shf with data from external file LSF_DATA. Will
2094	!-- be done directly in ls_foring_surf
2095	!-- Attention: Just a workaround, need to be revised!!!
2096	IF ( large_scale_forcing .AND. lsf_surf ) THEN
2097	! CALL ls_forcing_surf ( simulated_time )
2098	! surf%shf(num_h) = shf(j,i)
2099	ENDIF
2100	ENDIF
2101	ELSE
2102	surf%shf(num_h) = 0.0_wp
2103	ENDIF
2104	!
2105	!-- Set heat-flux at downward-facing surfaces
2106	ELSE
2107	surf%shf(num_h) = wall_heatflux(5) * &
2108	heatflux_input_conversion(k)
2109	ENDIF
2110
2111	IF ( humidity ) THEN
2112	IF ( upward_facing ) THEN
2113	IF ( constant_waterflux ) THEN
2114	surf%qsws(num_h) = surface_waterflux * &
2115	waterflux_input_conversion(k-1)
2116	IF ( k-1 /= 0 ) THEN
2117	surf%qsws(num_h) = wall_humidityflux(0) * &
2118	waterflux_input_conversion(k-1)
2119	ENDIF
2120	ELSE
2121	surf%qsws(num_h) = 0.0_wp
2122	ENDIF
2123	ELSE
2124	surf%qsws(num_h) = wall_humidityflux(5) * &
2125	heatflux_input_conversion(k)
2126	ENDIF
2127	ENDIF
2128
2129	IF ( passive_scalar ) THEN
2130	IF ( upward_facing ) THEN
2131	IF ( constant_scalarflux ) THEN
2132	surf%ssws(num_h) = surface_scalarflux * rho_air_zw(k-1)
2133
2134	IF ( k-1 /= 0 ) &
2135	surf%ssws(num_h) = wall_scalarflux(0) * &
2136	rho_air_zw(k-1)
2137
2138	ELSE
2139	surf%ssws(num_h) = 0.0_wp
2140	ENDIF
2141	ELSE
2142	surf%ssws(num_h) = wall_scalarflux(5) * rho_air_zw(k)
2143	ENDIF
2144	ENDIF
2145
2146	IF ( air_chemistry ) THEN
2147	lsp_pr = 1
2148	DO WHILE ( TRIM( surface_csflux_name( lsp_pr ) ) /= 'novalue' ) !<'novalue' is the default
2149	DO lsp = 1, nvar
2150	!
2151	!-- Assign surface flux for each variable species
2152	IF ( TRIM( spc_names(lsp) ) == TRIM( surface_csflux_name(lsp_pr) ) ) THEN
2153	IF ( upward_facing ) THEN
2154	IF ( constant_csflux(lsp_pr) ) THEN
2155	surf%cssws(lsp,num_h) = &
2156	surface_csflux(lsp_pr) *&
2157	rho_air_zw(k-1)
2158
2159	IF ( k-1 /= 0 ) &
2160	surf%cssws(lsp,num_h) = &
2161	wall_csflux(lsp,0) * &
2162	rho_air_zw(k-1)
2163	ELSE
2164	surf%cssws(lsp,num_h) = 0.0_wp
2165	ENDIF
2166	ELSE
2167	surf%cssws(lsp,num_h) = wall_csflux(lsp,5) * &
2168	rho_air_zw(k)
2169	ENDIF
2170	ENDIF
2171	ENDDO
2172	lsp_pr = lsp_pr + 1
2173	ENDDO
2174	ENDIF
2175
2176	IF ( ocean ) THEN
2177	IF ( upward_facing ) THEN
2178	surf%sasws(num_h) = bottom_salinityflux * rho_air_zw(k-1)
2179	ELSE
2180	surf%sasws(num_h) = 0.0_wp
2181	ENDIF
2182	ENDIF
2183	ENDIF
2184	!
2185	!-- Increment surface indices
2186	num_h = num_h + 1
2187	num_h_kji = num_h_kji + 1
2188
2189
2190	END SUBROUTINE initialize_horizontal_surfaces
2191
2192
2193	!------------------------------------------------------------------------------!
2194	! Description:
2195	! ------------
2196	!> Initialize model-top fluxes. Currently, only the heatflux and salinity flux
2197	!> can be prescribed, latent flux is zero in this case!
2198	!------------------------------------------------------------------------------!
2199	SUBROUTINE initialize_top( k, j, i, surf, num_h, num_h_kji )
2200
2201	IMPLICIT NONE
2202
2203	INTEGER(iwp) :: i !< running index x-direction
2204	INTEGER(iwp) :: j !< running index y-direction
2205	INTEGER(iwp) :: k !< running index z-direction
2206	INTEGER(iwp) :: num_h !< current number of surface element
2207	INTEGER(iwp) :: num_h_kji !< dummy increment
2208	INTEGER(iwp) :: lsp !< running index for chemical species
2209
2210	TYPE( surf_type ) :: surf !< respective surface type
2211	!
2212	!-- Store indices of respective surface element
2213	surf%i(num_h) = i
2214	surf%j(num_h) = j
2215	surf%k(num_h) = k
2216	!
2217	!-- Initialize top heat flux
2218	IF ( constant_top_heatflux ) &
2219	surf%shf(num_h) = top_heatflux * heatflux_input_conversion(nzt+1)
2220	!
2221	!-- Initialization in case of a coupled model run
2222	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
2223	surf%shf(num_h) = 0.0_wp
2224	surf%qsws(num_h) = 0.0_wp
2225	ENDIF
2226	!
2227	!-- Prescribe latent heat flux at the top
2228	IF ( humidity ) THEN
2229	surf%qsws(num_h) = 0.0_wp
2230	IF ( cloud_physics .AND. microphysics_morrison ) THEN
2231	surf%ncsws(num_h) = 0.0_wp
2232	surf%qcsws(num_h) = 0.0_wp
2233	ENDIF
2234	IF ( cloud_physics .AND. microphysics_seifert ) THEN
2235	surf%nrsws(num_h) = 0.0_wp
2236	surf%qrsws(num_h) = 0.0_wp
2237	ENDIF
2238	ENDIF
2239	!
2240	!-- Prescribe top scalar flux
2241	IF ( passive_scalar .AND. constant_top_scalarflux ) &
2242	surf%ssws(num_h) = top_scalarflux * rho_air_zw(nzt+1)
2243	!
2244	!-- Prescribe top chemical species' flux
2245	DO lsp = 1, nvar
2246	IF ( air_chemistry .AND. constant_top_csflux(lsp) ) THEN
2247	surf%cssws(lsp,num_h) = top_csflux(lsp) * rho_air_zw(nzt+1)
2248	ENDIF
2249	ENDDO
2250	!
2251	!-- Prescribe top salinity flux
2252	IF ( ocean .AND. constant_top_salinityflux) &
2253	surf%sasws(num_h) = top_salinityflux * rho_air_zw(nzt+1)
2254	!
2255	!-- Top momentum fluxes
2256	IF ( constant_top_momentumflux ) THEN
2257	surf%usws(num_h) = top_momentumflux_u * &
2258	momentumflux_input_conversion(nzt+1)
2259	surf%vsws(num_h) = top_momentumflux_v * &
2260	momentumflux_input_conversion(nzt+1)
2261	ENDIF
2262	!
2263	!-- Increment surface indices
2264	num_h = num_h + 1
2265	num_h_kji = num_h_kji + 1
2266
2267
2268	END SUBROUTINE initialize_top
2269
2270
2271	!------------------------------------------------------------------------------!
2272	! Description:
2273	! ------------
2274	!> Initialize vertical surface elements.
2275	!------------------------------------------------------------------------------!
2276	SUBROUTINE initialize_vertical_surfaces( l, k, j, i, surf, num_v, &
2277	num_v_kji, east_facing, &
2278	west_facing, south_facing, &
2279	north_facing )
2280
2281	IMPLICIT NONE
2282
2283	INTEGER(iwp) :: component !< index of wall_fluxes_ array for respective orientation
2284	INTEGER(iwp) :: i !< running index x-direction
2285	INTEGER(iwp) :: j !< running index x-direction
2286	INTEGER(iwp) :: k !< running index x-direction
2287	INTEGER(iwp) :: l !< index variable for the surface type, indicating the facing
2288	INTEGER(iwp) :: num_v !< current number of surface element
2289	INTEGER(iwp) :: num_v_kji !< current number of surface element at (j,i)
2290	INTEGER(iwp) :: lsp !< running index for chemical species
2291
2292
2293	LOGICAL :: east_facing !< flag indicating east-facing surfaces
2294	LOGICAL :: north_facing !< flag indicating north-facing surfaces
2295	LOGICAL :: south_facing !< flag indicating south-facing surfaces
2296	LOGICAL :: west_facing !< flag indicating west-facing surfaces
2297
2298	TYPE( surf_type ) :: surf !< respective surface type
2299
2300	!
2301	!-- Store indices of respective wall element
2302	surf%i(num_v) = i
2303	surf%j(num_v) = j
2304	surf%k(num_v) = k
2305	!
2306	!-- Initialize surface-layer height, or more precisely, distance to surface
2307	IF ( north_facing .OR. south_facing ) THEN
2308	surf%z_mo(num_v) = 0.5_wp * dy
2309	ELSE
2310	surf%z_mo(num_v) = 0.5_wp * dx
2311	ENDIF
2312
2313	surf%facing(num_v) = 0
2314	!
2315	!-- Surface orientation. Moreover, set component id to map wall_heatflux,
2316	!-- etc., on surface type (further below)
2317	IF ( north_facing ) THEN
2318	surf%facing(num_v) = 5 !IBSET( surf%facing(num_v), 0 )
2319	component = 4
2320	ENDIF
2321
2322	IF ( south_facing ) THEN
2323	surf%facing(num_v) = 6 !IBSET( surf%facing(num_v), 1 )
2324	component = 3
2325	ENDIF
2326
2327	IF ( east_facing ) THEN
2328	surf%facing(num_v) = 7 !IBSET( surf%facing(num_v), 2 )
2329	component = 2
2330	ENDIF
2331
2332	IF ( west_facing ) THEN
2333	surf%facing(num_v) = 8 !IBSET( surf%facing(num_v), 3 )
2334	component = 1
2335	ENDIF
2336
2337
2338	surf%z0(num_v) = roughness_length
2339	surf%z0h(num_v) = z0h_factor * roughness_length
2340	surf%z0q(num_v) = z0h_factor * roughness_length
2341
2342	surf%us(num_v) = 0.0_wp
2343	!
2344	!-- If required, initialize Obukhov length
2345	IF ( ALLOCATED( surf%ol ) ) &
2346	surf%ol(num_v) = surf%z_mo(num_v) / zeta_min
2347
2348	surf%uvw_abs(num_v) = 0.0_wp
2349
2350	surf%mom_flux_uv(num_v) = 0.0_wp
2351	surf%mom_flux_w(num_v) = 0.0_wp
2352	surf%mom_flux_tke(0:1,num_v) = 0.0_wp
2353
2354	surf%ts(num_v) = 0.0_wp
2355	surf%shf(num_v) = wall_heatflux(component)
2356	!
2357	!-- Set initial value for surface temperature
2358	surf%pt_surface(num_v) = pt_surface
2359
2360	IF ( humidity ) THEN
2361	surf%qs(num_v) = 0.0_wp
2362	surf%qsws(num_v) = wall_humidityflux(component)
2363	!
2364	!-- Following wall fluxes are assumed to be zero
2365	IF ( cloud_physics .AND. microphysics_morrison) THEN
2366	surf%qcs(num_v) = 0.0_wp
2367	surf%ncs(num_v) = 0.0_wp
2368
2369	surf%qcsws(num_v) = 0.0_wp
2370	surf%ncsws(num_v) = 0.0_wp
2371	ENDIF
2372	IF ( cloud_physics .AND. microphysics_seifert) THEN
2373	surf%qrs(num_v) = 0.0_wp
2374	surf%nrs(num_v) = 0.0_wp
2375
2376	surf%qrsws(num_v) = 0.0_wp
2377	surf%nrsws(num_v) = 0.0_wp
2378	ENDIF
2379	ENDIF
2380
2381	IF ( passive_scalar ) THEN
2382	surf%ss(num_v) = 0.0_wp
2383	surf%ssws(num_v) = wall_scalarflux(component)
2384	ENDIF
2385
2386	IF ( air_chemistry ) THEN
2387	DO lsp = 1, nvar
2388	surf%css(lsp,num_v) = 0.0_wp
2389	surf%cssws(lsp,num_v) = wall_csflux(lsp,component)
2390	ENDDO
2391	ENDIF
2392
2393	!
2394	!-- So far, salinityflux at vertical surfaces is simply zero
2395	!-- at the moment
2396	IF ( ocean ) surf%sasws(num_v) = wall_salinityflux(component)
2397	!
2398	!-- Increment wall indices
2399	num_v = num_v + 1
2400	num_v_kji = num_v_kji + 1
2401
2402	END SUBROUTINE initialize_vertical_surfaces
2403
2404	END SUBROUTINE init_surfaces
2405
2406
2407	!------------------------------------------------------------------------------!
2408	! Description:
2409	! ------------
2410	!> Determines topography-top index at given (j,i)-position.
2411	!------------------------------------------------------------------------------!
2412	FUNCTION get_topography_top_index_ji( j, i, grid )
2413
2414	IMPLICIT NONE
2415
2416	CHARACTER(LEN=*) :: grid !< flag to distinquish between staggered grids
2417	INTEGER(iwp) :: i !< grid index in x-dimension
2418	INTEGER(iwp) :: ibit !< bit position where topography information is stored on respective grid
2419	INTEGER(iwp) :: j !< grid index in y-dimension
2420	INTEGER(iwp) :: get_topography_top_index_ji !< topography top index
2421
2422	SELECT CASE ( TRIM( grid ) )
2423
2424	CASE ( 's' )
2425	ibit = 12
2426	CASE ( 'u' )
2427	ibit = 14
2428	CASE ( 'v' )
2429	ibit = 16
2430	CASE ( 'w' )
2431	ibit = 18
2432	CASE ( 's_out' )
2433	ibit = 24
2434	CASE DEFAULT
2435	!
2436	!-- Set default to scalar grid
2437	ibit = 12
2438
2439	END SELECT
2440
2441	get_topography_top_index_ji = MAXLOC( &
2442	MERGE( 1, 0, &
2443	BTEST( wall_flags_0(:,j,i), ibit ) &
2444	), DIM = 1 &
2445	) - 1
2446
2447	RETURN
2448
2449	END FUNCTION get_topography_top_index_ji
2450
2451	!------------------------------------------------------------------------------!
2452	! Description:
2453	! ------------
2454	!> Determines topography-top index at each (j,i)-position.
2455	!------------------------------------------------------------------------------!
2456	FUNCTION get_topography_top_index( grid )
2457
2458	IMPLICIT NONE
2459
2460	CHARACTER(LEN=*) :: grid !< flag to distinquish between staggered grids
2461	INTEGER(iwp) :: ibit !< bit position where topography information is stored on respective grid
2462	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: get_topography_top_index !< topography top index
2463
2464	SELECT CASE ( TRIM( grid ) )
2465
2466	CASE ( 's' )
2467	ibit = 12
2468	CASE ( 'u' )
2469	ibit = 14
2470	CASE ( 'v' )
2471	ibit = 16
2472	CASE ( 'w' )
2473	ibit = 18
2474	CASE ( 's_out' )
2475	ibit = 24
2476	CASE DEFAULT
2477	!
2478	!-- Set default to scalar grid
2479	ibit = 12
2480
2481	END SELECT
2482
2483	get_topography_top_index(nys:nyn,nxl:nxr) = MAXLOC( &
2484	MERGE( 1, 0, &
2485	BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), ibit )&
2486	), DIM = 1 &
2487	) - 1
2488
2489	RETURN
2490
2491	END FUNCTION get_topography_top_index
2492
2493	!------------------------------------------------------------------------------!
2494	! Description:
2495	! ------------
2496	!> Gathers all surface elements with the same facing (but possibly different
2497	!> type) onto a surface type, and writes binary data into restart files.
2498	!------------------------------------------------------------------------------!
2499	SUBROUTINE surface_wrd_local
2500
2501
2502	IMPLICIT NONE
2503
2504	CHARACTER(LEN=1) :: dum !< dummy string to create output-variable name
2505
2506	INTEGER(iwp) :: i !< running index x-direction
2507	INTEGER(iwp) :: j !< running index y-direction
2508	INTEGER(iwp) :: l !< index surface type orientation
2509	INTEGER(iwp) :: lsp !< running index chemical species
2510	INTEGER(iwp) :: m !< running index for surface elements on individual surface array
2511	INTEGER(iwp), DIMENSION(0:2) :: start_index_h !< start index for horizontal surface elements on gathered surface array
2512	INTEGER(iwp), DIMENSION(0:3) :: mm !< running index for surface elements on gathered surface array
2513	INTEGER(iwp), DIMENSION(0:3) :: start_index_v !< start index for vertical surface elements on gathered surface array
2514
2515	TYPE(surf_type), DIMENSION(0:2) :: surf_h !< gathered horizontal surfaces, contains all surface types
2516	TYPE(surf_type), DIMENSION(0:3) :: surf_v !< gathered vertical surfaces, contains all surface types
2517
2518	!
2519	!-- Determine total number of horizontal and vertical surface elements before
2520	!-- writing var_list
2521	CALL surface_last_actions
2522	!
2523	!-- Count number of grid points with same facing and allocate attributes respectively
2524	!-- Horizontal upward facing
2525	surf_h(0)%ns = ns_h_on_file(0)
2526	CALL allocate_surface_attributes_h( surf_h(0), nys, nyn, nxl, nxr )
2527	!
2528	!-- Horizontal downward facing
2529	surf_h(1)%ns = ns_h_on_file(1)
2530	CALL allocate_surface_attributes_h( surf_h(1), nys, nyn, nxl, nxr )
2531	!
2532	!-- Model top
2533	surf_h(2)%ns = ns_h_on_file(2)
2534	CALL allocate_surface_attributes_h_top( surf_h(2), nys, nyn, nxl, nxr )
2535	!
2536	!-- Vertical surfaces
2537	DO l = 0, 3
2538	surf_v(l)%ns = ns_v_on_file(l)
2539	CALL allocate_surface_attributes_v( surf_v(l), &
2540	nys, nyn, nxl, nxr )
2541	ENDDO
2542	!
2543	!-- In the following, gather data from surfaces elements with the same
2544	!-- facing (but possibly differt type) on 1 data-type array.
2545	mm(0:2) = 1
2546	DO l = 0, 2
2547	DO i = nxl, nxr
2548	DO j = nys, nyn
2549	DO m = surf_def_h(l)%start_index(j,i), &
2550	surf_def_h(l)%end_index(j,i)
2551	IF ( ALLOCATED( surf_def_h(l)%us ) ) &
2552	surf_h(l)%us(mm(l)) = surf_def_h(l)%us(m)
2553	IF ( ALLOCATED( surf_def_h(l)%ts ) ) &
2554	surf_h(l)%ts(mm(l)) = surf_def_h(l)%ts(m)
2555	IF ( ALLOCATED( surf_def_h(l)%qs ) ) &
2556	surf_h(l)%qs(mm(l)) = surf_def_h(l)%qs(m)
2557	IF ( ALLOCATED( surf_def_h(l)%ss ) ) &
2558	surf_h(l)%ss(mm(l)) = surf_def_h(l)%ss(m)
2559	IF ( ALLOCATED( surf_def_h(l)%qcs ) ) &
2560	surf_h(l)%qcs(mm(l)) = surf_def_h(l)%qcs(m)
2561	IF ( ALLOCATED( surf_def_h(l)%ncs ) ) &
2562	surf_h(l)%ncs(mm(l)) = surf_def_h(l)%ncs(m)
2563	IF ( ALLOCATED( surf_def_h(l)%qrs ) ) &
2564	surf_h(l)%qrs(mm(l)) = surf_def_h(l)%qrs(m)
2565	IF ( ALLOCATED( surf_def_h(l)%nrs ) ) &
2566	surf_h(l)%nrs(mm(l)) = surf_def_h(l)%nrs(m)
2567	IF ( ALLOCATED( surf_def_h(l)%ol ) ) &
2568	surf_h(l)%ol(mm(l)) = surf_def_h(l)%ol(m)
2569	IF ( ALLOCATED( surf_def_h(l)%rib ) ) &
2570	surf_h(l)%rib(mm(l)) = surf_def_h(l)%rib(m)
2571	IF ( ALLOCATED( surf_def_h(l)%pt_surface ) ) &
2572	surf_h(l)%pt_surface(mm(l)) = surf_def_h(l)%pt_surface(m)
2573	IF ( ALLOCATED( surf_def_h(l)%usws ) ) &
2574	surf_h(l)%usws(mm(l)) = surf_def_h(l)%usws(m)
2575	IF ( ALLOCATED( surf_def_h(l)%vsws ) ) &
2576	surf_h(l)%vsws(mm(l)) = surf_def_h(l)%vsws(m)
2577	IF ( ALLOCATED( surf_def_h(l)%shf ) ) &
2578	surf_h(l)%shf(mm(l)) = surf_def_h(l)%shf(m)
2579	IF ( ALLOCATED( surf_def_h(l)%qsws ) ) &
2580	surf_h(l)%qsws(mm(l)) = surf_def_h(l)%qsws(m)
2581	IF ( ALLOCATED( surf_def_h(l)%ssws ) ) &
2582	surf_h(l)%ssws(mm(l)) = surf_def_h(l)%ssws(m)
2583	IF ( ALLOCATED( surf_def_h(l)%css ) ) THEN
2584	DO lsp = 1,nvar
2585	surf_h(l)%css(lsp,mm(l)) = surf_def_h(l)%css(lsp,m)
2586	ENDDO
2587	ENDIF
2588	IF ( ALLOCATED( surf_def_h(l)%cssws ) ) THEN
2589	DO lsp = 1,nvar
2590	surf_h(l)%cssws(lsp,mm(l)) = surf_def_h(l)%cssws(lsp,m)
2591	ENDDO
2592	ENDIF
2593	IF ( ALLOCATED( surf_def_h(l)%ncsws ) ) &
2594	surf_h(l)%ncsws(mm(l)) = surf_def_h(l)%ncsws(m)
2595	IF ( ALLOCATED( surf_def_h(l)%nrsws ) ) &
2596	surf_h(l)%nrsws(mm(l)) = surf_def_h(l)%nrsws(m)
2597	IF ( ALLOCATED( surf_def_h(l)%sasws ) ) &
2598	surf_h(l)%sasws(mm(l)) = surf_def_h(l)%sasws(m)
2599
2600	mm(l) = mm(l) + 1
2601	ENDDO
2602
2603	IF ( l == 0 ) THEN
2604	DO m = surf_lsm_h%start_index(j,i), &
2605	surf_lsm_h%end_index(j,i)
2606	IF ( ALLOCATED( surf_lsm_h%us ) ) &
2607	surf_h(0)%us(mm(0)) = surf_lsm_h%us(m)
2608	IF ( ALLOCATED( surf_lsm_h%ts ) ) &
2609	surf_h(0)%ts(mm(0)) = surf_lsm_h%ts(m)
2610	IF ( ALLOCATED( surf_lsm_h%qs ) ) &
2611	surf_h(0)%qs(mm(0)) = surf_lsm_h%qs(m)
2612	IF ( ALLOCATED( surf_lsm_h%ss ) ) &
2613	surf_h(0)%ss(mm(0)) = surf_lsm_h%ss(m)
2614	IF ( ALLOCATED( surf_lsm_h%qcs ) ) &
2615	surf_h(0)%qcs(mm(0)) = surf_lsm_h%qcs(m)
2616	IF ( ALLOCATED( surf_lsm_h%ncs ) ) &
2617	surf_h(0)%ncs(mm(0)) = surf_lsm_h%ncs(m)
2618	IF ( ALLOCATED( surf_lsm_h%qrs ) ) &
2619	surf_h(0)%qrs(mm(0)) = surf_lsm_h%qrs(m)
2620	IF ( ALLOCATED( surf_lsm_h%nrs ) ) &
2621	surf_h(0)%nrs(mm(0)) = surf_lsm_h%nrs(m)
2622	IF ( ALLOCATED( surf_lsm_h%ol ) ) &
2623	surf_h(0)%ol(mm(0)) = surf_lsm_h%ol(m)
2624	IF ( ALLOCATED( surf_lsm_h%rib ) ) &
2625	surf_h(0)%rib(mm(0)) = surf_lsm_h%rib(m)
2626	IF ( ALLOCATED( surf_lsm_h%pt_surface ) ) &
2627	surf_h(l)%pt_surface(mm(l)) = surf_lsm_h%pt_surface(m)
2628	IF ( ALLOCATED( surf_lsm_h%usws ) ) &
2629	surf_h(0)%usws(mm(0)) = surf_lsm_h%usws(m)
2630	IF ( ALLOCATED( surf_lsm_h%vsws ) ) &
2631	surf_h(0)%vsws(mm(0)) = surf_lsm_h%vsws(m)
2632	IF ( ALLOCATED( surf_lsm_h%shf ) ) &
2633	surf_h(0)%shf(mm(0)) = surf_lsm_h%shf(m)
2634	IF ( ALLOCATED( surf_lsm_h%qsws ) ) &
2635	surf_h(0)%qsws(mm(0)) = surf_lsm_h%qsws(m)
2636	IF ( ALLOCATED( surf_lsm_h%ssws ) ) &
2637	surf_h(0)%ssws(mm(0)) = surf_lsm_h%ssws(m)
2638	IF ( ALLOCATED( surf_lsm_h%css ) ) THEN
2639	DO lsp = 1, nvar
2640	surf_h(0)%css(lsp,mm(0)) = surf_lsm_h%css(lsp,m)
2641	ENDDO
2642	ENDIF
2643	IF ( ALLOCATED( surf_lsm_h%cssws ) ) THEN
2644	DO lsp = 1, nvar
2645	surf_h(0)%cssws(lsp,mm(0)) = surf_lsm_h%cssws(lsp,m)
2646	ENDDO
2647	ENDIF
2648	IF ( ALLOCATED( surf_lsm_h%ncsws ) ) &
2649	surf_h(0)%ncsws(mm(0)) = surf_lsm_h%ncsws(m)
2650	IF ( ALLOCATED( surf_lsm_h%nrsws ) ) &
2651	surf_h(0)%nrsws(mm(0)) = surf_lsm_h%nrsws(m)
2652	IF ( ALLOCATED( surf_lsm_h%sasws ) ) &
2653	surf_h(0)%sasws(mm(0)) = surf_lsm_h%sasws(m)
2654
2655	mm(0) = mm(0) + 1
2656
2657	ENDDO
2658
2659	DO m = surf_usm_h%start_index(j,i), &
2660	surf_usm_h%end_index(j,i)
2661	IF ( ALLOCATED( surf_usm_h%us ) ) &
2662	surf_h(0)%us(mm(0)) = surf_usm_h%us(m)
2663	IF ( ALLOCATED( surf_usm_h%ts ) ) &
2664	surf_h(0)%ts(mm(0)) = surf_usm_h%ts(m)
2665	IF ( ALLOCATED( surf_usm_h%qs ) ) &
2666	surf_h(0)%qs(mm(0)) = surf_usm_h%qs(m)
2667	IF ( ALLOCATED( surf_usm_h%ss ) ) &
2668	surf_h(0)%ss(mm(0)) = surf_usm_h%ss(m)
2669	IF ( ALLOCATED( surf_usm_h%qcs ) ) &
2670	surf_h(0)%qcs(mm(0)) = surf_usm_h%qcs(m)
2671	IF ( ALLOCATED( surf_usm_h%ncs ) ) &
2672	surf_h(0)%ncs(mm(0)) = surf_usm_h%ncs(m)
2673	IF ( ALLOCATED( surf_usm_h%qrs ) ) &
2674	surf_h(0)%qrs(mm(0)) = surf_usm_h%qrs(m)
2675	IF ( ALLOCATED( surf_usm_h%nrs ) ) &
2676	surf_h(0)%nrs(mm(0)) = surf_usm_h%nrs(m)
2677	IF ( ALLOCATED( surf_usm_h%ol ) ) &
2678	surf_h(0)%ol(mm(0)) = surf_usm_h%ol(m)
2679	IF ( ALLOCATED( surf_usm_h%rib ) ) &
2680	surf_h(0)%rib(mm(0)) = surf_usm_h%rib(m)
2681	IF ( ALLOCATED( surf_usm_h%pt_surface ) ) &
2682	surf_h(l)%pt_surface(mm(l)) = surf_usm_h%pt_surface(m)
2683	IF ( ALLOCATED( surf_usm_h%usws ) ) &
2684	surf_h(0)%usws(mm(0)) = surf_usm_h%usws(m)
2685	IF ( ALLOCATED( surf_usm_h%vsws ) ) &
2686	surf_h(0)%vsws(mm(0)) = surf_usm_h%vsws(m)
2687	IF ( ALLOCATED( surf_usm_h%shf ) ) &
2688	surf_h(0)%shf(mm(0)) = surf_usm_h%shf(m)
2689	IF ( ALLOCATED( surf_usm_h%qsws ) ) &
2690	surf_h(0)%qsws(mm(0)) = surf_usm_h%qsws(m)
2691	IF ( ALLOCATED( surf_usm_h%ssws ) ) &
2692	surf_h(0)%ssws(mm(0)) = surf_usm_h%ssws(m)
2693	IF ( ALLOCATED( surf_usm_h%css ) ) THEN
2694	DO lsp = 1, nvar
2695	surf_h(0)%css(lsp,mm(0)) = surf_usm_h%css(lsp,m)
2696	ENDDO
2697	ENDIF
2698	IF ( ALLOCATED( surf_usm_h%cssws ) ) THEN
2699	DO lsp = 1, nvar
2700	surf_h(0)%cssws(lsp,mm(0)) = surf_usm_h%cssws(lsp,m)
2701	ENDDO
2702	ENDIF
2703	IF ( ALLOCATED( surf_usm_h%ncsws ) ) &
2704	surf_h(0)%ncsws(mm(0)) = surf_usm_h%ncsws(m)
2705	IF ( ALLOCATED( surf_usm_h%nrsws ) ) &
2706	surf_h(0)%nrsws(mm(0)) = surf_usm_h%nrsws(m)
2707	IF ( ALLOCATED( surf_usm_h%sasws ) ) &
2708	surf_h(0)%sasws(mm(0)) = surf_usm_h%sasws(m)
2709
2710	mm(0) = mm(0) + 1
2711
2712	ENDDO
2713
2714
2715	ENDIF
2716
2717	ENDDO
2718
2719	ENDDO
2720	!
2721	!-- Gather start- and end indices
2722	start_index_h(l) = 1
2723	DO i = nxl, nxr
2724	DO j = nys, nyn
2725
2726	surf_h(l)%start_index(j,i) = start_index_h(l)
2727	surf_h(l)%end_index(j,i) = surf_h(l)%start_index(j,i) -1
2728
2729	DO m = surf_def_h(l)%start_index(j,i), &
2730	surf_def_h(l)%end_index(j,i)
2731	surf_h(l)%end_index(j,i) = surf_h(l)%end_index(j,i) + 1
2732	ENDDO
2733	IF ( l == 0 ) THEN
2734	DO m = surf_lsm_h%start_index(j,i), &
2735	surf_lsm_h%end_index(j,i)
2736	surf_h(l)%end_index(j,i) = surf_h(l)%end_index(j,i) + 1
2737	ENDDO
2738	DO m = surf_usm_h%start_index(j,i), &
2739	surf_usm_h%end_index(j,i)
2740	surf_h(l)%end_index(j,i) = surf_h(l)%end_index(j,i) + 1
2741	ENDDO
2742	ENDIF
2743
2744	start_index_h(l) = surf_h(l)%end_index(j,i) + 1
2745
2746	ENDDO
2747	ENDDO
2748	ENDDO
2749
2750
2751	mm(0:3) = 1
2752	DO l = 0, 3
2753	DO i = nxl, nxr
2754	DO j = nys, nyn
2755	DO m = surf_def_v(l)%start_index(j,i), &
2756	surf_def_v(l)%end_index(j,i)
2757	IF ( ALLOCATED( surf_def_v(l)%us ) ) &
2758	surf_v(l)%us(mm(l)) = surf_def_v(l)%us(m)
2759	IF ( ALLOCATED( surf_def_v(l)%ts ) ) &
2760	surf_v(l)%ts(mm(l)) = surf_def_v(l)%ts(m)
2761	IF ( ALLOCATED( surf_def_v(l)%qs ) ) &
2762	surf_v(l)%qs(mm(l)) = surf_def_v(l)%qs(m)
2763	IF ( ALLOCATED( surf_def_v(l)%ss ) ) &
2764	surf_v(l)%ss(mm(l)) = surf_def_v(l)%ss(m)
2765	IF ( ALLOCATED( surf_def_v(l)%qcs ) ) &
2766	surf_v(l)%qcs(mm(l)) = surf_def_v(l)%qcs(m)
2767	IF ( ALLOCATED( surf_def_v(l)%ncs ) ) &
2768	surf_v(l)%ncs(mm(l)) = surf_def_v(l)%ncs(m)
2769	IF ( ALLOCATED( surf_def_v(l)%qrs ) ) &
2770	surf_v(l)%qrs(mm(l)) = surf_def_v(l)%qrs(m)
2771	IF ( ALLOCATED( surf_def_v(l)%nrs ) ) &
2772	surf_v(l)%nrs(mm(l)) = surf_def_v(l)%nrs(m)
2773	IF ( ALLOCATED( surf_def_v(l)%ol ) ) &
2774	surf_v(l)%ol(mm(l)) = surf_def_v(l)%ol(m)
2775	IF ( ALLOCATED( surf_def_v(l)%rib ) ) &
2776	surf_v(l)%rib(mm(l)) = surf_def_v(l)%rib(m)
2777	IF ( ALLOCATED( surf_def_v(l)%pt_surface ) ) &
2778	surf_v(l)%pt_surface(mm(l)) = surf_def_v(l)%pt_surface(m)
2779	IF ( ALLOCATED( surf_def_v(l)%shf ) ) &
2780	surf_v(l)%shf(mm(l)) = surf_def_v(l)%shf(m)
2781	IF ( ALLOCATED( surf_def_v(l)%qsws ) ) &
2782	surf_v(l)%qsws(mm(l)) = surf_def_v(l)%qsws(m)
2783	IF ( ALLOCATED( surf_def_v(l)%ssws ) ) &
2784	surf_v(l)%ssws(mm(l)) = surf_def_v(l)%ssws(m)
2785	IF ( ALLOCATED( surf_def_v(l)%css ) ) THEN
2786	DO lsp = 1, nvar
2787	surf_v(l)%css(lsp,mm(l)) = surf_def_v(l)%css(lsp,m)
2788	ENDDO
2789	ENDIF
2790	IF ( ALLOCATED( surf_def_v(l)%cssws ) ) THEN
2791	DO lsp = 1, nvar
2792	surf_v(l)%cssws(lsp,mm(l)) = surf_def_v(l)%cssws(lsp,m)
2793	ENDDO
2794	ENDIF
2795	IF ( ALLOCATED( surf_def_v(l)%ncsws ) ) &
2796	surf_v(l)%ncsws(mm(l)) = surf_def_v(l)%ncsws(m)
2797	IF ( ALLOCATED( surf_def_v(l)%nrsws ) ) &
2798	surf_v(l)%nrsws(mm(l)) = surf_def_v(l)%nrsws(m)
2799	IF ( ALLOCATED( surf_def_v(l)%sasws ) ) &
2800	surf_v(l)%sasws(mm(l)) = surf_def_v(l)%sasws(m)
2801	IF ( ALLOCATED( surf_def_v(l)%mom_flux_uv) ) &
2802	surf_v(l)%mom_flux_uv(mm(l)) = surf_def_v(l)%mom_flux_uv(m)
2803	IF ( ALLOCATED( surf_def_v(l)%mom_flux_w) ) &
2804	surf_v(l)%mom_flux_w(mm(l)) = surf_def_v(l)%mom_flux_w(m)
2805	IF ( ALLOCATED( surf_def_v(l)%mom_flux_tke) ) &
2806	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_def_v(l)%mom_flux_tke(0:1,m)
2807
2808	mm(l) = mm(l) + 1
2809	ENDDO
2810
2811	DO m = surf_lsm_v(l)%start_index(j,i), &
2812	surf_lsm_v(l)%end_index(j,i)
2813	IF ( ALLOCATED( surf_lsm_v(l)%us ) ) &
2814	surf_v(l)%us(mm(l)) = surf_lsm_v(l)%us(m)
2815	IF ( ALLOCATED( surf_lsm_v(l)%ts ) ) &
2816	surf_v(l)%ts(mm(l)) = surf_lsm_v(l)%ts(m)
2817	IF ( ALLOCATED( surf_lsm_v(l)%qs ) ) &
2818	surf_v(l)%qs(mm(l)) = surf_lsm_v(l)%qs(m)
2819	IF ( ALLOCATED( surf_lsm_v(l)%ss ) ) &
2820	surf_v(l)%ss(mm(l)) = surf_lsm_v(l)%ss(m)
2821	IF ( ALLOCATED( surf_lsm_v(l)%qcs ) ) &
2822	surf_v(l)%qcs(mm(l)) = surf_lsm_v(l)%qcs(m)
2823	IF ( ALLOCATED( surf_lsm_v(l)%ncs ) ) &
2824	surf_v(l)%ncs(mm(l)) = surf_lsm_v(l)%ncs(m)
2825	IF ( ALLOCATED( surf_lsm_v(l)%qrs ) ) &
2826	surf_v(l)%qrs(mm(l)) = surf_lsm_v(l)%qrs(m)
2827	IF ( ALLOCATED( surf_lsm_v(l)%nrs ) ) &
2828	surf_v(l)%nrs(mm(l)) = surf_lsm_v(l)%nrs(m)
2829	IF ( ALLOCATED( surf_lsm_v(l)%ol ) ) &
2830	surf_v(l)%ol(mm(l)) = surf_lsm_v(l)%ol(m)
2831	IF ( ALLOCATED( surf_lsm_v(l)%rib ) ) &
2832	surf_v(l)%rib(mm(l)) = surf_lsm_v(l)%rib(m)
2833	IF ( ALLOCATED( surf_lsm_v(l)%pt_surface ) ) &
2834	surf_v(l)%pt_surface(mm(l)) = surf_lsm_v(l)%pt_surface(m)
2835	IF ( ALLOCATED( surf_lsm_v(l)%usws ) ) &
2836	surf_v(l)%usws(mm(l)) = surf_lsm_v(l)%usws(m)
2837	IF ( ALLOCATED( surf_lsm_v(l)%vsws ) ) &
2838	surf_v(l)%vsws(mm(l)) = surf_lsm_v(l)%vsws(m)
2839	IF ( ALLOCATED( surf_lsm_v(l)%shf ) ) &
2840	surf_v(l)%shf(mm(l)) = surf_lsm_v(l)%shf(m)
2841	IF ( ALLOCATED( surf_lsm_v(l)%qsws ) ) &
2842	surf_v(l)%qsws(mm(l)) = surf_lsm_v(l)%qsws(m)
2843	IF ( ALLOCATED( surf_lsm_v(l)%ssws ) ) &
2844	surf_v(l)%ssws(mm(l)) = surf_lsm_v(l)%ssws(m)
2845	IF ( ALLOCATED( surf_lsm_v(l)%css ) ) THEN
2846	DO lsp = 1, nvar
2847	surf_v(l)%css(lsp,mm(l)) = surf_lsm_v(l)%css(lsp,m)
2848	ENDDO
2849	ENDIF
2850	IF ( ALLOCATED( surf_lsm_v(l)%cssws ) ) THEN
2851	DO lsp = 1, nvar
2852	surf_v(l)%cssws(lsp,mm(l)) = surf_lsm_v(l)%cssws(lsp,m)
2853	ENDDO
2854	ENDIF
2855	IF ( ALLOCATED( surf_lsm_v(l)%ncsws ) ) &
2856	surf_v(l)%ncsws(mm(l)) = surf_lsm_v(l)%ncsws(m)
2857	IF ( ALLOCATED( surf_lsm_v(l)%nrsws ) ) &
2858	surf_v(l)%nrsws(mm(l)) = surf_lsm_v(l)%nrsws(m)
2859	IF ( ALLOCATED( surf_lsm_v(l)%sasws ) ) &
2860	surf_v(l)%sasws(mm(l)) = surf_lsm_v(l)%sasws(m)
2861	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_uv) ) &
2862	surf_v(l)%mom_flux_uv(mm(l)) = surf_lsm_v(l)%mom_flux_uv(m)
2863	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_w) ) &
2864	surf_v(l)%mom_flux_w(mm(l)) = surf_lsm_v(l)%mom_flux_w(m)
2865	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_tke) ) &
2866	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_lsm_v(l)%mom_flux_tke(0:1,m)
2867
2868	mm(l) = mm(l) + 1
2869	ENDDO
2870
2871	DO m = surf_usm_v(l)%start_index(j,i), &
2872	surf_usm_v(l)%end_index(j,i)
2873	IF ( ALLOCATED( surf_usm_v(l)%us ) ) &
2874	surf_v(l)%us(mm(l)) = surf_usm_v(l)%us(m)
2875	IF ( ALLOCATED( surf_usm_v(l)%ts ) ) &
2876	surf_v(l)%ts(mm(l)) = surf_usm_v(l)%ts(m)
2877	IF ( ALLOCATED( surf_usm_v(l)%qs ) ) &
2878	surf_v(l)%qs(mm(l)) = surf_usm_v(l)%qs(m)
2879	IF ( ALLOCATED( surf_usm_v(l)%ss ) ) &
2880	surf_v(l)%ss(mm(l)) = surf_usm_v(l)%ss(m)
2881	IF ( ALLOCATED( surf_usm_v(l)%qcs ) ) &
2882	surf_v(l)%qcs(mm(l)) = surf_usm_v(l)%qcs(m)
2883	IF ( ALLOCATED( surf_usm_v(l)%ncs ) ) &
2884	surf_v(l)%ncs(mm(l)) = surf_usm_v(l)%ncs(m)
2885	IF ( ALLOCATED( surf_usm_v(l)%qrs ) ) &
2886	surf_v(l)%qrs(mm(l)) = surf_usm_v(l)%qrs(m)
2887	IF ( ALLOCATED( surf_usm_v(l)%nrs ) ) &
2888	surf_v(l)%nrs(mm(l)) = surf_usm_v(l)%nrs(m)
2889	IF ( ALLOCATED( surf_usm_v(l)%ol ) ) &
2890	surf_v(l)%ol(mm(l)) = surf_usm_v(l)%ol(m)
2891	IF ( ALLOCATED( surf_usm_v(l)%rib ) ) &
2892	surf_v(l)%rib(mm(l)) = surf_usm_v(l)%rib(m)
2893	IF ( ALLOCATED( surf_usm_v(l)%pt_surface ) ) &
2894	surf_v(l)%pt_surface(mm(l)) = surf_usm_v(l)%pt_surface(m)
2895	IF ( ALLOCATED( surf_usm_v(l)%usws ) ) &
2896	surf_v(l)%usws(mm(l)) = surf_usm_v(l)%usws(m)
2897	IF ( ALLOCATED( surf_usm_v(l)%vsws ) ) &
2898	surf_v(l)%vsws(mm(l)) = surf_usm_v(l)%vsws(m)
2899	IF ( ALLOCATED( surf_usm_v(l)%shf ) ) &
2900	surf_v(l)%shf(mm(l)) = surf_usm_v(l)%shf(m)
2901	IF ( ALLOCATED( surf_usm_v(l)%qsws ) ) &
2902	surf_v(l)%qsws(mm(l)) = surf_usm_v(l)%qsws(m)
2903	IF ( ALLOCATED( surf_usm_v(l)%ssws ) ) &
2904	surf_v(l)%ssws(mm(l)) = surf_usm_v(l)%ssws(m)
2905	IF ( ALLOCATED( surf_usm_v(l)%css ) ) THEN
2906	DO lsp = 1, nvar
2907	surf_v(l)%css(lsp,mm(l)) = surf_usm_v(l)%css(lsp,m)
2908	ENDDO
2909	ENDIF
2910	IF ( ALLOCATED( surf_usm_v(l)%cssws ) ) THEN
2911	DO lsp = 1, nvar
2912	surf_v(l)%cssws(lsp,mm(l)) = surf_usm_v(l)%cssws(lsp,m)
2913	ENDDO
2914	ENDIF
2915	IF ( ALLOCATED( surf_usm_v(l)%ncsws ) ) &
2916	surf_v(l)%ncsws(mm(l)) = surf_usm_v(l)%ncsws(m)
2917	IF ( ALLOCATED( surf_usm_v(l)%nrsws ) ) &
2918	surf_v(l)%nrsws(mm(l)) = surf_usm_v(l)%nrsws(m)
2919	IF ( ALLOCATED( surf_usm_v(l)%sasws ) ) &
2920	surf_v(l)%sasws(mm(l)) = surf_usm_v(l)%sasws(m)
2921	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_uv) ) &
2922	surf_v(l)%mom_flux_uv(mm(l)) = surf_usm_v(l)%mom_flux_uv(m)
2923	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_w) ) &
2924	surf_v(l)%mom_flux_w(mm(l)) = surf_usm_v(l)%mom_flux_w(m)
2925	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_tke) ) &
2926	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_usm_v(l)%mom_flux_tke(0:1,m)
2927
2928	mm(l) = mm(l) + 1
2929	ENDDO
2930
2931	ENDDO
2932	ENDDO
2933	!
2934	!-- Gather start- and end indices
2935	start_index_v(l) = 1
2936	DO i = nxl, nxr
2937	DO j = nys, nyn
2938
2939	surf_v(l)%start_index(j,i) = start_index_v(l)
2940	surf_v(l)%end_index(j,i) = surf_v(l)%start_index(j,i) -1
2941
2942	DO m = surf_def_v(l)%start_index(j,i), &
2943	surf_def_v(l)%end_index(j,i)
2944	surf_v(l)%end_index(j,i) = surf_v(l)%end_index(j,i) + 1
2945	ENDDO
2946	DO m = surf_lsm_v(l)%start_index(j,i), &
2947	surf_lsm_v(l)%end_index(j,i)
2948	surf_v(l)%end_index(j,i) = surf_v(l)%end_index(j,i) + 1
2949	ENDDO
2950	DO m = surf_usm_v(l)%start_index(j,i), &
2951	surf_usm_v(l)%end_index(j,i)
2952	surf_v(l)%end_index(j,i) = surf_v(l)%end_index(j,i) + 1
2953	ENDDO
2954
2955	start_index_v(l) = surf_v(l)%end_index(j,i) + 1
2956	ENDDO
2957	ENDDO
2958
2959	ENDDO
2960
2961
2962	CALL wrd_write_string( 'ns_h_on_file' )
2963	WRITE ( 14 ) ns_h_on_file
2964
2965	CALL wrd_write_string( 'ns_v_on_file' )
2966	WRITE ( 14 ) ns_v_on_file
2967
2968	!
2969	!-- Write required restart data.
2970	!-- Start with horizontal surfaces (upward-, downward-facing, and model top)
2971	DO l = 0, 2
2972	WRITE( dum, '(I1)') l
2973
2974	CALL wrd_write_string( 'surf_h(' // dum // ')%start_index' )
2975	WRITE ( 14 ) surf_h(l)%start_index
2976
2977	CALL wrd_write_string( 'surf_h(' // dum // ')%end_index' )
2978	WRITE ( 14 ) surf_h(l)%end_index
2979
2980	IF ( ALLOCATED ( surf_h(l)%us ) ) THEN
2981	CALL wrd_write_string( 'surf_h(' // dum // ')%us' )
2982	WRITE ( 14 ) surf_h(l)%us
2983	ENDIF
2984
2985	IF ( ALLOCATED ( surf_h(l)%ts ) ) THEN
2986	CALL wrd_write_string( 'surf_h(' // dum // ')%ts' )
2987	WRITE ( 14 ) surf_h(l)%ts
2988	ENDIF
2989
2990	IF ( ALLOCATED ( surf_h(l)%qs ) ) THEN
2991	CALL wrd_write_string( 'surf_h(' // dum // ')%qs' )
2992	WRITE ( 14 ) surf_h(l)%qs
2993	ENDIF
2994
2995	IF ( ALLOCATED ( surf_h(l)%ss ) ) THEN
2996	CALL wrd_write_string( 'surf_h(' // dum // ')%ss' )
2997	WRITE ( 14 ) surf_h(l)%ss
2998	ENDIF
2999
3000	IF ( ALLOCATED ( surf_h(l)%qcs ) ) THEN
3001	CALL wrd_write_string( 'surf_h(' // dum // ')%qcs' )
3002	WRITE ( 14 ) surf_h(l)%qcs
3003	ENDIF
3004
3005	IF ( ALLOCATED ( surf_h(l)%ncs ) ) THEN
3006	CALL wrd_write_string( 'surf_h(' // dum // ')%ncs' )
3007	WRITE ( 14 ) surf_h(l)%ncs
3008	ENDIF
3009
3010	IF ( ALLOCATED ( surf_h(l)%qrs ) ) THEN
3011	CALL wrd_write_string( 'surf_h(' // dum // ')%qrs' )
3012	WRITE ( 14 ) surf_h(l)%qrs
3013	ENDIF
3014
3015	IF ( ALLOCATED ( surf_h(l)%nrs ) ) THEN
3016	CALL wrd_write_string( 'surf_h(' // dum // ')%nrs' )
3017	WRITE ( 14 ) surf_h(l)%nrs
3018	ENDIF
3019
3020	IF ( ALLOCATED ( surf_h(l)%ol ) ) THEN
3021	CALL wrd_write_string( 'surf_h(' // dum // ')%ol' )
3022	WRITE ( 14 ) surf_h(l)%ol
3023	ENDIF
3024
3025	IF ( ALLOCATED ( surf_h(l)%rib ) ) THEN
3026	CALL wrd_write_string( 'surf_h(' // dum // ')%rib' )
3027	WRITE ( 14 ) surf_h(l)%rib
3028	ENDIF
3029
3030	IF ( ALLOCATED ( surf_h(l)%pt_surface ) ) THEN
3031	CALL wrd_write_string( 'surf_h(' // dum // ')%pt_surface' )
3032	WRITE ( 14 ) surf_h(l)%pt_surface
3033	ENDIF
3034
3035	IF ( ALLOCATED ( surf_h(l)%usws ) ) THEN
3036	CALL wrd_write_string( 'surf_h(' // dum // ')%usws' )
3037	WRITE ( 14 ) surf_h(l)%usws
3038	ENDIF
3039
3040	IF ( ALLOCATED ( surf_h(l)%vsws ) ) THEN
3041	CALL wrd_write_string( 'surf_h(' // dum // ')%vsws' )
3042	WRITE ( 14 ) surf_h(l)%vsws
3043	ENDIF
3044
3045	IF ( ALLOCATED ( surf_h(l)%shf ) ) THEN
3046	CALL wrd_write_string( 'surf_h(' // dum // ')%shf' )
3047	WRITE ( 14 ) surf_h(l)%shf
3048	ENDIF
3049
3050	IF ( ALLOCATED ( surf_h(l)%qsws ) ) THEN
3051	CALL wrd_write_string( 'surf_h(' // dum // ')%qsws' )
3052	WRITE ( 14 ) surf_h(l)%qsws
3053	ENDIF
3054
3055	IF ( ALLOCATED ( surf_h(l)%ssws ) ) THEN
3056	CALL wrd_write_string( 'surf_h(' // dum // ')%ssws' )
3057	WRITE ( 14 ) surf_h(l)%ssws
3058	ENDIF
3059
3060	IF ( ALLOCATED ( surf_h(l)%css ) ) THEN
3061	CALL wrd_write_string( 'surf_h(' // dum // ')%css' )
3062	WRITE ( 14 ) surf_h(l)%css
3063	ENDIF
3064
3065	IF ( ALLOCATED ( surf_h(l)%cssws ) ) THEN
3066	CALL wrd_write_string( 'surf_h(' // dum // ')%cssws' )
3067	WRITE ( 14 ) surf_h(l)%cssws
3068	ENDIF
3069
3070	IF ( ALLOCATED ( surf_h(l)%qcsws ) ) THEN
3071	CALL wrd_write_string( 'surf_h(' // dum // ')%qcsws' )
3072	WRITE ( 14 ) surf_h(l)%qcsws
3073	ENDIF
3074
3075	IF ( ALLOCATED ( surf_h(l)%ncsws ) ) THEN
3076	CALL wrd_write_string( 'surf_h(' // dum // ')%ncsws' )
3077	WRITE ( 14 ) surf_h(l)%ncsws
3078	ENDIF
3079
3080	IF ( ALLOCATED ( surf_h(l)%qrsws ) ) THEN
3081	CALL wrd_write_string( 'surf_h(' // dum // ')%qrsws' )
3082	WRITE ( 14 ) surf_h(l)%qrsws
3083	ENDIF
3084
3085	IF ( ALLOCATED ( surf_h(l)%nrsws ) ) THEN
3086	CALL wrd_write_string( 'surf_h(' // dum // ')%nrsws' )
3087	WRITE ( 14 ) surf_h(l)%nrsws
3088	ENDIF
3089
3090	IF ( ALLOCATED ( surf_h(l)%sasws ) ) THEN
3091	CALL wrd_write_string( 'surf_h(' // dum // ')%sasws' )
3092	WRITE ( 14 ) surf_h(l)%sasws
3093	ENDIF
3094
3095	ENDDO
3096	!
3097	!-- Write vertical surfaces
3098	DO l = 0, 3
3099	WRITE( dum, '(I1)') l
3100
3101	CALL wrd_write_string( 'surf_v(' // dum // ')%start_index' )
3102	WRITE ( 14 ) surf_v(l)%start_index
3103
3104	CALL wrd_write_string( 'surf_v(' // dum // ')%end_index' )
3105	WRITE ( 14 ) surf_v(l)%end_index
3106
3107	IF ( ALLOCATED ( surf_v(l)%us ) ) THEN
3108	CALL wrd_write_string( 'surf_v(' // dum // ')%us' )
3109	WRITE ( 14 ) surf_v(l)%us
3110	ENDIF
3111
3112	IF ( ALLOCATED ( surf_v(l)%ts ) ) THEN
3113	CALL wrd_write_string( 'surf_v(' // dum // ')%ts' )
3114	WRITE ( 14 ) surf_v(l)%ts
3115	ENDIF
3116
3117	IF ( ALLOCATED ( surf_v(l)%qs ) ) THEN
3118	CALL wrd_write_string( 'surf_v(' // dum // ')%qs' )
3119	WRITE ( 14 ) surf_v(l)%qs
3120	ENDIF
3121
3122	IF ( ALLOCATED ( surf_v(l)%ss ) ) THEN
3123	CALL wrd_write_string( 'surf_v(' // dum // ')%ss' )
3124	WRITE ( 14 ) surf_v(l)%ss
3125	ENDIF
3126
3127	IF ( ALLOCATED ( surf_v(l)%qcs ) ) THEN
3128	CALL wrd_write_string( 'surf_v(' // dum // ')%qcs' )
3129	WRITE ( 14 ) surf_v(l)%qcs
3130	ENDIF
3131
3132	IF ( ALLOCATED ( surf_v(l)%ncs ) ) THEN
3133	CALL wrd_write_string( 'surf_v(' // dum // ')%ncs' )
3134	WRITE ( 14 ) surf_v(l)%ncs
3135	ENDIF
3136
3137	IF ( ALLOCATED ( surf_v(l)%qrs ) ) THEN
3138	CALL wrd_write_string( 'surf_v(' // dum // ')%qrs' )
3139	WRITE ( 14 ) surf_v(l)%qrs
3140	ENDIF
3141
3142	IF ( ALLOCATED ( surf_v(l)%nrs ) ) THEN
3143	CALL wrd_write_string( 'surf_v(' // dum // ')%nrs' )
3144	WRITE ( 14 ) surf_v(l)%nrs
3145	ENDIF
3146
3147	IF ( ALLOCATED ( surf_v(l)%ol ) ) THEN
3148	CALL wrd_write_string( 'surf_v(' // dum // ')%ol' )
3149	WRITE ( 14 ) surf_v(l)%ol
3150	ENDIF
3151
3152	IF ( ALLOCATED ( surf_v(l)%rib ) ) THEN
3153	CALL wrd_write_string( 'surf_v(' // dum // ')%rib' )
3154	WRITE ( 14 ) surf_v(l)%rib
3155	ENDIF
3156
3157	IF ( ALLOCATED ( surf_v(l)%pt_surface ) ) THEN
3158	CALL wrd_write_string( 'surf_v(' // dum // ')%pt_surface' )
3159	WRITE ( 14 ) surf_v(l)%pt_surface
3160	ENDIF
3161
3162	IF ( ALLOCATED ( surf_v(l)%shf ) ) THEN
3163	CALL wrd_write_string( 'surf_v(' // dum // ')%shf' )
3164	WRITE ( 14 ) surf_v(l)%shf
3165	ENDIF
3166
3167	IF ( ALLOCATED ( surf_v(l)%qsws ) ) THEN
3168	CALL wrd_write_string( 'surf_v(' // dum // ')%qsws' )
3169	WRITE ( 14 ) surf_v(l)%qsws
3170	ENDIF
3171
3172	IF ( ALLOCATED ( surf_v(l)%ssws ) ) THEN
3173	CALL wrd_write_string( 'surf_v(' // dum // ')%ssws' )
3174	WRITE ( 14 ) surf_v(l)%ssws
3175	ENDIF
3176
3177	IF ( ALLOCATED ( surf_v(l)%css ) ) THEN
3178	CALL wrd_write_string( 'surf_v(' // dum // ')%css' )
3179	WRITE ( 14 ) surf_v(l)%css
3180	ENDIF
3181
3182	IF ( ALLOCATED ( surf_v(l)%cssws ) ) THEN
3183	CALL wrd_write_string( 'surf_v(' // dum // ')%cssws' )
3184	WRITE ( 14 ) surf_v(l)%cssws
3185	ENDIF
3186
3187	IF ( ALLOCATED ( surf_v(l)%qcsws ) ) THEN
3188	CALL wrd_write_string( 'surf_v(' // dum // ')%qcsws' )
3189	WRITE ( 14 ) surf_v(l)%qcsws
3190	ENDIF
3191
3192	IF ( ALLOCATED ( surf_v(l)%ncsws ) ) THEN
3193	CALL wrd_write_string( 'surf_v(' // dum // ')%ncsws' )
3194	WRITE ( 14 ) surf_v(l)%ncsws
3195	ENDIF
3196
3197	IF ( ALLOCATED ( surf_v(l)%qrsws ) ) THEN
3198	CALL wrd_write_string( 'surf_v(' // dum // ')%qrsws' )
3199	WRITE ( 14 ) surf_v(l)%qrsws
3200	ENDIF
3201
3202	IF ( ALLOCATED ( surf_v(l)%nrsws ) ) THEN
3203	CALL wrd_write_string( 'surf_v(' // dum // ')%nrsws' )
3204	WRITE ( 14 ) surf_v(l)%nrsws
3205	ENDIF
3206
3207	IF ( ALLOCATED ( surf_v(l)%sasws ) ) THEN
3208	CALL wrd_write_string( 'surf_v(' // dum // ')%sasws' )
3209	WRITE ( 14 ) surf_v(l)%sasws
3210	ENDIF
3211
3212	IF ( ALLOCATED ( surf_v(l)%mom_flux_uv ) ) THEN
3213	CALL wrd_write_string( 'surf_v(' // dum // ')%mom_uv' )
3214	WRITE ( 14 ) surf_v(l)%mom_flux_uv
3215	ENDIF
3216
3217	IF ( ALLOCATED ( surf_v(l)%mom_flux_w ) ) THEN
3218	CALL wrd_write_string( 'surf_v(' // dum // ')%mom_w' )
3219	WRITE ( 14 ) surf_v(l)%mom_flux_w
3220	ENDIF
3221
3222	IF ( ALLOCATED ( surf_v(l)%mom_flux_tke ) ) THEN
3223	CALL wrd_write_string( 'surf_v(' // dum // ')%mom_tke' )
3224	WRITE ( 14 ) surf_v(l)%mom_flux_tke
3225	ENDIF
3226
3227	ENDDO
3228
3229
3230	END SUBROUTINE surface_wrd_local
3231
3232
3233	!------------------------------------------------------------------------------!
3234	! Description:
3235	! ------------
3236	!> Reads surface-related restart data. Please note, restart data for a certain
3237	!> surface orientation (e.g. horizontal upward-facing) is stored in one
3238	!> array, even if surface elements may belong to different surface types
3239	!> natural or urban for example). Surface elements are redistributed into its
3240	!> respective surface types within this routine. This allows e.g. changing the
3241	!> surface type after reading the restart data, which might be required in case
3242	!> of cyclic_fill mode.
3243	!------------------------------------------------------------------------------!
3244	SUBROUTINE surface_rrd_local( ii, kk, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
3245	nxr_on_file, nynf, nync, nyn_on_file, nysf, &
3246	nysc, nys_on_file, found )
3247
3248	USE pegrid, &
3249	ONLY: myid, numprocs_previous_run
3250
3251	IMPLICIT NONE
3252
3253	CHARACTER (LEN=1) :: dum !< dummy to create correct string for reading input variable
3254
3255	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
3256	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
3257	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
3258	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
3259	INTEGER(iwp) :: k !< running index along z-direction
3260	INTEGER(iwp) :: m !< running index for surface elements, refers to gathered array encompassing all surface types
3261	INTEGER(iwp) :: mm !< running index for surface elements, refers to individual surface types
3262	INTEGER(iwp) :: ii !< running index over input files
3263	INTEGER(iwp) :: kk !< running index over previous input files covering current local domain
3264	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
3265	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
3266	INTEGER(iwp) :: nxl_on_file !< index of left boundary on former local domain
3267	INTEGER(iwp) :: nxrc !< index of right boundary on current subdomain
3268	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
3269	INTEGER(iwp) :: nxr_on_file !< index of right boundary on former local domain
3270	INTEGER(iwp) :: nync !< index of north boundary on current subdomain
3271	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
3272	INTEGER(iwp) :: nyn_on_file !< index of norht boundary on former local domain
3273	INTEGER(iwp) :: nysc !< index of south boundary on current subdomain
3274	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
3275	INTEGER(iwp) :: nys_on_file !< index of south boundary on former local domain
3276
3277	INTEGER(iwp), SAVE :: l !< index variable for surface type
3278
3279	LOGICAL :: surf_match_def !< flag indicating that surface element is of default type
3280	LOGICAL :: surf_match_lsm !< flag indicating that surface element is of natural type
3281	LOGICAL :: surf_match_usm !< flag indicating that surface element is of urban type
3282
3283	LOGICAL, INTENT(OUT) :: found
3284
3285	LOGICAL, SAVE :: horizontal_surface !< flag indicating horizontal surfaces
3286	LOGICAL, SAVE :: vertical_surface !< flag indicating vertical surfaces
3287
3288	TYPE(surf_type), DIMENSION(0:2), SAVE :: surf_h !< horizontal surface type on file
3289	TYPE(surf_type), DIMENSION(0:3), SAVE :: surf_v !< vertical surface type on file
3290
3291
3292	found = .TRUE.
3293
3294
3295	SELECT CASE ( restart_string(1:length) )
3296
3297	CASE ( 'ns_h_on_file' )
3298	IF ( kk == 1 ) THEN
3299	READ ( 13 ) ns_h_on_file
3300
3301	IF ( ALLOCATED( surf_h(0)%start_index ) ) &
3302	CALL deallocate_surface_attributes_h( surf_h(0) )
3303	IF ( ALLOCATED( surf_h(1)%start_index ) ) &
3304	CALL deallocate_surface_attributes_h( surf_h(1) )
3305	IF ( ALLOCATED( surf_h(2)%start_index ) ) &
3306	CALL deallocate_surface_attributes_h_top( surf_h(2) )
3307
3308	!-- Allocate memory for number of surface elements on file.
3309	!-- Please note, these number is not necessarily the same as
3310	!-- the final number of surface elements on local domain,
3311	!-- which is the case if processor topology changes during
3312	!-- restart runs.
3313	!-- Horizontal upward facing
3314	surf_h(0)%ns = ns_h_on_file(0)
3315	CALL allocate_surface_attributes_h( surf_h(0), &
3316	nys_on_file, nyn_on_file, &
3317	nxl_on_file, nxr_on_file )
3318
3319	!-- Horizontal downward facing
3320	surf_h(1)%ns = ns_h_on_file(1)
3321	CALL allocate_surface_attributes_h( surf_h(1), &
3322	nys_on_file, nyn_on_file, &
3323	nxl_on_file, nxr_on_file )
3324
3325	!-- Model top
3326	surf_h(2)%ns = ns_h_on_file(2)
3327	CALL allocate_surface_attributes_h_top( surf_h(2), &
3328	nys_on_file, nyn_on_file, &
3329	nxl_on_file, nxr_on_file )
3330
3331	!
3332	!-- Initial setting of flags for horizontal and vertical surfaces,
3333	!-- will be set after start- and end-indices are read.
3334	horizontal_surface = .FALSE.
3335	vertical_surface = .FALSE.
3336
3337	ENDIF
3338
3339	CASE ( 'ns_v_on_file' )
3340	IF ( kk == 1 ) THEN
3341	READ ( 13 ) ns_v_on_file
3342
3343	DO l = 0, 3
3344	IF ( ALLOCATED( surf_v(l)%start_index ) ) &
3345	CALL deallocate_surface_attributes_v( surf_v(l) )
3346	ENDDO
3347
3348	!-- Vertical surfaces
3349	DO l = 0, 3
3350	surf_v(l)%ns = ns_v_on_file(l)
3351	CALL allocate_surface_attributes_v( surf_v(l), &
3352	nys_on_file, nyn_on_file, &
3353	nxl_on_file, nxr_on_file )
3354	ENDDO
3355
3356	ENDIF
3357
3358	CASE ( 'surf_h(0)%start_index' )
3359	IF ( kk == 1 ) &
3360	READ ( 13 ) surf_h(0)%start_index
3361	l = 0
3362	CASE ( 'surf_h(0)%end_index' )
3363	IF ( kk == 1 ) &
3364	READ ( 13 ) surf_h(0)%end_index
3365	horizontal_surface = .TRUE.
3366	vertical_surface = .FALSE.
3367	CASE ( 'surf_h(0)%us' )
3368	IF ( ALLOCATED( surf_h(0)%us ) .AND. kk == 1 ) &
3369	READ ( 13 ) surf_h(0)%us
3370	CASE ( 'surf_h(0)%ts' )
3371	IF ( ALLOCATED( surf_h(0)%ts ) .AND. kk == 1 ) &
3372	READ ( 13 ) surf_h(0)%ts
3373	CASE ( 'surf_h(0)%qs' )
3374	IF ( ALLOCATED( surf_h(0)%qs ) .AND. kk == 1 ) &
3375	READ ( 13 ) surf_h(0)%qs
3376	CASE ( 'surf_h(0)%ss' )
3377	IF ( ALLOCATED( surf_h(0)%ss ) .AND. kk == 1 ) &
3378	READ ( 13 ) surf_h(0)%ss
3379	CASE ( 'surf_h(0)%qcs' )
3380	IF ( ALLOCATED( surf_h(0)%qcs ) .AND. kk == 1 ) &
3381	READ ( 13 ) surf_h(0)%qcs
3382	CASE ( 'surf_h(0)%ncs' )
3383	IF ( ALLOCATED( surf_h(0)%ncs ) .AND. kk == 1 ) &
3384	READ ( 13 ) surf_h(0)%ncs
3385	CASE ( 'surf_h(0)%qrs' )
3386	IF ( ALLOCATED( surf_h(0)%qrs ) .AND. kk == 1 ) &
3387	READ ( 13 ) surf_h(0)%qrs
3388	CASE ( 'surf_h(0)%nrs' )
3389	IF ( ALLOCATED( surf_h(0)%nrs ) .AND. kk == 1 ) &
3390	READ ( 13 ) surf_h(0)%nrs
3391	CASE ( 'surf_h(0)%ol' )
3392	IF ( ALLOCATED( surf_h(0)%ol ) .AND. kk == 1 ) &
3393	READ ( 13 ) surf_h(0)%ol
3394	CASE ( 'surf_h(0)%rib' )
3395	IF ( ALLOCATED( surf_h(0)%rib ) .AND. kk == 1 ) &
3396	READ ( 13 ) surf_h(0)%rib
3397	CASE ( 'surf_h(0)%pt_surface' )
3398	IF ( ALLOCATED( surf_h(0)%pt_surface ) .AND. kk == 1 ) &
3399	READ ( 13 ) surf_h(0)%pt_surface
3400	CASE ( 'surf_h(0)%usws' )
3401	IF ( ALLOCATED( surf_h(0)%usws ) .AND. kk == 1 ) &
3402	READ ( 13 ) surf_h(0)%usws
3403	CASE ( 'surf_h(0)%vsws' )
3404	IF ( ALLOCATED( surf_h(0)%vsws ) .AND. kk == 1 ) &
3405	READ ( 13 ) surf_h(0)%vsws
3406	CASE ( 'surf_h(0)%shf' )
3407	IF ( ALLOCATED( surf_h(0)%shf ) .AND. kk == 1 ) &
3408	READ ( 13 ) surf_h(0)%shf
3409	CASE ( 'surf_h(0)%qsws' )
3410	IF ( ALLOCATED( surf_h(0)%qsws ) .AND. kk == 1 ) &
3411	READ ( 13 ) surf_h(0)%qsws
3412	CASE ( 'surf_h(0)%ssws' )
3413	IF ( ALLOCATED( surf_h(0)%ssws ) .AND. kk == 1 ) &
3414	READ ( 13 ) surf_h(0)%ssws
3415	CASE ( 'surf_h(0)%css' )
3416	IF ( ALLOCATED( surf_h(0)%css ) .AND. kk == 1 ) &
3417	READ ( 13 ) surf_h(0)%css
3418	CASE ( 'surf_h(0)%cssws' )
3419	IF ( ALLOCATED( surf_h(0)%cssws ) .AND. kk == 1 ) &
3420	READ ( 13 ) surf_h(0)%cssws
3421	CASE ( 'surf_h(0)%qcsws' )
3422	IF ( ALLOCATED( surf_h(0)%qcsws ) .AND. kk == 1 ) &
3423	READ ( 13 ) surf_h(0)%qcsws
3424	CASE ( 'surf_h(0)%ncsws' )
3425	IF ( ALLOCATED( surf_h(0)%ncsws ) .AND. kk == 1 ) &
3426	READ ( 13 ) surf_h(0)%ncsws
3427	CASE ( 'surf_h(0)%qrsws' )
3428	IF ( ALLOCATED( surf_h(0)%qrsws ) .AND. kk == 1 ) &
3429	READ ( 13 ) surf_h(0)%qrsws
3430	CASE ( 'surf_h(0)%nrsws' )
3431	IF ( ALLOCATED( surf_h(0)%nrsws ) .AND. kk == 1 ) &
3432	READ ( 13 ) surf_h(0)%nrsws
3433	CASE ( 'surf_h(0)%sasws' )
3434	IF ( ALLOCATED( surf_h(0)%sasws ) .AND. kk == 1 ) &
3435	READ ( 13 ) surf_h(0)%sasws
3436
3437	CASE ( 'surf_h(1)%start_index' )
3438	IF ( kk == 1 ) &
3439	READ ( 13 ) surf_h(1)%start_index
3440	l = 1
3441	CASE ( 'surf_h(1)%end_index' )
3442	IF ( kk == 1 ) &
3443	READ ( 13 ) surf_h(1)%end_index
3444	CASE ( 'surf_h(1)%us' )
3445	IF ( ALLOCATED( surf_h(1)%us ) .AND. kk == 1 ) &
3446	READ ( 13 ) surf_h(1)%us
3447	CASE ( 'surf_h(1)%ts' )
3448	IF ( ALLOCATED( surf_h(1)%ts ) .AND. kk == 1 ) &
3449	READ ( 13 ) surf_h(1)%ts
3450	CASE ( 'surf_h(1)%qs' )
3451	IF ( ALLOCATED( surf_h(1)%qs ) .AND. kk == 1 ) &
3452	READ ( 13 ) surf_h(1)%qs
3453	CASE ( 'surf_h(1)%ss' )
3454	IF ( ALLOCATED( surf_h(1)%ss ) .AND. kk == 1 ) &
3455	READ ( 13 ) surf_h(1)%ss
3456	CASE ( 'surf_h(1)%qcs' )
3457	IF ( ALLOCATED( surf_h(1)%qcs ) .AND. kk == 1 ) &
3458	READ ( 13 ) surf_h(1)%qcs
3459	CASE ( 'surf_h(1)%ncs' )
3460	IF ( ALLOCATED( surf_h(1)%ncs ) .AND. kk == 1 ) &
3461	READ ( 13 ) surf_h(1)%ncs
3462	CASE ( 'surf_h(1)%qrs' )
3463	IF ( ALLOCATED( surf_h(1)%qrs ) .AND. kk == 1 ) &
3464	READ ( 13 ) surf_h(1)%qrs
3465	CASE ( 'surf_h(1)%nrs' )
3466	IF ( ALLOCATED( surf_h(1)%nrs ) .AND. kk == 1 ) &
3467	READ ( 13 ) surf_h(1)%nrs
3468	CASE ( 'surf_h(1)%ol' )
3469	IF ( ALLOCATED( surf_h(1)%ol ) .AND. kk == 1 ) &
3470	READ ( 13 ) surf_h(1)%ol
3471	CASE ( 'surf_h(1)%rib' )
3472	IF ( ALLOCATED( surf_h(1)%rib ) .AND. kk == 1 ) &
3473	READ ( 13 ) surf_h(1)%rib
3474	CASE ( 'surf_h(1)%pt_surface' )
3475	IF ( ALLOCATED( surf_h(1)%pt_surface ) .AND. kk == 1 ) &
3476	READ ( 13 ) surf_h(1)%pt_surface
3477	CASE ( 'surf_h(1)%usws' )
3478	IF ( ALLOCATED( surf_h(1)%usws ) .AND. kk == 1 ) &
3479	READ ( 13 ) surf_h(1)%usws
3480	CASE ( 'surf_h(1)%vsws' )
3481	IF ( ALLOCATED( surf_h(1)%vsws ) .AND. kk == 1 ) &
3482	READ ( 13 ) surf_h(1)%vsws
3483	CASE ( 'surf_h(1)%shf' )
3484	IF ( ALLOCATED( surf_h(1)%shf ) .AND. kk == 1 ) &
3485	READ ( 13 ) surf_h(1)%shf
3486	CASE ( 'surf_h(1)%qsws' )
3487	IF ( ALLOCATED( surf_h(1)%qsws ) .AND. kk == 1 ) &
3488	READ ( 13 ) surf_h(1)%qsws
3489	CASE ( 'surf_h(1)%ssws' )
3490	IF ( ALLOCATED( surf_h(1)%ssws ) .AND. kk == 1 ) &
3491	READ ( 13 ) surf_h(1)%ssws
3492	CASE ( 'surf_h(1)%css' )
3493	IF ( ALLOCATED( surf_h(1)%css ) .AND. kk == 1 ) &
3494	READ ( 13 ) surf_h(1)%css
3495	CASE ( 'surf_h(1)%cssws' )
3496	IF ( ALLOCATED( surf_h(1)%cssws ) .AND. kk == 1 ) &
3497	READ ( 13 ) surf_h(1)%cssws
3498	CASE ( 'surf_h(1)%qcsws' )
3499	IF ( ALLOCATED( surf_h(1)%qcsws ) .AND. kk == 1 ) &
3500	READ ( 13 ) surf_h(1)%qcsws
3501	CASE ( 'surf_h(1)%ncsws' )
3502	IF ( ALLOCATED( surf_h(1)%ncsws ) .AND. kk == 1 ) &
3503	READ ( 13 ) surf_h(1)%ncsws
3504	CASE ( 'surf_h(1)%qrsws' )
3505	IF ( ALLOCATED( surf_h(1)%qrsws ) .AND. kk == 1 ) &
3506	READ ( 13 ) surf_h(1)%qrsws
3507	CASE ( 'surf_h(1)%nrsws' )
3508	IF ( ALLOCATED( surf_h(1)%nrsws ) .AND. kk == 1 ) &
3509	READ ( 13 ) surf_h(1)%nrsws
3510	CASE ( 'surf_h(1)%sasws' )
3511	IF ( ALLOCATED( surf_h(1)%sasws ) .AND. kk == 1 ) &
3512	READ ( 13 ) surf_h(1)%sasws
3513
3514	CASE ( 'surf_h(2)%start_index' )
3515	IF ( kk == 1 ) &
3516	READ ( 13 ) surf_h(2)%start_index
3517	l = 2
3518	CASE ( 'surf_h(2)%end_index' )
3519	IF ( kk == 1 ) &
3520	READ ( 13 ) surf_h(2)%end_index
3521	CASE ( 'surf_h(2)%us' )
3522	IF ( ALLOCATED( surf_h(2)%us ) .AND. kk == 1 ) &
3523	READ ( 13 ) surf_h(2)%us
3524	CASE ( 'surf_h(2)%ts' )
3525	IF ( ALLOCATED( surf_h(2)%ts ) .AND. kk == 1 ) &
3526	READ ( 13 ) surf_h(2)%ts
3527	CASE ( 'surf_h(2)%qs' )
3528	IF ( ALLOCATED( surf_h(2)%qs ) .AND. kk == 1 ) &
3529	READ ( 13 ) surf_h(2)%qs
3530	CASE ( 'surf_h(2)%ss' )
3531	IF ( ALLOCATED( surf_h(2)%ss ) .AND. kk == 1 ) &
3532	READ ( 13 ) surf_h(2)%ss
3533	CASE ( 'surf_h(2)%qcs' )
3534	IF ( ALLOCATED( surf_h(2)%qcs ) .AND. kk == 1 ) &
3535	READ ( 13 ) surf_h(2)%qcs
3536	CASE ( 'surf_h(2)%ncs' )
3537	IF ( ALLOCATED( surf_h(2)%ncs ) .AND. kk == 1 ) &
3538	READ ( 13 ) surf_h(2)%ncs
3539	CASE ( 'surf_h(2)%qrs' )
3540	IF ( ALLOCATED( surf_h(2)%qrs ) .AND. kk == 1 ) &
3541	READ ( 13 ) surf_h(2)%qrs
3542	CASE ( 'surf_h(2)%nrs' )
3543	IF ( ALLOCATED( surf_h(2)%nrs ) .AND. kk == 1 ) &
3544	READ ( 13 ) surf_h(2)%nrs
3545	CASE ( 'surf_h(2)%ol' )
3546	IF ( ALLOCATED( surf_h(2)%ol ) .AND. kk == 1 ) &
3547	READ ( 13 ) surf_h(2)%ol
3548	CASE ( 'surf_h(2)%rib' )
3549	IF ( ALLOCATED( surf_h(2)%rib ) .AND. kk == 1 ) &
3550	READ ( 13 ) surf_h(2)%rib
3551	CASE ( 'surf_h(2)%pt_surface' )
3552	IF ( ALLOCATED( surf_h(2)%pt_surface ) .AND. kk == 1 ) &
3553	READ ( 13 ) surf_h(2)%pt_surface
3554	CASE ( 'surf_h(2)%usws' )
3555	IF ( ALLOCATED( surf_h(2)%usws ) .AND. kk == 1 ) &
3556	READ ( 13 ) surf_h(2)%usws
3557	CASE ( 'surf_h(2)%vsws' )
3558	IF ( ALLOCATED( surf_h(2)%vsws ) .AND. kk == 1 ) &
3559	READ ( 13 ) surf_h(2)%vsws
3560	CASE ( 'surf_h(2)%shf' )
3561	IF ( ALLOCATED( surf_h(2)%shf ) .AND. kk == 1 ) &
3562	READ ( 13 ) surf_h(2)%shf
3563	CASE ( 'surf_h(2)%qsws' )
3564	IF ( ALLOCATED( surf_h(2)%qsws ) .AND. kk == 1 ) &
3565	READ ( 13 ) surf_h(2)%qsws
3566	CASE ( 'surf_h(2)%ssws' )
3567	IF ( ALLOCATED( surf_h(2)%ssws ) .AND. kk == 1 ) &
3568	READ ( 13 ) surf_h(2)%ssws
3569	CASE ( 'surf_h(2)%css' )
3570	IF ( ALLOCATED( surf_h(2)%css ) .AND. kk == 1 ) &
3571	READ ( 13 ) surf_h(2)%css
3572	CASE ( 'surf_h(2)%cssws' )
3573	IF ( ALLOCATED( surf_h(2)%cssws ) .AND. kk == 1 ) &
3574	READ ( 13 ) surf_h(2)%cssws
3575	CASE ( 'surf_h(2)%qcsws' )
3576	IF ( ALLOCATED( surf_h(2)%qcsws ) .AND. kk == 1 ) &
3577	READ ( 13 ) surf_h(2)%qcsws
3578	CASE ( 'surf_h(2)%ncsws' )
3579	IF ( ALLOCATED( surf_h(2)%ncsws ) .AND. kk == 1 ) &
3580	READ ( 13 ) surf_h(2)%ncsws
3581	CASE ( 'surf_h(2)%qrsws' )
3582	IF ( ALLOCATED( surf_h(2)%qrsws ) .AND. kk == 1 ) &
3583	READ ( 13 ) surf_h(2)%qrsws
3584	CASE ( 'surf_h(2)%nrsws' )
3585	IF ( ALLOCATED( surf_h(2)%nrsws ) .AND. kk == 1 ) &
3586	READ ( 13 ) surf_h(2)%nrsws
3587	CASE ( 'surf_h(2)%sasws' )
3588	IF ( ALLOCATED( surf_h(2)%sasws ) .AND. kk == 1 ) &
3589	READ ( 13 ) surf_h(2)%sasws
3590
3591	CASE ( 'surf_v(0)%start_index' )
3592	IF ( kk == 1 ) &
3593	READ ( 13 ) surf_v(0)%start_index
3594	l = 0
3595	horizontal_surface = .FALSE.
3596	vertical_surface = .TRUE.
3597	CASE ( 'surf_v(0)%end_index' )
3598	IF ( kk == 1 ) &
3599	READ ( 13 ) surf_v(0)%end_index
3600	CASE ( 'surf_v(0)%us' )
3601	IF ( ALLOCATED( surf_v(0)%us ) .AND. kk == 1 ) &
3602	READ ( 13 ) surf_v(0)%us
3603	CASE ( 'surf_v(0)%ts' )
3604	IF ( ALLOCATED( surf_v(0)%ts ) .AND. kk == 1 ) &
3605	READ ( 13 ) surf_v(0)%ts
3606	CASE ( 'surf_v(0)%qs' )
3607	IF ( ALLOCATED( surf_v(0)%qs ) .AND. kk == 1 ) &
3608	READ ( 13 ) surf_v(0)%qs
3609	CASE ( 'surf_v(0)%ss' )
3610	IF ( ALLOCATED( surf_v(0)%ss ) .AND. kk == 1 ) &
3611	READ ( 13 ) surf_v(0)%ss
3612	CASE ( 'surf_v(0)%qcs' )
3613	IF ( ALLOCATED( surf_v(0)%qcs ) .AND. kk == 1 ) &
3614	READ ( 13 ) surf_v(0)%qcs
3615	CASE ( 'surf_v(0)%ncs' )
3616	IF ( ALLOCATED( surf_v(0)%ncs ) .AND. kk == 1 ) &
3617	READ ( 13 ) surf_v(0)%ncs
3618	CASE ( 'surf_v(0)%qrs' )
3619	IF ( ALLOCATED( surf_v(0)%qrs ) .AND. kk == 1 ) &
3620	READ ( 13 ) surf_v(0)%qrs
3621	CASE ( 'surf_v(0)%nrs' )
3622	IF ( ALLOCATED( surf_v(0)%nrs ) .AND. kk == 1 ) &
3623	READ ( 13 ) surf_v(0)%nrs
3624	CASE ( 'surf_v(0)%ol' )
3625	IF ( ALLOCATED( surf_v(0)%ol ) .AND. kk == 1 ) &
3626	READ ( 13 ) surf_v(0)%ol
3627	CASE ( 'surf_v(0)%rib' )
3628	IF ( ALLOCATED( surf_v(0)%rib ) .AND. kk == 1 ) &
3629	READ ( 13 ) surf_v(0)%rib
3630	CASE ( 'surf_v(0)%pt_surface' )
3631	IF ( ALLOCATED( surf_v(0)%pt_surface ) .AND. kk == 1 ) &
3632	READ ( 13 ) surf_v(0)%pt_surface
3633	CASE ( 'surf_v(0)%shf' )
3634	IF ( ALLOCATED( surf_v(0)%shf ) .AND. kk == 1 ) &
3635	READ ( 13 ) surf_v(0)%shf
3636	CASE ( 'surf_v(0)%qsws' )
3637	IF ( ALLOCATED( surf_v(0)%qsws ) .AND. kk == 1 ) &
3638	READ ( 13 ) surf_v(0)%qsws
3639	CASE ( 'surf_v(0)%ssws' )
3640	IF ( ALLOCATED( surf_v(0)%ssws ) .AND. kk == 1 ) &
3641	READ ( 13 ) surf_v(0)%ssws
3642	CASE ( 'surf_v(0)%css' )
3643	IF ( ALLOCATED( surf_v(0)%css ) .AND. kk == 1 ) &
3644	READ ( 13 ) surf_v(0)%css
3645	CASE ( 'surf_v(0)%cssws' )
3646	IF ( ALLOCATED( surf_v(0)%cssws ) .AND. kk == 1 ) &
3647	READ ( 13 ) surf_v(0)%cssws
3648	CASE ( 'surf_v(0)%qcsws' )
3649	IF ( ALLOCATED( surf_v(0)%qcsws ) .AND. kk == 1 ) &
3650	READ ( 13 ) surf_v(0)%qcsws
3651	CASE ( 'surf_v(0)%ncsws' )
3652	IF ( ALLOCATED( surf_v(0)%ncsws ) .AND. kk == 1 ) &
3653	READ ( 13 ) surf_v(0)%ncsws
3654	CASE ( 'surf_v(0)%qrsws' )
3655	IF ( ALLOCATED( surf_v(0)%qrsws ) .AND. kk == 1 ) &
3656	READ ( 13 ) surf_v(0)%qrsws
3657	CASE ( 'surf_v(0)%nrsws' )
3658	IF ( ALLOCATED( surf_v(0)%nrsws ) .AND. kk == 1 ) &
3659	READ ( 13 ) surf_v(0)%nrsws
3660	CASE ( 'surf_v(0)%sasws' )
3661	IF ( ALLOCATED( surf_v(0)%sasws ) .AND. kk == 1 ) &
3662	READ ( 13 ) surf_v(0)%sasws
3663	CASE ( 'surf_v(0)%mom_uv' )
3664	IF ( ALLOCATED( surf_v(0)%mom_flux_uv ) .AND. kk == 1 ) &
3665	READ ( 13 ) surf_v(0)%mom_flux_uv
3666	CASE ( 'surf_v(0)%mom_w' )
3667	IF ( ALLOCATED( surf_v(0)%mom_flux_w ) .AND. kk == 1 ) &
3668	READ ( 13 ) surf_v(0)%mom_flux_w
3669	CASE ( 'surf_v(0)%mom_tke' )
3670	IF ( ALLOCATED( surf_v(0)%mom_flux_tke ) .AND. kk == 1 ) &
3671	READ ( 13 ) surf_v(0)%mom_flux_tke
3672
3673	CASE ( 'surf_v(1)%start_index' )
3674	IF ( kk == 1 ) &
3675	READ ( 13 ) surf_v(1)%start_index
3676	l = 1
3677	CASE ( 'surf_v(1)%end_index' )
3678	IF ( kk == 1 ) &
3679	READ ( 13 ) surf_v(1)%end_index
3680	CASE ( 'surf_v(1)%us' )
3681	IF ( ALLOCATED( surf_v(1)%us ) .AND. kk == 1 ) &
3682	READ ( 13 ) surf_v(1)%us
3683	CASE ( 'surf_v(1)%ts' )
3684	IF ( ALLOCATED( surf_v(1)%ts ) .AND. kk == 1 ) &
3685	READ ( 13 ) surf_v(1)%ts
3686	CASE ( 'surf_v(1)%qs' )
3687	IF ( ALLOCATED( surf_v(1)%qs ) .AND. kk == 1 ) &
3688	READ ( 13 ) surf_v(1)%qs
3689	CASE ( 'surf_v(1)%ss' )
3690	IF ( ALLOCATED( surf_v(1)%ss ) .AND. kk == 1 ) &
3691	READ ( 13 ) surf_v(1)%ss
3692	CASE ( 'surf_v(1)%qcs' )
3693	IF ( ALLOCATED( surf_v(1)%qcs ) .AND. kk == 1 ) &
3694	READ ( 13 ) surf_v(1)%qcs
3695	CASE ( 'surf_v(1)%ncs' )
3696	IF ( ALLOCATED( surf_v(1)%ncs ) .AND. kk == 1 ) &
3697	READ ( 13 ) surf_v(1)%ncs
3698	CASE ( 'surf_v(1)%qrs' )
3699	IF ( ALLOCATED( surf_v(1)%qrs ) .AND. kk == 1 ) &
3700	READ ( 13 ) surf_v(1)%qrs
3701	CASE ( 'surf_v(1)%nrs' )
3702	IF ( ALLOCATED( surf_v(1)%nrs ) .AND. kk == 1 ) &
3703	READ ( 13 ) surf_v(1)%nrs
3704	CASE ( 'surf_v(1)%ol' )
3705	IF ( ALLOCATED( surf_v(1)%ol ) .AND. kk == 1 ) &
3706	READ ( 13 ) surf_v(1)%ol
3707	CASE ( 'surf_v(1)%rib' )
3708	IF ( ALLOCATED( surf_v(1)%rib ) .AND. kk == 1 ) &
3709	READ ( 13 ) surf_v(1)%rib
3710	CASE ( 'surf_v(1)%pt_surface' )
3711	IF ( ALLOCATED( surf_v(1)%pt_surface ) .AND. kk == 1 ) &
3712	READ ( 13 ) surf_v(1)%pt_surface
3713	CASE ( 'surf_v(1)%shf' )
3714	IF ( ALLOCATED( surf_v(1)%shf ) .AND. kk == 1 ) &
3715	READ ( 13 ) surf_v(1)%shf
3716	CASE ( 'surf_v(1)%qsws' )
3717	IF ( ALLOCATED( surf_v(1)%qsws ) .AND. kk == 1 ) &
3718	READ ( 13 ) surf_v(1)%qsws
3719	CASE ( 'surf_v(1)%ssws' )
3720	IF ( ALLOCATED( surf_v(1)%ssws ) .AND. kk == 1 ) &
3721	READ ( 13 ) surf_v(1)%ssws
3722	CASE ( 'surf_v(1)%css' )
3723	IF ( ALLOCATED( surf_v(1)%css ) .AND. kk == 1 ) &
3724	READ ( 13 ) surf_v(1)%css
3725	CASE ( 'surf_v(1)%cssws' )
3726	IF ( ALLOCATED( surf_v(1)%cssws ) .AND. kk == 1 ) &
3727	READ ( 13 ) surf_v(1)%cssws
3728	CASE ( 'surf_v(1)%qcsws' )
3729	IF ( ALLOCATED( surf_v(1)%qcsws ) .AND. kk == 1 ) &
3730	READ ( 13 ) surf_v(1)%qcsws
3731	CASE ( 'surf_v(1)%ncsws' )
3732	IF ( ALLOCATED( surf_v(1)%ncsws ) .AND. kk == 1 ) &
3733	READ ( 13 ) surf_v(1)%ncsws
3734	CASE ( 'surf_v(1)%qrsws' )
3735	IF ( ALLOCATED( surf_v(1)%qrsws ) .AND. kk == 1 ) &
3736	READ ( 13 ) surf_v(1)%qrsws
3737	CASE ( 'surf_v(1)%nrsws' )
3738	IF ( ALLOCATED( surf_v(1)%nrsws ) .AND. kk == 1 ) &
3739	READ ( 13 ) surf_v(1)%nrsws
3740	CASE ( 'surf_v(1)%sasws' )
3741	IF ( ALLOCATED( surf_v(1)%sasws ) .AND. kk == 1 ) &
3742	READ ( 13 ) surf_v(1)%sasws
3743	CASE ( 'surf_v(1)%mom_uv' )
3744	IF ( ALLOCATED( surf_v(1)%mom_flux_uv ) .AND. kk == 1 ) &
3745	READ ( 13 ) surf_v(1)%mom_flux_uv
3746	CASE ( 'surf_v(1)%mom_w' )
3747	IF ( ALLOCATED( surf_v(1)%mom_flux_w ) .AND. kk == 1 ) &
3748	READ ( 13 ) surf_v(1)%mom_flux_w
3749	CASE ( 'surf_v(1)%mom_tke' )
3750	IF ( ALLOCATED( surf_v(1)%mom_flux_tke ) .AND. kk == 1 ) &
3751	READ ( 13 ) surf_v(1)%mom_flux_tke
3752
3753	CASE ( 'surf_v(2)%start_index' )
3754	IF ( kk == 1 ) &
3755	READ ( 13 ) surf_v(2)%start_index
3756	l = 2
3757	CASE ( 'surf_v(2)%end_index' )
3758	IF ( kk == 1 ) &
3759	READ ( 13 ) surf_v(2)%end_index
3760	CASE ( 'surf_v(2)%us' )
3761	IF ( ALLOCATED( surf_v(2)%us ) .AND. kk == 1 ) &
3762	READ ( 13 ) surf_v(2)%us
3763	CASE ( 'surf_v(2)%ts' )
3764	IF ( ALLOCATED( surf_v(2)%ts ) .AND. kk == 1 ) &
3765	READ ( 13 ) surf_v(2)%ts
3766	CASE ( 'surf_v(2)%qs' )
3767	IF ( ALLOCATED( surf_v(2)%qs ) .AND. kk == 1 ) &
3768	READ ( 13 ) surf_v(2)%qs
3769	CASE ( 'surf_v(2)%ss' )
3770	IF ( ALLOCATED( surf_v(2)%ss ) .AND. kk == 1 ) &
3771	READ ( 13 ) surf_v(2)%ss
3772	CASE ( 'surf_v(2)%qcs' )
3773	IF ( ALLOCATED( surf_v(2)%qcs ) .AND. kk == 1 ) &
3774	READ ( 13 ) surf_v(2)%qcs
3775	CASE ( 'surf_v(2)%ncs' )
3776	IF ( ALLOCATED( surf_v(2)%ncs ) .AND. kk == 1 ) &
3777	READ ( 13 ) surf_v(2)%ncs
3778	CASE ( 'surf_v(2)%qrs' )
3779	IF ( ALLOCATED( surf_v(2)%qrs ) .AND. kk == 1 ) &
3780	READ ( 13 ) surf_v(2)%qrs
3781	CASE ( 'surf_v(2)%nrs' )
3782	IF ( ALLOCATED( surf_v(2)%nrs ) .AND. kk == 1 ) &
3783	READ ( 13 ) surf_v(2)%nrs
3784	CASE ( 'surf_v(2)%ol' )
3785	IF ( ALLOCATED( surf_v(2)%ol ) .AND. kk == 1 ) &
3786	READ ( 13 ) surf_v(2)%ol
3787	CASE ( 'surf_v(2)%rib' )
3788	IF ( ALLOCATED( surf_v(2)%rib ) .AND. kk == 1 ) &
3789	READ ( 13 ) surf_v(2)%rib
3790	CASE ( 'surf_v(2)%pt_surface' )
3791	IF ( ALLOCATED( surf_v(2)%pt_surface ) .AND. kk == 1 ) &
3792	READ ( 13 ) surf_v(2)%pt_surface
3793	CASE ( 'surf_v(2)%shf' )
3794	IF ( ALLOCATED( surf_v(2)%shf ) .AND. kk == 1 ) &
3795	READ ( 13 ) surf_v(2)%shf
3796	CASE ( 'surf_v(2)%qsws' )
3797	IF ( ALLOCATED( surf_v(2)%qsws ) .AND. kk == 1 ) &
3798	READ ( 13 ) surf_v(2)%qsws
3799	CASE ( 'surf_v(2)%ssws' )
3800	IF ( ALLOCATED( surf_v(2)%ssws ) .AND. kk == 1 ) &
3801	READ ( 13 ) surf_v(2)%ssws
3802	CASE ( 'surf_v(2)%css' )
3803	IF ( ALLOCATED( surf_v(2)%css ) .AND. kk == 1 ) &
3804	READ ( 13 ) surf_v(2)%css
3805	CASE ( 'surf_v(2)%cssws' )
3806	IF ( ALLOCATED( surf_v(2)%cssws ) .AND. kk == 1 ) &
3807	READ ( 13 ) surf_v(2)%cssws
3808	CASE ( 'surf_v(2)%qcsws' )
3809	IF ( ALLOCATED( surf_v(2)%qcsws ) .AND. kk == 1 ) &
3810	READ ( 13 ) surf_v(2)%qcsws
3811	CASE ( 'surf_v(2)%ncsws' )
3812	IF ( ALLOCATED( surf_v(2)%ncsws ) .AND. kk == 1 ) &
3813	READ ( 13 ) surf_v(2)%ncsws
3814	CASE ( 'surf_v(2)%qrsws' )
3815	IF ( ALLOCATED( surf_v(2)%qrsws ) .AND. kk == 1 ) &
3816	READ ( 13 ) surf_v(2)%qrsws
3817	CASE ( 'surf_v(2)%nrsws' )
3818	IF ( ALLOCATED( surf_v(2)%nrsws ) .AND. kk == 1 ) &
3819	READ ( 13 ) surf_v(2)%nrsws
3820	CASE ( 'surf_v(2)%sasws' )
3821	IF ( ALLOCATED( surf_v(2)%sasws ) .AND. kk == 1 ) &
3822	READ ( 13 ) surf_v(2)%sasws
3823	CASE ( 'surf_v(2)%mom_uv' )
3824	IF ( ALLOCATED( surf_v(2)%mom_flux_uv ) .AND. kk == 1 ) &
3825	READ ( 13 ) surf_v(2)%mom_flux_uv
3826	CASE ( 'surf_v(2)%mom_w' )
3827	IF ( ALLOCATED( surf_v(2)%mom_flux_w ) .AND. kk == 1 ) &
3828	READ ( 13 ) surf_v(2)%mom_flux_w
3829	CASE ( 'surf_v(2)%mom_tke' )
3830	IF ( ALLOCATED( surf_v(2)%mom_flux_tke ) .AND. kk == 1 ) &
3831	READ ( 13 ) surf_v(2)%mom_flux_tke
3832
3833	CASE ( 'surf_v(3)%start_index' )
3834	IF ( kk == 1 ) &
3835	READ ( 13 ) surf_v(3)%start_index
3836	l = 3
3837	CASE ( 'surf_v(3)%end_index' )
3838	IF ( kk == 1 ) &
3839	READ ( 13 ) surf_v(3)%end_index
3840	CASE ( 'surf_v(3)%us' )
3841	IF ( ALLOCATED( surf_v(3)%us ) .AND. kk == 1 ) &
3842	READ ( 13 ) surf_v(3)%us
3843	CASE ( 'surf_v(3)%ts' )
3844	IF ( ALLOCATED( surf_v(3)%ts ) .AND. kk == 1 ) &
3845	READ ( 13 ) surf_v(3)%ts
3846	CASE ( 'surf_v(3)%qs' )
3847	IF ( ALLOCATED( surf_v(3)%qs ) .AND. kk == 1 ) &
3848	READ ( 13 ) surf_v(3)%qs
3849	CASE ( 'surf_v(3)%ss' )
3850	IF ( ALLOCATED( surf_v(3)%ss ) .AND. kk == 1 ) &
3851	READ ( 13 ) surf_v(3)%ss
3852	CASE ( 'surf_v(3)%qcs' )
3853	IF ( ALLOCATED( surf_v(3)%qcs ) .AND. kk == 1 ) &
3854	READ ( 13 ) surf_v(3)%qcs
3855	CASE ( 'surf_v(3)%ncs' )
3856	IF ( ALLOCATED( surf_v(3)%ncs ) .AND. kk == 1 ) &
3857	READ ( 13 ) surf_v(3)%ncs
3858	CASE ( 'surf_v(3)%qrs' )
3859	IF ( ALLOCATED( surf_v(3)%qrs ) .AND. kk == 1 ) &
3860	READ ( 13 ) surf_v(3)%qrs
3861	CASE ( 'surf_v(3)%nrs' )
3862	IF ( ALLOCATED( surf_v(3)%nrs ) .AND. kk == 1 ) &
3863	READ ( 13 ) surf_v(3)%nrs
3864	CASE ( 'surf_v(3)%ol' )
3865	IF ( ALLOCATED( surf_v(3)%ol ) .AND. kk == 1 ) &
3866	READ ( 13 ) surf_v(3)%ol
3867	CASE ( 'surf_v(3)%rib' )
3868	IF ( ALLOCATED( surf_v(3)%rib ) .AND. kk == 1 ) &
3869	READ ( 13 ) surf_v(3)%rib
3870	CASE ( 'surf_v(3)%pt_surface' )
3871	IF ( ALLOCATED( surf_v(3)%pt_surface ) .AND. kk == 1 ) &
3872	READ ( 13 ) surf_v(3)%pt_surface
3873	CASE ( 'surf_v(3)%shf' )
3874	IF ( ALLOCATED( surf_v(3)%shf ) .AND. kk == 1 ) &
3875	READ ( 13 ) surf_v(3)%shf
3876	CASE ( 'surf_v(3)%qsws' )
3877	IF ( ALLOCATED( surf_v(3)%qsws ) .AND. kk == 1 ) &
3878	READ ( 13 ) surf_v(3)%qsws
3879	CASE ( 'surf_v(3)%ssws' )
3880	IF ( ALLOCATED( surf_v(3)%ssws ) .AND. kk == 1 ) &
3881	READ ( 13 ) surf_v(3)%ssws
3882	CASE ( 'surf_v(3)%css' )
3883	IF ( ALLOCATED( surf_v(3)%css ) .AND. kk == 1 ) &
3884	READ ( 13 ) surf_v(3)%css
3885	CASE ( 'surf_v(3)%cssws' )
3886	IF ( ALLOCATED( surf_v(3)%cssws ) .AND. kk == 1 ) &
3887	READ ( 13 ) surf_v(3)%cssws
3888	CASE ( 'surf_v(3)%qcsws' )
3889	IF ( ALLOCATED( surf_v(3)%qcsws ) .AND. kk == 1 ) &
3890	READ ( 13 ) surf_v(3)%qcsws
3891	CASE ( 'surf_v(3)%ncsws' )
3892	IF ( ALLOCATED( surf_v(3)%ncsws ) .AND. kk == 1 ) &
3893	READ ( 13 ) surf_v(3)%ncsws
3894	CASE ( 'surf_v(3)%qrsws' )
3895	IF ( ALLOCATED( surf_v(3)%qrsws ) .AND. kk == 1 ) &
3896	READ ( 13 ) surf_v(3)%qrsws
3897	CASE ( 'surf_v(3)%nrsws' )
3898	IF ( ALLOCATED( surf_v(3)%nrsws ) .AND. kk == 1 ) &
3899	READ ( 13 ) surf_v(3)%nrsws
3900	CASE ( 'surf_v(3)%sasws' )
3901	IF ( ALLOCATED( surf_v(3)%sasws ) .AND. kk == 1 ) &
3902	READ ( 13 ) surf_v(3)%sasws
3903	CASE ( 'surf_v(3)%mom_uv' )
3904	IF ( ALLOCATED( surf_v(3)%mom_flux_uv ) .AND. kk == 1 ) &
3905	READ ( 13 ) surf_v(3)%mom_flux_uv
3906	CASE ( 'surf_v(3)%mom_w' )
3907	IF ( ALLOCATED( surf_v(3)%mom_flux_w ) .AND. kk == 1 ) &
3908	READ ( 13 ) surf_v(3)%mom_flux_w
3909	CASE ( 'surf_v(3)%mom_tke' )
3910	IF ( ALLOCATED( surf_v(3)%mom_flux_tke ) .AND. kk == 1 ) &
3911	READ ( 13 ) surf_v(3)%mom_flux_tke
3912
3913	CASE DEFAULT
3914
3915	found = .FALSE.
3916
3917	END SELECT
3918
3919	!
3920	!-- Redistribute surface elements on its respective type.
3921	IF ( horizontal_surface .AND. &
3922	.NOT. INDEX( restart_string(1:length), '%start_index' ) /= 0 ) &
3923	THEN
3924	ic = nxlc
3925	DO i = nxlf, nxrf
3926	jc = nysc
3927	DO j = nysf, nynf
3928
3929	surf_match_def = surf_def_h(l)%end_index(jc,ic) >= &
3930	surf_def_h(l)%start_index(jc,ic)
3931	surf_match_lsm = ( surf_lsm_h%end_index(jc,ic) >= &
3932	surf_lsm_h%start_index(jc,ic) ) &
3933	.AND. l == 0
3934	surf_match_usm = ( surf_usm_h%end_index(jc,ic) >= &
3935	surf_usm_h%start_index(jc,ic) ) &
3936	.AND. l == 0
3937
3938	IF ( surf_match_def ) THEN
3939	mm = surf_def_h(l)%start_index(jc,ic)
3940	DO m = surf_h(l)%start_index(j,i), &
3941	surf_h(l)%end_index(j,i)
3942	IF ( surf_def_h(l)%end_index(jc,ic) >= mm ) &
3943	CALL restore_surface_elements( surf_def_h(l), &
3944	mm, surf_h(l), m )
3945	mm = mm + 1
3946	ENDDO
3947	ENDIF
3948
3949	IF ( surf_match_lsm ) THEN
3950	mm = surf_lsm_h%start_index(jc,ic)
3951	DO m = surf_h(l)%start_index(j,i), &
3952	surf_h(l)%end_index(j,i)
3953	IF ( surf_lsm_h%end_index(jc,ic) >= mm ) &
3954	CALL restore_surface_elements( surf_lsm_h, &
3955	mm, surf_h(l), m )
3956	mm = mm + 1
3957	ENDDO
3958	ENDIF
3959
3960	IF ( surf_match_usm ) THEN
3961	mm = surf_usm_h%start_index(jc,ic)
3962	DO m = surf_h(l)%start_index(j,i), &
3963	surf_h(l)%end_index(j,i)
3964	IF ( surf_usm_h%end_index(jc,ic) >= mm ) &
3965	CALL restore_surface_elements( surf_usm_h, &
3966	mm, surf_h(l), m )
3967	mm = mm + 1
3968	ENDDO
3969	ENDIF
3970
3971	jc = jc + 1
3972	ENDDO
3973	ic = ic + 1
3974	ENDDO
3975	ELSEIF ( vertical_surface .AND. &
3976	.NOT. INDEX( restart_string(1:length), '%start_index' ) /= 0 ) &
3977	THEN
3978	ic = nxlc
3979	DO i = nxlf, nxrf
3980	jc = nysc
3981	DO j = nysf, nynf
3982
3983	surf_match_def = surf_def_v(l)%end_index(jc,ic) >= &
3984	surf_def_v(l)%start_index(jc,ic)
3985	surf_match_lsm = surf_lsm_v(l)%end_index(jc,ic) >= &
3986	surf_lsm_v(l)%start_index(jc,ic)
3987	surf_match_usm = surf_usm_v(l)%end_index(jc,ic) >= &
3988	surf_usm_v(l)%start_index(jc,ic)
3989
3990	IF ( surf_match_def ) THEN
3991	mm = surf_def_v(l)%start_index(jc,ic)
3992	DO m = surf_v(l)%start_index(j,i), &
3993	surf_v(l)%end_index(j,i)
3994	IF ( surf_def_v(l)%end_index(jc,ic) >= mm ) &
3995	CALL restore_surface_elements( surf_def_v(l), &
3996	mm, surf_v(l), m )
3997	mm = mm + 1
3998	ENDDO
3999	ENDIF
4000
4001	IF ( surf_match_lsm ) THEN
4002	mm = surf_lsm_v(l)%start_index(jc,ic)
4003	DO m = surf_v(l)%start_index(j,i), &
4004	surf_v(l)%end_index(j,i)
4005	IF ( surf_lsm_v(l)%end_index(jc,ic) >= mm ) &
4006	CALL restore_surface_elements( surf_lsm_v(l), &
4007	mm, surf_v(l), m )
4008	mm = mm + 1
4009	ENDDO
4010	ENDIF
4011
4012	IF ( surf_match_usm ) THEN
4013	mm = surf_usm_v(l)%start_index(jc,ic)
4014	DO m = surf_v(l)%start_index(j,i), &
4015	surf_v(l)%end_index(j,i)
4016	IF ( surf_usm_v(l)%end_index(jc,ic) >= mm ) &
4017	CALL restore_surface_elements( surf_usm_v(l), &
4018	mm, surf_v(l), m )
4019	mm = mm + 1
4020	ENDDO
4021	ENDIF
4022
4023	jc = jc + 1
4024	ENDDO
4025	ic = ic + 1
4026	ENDDO
4027	ENDIF
4028
4029
4030	CONTAINS
4031	!------------------------------------------------------------------------------!
4032	! Description:
4033	! ------------
4034	!> Restores surfacle elements back on its respective type.
4035	!------------------------------------------------------------------------------!
4036	SUBROUTINE restore_surface_elements( surf_target, m_target, &
4037	surf_file, m_file )
4038
4039	IMPLICIT NONE
4040
4041	INTEGER(iwp) :: m_file !< respective surface-element index of current surface array
4042	INTEGER(iwp) :: m_target !< respecitve surface-element index of surface array on file
4043	INTEGER(iwp) :: lsp !< running index chemical species
4044
4045	TYPE( surf_type ) :: surf_target !< target surface type
4046	TYPE( surf_type ) :: surf_file !< surface type on file
4047
4048
4049	IF ( INDEX( restart_string(1:length), '%us' ) /= 0 ) THEN
4050	IF ( ALLOCATED( surf_target%us ) .AND. &
4051	ALLOCATED( surf_file%us ) ) &
4052	surf_target%us(m_target) = surf_file%us(m_file)
4053	ENDIF
4054
4055	IF ( INDEX( restart_string(1:length), '%ol' ) /= 0 ) THEN
4056	IF ( ALLOCATED( surf_target%ol ) .AND. &
4057	ALLOCATED( surf_file%ol ) ) &
4058	surf_target%ol(m_target) = surf_file%ol(m_file)
4059	ENDIF
4060
4061	IF ( INDEX( restart_string(1:length), '%pt_surface' ) /= 0 ) THEN
4062	IF ( ALLOCATED( surf_target%pt_surface ) .AND. &
4063	ALLOCATED( surf_file%pt_surface ) ) &
4064	surf_target%pt_surface(m_target) = surf_file%pt_surface(m_file)
4065	ENDIF
4066
4067	IF ( INDEX( restart_string(1:length), '%usws' ) /= 0 ) THEN
4068	IF ( ALLOCATED( surf_target%usws ) .AND. &
4069	ALLOCATED( surf_file%usws ) ) &
4070	surf_target%usws(m_target) = surf_file%usws(m_file)
4071	ENDIF
4072
4073	IF ( INDEX( restart_string(1:length), '%vsws' ) /= 0 ) THEN
4074	IF ( ALLOCATED( surf_target%vsws ) .AND. &
4075	ALLOCATED( surf_file%vsws ) ) &
4076	surf_target%vsws(m_target) = surf_file%vsws(m_file)
4077	ENDIF
4078
4079	IF ( INDEX( restart_string(1:length), '%ts' ) /= 0 ) THEN
4080	IF ( ALLOCATED( surf_target%ts ) .AND. &
4081	ALLOCATED( surf_file%ts ) ) &
4082	surf_target%ts(m_target) = surf_file%ts(m_file)
4083	ENDIF
4084
4085	IF ( INDEX( restart_string(1:length), '%shf' ) /= 0 ) THEN
4086	IF ( ALLOCATED( surf_target%shf ) .AND. &
4087	ALLOCATED( surf_file%shf ) ) &
4088	surf_target%shf(m_target) = surf_file%shf(m_file)
4089	ENDIF
4090
4091	IF ( INDEX( restart_string(1:length), '%qs' ) /= 0 ) THEN
4092	IF ( ALLOCATED( surf_target%qs ) .AND. &
4093	ALLOCATED( surf_file%qs ) ) &
4094	surf_target%qs(m_target) = surf_file%qs(m_file)
4095	ENDIF
4096
4097	IF ( INDEX( restart_string(1:length), '%qsws' ) /= 0 ) THEN
4098	IF ( ALLOCATED( surf_target%qsws ) .AND. &
4099	ALLOCATED( surf_file%qsws ) ) &
4100	surf_target%qsws(m_target) = surf_file%qsws(m_file)
4101	ENDIF
4102
4103	IF ( INDEX( restart_string(1:length), '%ss' ) /= 0 ) THEN
4104	IF ( ALLOCATED( surf_target%ss ) .AND. &
4105	ALLOCATED( surf_file%ss ) ) &
4106	surf_target%ss(m_target) = surf_file%ss(m_file)
4107	ENDIF
4108
4109	IF ( INDEX( restart_string(1:length), '%ssws' ) /= 0 ) THEN
4110	IF ( ALLOCATED( surf_target%ssws ) .AND. &
4111	ALLOCATED( surf_file%ssws ) ) &
4112	surf_target%ssws(m_target) = surf_file%ssws(m_file)
4113	ENDIF
4114
4115	IF ( INDEX( restart_string(1:length), '%css' ) /= 0 ) THEN
4116	IF ( ALLOCATED( surf_target%css ) .AND. &
4117	ALLOCATED( surf_file%css ) ) THEN
4118	DO lsp = 1, nvar
4119	surf_target%css(lsp,m_target) = surf_file%css(lsp,m_file)
4120	ENDDO
4121	ENDIF
4122	ENDIF
4123	IF ( INDEX( restart_string(1:length), '%cssws' ) /= 0 ) THEN
4124	IF ( ALLOCATED( surf_target%cssws ) .AND. &
4125	ALLOCATED( surf_file%cssws ) ) THEN
4126	DO lsp = 1, nvar
4127	surf_target%cssws(lsp,m_target) = surf_file%cssws(lsp,m_file)
4128	ENDDO
4129	ENDIF
4130	ENDIF
4131
4132	IF ( INDEX( restart_string(1:length), '%qcs' ) /= 0 ) THEN
4133	IF ( ALLOCATED( surf_target%qcs ) .AND. &
4134	ALLOCATED( surf_file%qcs ) ) &
4135	surf_target%qcs(m_target) = surf_file%qcs(m_file)
4136	ENDIF
4137
4138	IF ( INDEX( restart_string(1:length), '%qcsws' ) /= 0 ) THEN
4139	IF ( ALLOCATED( surf_target%qcsws ) .AND. &
4140	ALLOCATED( surf_file%qcsws ) ) &
4141	surf_target%qcsws(m_target) = surf_file%qcsws(m_file)
4142	ENDIF
4143
4144	IF ( INDEX( restart_string(1:length), '%ncs' ) /= 0 ) THEN
4145	IF ( ALLOCATED( surf_target%ncs ) .AND. &
4146	ALLOCATED( surf_file%ncs ) ) &
4147	surf_target%ncs(m_target) = surf_file%ncs(m_file)
4148	ENDIF
4149
4150	IF ( INDEX( restart_string(1:length), '%ncsws' ) /= 0 ) THEN
4151	IF ( ALLOCATED( surf_target%ncsws ) .AND. &
4152	ALLOCATED( surf_file%ncsws ) ) &
4153	surf_target%ncsws(m_target) = surf_file%ncsws(m_file)
4154	ENDIF
4155
4156	IF ( INDEX( restart_string(1:length), '%qrs' ) /= 0 ) THEN
4157	IF ( ALLOCATED( surf_target%qrs ) .AND. &
4158	ALLOCATED( surf_file%qrs ) ) &
4159	surf_target%qrs(m_target) = surf_file%qrs(m_file)
4160	ENDIF
4161
4162	IF ( INDEX( restart_string(1:length), '%qrsws' ) /= 0 ) THEN
4163	IF ( ALLOCATED( surf_target%qrsws ) .AND. &
4164	ALLOCATED( surf_file%qrsws ) ) &
4165	surf_target%qrsws(m_target) = surf_file%qrsws(m_file)
4166	ENDIF
4167
4168	IF ( INDEX( restart_string(1:length), '%nrs' ) /= 0 ) THEN
4169	IF ( ALLOCATED( surf_target%nrs ) .AND. &
4170	ALLOCATED( surf_file%nrs ) ) &
4171	surf_target%nrs(m_target) = surf_file%nrs(m_file)
4172	ENDIF
4173
4174	IF ( INDEX( restart_string(1:length), '%nrsws' ) /= 0 ) THEN
4175	IF ( ALLOCATED( surf_target%nrsws ) .AND. &
4176	ALLOCATED( surf_file%nrsws ) ) &
4177	surf_target%nrsws(m_target) = surf_file%nrsws(m_file)
4178	ENDIF
4179
4180	IF ( INDEX( restart_string(1:length), '%sasws' ) /= 0 ) THEN
4181	IF ( ALLOCATED( surf_target%sasws ) .AND. &
4182	ALLOCATED( surf_file%sasws ) ) &
4183	surf_target%sasws(m_target) = surf_file%sasws(m_file)
4184	ENDIF
4185
4186	IF ( INDEX( restart_string(1:length), '%mom_uv' ) /= 0 ) THEN
4187	IF ( ALLOCATED( surf_target%mom_flux_uv ) .AND. &
4188	ALLOCATED( surf_file%mom_flux_uv ) ) &
4189	surf_target%mom_flux_uv(m_target) = &
4190	surf_file%mom_flux_uv(m_file)
4191	ENDIF
4192
4193	IF ( INDEX( restart_string(1:length), '%mom_w' ) /= 0 ) THEN
4194	IF ( ALLOCATED( surf_target%mom_flux_w ) .AND. &
4195	ALLOCATED( surf_file%mom_flux_w ) ) &
4196	surf_target%mom_flux_w(m_target) = &
4197	surf_file%mom_flux_w(m_file)
4198	ENDIF
4199
4200	IF ( INDEX( restart_string(1:length), '%mom_tke' ) /= 0 ) THEN
4201	IF ( ALLOCATED( surf_target%mom_flux_tke ) .AND. &
4202	ALLOCATED( surf_file%mom_flux_tke ) ) &
4203	surf_target%mom_flux_tke(0:1,m_target) = &
4204	surf_file%mom_flux_tke(0:1,m_file)
4205	ENDIF
4206
4207
4208	END SUBROUTINE restore_surface_elements
4209
4210
4211	END SUBROUTINE surface_rrd_local
4212
4213
4214	!------------------------------------------------------------------------------!
4215	! Description:
4216	! ------------
4217	!> Counts the number of surface elements with the same facing, required for
4218	!> reading and writing restart data.
4219	!------------------------------------------------------------------------------!
4220	SUBROUTINE surface_last_actions
4221
4222	IMPLICIT NONE
4223	!
4224	!-- Horizontal surfaces
4225	ns_h_on_file(0) = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns
4226	ns_h_on_file(1) = surf_def_h(1)%ns
4227	ns_h_on_file(2) = surf_def_h(2)%ns
4228	!
4229	!-- Vertical surfaces
4230	ns_v_on_file(0) = surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns
4231	ns_v_on_file(1) = surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns
4232	ns_v_on_file(2) = surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns
4233	ns_v_on_file(3) = surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns
4234
4235	END SUBROUTINE surface_last_actions
4236
4237	!------------------------------------------------------------------------------!
4238	! Description:
4239	! ------------
4240	!> Routine maps surface data read from file after restart - 1D arrays
4241	!------------------------------------------------------------------------------!
4242	SUBROUTINE surface_restore_elements_1d( surf_target, surf_file, &
4243	start_index_c, &
4244	start_index_on_file, &
4245	end_index_on_file, &
4246	nxlc, nysc, nxlf, nxrf, nysf, nynf,&
4247	nys_on_file, nyn_on_file, &
4248	nxl_on_file,nxr_on_file )
4249
4250	IMPLICIT NONE
4251
4252	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
4253	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
4254	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
4255	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
4256	INTEGER(iwp) :: m !< surface-element index on file
4257	INTEGER(iwp) :: mm !< surface-element index on current subdomain
4258	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
4259	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
4260	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
4261	INTEGER(iwp) :: nysc !< index of north boundary on current subdomain
4262	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
4263	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
4264
4265	INTEGER(iwp) :: nxl_on_file !< leftmost index on file
4266	INTEGER(iwp) :: nxr_on_file !< rightmost index on file
4267	INTEGER(iwp) :: nyn_on_file !< northmost index on file
4268	INTEGER(iwp) :: nys_on_file !< southmost index on file
4269
4270	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: start_index_c
4271	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
4272	start_index_on_file !< start index of surface elements on file
4273	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
4274	end_index_on_file !< end index of surface elements on file
4275
4276	REAL(wp), DIMENSION(:) :: surf_target !< target surface type
4277	REAL(wp), DIMENSION(:) :: surf_file !< surface type on file
4278
4279	ic = nxlc
4280	DO i = nxlf, nxrf
4281	jc = nysc
4282	DO j = nysf, nynf
4283
4284	mm = start_index_c(jc,ic)
4285	DO m = start_index_on_file(j,i), end_index_on_file(j,i)
4286	surf_target(mm) = surf_file(m)
4287	mm = mm + 1
4288	ENDDO
4289
4290	jc = jc + 1
4291	ENDDO
4292	ic = ic + 1
4293	ENDDO
4294
4295
4296	END SUBROUTINE surface_restore_elements_1d
4297
4298	!------------------------------------------------------------------------------!
4299	! Description:
4300	! ------------
4301	!> Routine maps surface data read from file after restart - 2D arrays
4302	!------------------------------------------------------------------------------!
4303	SUBROUTINE surface_restore_elements_2d( surf_target, surf_file, &
4304	start_index_c, &
4305	start_index_on_file, &
4306	end_index_on_file, &
4307	nxlc, nysc, nxlf, nxrf, nysf, nynf,&
4308	nys_on_file, nyn_on_file, &
4309	nxl_on_file,nxr_on_file )
4310
4311	IMPLICIT NONE
4312
4313	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
4314	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
4315	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
4316	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
4317	INTEGER(iwp) :: m !< surface-element index on file
4318	INTEGER(iwp) :: mm !< surface-element index on current subdomain
4319	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
4320	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
4321	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
4322	INTEGER(iwp) :: nysc !< index of north boundary on current subdomain
4323	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
4324	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
4325
4326	INTEGER(iwp) :: nxl_on_file !< leftmost index on file
4327	INTEGER(iwp) :: nxr_on_file !< rightmost index on file
4328	INTEGER(iwp) :: nyn_on_file !< northmost index on file
4329	INTEGER(iwp) :: nys_on_file !< southmost index on file
4330
4331	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: start_index_c !< start index of surface type
4332	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
4333	start_index_on_file !< start index of surface elements on file
4334	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
4335	end_index_on_file !< end index of surface elements on file
4336
4337	REAL(wp), DIMENSION(:,:) :: surf_target !< target surface type
4338	REAL(wp), DIMENSION(:,:) :: surf_file !< surface type on file
4339
4340	ic = nxlc
4341	DO i = nxlf, nxrf
4342	jc = nysc
4343	DO j = nysf, nynf
4344	mm = start_index_c(jc,ic)
4345	DO m = start_index_on_file(j,i), end_index_on_file(j,i)
4346	surf_target(:,mm) = surf_file(:,m)
4347	mm = mm + 1
4348	ENDDO
4349
4350	jc = jc + 1
4351	ENDDO
4352	ic = ic + 1
4353	ENDDO
4354
4355	END SUBROUTINE surface_restore_elements_2d
4356
4357
4358	END MODULE surface_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |