Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |