Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |