Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |