Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |