Home

Context Navigation

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

Last change on this file since 2980 was 2977, checked in by kanani, 6 years ago
Fixes for radiative transfer model
Property svn:keywords set to `Id`
File size: 210.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |