Home

Context Navigation

surface_mod.f90 @ 4156

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