Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |