Home

Context Navigation

surface_mod.f90 @ 2942

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