Home

Context Navigation

surface_mod.f90 @ 3933

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