Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |