Home

Context Navigation

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

Last change on this file since 2759 was 2759, checked in by suehring, 6 years ago
Major bugfix in horizontal diffusion of all scalar quantities at vertical surfaces; further bugfix: density is considered for vertical fluxes of passive scalar and salinity
Property svn:keywords set to `Id`
File size: 202.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |