Home

Context Navigation

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

Last change on this file since 3299 was 3298, checked in by kanani, 6 years ago
Merge chemistry branch at r3297 to trunk
Property svn:keywords set to `Id`
File size: 220.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |