Home

Context Navigation

surface_mod.f90 @ 2812

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