Home

Context Navigation

surface_mod.f90 @ 3943

Last change on this file since 3943 was 3943, checked in by maronga, 5 years ago
bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements
Property svn:keywords set to `Id`
File size: 241.3 KB

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