Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |