Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |