Home

Context Navigation

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

Last change on this file since 3255 was 3254, checked in by suehring, 6 years ago
Additional checks for surface_fractions and building_id; Remove redundant subroutine argument in surface_mod
Property svn:keywords set to `Id`
File size: 219.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |