Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |