Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |