Home

Context Navigation

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

Last change on this file since 4355 was 4355, checked in by suehring, 5 years ago
surface_mod: bugfix in message call and define appropriate error number
Property svn:keywords set to `Id`
File size: 242.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |