Home

Context Navigation

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

Last change on this file since 2712 was 2707, checked in by suehring, 6 years ago
changes documented
Property svn:keywords set to `Id`
File size: 201.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |