Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |