Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |