Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |