Home

Context Navigation

surface_mod.f90 @ 4102

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