Home

Context Navigation

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

Last change on this file since 2269 was 2269, checked in by suehring, 7 years ago
Enable restarts with USM with different number of PEs; some bugfixes in new surface structure in USM; formatting adjustments and descriptions in surface_mod
Property svn:keywords set to `Id`
File size: 146.9 KB

Line
1	!> @file surface_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
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-2017 Leibniz Universitaet Hannover
18	!
19	!------------------------------------------------------------------------------!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: surface_mod.f90 2269 2017-06-09 11:57:32Z suehring $
27	! Formatting and description adjustments
28	!
29	! 2256 2017-06-07 13:58:08Z suehring
30	! Enable heating at downward-facing surfaces
31	!
32	! 2233 2017-05-30 18:08:54Z suehring
33	! Initial revision
34	!
35	!
36	! Description:
37	! ------------
38	!> Surface module defines derived data structures to treat surface-
39	!> bounded grid cells. Three different types of surfaces are defined:
40	!> default surfaces, natural surfaces, and urban surfaces. The module
41	!> encompasses the allocation and initialization of surface arrays, and handles
42	!> reading and writing restart data.
43	!> In addition, a further derived data structure is defined, in order to set
44	!> boundary conditions at surfaces.
45	!------------------------------------------------------------------------------!
46	MODULE surface_mod
47
48	USE arrays_3d, &
49	ONLY: zu, zw, heatflux_input_conversion, waterflux_input_conversion, &
50	momentumflux_input_conversion
51
52	USE control_parameters
53
54	USE indices, &
55	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, wall_flags_0
56
57	USE grid_variables, &
58	ONLY: dx, dy
59
60	USE kinds
61
62	USE model_1d, &
63	ONLY: rif1d, us1d, usws1d, vsws1d
64
65
66	IMPLICIT NONE
67
68	!
69	!-- Data type used to identify grid-points where horizontal boundary conditions
70	!-- are applied
71	TYPE bc_type
72
73	INTEGER(iwp) :: ns !< number of surface elements on the PE
74
75	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid
76	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid
77	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid
78
79	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
80	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< end index within surface data type for given (j,i)
81
82	END TYPE bc_type
83	!
84	!-- Data type used to identify and treat surface-bounded grid points
85	TYPE surf_type
86
87	INTEGER(iwp) :: ns !< number of surface elements on the PE
88	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid
89	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid
90	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid
91
92	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: facing !< Bit indicating surface orientation
93
94	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i)
95	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< End index within surface data type for given (j,i)
96
97	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_mo !< surface-layer height
98	REAL(wp), DIMENSION(:), ALLOCATABLE :: uvw_abs !< absolute surface-parallel velocity
99	REAL(wp), DIMENSION(:), ALLOCATABLE :: us !< friction velocity
100	REAL(wp), DIMENSION(:), ALLOCATABLE :: ts !< scaling parameter temerature
101	REAL(wp), DIMENSION(:), ALLOCATABLE :: qs !< scaling parameter humidity
102	REAL(wp), DIMENSION(:), ALLOCATABLE :: ss !< scaling parameter passive scalar
103	REAL(wp), DIMENSION(:), ALLOCATABLE :: qrs !< scaling parameter qr
104	REAL(wp), DIMENSION(:), ALLOCATABLE :: nrs !< scaling parameter nr
105
106	REAL(wp), DIMENSION(:), ALLOCATABLE :: ol !< Obukhov length
107	REAL(wp), DIMENSION(:), ALLOCATABLE :: rib !< Richardson bulk number
108
109	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0 !< roughness length for momentum
110	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0h !< roughness length for heat
111	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0q !< roughness length for humidity
112
113	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt1 !< Specific humidity at first grid level (required for cloud_physics = .T.)
114	REAL(wp), DIMENSION(:), ALLOCATABLE :: qv1 !< Potential temperature at first grid level (required for cloud_physics = .T.)
115	!
116	!-- Define arrays for surface fluxes
117	REAL(wp), DIMENSION(:), ALLOCATABLE :: usws !< vertical momentum flux for u-component at horizontal surfaces
118	REAL(wp), DIMENSION(:), ALLOCATABLE :: vsws !< vertical momentum flux for v-component at horizontal surfaces
119
120	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf !< surface flux sensible heat
121	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws !< surface flux latent heat
122	REAL(wp), DIMENSION(:), ALLOCATABLE :: ssws !< surface flux passive scalar
123	REAL(wp), DIMENSION(:), ALLOCATABLE :: qrsws !< surface flux qr
124	REAL(wp), DIMENSION(:), ALLOCATABLE :: nrsws !< surface flux nr
125	REAL(wp), DIMENSION(:), ALLOCATABLE :: sasws !< surface flux salinity
126	!
127	!-- Required for horizontal walls in production_e
128	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_0 !< virtual velocity component (see production_e_init for further explanation)
129	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_0 !< virtual velocity component (see production_e_init for further explanation)
130
131	REAL(wp), DIMENSION(:), ALLOCATABLE :: mom_flux_uv !< momentum flux usvs and vsus at vertical surfaces (used in diffusion_u and diffusion_v)
132	REAL(wp), DIMENSION(:), ALLOCATABLE :: mom_flux_w !< momentum flux wsus and wsvs at vertical surfaces (used in diffusion_w)
133	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mom_flux_tke !< momentum flux usvs, vsus, wsus, wsvs at vertical surfaces at grid center (used in production_e)
134	!
135	!-- Variables required for LSM as well as for USM
136	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_surface !< skin-surface temperature
137	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< net radiation
138	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h !< heat conductivity of soil (W/m/K)
139
140	!
141	!-- Variables required for land-surface model
142	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: soil_type_2d !< soil type, 0: user-defined, 1-7: generic (see list)
143	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: veg_type_2d !< vegetation type, 0: user-defined, 1-19: generic (see list)
144
145	LOGICAL, DIMENSION(:), ALLOCATABLE :: building_surface !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
146	LOGICAL, DIMENSION(:), ALLOCATABLE :: pave_surface !< flag parameter for pavements (asphalt etc.) (class 20)
147	LOGICAL, DIMENSION(:), ALLOCATABLE :: water_surface !< flag parameter for water surfaces (classes 14+15)
148
149	REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_vg !< coef. of Van Genuchten
150	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_liq !< liquid water coverage (of vegetated area)
151	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_veg !< vegetation coverage
152	REAL(wp), DIMENSION(:), ALLOCATABLE :: f_sw_in !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
153	REAL(wp), DIMENSION(:), ALLOCATABLE :: ghf_eb !< ground heat flux
154	REAL(wp), DIMENSION(:), ALLOCATABLE :: gamma_w_sat !< hydraulic conductivity at saturation
155	REAL(wp), DIMENSION(:), ALLOCATABLE :: g_d !< coefficient for dependence of r_canopy on water vapour pressure deficit
156	REAL(wp), DIMENSION(:), ALLOCATABLE :: l_vg !< coef. of Van Genuchten
157	REAL(wp), DIMENSION(:), ALLOCATABLE :: lai !< leaf area index
158	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surface_u !< coupling between surface and soil (depends on vegetation type)
159	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surface_s !< coupling between surface and soil (depends on vegetation type)
160	REAL(wp), DIMENSION(:), ALLOCATABLE :: m_fc !< soil moisture at field capacity (m3/m3)
161	REAL(wp), DIMENSION(:), ALLOCATABLE :: m_res !< residual soil moisture
162	REAL(wp), DIMENSION(:), ALLOCATABLE :: m_sat !< saturation soil moisture (m3/m3)
163	REAL(wp), DIMENSION(:), ALLOCATABLE :: m_wilt !< soil moisture at permanent wilting point (m3/m3)
164	REAL(wp), DIMENSION(:), ALLOCATABLE :: n_vg !< coef. Van Genuchten
165	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_eb !< surface flux of latent heat (total)
166	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_liq_eb !< surface flux of latent heat (liquid water portion)
167	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_soil_eb !< surface flux of latent heat (soil portion)
168	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_veg_eb !< surface flux of latent heat (vegetation portion)
169	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_a !< aerodynamic resistance
170	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_canopy !< canopy resistance
171	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_soil !< soil resistance
172	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_soil_min !< minimum soil resistance
173	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_s !< total surface resistance (combination of r_soil and r_canopy)
174	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_canopy_min !< minimum canopy (stomatal) resistance
175	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf_eb !< surface flux of sensible heat
176
177	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_w !< hydraulic diffusivity of soil (?)
178	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gamma_w !< hydraulic conductivity of soil (W/m/K)
179	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_total !< volumetric heat capacity of the actual soil matrix (?)
180	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: root_fr !< root fraction within the soil layers
181	!
182	!-- Urban surface variables
183	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_types !< array of types of wall parameters
184
185	LOGICAL, DIMENSION(:), ALLOCATABLE :: isroof_surf !< flag indication roof surfaces
186
187	REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface
188	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface !< heat capacity of the wall surface skin ( J mâ2 Kâ1 )
189	REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface
190	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 )
191	REAL(wp), DIMENSION(:), ALLOCATABLE :: roughness_wall !< roughness relative to concrete
192	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_wall !< thickness of the wall, roof and soil layers
193
194	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< reflected shortwave radiation for local surface in i-th reflection
195	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< reflected + emitted longwave radiation for local surface in i-th reflection
196	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf !< total radiation flux incoming to minus outgoing from local surface
197
198	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K)
199	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf !< kinematic wall heat flux of sensible heat (actually no longer needed)
200	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb !< wall heat flux of sensible heat in wall normal direction
201
202
203	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb !< wall ground heat flux
204
205	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_in_sw !< incoming shortwave radiation
206	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_out_sw !< emitted shortwave radiation
207	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_in_lw !< incoming longwave radiation
208	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_out_lw !< emitted longwave radiation
209
210	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< shortwave radiation falling to local surface including radiation from reflections
211	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< total shortwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
212	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< longwave radiation falling to local surface including radiation from reflections
213	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< total longwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
214
215
216	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_wall !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)
217	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall !< wall grid spacing (center-center)
218	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall !< 1/dz_wall
219	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall_stag !< wall grid spacing (edge-edge)
220	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall_stag !< 1/dz_wall_stag
221	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array
222	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw !< wall layer depths (m)
223
224
225	!-- arrays for time averages
226	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l
227	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
228	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
229	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
230	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
231	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
232	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
233	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
234	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
235	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
236	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
237	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
238	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface
239	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_av !< average of wghf_eb
240	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb_av !< average of wshf_eb
241	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_av !< average of wall surface temperature (K)
242
243	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_wall_av !< Average of t_wall
244
245	END TYPE surf_type
246
247	TYPE (bc_type), DIMENSION(0:1) :: bc_h !< boundary condition data type, horizontal upward- and downward facing surfaces
248
249	TYPE (surf_type), DIMENSION(0:2), TARGET :: surf_def_h !< horizontal default surfaces (Up, Down, and Top)
250	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_def_v !< vertical default surfaces (North, South, West, East)
251	TYPE (surf_type) , TARGET :: surf_lsm_h !< horizontal natural land surfaces, so far only upward-facing
252	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_lsm_v !< vertical land surfaces (North, South, West, East)
253	TYPE (surf_type) , TARGET :: surf_usm_h !< horizontal urban surfaces, so far only upward-facing
254	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_usm_v !< vertical urban surfaces (North, South, West, East)
255
256	INTEGER(iwp) :: ns_h_on_file(0:2) !< total number of horizontal surfaces with the same facing, required for writing restart data
257	INTEGER(iwp) :: ns_v_on_file(0:3) !< total number of vertical surfaces with the same facing, required for writing restart data
258
259
260	SAVE
261
262	PRIVATE
263	!
264	!-- Public variables
265	PUBLIC bc_h, ns_h_on_file, ns_v_on_file, surf_def_h, surf_def_v, &
266	surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type
267	!
268	!-- Public subroutines
269	PUBLIC init_bc, init_surfaces, init_surface_arrays, &
270	surface_read_restart_data, surface_write_restart_data, &
271	surface_last_actions
272
273	INTERFACE init_bc
274	MODULE PROCEDURE init_bc
275	END INTERFACE init_bc
276
277	INTERFACE init_surfaces
278	MODULE PROCEDURE init_surfaces
279	END INTERFACE init_surfaces
280
281	INTERFACE init_surface_arrays
282	MODULE PROCEDURE init_surface_arrays
283	END INTERFACE init_surface_arrays
284
285	INTERFACE surface_read_restart_data
286	MODULE PROCEDURE surface_read_restart_data
287	END INTERFACE surface_read_restart_data
288
289	INTERFACE surface_write_restart_data
290	MODULE PROCEDURE surface_write_restart_data
291	END INTERFACE surface_write_restart_data
292
293	INTERFACE surface_last_actions
294	MODULE PROCEDURE surface_last_actions
295	END INTERFACE surface_last_actions
296
297
298	CONTAINS
299
300	!------------------------------------------------------------------------------!
301	! Description:
302	! ------------
303	!> Initialize data type for setting boundary conditions at horizontal surfaces.
304	!------------------------------------------------------------------------------!
305	SUBROUTINE init_bc
306
307	IMPLICIT NONE
308
309	INTEGER(iwp) :: i !<
310	INTEGER(iwp) :: j !<
311	INTEGER(iwp) :: k !<
312
313	INTEGER(iwp), DIMENSION(0:1) :: num_h !<
314	INTEGER(iwp), DIMENSION(0:1) :: num_h_kji !<
315	INTEGER(iwp), DIMENSION(0:1) :: start_index_h !<
316
317	!
318	!-- First of all, count the number of upward- and downward-facing surfaces
319	num_h = 0
320	DO i = nxlg, nxrg
321	DO j = nysg, nyng
322	DO k = nzb+1, nzt
323	!
324	!-- Check if current gridpoint belongs to the atmosphere
325	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
326	!
327	!-- Upward-facing
328	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) &
329	num_h(0) = num_h(0) + 1
330	!
331	!-- Downward-facing
332	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
333	num_h(1) = num_h(1) + 1
334	ENDIF
335	ENDDO
336	ENDDO
337	ENDDO
338	!
339	!-- Save the number of surface elements
340	bc_h(0)%ns = num_h(0)
341	bc_h(1)%ns = num_h(1)
342	!
343	!-- ALLOCATE data type variables
344	!-- Upward facing
345	ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
346	ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
347	ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
348	ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
349	ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg) )
350	bc_h(0)%start_index = 1
351	bc_h(0)%end_index = 0
352	!
353	!-- Downward facing
354	ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
355	ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
356	ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
357	ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
358	ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg) )
359	bc_h(1)%start_index = 1
360	bc_h(1)%end_index = 0
361	!
362	!-- Store the respective indices on data type
363	num_h(0:1) = 1
364	start_index_h(0:1) = 1
365	DO i = nxlg, nxrg
366	DO j = nysg, nyng
367
368	num_h_kji(0:1) = 0
369	DO k = nzb+1, nzt
370	!
371	!-- Check if current gridpoint belongs to the atmosphere
372	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
373	!
374	!-- Upward-facing
375	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
376	bc_h(0)%i(num_h(0)) = i
377	bc_h(0)%j(num_h(0)) = j
378	bc_h(0)%k(num_h(0)) = k
379	num_h_kji(0) = num_h_kji(0) + 1
380	num_h(0) = num_h(0) + 1
381	ENDIF
382	!
383	!-- Downward-facing
384	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
385	bc_h(1)%i(num_h(1)) = i
386	bc_h(1)%j(num_h(1)) = j
387	bc_h(1)%k(num_h(1)) = k
388	num_h_kji(1) = num_h_kji(1) + 1
389	num_h(1) = num_h(1) + 1
390	ENDIF
391	ENDIF
392	ENDDO
393	bc_h(0)%start_index(j,i) = start_index_h(0)
394	bc_h(0)%end_index(j,i) = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
395	start_index_h(0) = bc_h(0)%end_index(j,i) + 1
396
397	bc_h(1)%start_index(j,i) = start_index_h(1)
398	bc_h(1)%end_index(j,i) = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
399	start_index_h(1) = bc_h(1)%end_index(j,i) + 1
400	ENDDO
401	ENDDO
402
403
404	END SUBROUTINE init_bc
405
406
407	!------------------------------------------------------------------------------!
408	! Description:
409	! ------------
410	!> Initialize horizontal and vertical surfaces. Counts the number of default-,
411	!> natural and urban surfaces and allocates memory, respectively.
412	!------------------------------------------------------------------------------!
413	SUBROUTINE init_surface_arrays
414
415	IMPLICIT NONE
416
417	INTEGER(iwp) :: i !< running index x-direction
418	INTEGER(iwp) :: j !< running index y-direction
419	INTEGER(iwp) :: k !< running index z-direction
420	INTEGER(iwp) :: l !< index variable for surface facing
421	INTEGER(iwp) :: num_lsm_h !< number of horizontally-aligned natural surfaces
422	INTEGER(iwp) :: num_usm_h !< number of horizontally-aligned urban surfaces
423
424	INTEGER(iwp), DIMENSION(0:2) :: num_def_h !< number of horizontally-aligned default surfaces
425	INTEGER(iwp), DIMENSION(0:3) :: num_def_v !< number of vertically-aligned default surfaces
426	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< number of vertically-aligned natural surfaces
427	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< number of vertically-aligned urban surfaces
428
429
430	num_def_h = 0
431	num_def_v = 0
432	num_lsm_h = 0
433	num_lsm_v = 0
434	num_usm_h = 0
435	num_usm_v = 0
436	!
437	!-- Count number of horizontal surfaces on local domain
438	DO i = nxl, nxr
439	DO j = nys, nyn
440	DO k = nzb+1, nzt
441	!
442	!-- Check if current gridpoint belongs to the atmosphere
443	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
444	!
445	!-- Check if grid point adjoins to any upward-facing horizontal
446	!-- surface, e.g. the Earth surface, plane roofs, or ceilings.
447	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
448	!
449	!-- Land-surface type
450	IF ( land_surface ) THEN
451	num_lsm_h = num_lsm_h + 1
452	!
453	!-- Urban surface tpye
454	ELSEIF ( urban_surface ) THEN
455	num_usm_h = num_usm_h + 1
456	!
457	!-- Default-surface type
458	ELSE
459	num_def_h(0) = num_def_h(0) + 1
460	ENDIF
461
462	ENDIF
463	!
464	!-- Check for top-fluxes
465	IF ( k == nzt .AND. use_top_fluxes ) THEN
466	num_def_h(2) = num_def_h(2) + 1
467	!
468	!-- Check for any other downward-facing surface. So far only for
469	!-- default surface type.
470	ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
471	num_def_h(1) = num_def_h(1) + 1
472	ENDIF
473
474	ENDIF
475	ENDDO
476	ENDDO
477	ENDDO
478	!
479	!-- Count number of vertical surfaces on local domain
480	DO i = nxl, nxr
481	DO j = nys, nyn
482	DO k = nzb+1, nzt
483	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
484	!
485	!-- Northward-facing
486	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) ) THEN
487	IF ( urban_surface ) THEN
488	num_usm_v(0) = num_usm_v(0) + 1
489	ELSEIF ( land_surface ) THEN
490	num_lsm_v(0) = num_lsm_v(0) + 1
491	ELSE
492	num_def_v(0) = num_def_v(0) + 1
493	ENDIF
494	ENDIF
495	!
496	!-- Southward-facing
497	IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) THEN
498	IF ( urban_surface ) THEN
499	num_usm_v(1) = num_usm_v(1) + 1
500	ELSEIF ( land_surface ) THEN
501	num_lsm_v(1) = num_lsm_v(1) + 1
502	ELSE
503	num_def_v(1) = num_def_v(1) + 1
504	ENDIF
505	ENDIF
506	!
507	!-- Eastward-facing
508	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) ) THEN
509	IF ( urban_surface ) THEN
510	num_usm_v(2) = num_usm_v(2) + 1
511	ELSEIF ( land_surface ) THEN
512	num_lsm_v(2) = num_lsm_v(2) + 1
513	ELSE
514	num_def_v(2) = num_def_v(2) + 1
515	ENDIF
516	ENDIF
517	!
518	!-- Westward-facing
519	IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) THEN
520	IF ( urban_surface ) THEN
521	num_usm_v(3) = num_usm_v(3) + 1
522	ELSEIF ( land_surface ) THEN
523	num_lsm_v(3) = num_lsm_v(3) + 1
524	ELSE
525	num_def_v(3) = num_def_v(3) + 1
526	ENDIF
527	ENDIF
528	ENDIF
529	ENDDO
530	ENDDO
531	ENDDO
532
533	!
534	!-- Store number of surfaces per core.
535	!-- Horizontal surface, default type, upward facing
536	surf_def_h(0)%ns = num_def_h(0)
537	!
538	!-- Horizontal surface, default type, downward facing
539	surf_def_h(1)%ns = num_def_h(1)
540	!
541	!-- Horizontal surface, default type, top downward facing
542	surf_def_h(2)%ns = num_def_h(2)
543	!
544	!-- Horizontal surface, natural type, so far only upward-facing
545	surf_lsm_h%ns = num_lsm_h
546	!
547	!-- Horizontal surface, urban type, so far only upward-facing
548	surf_usm_h%ns = num_usm_h
549	!
550	!-- Vertical surface, default type, northward facing
551	surf_def_v(0)%ns = num_def_v(0)
552	!
553	!-- Vertical surface, default type, southward facing
554	surf_def_v(1)%ns = num_def_v(1)
555	!
556	!-- Vertical surface, default type, eastward facing
557	surf_def_v(2)%ns = num_def_v(2)
558	!
559	!-- Vertical surface, default type, westward facing
560	surf_def_v(3)%ns = num_def_v(3)
561	!
562	!-- Vertical surface, natural type, northward facing
563	surf_lsm_v(0)%ns = num_lsm_v(0)
564	!
565	!-- Vertical surface, natural type, southward facing
566	surf_lsm_v(1)%ns = num_lsm_v(1)
567	!
568	!-- Vertical surface, natural type, eastward facing
569	surf_lsm_v(2)%ns = num_lsm_v(2)
570	!
571	!-- Vertical surface, natural type, westward facing
572	surf_lsm_v(3)%ns = num_lsm_v(3)
573	!
574	!-- Vertical surface, urban type, northward facing
575	surf_usm_v(0)%ns = num_usm_v(0)
576	!
577	!-- Vertical surface, urban type, southward facing
578	surf_usm_v(1)%ns = num_usm_v(1)
579	!
580	!-- Vertical surface, urban type, eastward facing
581	surf_usm_v(2)%ns = num_usm_v(2)
582	!
583	!-- Vertical surface, urban type, westward facing
584	surf_usm_v(3)%ns = num_usm_v(3)
585	!
586	!-- Allocate required attributes for horizontal surfaces - default type.
587	!-- Upward-facing (l=0) and downward-facing (l=1).
588	DO l = 0, 1
589	CALL allocate_surface_attributes_h ( surf_def_h(l), nys, nyn, nxl, nxr )
590	ENDDO
591	!
592	!-- Allocate required attributes for model top
593	CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
594	!
595	!-- Allocate required attributes for horizontal surfaces - natural type.
596	CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
597	!
598	!-- Allocate required attributes for horizontal surfaces - urban type.
599	CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
600
601	!
602	!-- Allocate required attributes for vertical surfaces.
603	!-- Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2)
604	!-- and westward-facing (l=3).
605	!-- Default type.
606	DO l = 0, 3
607	CALL allocate_surface_attributes_v ( surf_def_v(l), .FALSE., &
608	nys, nyn, nxl, nxr )
609	ENDDO
610	!
611	!-- Natural type
612	DO l = 0, 3
613	CALL allocate_surface_attributes_v ( surf_lsm_v(l), .TRUE., &
614	nys, nyn, nxl, nxr )
615	ENDDO
616	!
617	!-- Urban type
618	DO l = 0, 3
619	CALL allocate_surface_attributes_v ( surf_usm_v(l), .FALSE., &
620	nys, nyn, nxl, nxr )
621	ENDDO
622
623	END SUBROUTINE init_surface_arrays
624
625	!------------------------------------------------------------------------------!
626	! Description:
627	! ------------
628	!> Allocating memory for upward and downward-facing horizontal surface types,
629	!> except for top fluxes.
630	!------------------------------------------------------------------------------!
631	SUBROUTINE allocate_surface_attributes_h( surfaces, &
632	nys_l, nyn_l, nxl_l, nxr_l )
633
634	IMPLICIT NONE
635
636	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
637	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
638	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
639	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
640
641	TYPE(surf_type) :: surfaces !< respective surface type
642
643	!
644	!-- Allocate arrays for start and end index of horizontal surface type
645	!-- for each (j,i)-grid point. This is required e.g. in diffion_x, which is
646	!-- called for each (j,i). In order to find the location where the
647	!-- respective flux is store within the surface-type, start- and end-
648	!-- index are stored for each (j,i). For example, each (j,i) can have
649	!-- several entries where fluxes for horizontal surfaces might be stored,
650	!-- e.g. for overhanging structures where several upward-facing surfaces
651	!-- might exist for given (j,i).
652	!-- If no surface of respective type exist at current (j,i), set indicies
653	!-- such that loop in diffusion routines will not be entered.
654	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
655	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
656	surfaces%start_index = 0
657	surfaces%end_index = -1
658	!
659	!-- Indices to locate surface element
660	ALLOCATE ( surfaces%i(1:surfaces%ns) )
661	ALLOCATE ( surfaces%j(1:surfaces%ns) )
662	ALLOCATE ( surfaces%k(1:surfaces%ns) )
663	!
664	!-- Surface-layer height
665	ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
666	!
667	!-- Surface orientation
668	ALLOCATE ( surfaces%facing(1:surfaces%ns) )
669	!
670	!-- Surface-parallel wind velocity
671	ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
672	! ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
673	!
674	!-- Roughness
675	ALLOCATE ( surfaces%z0(1:surfaces%ns) )
676	ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
677	ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
678	!
679	!-- Friction velocity
680	ALLOCATE ( surfaces%us(1:surfaces%ns) )
681	!
682	!-- Stability parameter
683	ALLOCATE ( surfaces%ol(1:surfaces%ns) )
684	!
685	!-- Bulk Richardson number
686	ALLOCATE ( surfaces%rib(1:surfaces%ns) )
687	!
688	!-- Vertical momentum fluxes of u and v
689	ALLOCATE ( surfaces%usws(1:surfaces%ns) )
690	ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
691	!
692	!-- Required in production_e
693	IF ( .NOT. constant_diffusion ) THEN
694	ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
695	ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
696	ENDIF
697	!
698	!-- Characteristic temperature and surface flux of sensible heat
699	ALLOCATE ( surfaces%ts(1:surfaces%ns) )
700	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
701	!
702	!-- Characteristic humidity and surface flux of latent heat
703	IF ( humidity ) THEN
704	ALLOCATE ( surfaces%qs(1:surfaces%ns) )
705	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
706	ENDIF
707	!
708	!-- Characteristic scalar and surface flux of scalar
709	IF ( passive_scalar ) THEN
710	ALLOCATE ( surfaces%ss(1:surfaces%ns) )
711	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
712	ENDIF
713	!
714	!-- When cloud physics is used, arrays for storing potential temperature and
715	!-- specific humidity at first grid level are required.
716	IF ( cloud_physics ) THEN
717	ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
718	ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
719	ENDIF
720	!
721	!--
722	IF ( cloud_physics .AND. microphysics_seifert) THEN
723	ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
724	ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
725	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
726	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
727	ENDIF
728	!
729	!-- Salinity surface flux
730	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
731
732	END SUBROUTINE allocate_surface_attributes_h
733
734
735	!------------------------------------------------------------------------------!
736	! Description:
737	! ------------
738	!> Allocating memory for model-top fluxes
739	!------------------------------------------------------------------------------!
740	SUBROUTINE allocate_surface_attributes_h_top( surfaces, &
741	nys_l, nyn_l, nxl_l, nxr_l )
742
743	IMPLICIT NONE
744
745	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
746	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
747	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
748	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
749
750	TYPE(surf_type) :: surfaces !< respective surface type
751
752	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
753	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
754	surfaces%start_index = 0
755	surfaces%end_index = -1
756	!
757	!-- Indices to locate surface (model-top) element
758	ALLOCATE ( surfaces%i(1:surfaces%ns) )
759	ALLOCATE ( surfaces%j(1:surfaces%ns) )
760	ALLOCATE ( surfaces%k(1:surfaces%ns) )
761	!
762	!-- Vertical momentum fluxes of u and v
763	ALLOCATE ( surfaces%usws(1:surfaces%ns) )
764	ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
765	!
766	!-- Sensible heat flux
767	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
768	!
769	!-- Latent heat flux
770	IF ( humidity ) THEN
771	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
772	ENDIF
773	!
774	!-- Scalar flux
775	IF ( passive_scalar ) THEN
776	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
777	ENDIF
778	!
779	!--
780	IF ( cloud_physics .AND. microphysics_seifert) THEN
781	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
782	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
783	ENDIF
784	!
785	!-- Salinity flux
786	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
787
788	END SUBROUTINE allocate_surface_attributes_h_top
789
790	!------------------------------------------------------------------------------!
791	! Description:
792	! ------------
793	!> Allocating memory for vertical surface types.
794	!------------------------------------------------------------------------------!
795	SUBROUTINE allocate_surface_attributes_v( surfaces, lsm, &
796	nys_l, nyn_l, nxl_l, nxr_l )
797
798	IMPLICIT NONE
799
800	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
801	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
802	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
803	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
804
805	LOGICAL :: lsm !< flag indicating data type of natural land surface
806
807	TYPE(surf_type) :: surfaces !< respective surface type
808
809	!
810	!-- Allocate arrays for start and end index of vertical surface type
811	!-- for each (j,i)-grid point. This is required in diffion_x, which is
812	!-- called for each (j,i). In order to find the location where the
813	!-- respective flux is store within the surface-type, start- and end-
814	!-- index are stored for each (j,i). For example, each (j,i) can have
815	!-- several entries where fluxes for vertical surfaces might be stored.
816	!-- In the flat case, where no vertical walls exit, set indicies such
817	!-- that loop in diffusion routines will not be entered.
818	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
819	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
820	surfaces%start_index = 0
821	surfaces%end_index = -1
822	!
823	!-- Indices to locate surface element.
824	ALLOCATE ( surfaces%i(1:surfaces%ns) )
825	ALLOCATE ( surfaces%j(1:surfaces%ns) )
826	ALLOCATE ( surfaces%k(1:surfaces%ns) )
827	!
828	!-- Surface-layer height
829	ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
830	!
831	!-- Surface orientation
832	ALLOCATE ( surfaces%facing(1:surfaces%ns) )
833	!
834	!-- Surface parallel wind velocity
835	ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
836	!
837	!-- Roughness
838	ALLOCATE ( surfaces%z0(1:surfaces%ns) )
839	ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
840	ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
841
842	!
843	!-- Friction velocity
844	ALLOCATE ( surfaces%us(1:surfaces%ns) )
845	!
846	!-- Allocate Obukhov length and bulk Richardson number. Only required
847	!-- for natural land surfaces
848	IF ( lsm ) THEN
849	ALLOCATE( surfaces%ol(1:surfaces%ns) )
850	ALLOCATE( surfaces%rib(1:surfaces%ns) )
851	ENDIF
852	!
853	!-- Allocate arrays for surface momentum fluxes for u and v. For u at north-
854	!-- and south-facing surfaces, for v at east- and west-facing surfaces.
855	ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
856	!
857	!-- Allocate array for surface momentum flux for w - wsus and wsvs
858	ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) )
859	!
860	!-- Allocate array for surface momentum flux for subgrid-scale tke wsus and
861	!-- wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
862	!-- on surface.
863	ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )
864	!
865	!-- Characteristic temperature and surface flux of sensible heat
866	ALLOCATE ( surfaces%ts(1:surfaces%ns) )
867	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
868	!
869	!-- Characteristic humidity and surface flux of latent heat
870	IF ( humidity ) THEN
871	ALLOCATE ( surfaces%qs(1:surfaces%ns) )
872	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
873	ENDIF
874	!
875	!-- Characteristic scalar and surface flux of scalar
876	IF ( passive_scalar ) THEN
877	ALLOCATE ( surfaces%ss(1:surfaces%ns) )
878	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
879	ENDIF
880
881	IF ( cloud_physics .AND. microphysics_seifert) THEN
882	ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
883	ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
884	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
885	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
886	ENDIF
887	!
888	!-- Salinity surface flux
889	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
890
891	END SUBROUTINE allocate_surface_attributes_v
892
893	!------------------------------------------------------------------------------!
894	! Description:
895	! ------------
896	!> Initialize surface elements.
897	!------------------------------------------------------------------------------!
898	SUBROUTINE init_surfaces
899
900	IMPLICIT NONE
901
902	INTEGER(iwp) :: i !< running index x-direction
903	INTEGER(iwp) :: j !< running index y-direction
904	INTEGER(iwp) :: k !< running index z-direction
905	INTEGER(iwp) :: l !< index variable used to distinguish surface facing
906	INTEGER(iwp) :: m !< running index surface elements
907
908	INTEGER(iwp) :: start_index_lsm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal natural surfaces
909	INTEGER(iwp) :: start_index_usm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal urban surfaces
910
911	INTEGER(iwp) :: num_lsm_h !< current number of horizontal surface element, natural type
912	INTEGER(iwp) :: num_lsm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal natural surfaces
913	INTEGER(iwp) :: num_usm_h !< current number of horizontal surface element, urban type
914	INTEGER(iwp) :: num_usm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban surfaces
915
916	INTEGER(iwp), DIMENSION(0:2) :: num_def_h !< current number of horizontal surface element, default type
917	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
918	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
919
920	INTEGER(iwp), DIMENSION(0:3) :: num_def_v !< current number of vertical surface element, default type
921	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
922	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< current number of vertical surface element, natural type
923	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
924	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< current number of vertical surface element, urban type
925	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
926
927	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
928	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
929	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
930
931
932	!
933	!-- Initialize surface attributes, store indicies, surfaces orientation, etc.,
934	num_def_h(0:2) = 1
935	num_def_v(0:3) = 1
936
937	num_lsm_h = 1
938	num_lsm_v(0:3) = 1
939
940	num_usm_h = 1
941	num_usm_v(0:3) = 1
942
943	start_index_def_h(0:2) = 1
944	start_index_def_v(0:3) = 1
945
946	start_index_lsm_h = 1
947	start_index_lsm_v(0:3) = 1
948
949	start_index_usm_h = 1
950	start_index_usm_v(0:3) = 1
951
952	DO i = nxl, nxr
953	DO j = nys, nyn
954
955	num_def_h_kji = 0
956	num_def_v_kji = 0
957	num_lsm_h_kji = 0
958	num_lsm_v_kji = 0
959	num_usm_h_kji = 0
960	num_usm_v_kji = 0
961
962	DO k = nzb+1, nzt
963	!
964	!-- Check if current gridpoint belongs to the atmosphere
965	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
966	!
967	!-- Upward-facing surface. Distinguish between differet surface types.
968	!-- To do, think about method to flag natural and non-natural
969	!-- surfaces. Only to ask for land_surface or urban surface
970	!-- is just a work-around.
971	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
972	!
973	!-- Natural surface type
974	IF ( land_surface ) THEN
975	CALL initialize_horizontal_surfaces( k, j, i, &
976	surf_lsm_h, &
977	num_lsm_h, &
978	num_lsm_h_kji, &
979	.TRUE., .FALSE. )
980	!
981	!-- Urban surface tpye
982	ELSEIF ( urban_surface ) THEN
983	CALL initialize_horizontal_surfaces( k, j, i, &
984	surf_usm_h, &
985	num_usm_h, &
986	num_usm_h_kji, &
987	.TRUE., .FALSE. )
988	!
989	!-- Default surface type
990	ELSE
991	CALL initialize_horizontal_surfaces( k, j, i, &
992	surf_def_h(0), &
993	num_def_h(0), &
994	num_def_h_kji(0),&
995	.TRUE., .FALSE. )
996	ENDIF
997	ENDIF
998	!
999	!-- downward-facing surface, first, model top
1000	IF ( k == nzt .AND. use_top_fluxes ) THEN
1001	CALL initialize_top( k, j, i, surf_def_h(2), &
1002	num_def_h(2), num_def_h_kji(2) )
1003	!
1004	!-- Check for any other downward-facing surface. So far only for
1005	!-- default surface type.
1006	ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
1007	CALL initialize_horizontal_surfaces( k, j, i, &
1008	surf_def_h(1), &
1009	num_def_h(1), &
1010	num_def_h_kji(1), &
1011	.FALSE., .TRUE. )
1012	ENDIF
1013	!
1014	!-- Check for vertical walls and, if required, initialize it.
1015	! Start with northward-facing surface.
1016	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) ) THEN
1017	IF ( urban_surface ) THEN
1018	CALL initialize_vertical_surfaces( 0, k, j, i, &
1019	surf_usm_v(0), &
1020	num_usm_v(0), &
1021	num_usm_v_kji(0), &
1022	.FALSE., .FALSE., &
1023	.FALSE., .TRUE. )
1024	ELSEIF ( land_surface ) THEN
1025	CALL initialize_vertical_surfaces( 0, k, j, i, &
1026	surf_lsm_v(0), &
1027	num_lsm_v(0), &
1028	num_lsm_v_kji(0), &
1029	.FALSE., .FALSE., &
1030	.FALSE., .TRUE. )
1031	ELSE
1032	CALL initialize_vertical_surfaces( 0, k, j, i, &
1033	surf_def_v(0), &
1034	num_def_v(0), &
1035	num_def_v_kji(0), &
1036	.FALSE., .FALSE., &
1037	.FALSE., .TRUE. )
1038	ENDIF
1039	ENDIF
1040	!
1041	!-- southward-facing surface
1042	IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) THEN
1043	IF ( urban_surface ) THEN
1044	CALL initialize_vertical_surfaces( 1, k, j, i, &
1045	surf_usm_v(1), &
1046	num_usm_v(1), &
1047	num_usm_v_kji(1), &
1048	.FALSE., .FALSE., &
1049	.TRUE., .FALSE. )
1050	ELSEIF ( land_surface ) THEN
1051	CALL initialize_vertical_surfaces( 1, k, j, i, &
1052	surf_lsm_v(1), &
1053	num_lsm_v(1), &
1054	num_lsm_v_kji(1), &
1055	.FALSE., .FALSE., &
1056	.TRUE., .FALSE. )
1057	ELSE
1058	CALL initialize_vertical_surfaces( 1, k, j, i, &
1059	surf_def_v(1), &
1060	num_def_v(1), &
1061	num_def_v_kji(1), &
1062	.FALSE., .FALSE., &
1063	.TRUE., .FALSE. )
1064	ENDIF
1065	ENDIF
1066	!
1067	!-- eastward-facing surface
1068	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) ) THEN
1069	IF ( urban_surface ) THEN
1070	CALL initialize_vertical_surfaces( 2, k, j, i, &
1071	surf_usm_v(2), &
1072	num_usm_v(2), &
1073	num_usm_v_kji(2), &
1074	.TRUE., .FALSE., &
1075	.FALSE., .FALSE. )
1076	ELSEIF ( land_surface ) THEN
1077	CALL initialize_vertical_surfaces( 2, k, j, i, &
1078	surf_lsm_v(2), &
1079	num_lsm_v(2), &
1080	num_lsm_v_kji(2), &
1081	.TRUE., .FALSE., &
1082	.FALSE., .FALSE. )
1083	ELSE
1084	CALL initialize_vertical_surfaces( 2, k, j, i, &
1085	surf_def_v(2), &
1086	num_def_v(2), &
1087	num_def_v_kji(2), &
1088	.TRUE., .FALSE., &
1089	.FALSE., .FALSE. )
1090	ENDIF
1091	ENDIF
1092	!
1093	!-- westward-facing surface
1094	IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) THEN
1095	IF ( urban_surface ) THEN
1096	CALL initialize_vertical_surfaces( 3, k, j, i, &
1097	surf_usm_v(3), &
1098	num_usm_v(3), &
1099	num_usm_v_kji(3), &
1100	.FALSE., .TRUE., &
1101	.FALSE., .FALSE. )
1102	ELSEIF ( land_surface ) THEN
1103	CALL initialize_vertical_surfaces( 3, k, j, i, &
1104	surf_lsm_v(3), &
1105	num_lsm_v(3), &
1106	num_lsm_v_kji(3), &
1107	.FALSE., .TRUE., &
1108	.FALSE., .FALSE. )
1109	ELSE
1110	CALL initialize_vertical_surfaces( 3, k, j, i, &
1111	surf_def_v(3), &
1112	num_def_v(3), &
1113	num_def_v_kji(3), &
1114	.FALSE., .TRUE., &
1115	.FALSE., .FALSE. )
1116	ENDIF
1117	ENDIF
1118	ENDIF
1119
1120
1121	ENDDO
1122	!
1123	!-- Determine start- and end-index at grid point (j,i). Also, for
1124	!-- horizontal surfaces more than 1 horizontal surface element can
1125	!-- exist at grid point (j,i) if overhanging structures are present.
1126	!-- Upward-facing surfaces
1127	surf_def_h(0)%start_index(j,i) = start_index_def_h(0)
1128	surf_def_h(0)%end_index(j,i) = surf_def_h(0)%start_index(j,i) + &
1129	num_def_h_kji(0) - 1
1130	start_index_def_h(0) = surf_def_h(0)%end_index(j,i) + 1
1131	!
1132	!-- Downward-facing surfaces, except model top
1133	surf_def_h(1)%start_index(j,i) = start_index_def_h(1)
1134	surf_def_h(1)%end_index(j,i) = surf_def_h(1)%start_index(j,i) + &
1135	num_def_h_kji(1) - 1
1136	start_index_def_h(1) = surf_def_h(1)%end_index(j,i) + 1
1137	!
1138	!-- Downward-facing surfaces -- model top fluxes
1139	surf_def_h(2)%start_index(j,i) = start_index_def_h(2)
1140	surf_def_h(2)%end_index(j,i) = surf_def_h(2)%start_index(j,i) + &
1141	num_def_h_kji(2) - 1
1142	start_index_def_h(2) = surf_def_h(2)%end_index(j,i) + 1
1143	!
1144	!-- Horizontal natural land surfaces
1145	surf_lsm_h%start_index(j,i) = start_index_lsm_h
1146	surf_lsm_h%end_index(j,i) = surf_lsm_h%start_index(j,i) + &
1147	num_lsm_h_kji - 1
1148	start_index_lsm_h = surf_lsm_h%end_index(j,i) + 1
1149	!
1150	!-- Horizontal urban surfaces
1151	surf_usm_h%start_index(j,i) = start_index_usm_h
1152	surf_usm_h%end_index(j,i) = surf_usm_h%start_index(j,i) + &
1153	num_usm_h_kji - 1
1154	start_index_usm_h = surf_usm_h%end_index(j,i) + 1
1155
1156	!
1157	!-- Vertical surfaces - Default type
1158	surf_def_v(0)%start_index(j,i) = start_index_def_v(0)
1159	surf_def_v(1)%start_index(j,i) = start_index_def_v(1)
1160	surf_def_v(2)%start_index(j,i) = start_index_def_v(2)
1161	surf_def_v(3)%start_index(j,i) = start_index_def_v(3)
1162	surf_def_v(0)%end_index(j,i) = start_index_def_v(0) + &
1163	num_def_v_kji(0) - 1
1164	surf_def_v(1)%end_index(j,i) = start_index_def_v(1) + &
1165	num_def_v_kji(1) - 1
1166	surf_def_v(2)%end_index(j,i) = start_index_def_v(2) + &
1167	num_def_v_kji(2) - 1
1168	surf_def_v(3)%end_index(j,i) = start_index_def_v(3) + &
1169	num_def_v_kji(3) - 1
1170	start_index_def_v(0) = surf_def_v(0)%end_index(j,i) + 1
1171	start_index_def_v(1) = surf_def_v(1)%end_index(j,i) + 1
1172	start_index_def_v(2) = surf_def_v(2)%end_index(j,i) + 1
1173	start_index_def_v(3) = surf_def_v(3)%end_index(j,i) + 1
1174	!
1175	!-- Natural type
1176	surf_lsm_v(0)%start_index(j,i) = start_index_lsm_v(0)
1177	surf_lsm_v(1)%start_index(j,i) = start_index_lsm_v(1)
1178	surf_lsm_v(2)%start_index(j,i) = start_index_lsm_v(2)
1179	surf_lsm_v(3)%start_index(j,i) = start_index_lsm_v(3)
1180	surf_lsm_v(0)%end_index(j,i) = start_index_lsm_v(0) + &
1181	num_lsm_v_kji(0) - 1
1182	surf_lsm_v(1)%end_index(j,i) = start_index_lsm_v(1) + &
1183	num_lsm_v_kji(1) - 1
1184	surf_lsm_v(2)%end_index(j,i) = start_index_lsm_v(2) + &
1185	num_lsm_v_kji(2) - 1
1186	surf_lsm_v(3)%end_index(j,i) = start_index_lsm_v(3) + &
1187	num_lsm_v_kji(3) - 1
1188	start_index_lsm_v(0) = surf_lsm_v(0)%end_index(j,i) + 1
1189	start_index_lsm_v(1) = surf_lsm_v(1)%end_index(j,i) + 1
1190	start_index_lsm_v(2) = surf_lsm_v(2)%end_index(j,i) + 1
1191	start_index_lsm_v(3) = surf_lsm_v(3)%end_index(j,i) + 1
1192	!
1193	!-- Urban type
1194	surf_usm_v(0)%start_index(j,i) = start_index_usm_v(0)
1195	surf_usm_v(1)%start_index(j,i) = start_index_usm_v(1)
1196	surf_usm_v(2)%start_index(j,i) = start_index_usm_v(2)
1197	surf_usm_v(3)%start_index(j,i) = start_index_usm_v(3)
1198	surf_usm_v(0)%end_index(j,i) = start_index_usm_v(0) + &
1199	num_usm_v_kji(0) - 1
1200	surf_usm_v(1)%end_index(j,i) = start_index_usm_v(1) + &
1201	num_usm_v_kji(1) - 1
1202	surf_usm_v(2)%end_index(j,i) = start_index_usm_v(2) + &
1203	num_usm_v_kji(2) - 1
1204	surf_usm_v(3)%end_index(j,i) = start_index_usm_v(3) + &
1205	num_usm_v_kji(3) - 1
1206	start_index_usm_v(0) = surf_usm_v(0)%end_index(j,i) + 1
1207	start_index_usm_v(1) = surf_usm_v(1)%end_index(j,i) + 1
1208	start_index_usm_v(2) = surf_usm_v(2)%end_index(j,i) + 1
1209	start_index_usm_v(3) = surf_usm_v(3)%end_index(j,i) + 1
1210
1211
1212	ENDDO
1213	ENDDO
1214
1215	CONTAINS
1216
1217	!------------------------------------------------------------------------------!
1218	! Description:
1219	! ------------
1220	!> Initialize horizontal surface elements, upward- and downward-facing.
1221	!> Note, horizontal surface type alsw comprises model-top fluxes, which are,
1222	!> initialized in a different routine.
1223	!------------------------------------------------------------------------------!
1224	SUBROUTINE initialize_horizontal_surfaces( k, j, i, surf, num_h, &
1225	num_h_kji, upward_facing, &
1226	downward_facing )
1227
1228	IMPLICIT NONE
1229
1230	INTEGER(iwp) :: i !< running index x-direction
1231	INTEGER(iwp) :: j !< running index y-direction
1232	INTEGER(iwp) :: k !< running index z-direction
1233	INTEGER(iwp) :: num_h !< current number of surface element
1234	INTEGER(iwp) :: num_h_kji !< dummy increment
1235
1236	LOGICAL :: upward_facing !< flag indicating upward-facing surface
1237	LOGICAL :: downward_facing !< flag indicating downward-facing surface
1238
1239	TYPE( surf_type ) :: surf !< respective surface type
1240	!
1241	!-- Store indices of respective surface element
1242	surf%i(num_h) = i
1243	surf%j(num_h) = j
1244	surf%k(num_h) = k
1245	!
1246	!-- Surface orientation, bit 0 is set to 1 for upward-facing surfaces,
1247	!-- bit 1 is for downward-facing surfaces.
1248	IF ( upward_facing ) surf%facing(num_h) = IBSET( surf%facing(num_h), 0 )
1249	IF ( downward_facing ) surf%facing(num_h) = IBSET( surf%facing(num_h), 1 )
1250	!
1251	!-- Initialize surface-layer height
1252	IF ( upward_facing ) THEN
1253	surf%z_mo(num_h) = zu(k) - zw(k-1)
1254	ELSE
1255	surf%z_mo(num_h) = zw(k) - zu(k)
1256	ENDIF
1257
1258	surf%z0(num_h) = roughness_length
1259	surf%z0h(num_h) = z0h_factor * roughness_length
1260	surf%z0q(num_h) = z0h_factor * roughness_length
1261	!
1262	!-- Initialization in case of 1D pre-cursor run
1263	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )&
1264	THEN
1265	IF ( .NOT. constant_diffusion ) THEN
1266	IF ( constant_flux_layer ) THEN
1267	surf%ol(num_h) = surf%z_mo(num_h) / &
1268	( rif1d(nzb+1) + 1.0E-20_wp )
1269	surf%us(num_h) = us1d
1270	surf%usws(num_h) = usws1d
1271	surf%vsws(num_h) = vsws1d
1272	ELSE
1273	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
1274	surf%us(num_h) = 0.0_wp
1275	surf%usws(num_h) = 0.0_wp
1276	surf%vsws(num_h) = 0.0_wp
1277	ENDIF
1278	ELSE
1279	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
1280	surf%us(num_h) = 0.0_wp
1281	surf%usws(num_h) = 0.0_wp
1282	surf%vsws(num_h) = 0.0_wp
1283	ENDIF
1284	!
1285	!-- Initialization in case of constant profiles
1286	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )&
1287	THEN
1288
1289	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
1290	!
1291	!-- Very small number is required for calculation of Obukhov length
1292	!-- at first timestep
1293	surf%us(num_h) = 1E-30_wp
1294	surf%usws(num_h) = 0.0_wp
1295	surf%vsws(num_h) = 0.0_wp
1296
1297	ENDIF
1298
1299	surf%rib(num_h) = 0.0_wp
1300	surf%uvw_abs(num_h) = 0.0_wp
1301
1302	IF ( .NOT. constant_diffusion ) THEN
1303	surf%u_0(num_h) = 0.0_wp
1304	surf%v_0(num_h) = 0.0_wp
1305	ENDIF
1306
1307	surf%ts(num_h) = 0.0_wp
1308
1309	IF ( humidity ) THEN
1310	surf%qs(num_h) = 0.0_wp
1311	IF ( cloud_physics .AND. microphysics_seifert) THEN
1312	surf%qrs(num_h) = 0.0_wp
1313	surf%nrs(num_h) = 0.0_wp
1314
1315	surf%qrsws(num_h) = 0.0_wp
1316	surf%nrsws(num_h) = 0.0_wp
1317
1318	surf%pt1(num_h) = 0.0_wp
1319	surf%qv1(num_h) = 0.0_wp
1320
1321	ENDIF
1322	ENDIF
1323
1324	IF ( passive_scalar ) surf%ss(num_h) = 0.0_wp
1325	!
1326	!-- Inititalize surface fluxes of sensible and latent heat, as well as
1327	!-- passive scalar
1328	IF ( use_surface_fluxes ) THEN
1329
1330	IF ( upward_facing ) THEN
1331	IF ( constant_heatflux ) THEN
1332	!
1333	!-- Initialize surface heatflux. However, skip this for now if
1334	!-- if random_heatflux is set. This case, shf is initialized later.
1335	IF ( .NOT. random_heatflux ) THEN
1336	surf%shf(num_h) = surface_heatflux * &
1337	heatflux_input_conversion(nzb)
1338	!
1339	!-- Check if surface heat flux might be replaced by
1340	!-- prescribed wall heatflux
1341	IF ( k-1 /= 0 ) THEN
1342	surf%shf(num_h) = wall_heatflux(0) * &
1343	heatflux_input_conversion(k-1)
1344	ENDIF
1345	!
1346	!-- Initialize shf with data from external file LSF_DATA. Will
1347	!-- be done directly in ls_foring_surf
1348	!-- Attention: Just a workaround, need to be revised!!!
1349	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1350	! CALL ls_forcing_surf ( simulated_time )
1351	! surf%shf(num_h) = shf(j,i)
1352	ENDIF
1353	ENDIF
1354	ELSE
1355	surf%shf(num_h) = 0.0_wp
1356	ENDIF
1357	!
1358	!-- Set heat-flux at downward-facing surfaces
1359	ELSE
1360	surf%shf(num_h) = wall_heatflux(5) * &
1361	heatflux_input_conversion(k)
1362	ENDIF
1363
1364	IF ( humidity ) THEN
1365	IF ( upward_facing ) THEN
1366	IF ( constant_waterflux ) THEN
1367	surf%qsws(num_h) = surface_waterflux * &
1368	waterflux_input_conversion(nzb)
1369	IF ( k-1 /= 0 ) THEN
1370	surf%qsws(num_h) = wall_humidityflux(0) * &
1371	waterflux_input_conversion(k-1)
1372	ENDIF
1373	ELSE
1374	surf%qsws(num_h) = 0.0_wp
1375	ENDIF
1376	ELSE
1377	surf%qsws(num_h) = wall_humidityflux(5) * &
1378	heatflux_input_conversion(k)
1379	ENDIF
1380	ENDIF
1381
1382	IF ( passive_scalar ) THEN
1383	IF ( upward_facing ) THEN
1384	IF ( constant_scalarflux ) THEN
1385	surf%ssws(num_h) = surface_scalarflux
1386
1387	IF ( k-1 /= 0 ) &
1388	surf%ssws(num_h) = wall_scalarflux(0)
1389
1390	ELSE
1391	surf%ssws(num_h) = 0.0_wp
1392	ENDIF
1393	ELSE
1394	surf%ssws(num_h) = wall_scalarflux(5)
1395	ENDIF
1396	ENDIF
1397
1398	IF ( ocean ) THEN
1399	IF ( upward_facing ) THEN
1400	surf%sasws(num_h) = bottom_salinityflux
1401	ELSE
1402	surf%sasws(num_h) = 0.0_wp
1403	ENDIF
1404	ENDIF
1405	ENDIF
1406	!
1407	!-- Increment surface indices
1408	num_h = num_h + 1
1409	num_h_kji = num_h_kji + 1
1410
1411
1412	END SUBROUTINE initialize_horizontal_surfaces
1413
1414
1415	!------------------------------------------------------------------------------!
1416	! Description:
1417	! ------------
1418	!> Initialize model-top fluxes. Currently, only the heatflux and salinity flux
1419	!> can be prescribed, latent flux is zero in this case!
1420	!------------------------------------------------------------------------------!
1421	SUBROUTINE initialize_top( k, j, i, surf, num_h, num_h_kji )
1422
1423	IMPLICIT NONE
1424
1425	INTEGER(iwp) :: i !< running index x-direction
1426	INTEGER(iwp) :: j !< running index y-direction
1427	INTEGER(iwp) :: k !< running index z-direction
1428	INTEGER(iwp) :: num_h !< current number of surface element
1429	INTEGER(iwp) :: num_h_kji !< dummy increment
1430
1431	TYPE( surf_type ) :: surf !< respective surface type
1432	!
1433	!-- Store indices of respective surface element
1434	surf%i(num_h) = i
1435	surf%j(num_h) = j
1436	surf%k(num_h) = k
1437	!
1438	!-- Initialize top heat flux
1439	IF ( constant_top_heatflux ) &
1440	surf%shf = top_heatflux * heatflux_input_conversion(nzt+1)
1441	!
1442	!-- Initialization in case of a coupled model run
1443	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
1444	surf%shf = 0.0_wp
1445	ENDIF
1446	!
1447	!-- Prescribe latent heat flux at the top
1448	IF ( humidity ) THEN
1449	surf%qsws = 0.0_wp
1450	IF ( cloud_physics .AND. microphysics_seifert ) THEN
1451	surf%nrsws = 0.0_wp
1452	surf%qrsws = 0.0_wp
1453	ENDIF
1454	ENDIF
1455	!
1456	!-- Prescribe top scalar flux
1457	IF ( passive_scalar .AND. constant_top_scalarflux ) &
1458	surf%ssws = top_scalarflux
1459	!
1460	!-- Prescribe top salinity flux
1461	IF ( ocean .AND. constant_top_salinityflux) &
1462	surf%sasws = top_salinityflux
1463	!
1464	!-- Initialization in case of a coupled model run
1465	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
1466	surf%shf = 0.0_wp
1467	ENDIF
1468	!
1469	!-- Top momentum fluxes
1470	surf%usws = top_momentumflux_u * momentumflux_input_conversion(nzt+1)
1471	surf%vsws = top_momentumflux_v * momentumflux_input_conversion(nzt+1)
1472	!
1473	!-- Increment surface indices
1474	num_h = num_h + 1
1475	num_h_kji = num_h_kji + 1
1476
1477
1478	END SUBROUTINE initialize_top
1479
1480
1481	!------------------------------------------------------------------------------!
1482	! Description:
1483	! ------------
1484	!> Initialize vertical surface elements.
1485	!------------------------------------------------------------------------------!
1486	SUBROUTINE initialize_vertical_surfaces( l, k, j, i, surf, num_v, &
1487	num_v_kji, east_facing, &
1488	west_facing, south_facing, &
1489	north_facing )
1490
1491	IMPLICIT NONE
1492
1493	INTEGER(iwp) :: component !<
1494	INTEGER(iwp) :: i !< running index x-direction
1495	INTEGER(iwp) :: j !< running index x-direction
1496	INTEGER(iwp) :: k !< running index x-direction
1497	INTEGER(iwp) :: l !< index variable for the surface type, indicating the facing
1498	INTEGER(iwp) :: num_v !< current number of surface element
1499	INTEGER(iwp) :: num_v_kji !< current number of surface element at (j,i)
1500
1501
1502	LOGICAL :: east_facing !< flag indicating east-facing surfaces
1503	LOGICAL :: north_facing !< flag indicating north-facing surfaces
1504	LOGICAL :: south_facing !< flag indicating south-facing surfaces
1505	LOGICAL :: west_facing !< flag indicating west-facing surfaces
1506
1507	TYPE( surf_type ) :: surf !< respective surface type
1508
1509	!
1510	!-- Store indices of respective wall element
1511	surf%i(num_v) = i
1512	surf%j(num_v) = j
1513	surf%k(num_v) = k
1514	!
1515	!-- Initialize surface-layer height, or more precisely, distance to surface
1516	IF ( north_facing .OR. south_facing ) THEN
1517	surf%z_mo(num_v) = 0.5_wp * dy
1518	ELSE
1519	surf%z_mo(num_v) = 0.5_wp * dx
1520	ENDIF
1521
1522	surf%facing(num_v) = 0
1523	!
1524	!-- Surface orientation. Moreover, set component id to map wall_heatflux,
1525	!-- etc., on surface type (further below)
1526	IF ( north_facing ) THEN
1527	surf%facing(num_v) = IBSET( surf%facing(num_v), 0 )
1528	component = 4
1529	ENDIF
1530
1531	IF ( south_facing ) THEN
1532	surf%facing(num_v) = IBSET( surf%facing(num_v), 1 )
1533	component = 3
1534	ENDIF
1535
1536	IF ( east_facing ) THEN
1537	surf%facing(num_v) = IBSET( surf%facing(num_v), 2 )
1538	component = 2
1539	ENDIF
1540
1541	IF ( west_facing ) THEN
1542	surf%facing(num_v) = IBSET( surf%facing(num_v), 3 )
1543	component = 1
1544	ENDIF
1545
1546
1547	surf%z0(num_v) = roughness_length
1548	surf%z0h(num_v) = z0h_factor * roughness_length
1549	surf%z0q(num_v) = z0h_factor * roughness_length
1550
1551	surf%us(num_v) = 0.0_wp
1552	!
1553	!-- If required, initialize Obukhov length
1554	IF ( ALLOCATED( surf%ol ) ) &
1555	surf%ol(num_v) = surf%z_mo(num_v) / zeta_min
1556
1557	surf%uvw_abs(num_v) = 0.0_wp
1558
1559	surf%mom_flux_uv(num_v) = 0.0_wp
1560	surf%mom_flux_w(num_v) = 0.0_wp
1561	surf%mom_flux_tke(0:1,num_v) = 0.0_wp
1562
1563	surf%ts(num_v) = 0.0_wp
1564	surf%shf(num_v) = wall_heatflux(component)
1565
1566	IF ( humidity ) THEN
1567	surf%qs(num_v) = 0.0_wp
1568	surf%qsws(num_v) = wall_humidityflux(component)
1569	!
1570	!-- Following wall fluxes are assumed to be zero
1571	IF ( cloud_physics .AND. microphysics_seifert) THEN
1572	surf%qrs(num_v) = 0.0_wp
1573	surf%nrs(num_v) = 0.0_wp
1574
1575	surf%qrsws(num_v) = 0.0_wp
1576	surf%nrsws(num_v) = 0.0_wp
1577	ENDIF
1578	ENDIF
1579
1580	IF ( passive_scalar ) THEN
1581	surf%ss(num_v) = 0.0_wp
1582	surf%ssws(num_v) = wall_scalarflux(component)
1583	ENDIF
1584	!
1585	!-- So far, salinityflux at vertical surfaces is simply zero
1586	!-- at the moment
1587	IF ( ocean ) surf%sasws(num_v) = wall_salinityflux(component)
1588	!
1589	!-- Increment wall indices
1590	num_v = num_v + 1
1591	num_v_kji = num_v_kji + 1
1592
1593	END SUBROUTINE initialize_vertical_surfaces
1594
1595	END SUBROUTINE init_surfaces
1596
1597	!------------------------------------------------------------------------------!
1598	! Description:
1599	! ------------
1600	!> Gathers all surface elements with the same facing (but possibly different
1601	!> type) onto a surface type, and writes binary data into restart files.
1602	!------------------------------------------------------------------------------!
1603	SUBROUTINE surface_write_restart_data
1604
1605	IMPLICIT NONE
1606
1607	CHARACTER(LEN=1) :: dum !< dummy string to create output-variable name
1608
1609	INTEGER(iwp) :: i !< running index x-direction
1610	INTEGER(iwp) :: j !< running index y-direction
1611	INTEGER(iwp) :: l !< index surface type orientation
1612	INTEGER(iwp) :: m !< running index for surface elements on individual surface array
1613	INTEGER(iwp), DIMENSION(0:3) :: mm !< running index for surface elements on gathered surface array
1614
1615	TYPE(surf_type), DIMENSION(0:2) :: surf_h !< gathered horizontal surfaces, contains all surface types
1616	TYPE(surf_type), DIMENSION(0:3) :: surf_v !< gathered vertical surfaces, contains all surface types
1617
1618	!
1619	!-- Determine total number of horizontal and vertical surface elements before
1620	!-- writing var_list
1621	CALL surface_last_actions
1622	!
1623	!-- Count number of grid points with same facing and allocate attributes respectively
1624	!-- Horizontal upward facing
1625	surf_h(0)%ns = ns_h_on_file(0)
1626	CALL allocate_surface_attributes_h( surf_h(0), nys, nyn, nxl, nxr )
1627	!
1628	!-- Horizontal downward facing
1629	surf_h(1)%ns = ns_h_on_file(1)
1630	CALL allocate_surface_attributes_h( surf_h(1), nys, nyn, nxl, nxr )
1631	!
1632	!-- Model top
1633	surf_h(2)%ns = ns_h_on_file(2)
1634	CALL allocate_surface_attributes_h_top( surf_h(2), nys, nyn, nxl, nxr )
1635	!
1636	!-- Vertical surfaces
1637	DO l = 0, 3
1638	surf_v(l)%ns = ns_v_on_file(l)
1639	CALL allocate_surface_attributes_v( surf_v(l), .FALSE., &
1640	nys, nyn, nxl, nxr )
1641	ENDDO
1642	!
1643	!-- In the following, gather data from surfaces elements with the same
1644	!-- facing (but possibly differt type) on 1 data-type array.
1645	mm(0:2) = 1
1646	DO l = 0, 2
1647	DO i = nxl, nxr
1648	DO j = nys, nyn
1649	DO m = surf_def_h(l)%start_index(j,i), &
1650	surf_def_h(l)%end_index(j,i)
1651	IF ( ALLOCATED( surf_def_h(l)%us ) ) &
1652	surf_h(l)%us(mm(l)) = surf_def_h(l)%us(m)
1653	IF ( ALLOCATED( surf_def_h(l)%ts ) ) &
1654	surf_h(l)%ts(mm(l)) = surf_def_h(l)%ts(m)
1655	IF ( ALLOCATED( surf_def_h(l)%qs ) ) &
1656	surf_h(l)%qs(mm(l)) = surf_def_h(l)%qs(m)
1657	IF ( ALLOCATED( surf_def_h(l)%ss ) ) &
1658	surf_h(l)%ss(mm(l)) = surf_def_h(l)%ss(m)
1659	IF ( ALLOCATED( surf_def_h(l)%qrs ) ) &
1660	surf_h(l)%qrs(mm(l)) = surf_def_h(l)%qrs(m)
1661	IF ( ALLOCATED( surf_def_h(l)%nrs ) ) &
1662	surf_h(l)%nrs(mm(l)) = surf_def_h(l)%nrs(m)
1663	IF ( ALLOCATED( surf_def_h(l)%ol ) ) &
1664	surf_h(l)%ol(mm(l)) = surf_def_h(l)%ol(m)
1665	IF ( ALLOCATED( surf_def_h(l)%rib ) ) &
1666	surf_h(l)%rib(mm(l)) = surf_def_h(l)%rib(m)
1667	IF ( ALLOCATED( surf_def_h(l)%usws ) ) &
1668	surf_h(l)%usws(mm(l)) = surf_def_h(l)%usws(m)
1669	IF ( ALLOCATED( surf_def_h(l)%vsws ) ) &
1670	surf_h(l)%vsws(mm(l)) = surf_def_h(l)%vsws(m)
1671	IF ( ALLOCATED( surf_def_h(l)%shf ) ) &
1672	surf_h(l)%shf(mm(l)) = surf_def_h(l)%shf(m)
1673	IF ( ALLOCATED( surf_def_h(l)%qsws ) ) &
1674	surf_h(l)%qsws(mm(l)) = surf_def_h(l)%qsws(m)
1675	IF ( ALLOCATED( surf_def_h(l)%ssws ) ) &
1676	surf_h(l)%qsws(mm(l)) = surf_def_h(l)%ssws(m)
1677	IF ( ALLOCATED( surf_def_h(l)%nrsws ) ) &
1678	surf_h(l)%nrsws(mm(l)) = surf_def_h(l)%nrsws(m)
1679	IF ( ALLOCATED( surf_def_h(l)%sasws ) ) &
1680	surf_h(l)%sasws(mm(l)) = surf_def_h(l)%sasws(m)
1681
1682	mm(l) = mm(l) + 1
1683	ENDDO
1684
1685	IF ( l == 0 ) THEN
1686	DO m = surf_lsm_h%start_index(j,i), &
1687	surf_lsm_h%end_index(j,i)
1688	IF ( ALLOCATED( surf_lsm_h%us ) ) &
1689	surf_h(0)%us(mm(0)) = surf_lsm_h%us(m)
1690	IF ( ALLOCATED( surf_lsm_h%ts ) ) &
1691	surf_h(0)%ts(mm(0)) = surf_lsm_h%ts(m)
1692	IF ( ALLOCATED( surf_lsm_h%qs ) ) &
1693	surf_h(0)%qs(mm(0)) = surf_lsm_h%qs(m)
1694	IF ( ALLOCATED( surf_lsm_h%ss ) ) &
1695	surf_h(0)%ss(mm(0)) = surf_lsm_h%ss(m)
1696	IF ( ALLOCATED( surf_lsm_h%qrs ) ) &
1697	surf_h(0)%qrs(mm(0)) = surf_lsm_h%qrs(m)
1698	IF ( ALLOCATED( surf_lsm_h%nrs ) ) &
1699	surf_h(0)%nrs(mm(0)) = surf_lsm_h%nrs(m)
1700	IF ( ALLOCATED( surf_lsm_h%ol ) ) &
1701	surf_h(0)%ol(mm(0)) = surf_lsm_h%ol(m)
1702	IF ( ALLOCATED( surf_lsm_h%rib ) ) &
1703	surf_h(0)%rib(mm(0)) = surf_lsm_h%rib(m)
1704	IF ( ALLOCATED( surf_lsm_h%usws ) ) &
1705	surf_h(0)%usws(mm(0)) = surf_lsm_h%usws(m)
1706	IF ( ALLOCATED( surf_lsm_h%vsws ) ) &
1707	surf_h(0)%vsws(mm(0)) = surf_lsm_h%vsws(m)
1708	IF ( ALLOCATED( surf_lsm_h%shf ) ) &
1709	surf_h(0)%shf(mm(0)) = surf_lsm_h%shf(m)
1710	IF ( ALLOCATED( surf_lsm_h%qsws ) ) &
1711	surf_h(0)%qsws(mm(0)) = surf_lsm_h%qsws(m)
1712	IF ( ALLOCATED( surf_lsm_h%ssws ) ) &
1713	surf_h(0)%qsws(mm(0)) = surf_lsm_h%ssws(m)
1714	IF ( ALLOCATED( surf_lsm_h%nrsws ) ) &
1715	surf_h(0)%nrsws(mm(0)) = surf_lsm_h%nrsws(m)
1716	IF ( ALLOCATED( surf_lsm_h%sasws ) ) &
1717	surf_h(0)%sasws(mm(0)) = surf_lsm_h%sasws(m)
1718
1719	mm(0) = mm(0) + 1
1720
1721	ENDDO
1722
1723	DO m = surf_usm_h%start_index(j,i), &
1724	surf_usm_h%end_index(j,i)
1725	IF ( ALLOCATED( surf_usm_h%us ) ) &
1726	surf_h(0)%us(mm(0)) = surf_usm_h%us(m)
1727	IF ( ALLOCATED( surf_usm_h%ts ) ) &
1728	surf_h(0)%ts(mm(0)) = surf_usm_h%ts(m)
1729	IF ( ALLOCATED( surf_usm_h%qs ) ) &
1730	surf_h(0)%qs(mm(0)) = surf_usm_h%qs(m)
1731	IF ( ALLOCATED( surf_usm_h%ss ) ) &
1732	surf_h(0)%ss(mm(0)) = surf_usm_h%ss(m)
1733	IF ( ALLOCATED( surf_usm_h%qrs ) ) &
1734	surf_h(0)%qrs(mm(0)) = surf_usm_h%qrs(m)
1735	IF ( ALLOCATED( surf_usm_h%nrs ) ) &
1736	surf_h(0)%nrs(mm(0)) = surf_usm_h%nrs(m)
1737	IF ( ALLOCATED( surf_usm_h%ol ) ) &
1738	surf_h(0)%ol(mm(0)) = surf_usm_h%ol(m)
1739	IF ( ALLOCATED( surf_usm_h%rib ) ) &
1740	surf_h(0)%rib(mm(0)) = surf_usm_h%rib(m)
1741	IF ( ALLOCATED( surf_usm_h%usws ) ) &
1742	surf_h(0)%usws(mm(0)) = surf_usm_h%usws(m)
1743	IF ( ALLOCATED( surf_usm_h%vsws ) ) &
1744	surf_h(0)%vsws(mm(0)) = surf_usm_h%vsws(m)
1745	IF ( ALLOCATED( surf_usm_h%shf ) ) &
1746	surf_h(0)%shf(mm(0)) = surf_usm_h%shf(m)
1747	IF ( ALLOCATED( surf_usm_h%qsws ) ) &
1748	surf_h(0)%qsws(mm(0)) = surf_usm_h%qsws(m)
1749	IF ( ALLOCATED( surf_usm_h%ssws ) ) &
1750	surf_h(0)%qsws(mm(0)) = surf_usm_h%ssws(m)
1751	IF ( ALLOCATED( surf_usm_h%nrsws ) ) &
1752	surf_h(0)%nrsws(mm(0)) = surf_usm_h%nrsws(m)
1753	IF ( ALLOCATED( surf_usm_h%sasws ) ) &
1754	surf_h(0)%sasws(mm(0)) = surf_usm_h%sasws(m)
1755
1756	mm(0) = mm(0) + 1
1757
1758	ENDDO
1759
1760
1761	ENDIF
1762
1763	ENDDO
1764
1765	ENDDO
1766	IF ( l == 0 ) THEN
1767	surf_h(l)%start_index = MAX( surf_def_h(l)%start_index, &
1768	surf_lsm_h%start_index, &
1769	surf_usm_h%start_index )
1770	surf_h(l)%end_index = MAX( surf_def_h(l)%end_index, &
1771	surf_lsm_h%end_index, &
1772	surf_usm_h%end_index )
1773	ELSE
1774	surf_h(l)%start_index = surf_def_h(l)%start_index
1775	surf_h(l)%end_index = surf_def_h(l)%end_index
1776	ENDIF
1777	ENDDO
1778
1779
1780	mm(0:3) = 1
1781	DO l = 0, 3
1782	DO i = nxl, nxr
1783	DO j = nys, nyn
1784	DO m = surf_def_v(l)%start_index(j,i), &
1785	surf_def_v(l)%end_index(j,i)
1786	IF ( ALLOCATED( surf_def_v(l)%us ) ) &
1787	surf_v(l)%us(mm(l)) = surf_def_v(l)%us(m)
1788	IF ( ALLOCATED( surf_def_v(l)%ts ) ) &
1789	surf_v(l)%ts(mm(l)) = surf_def_v(l)%ts(m)
1790	IF ( ALLOCATED( surf_def_v(l)%qs ) ) &
1791	surf_v(l)%qs(mm(l)) = surf_def_v(l)%qs(m)
1792	IF ( ALLOCATED( surf_def_v(l)%ss ) ) &
1793	surf_v(l)%ss(mm(l)) = surf_def_v(l)%ss(m)
1794	IF ( ALLOCATED( surf_def_v(l)%qrs ) ) &
1795	surf_v(l)%qrs(mm(l)) = surf_def_v(l)%qrs(m)
1796	IF ( ALLOCATED( surf_def_v(l)%nrs ) ) &
1797	surf_v(l)%nrs(mm(l)) = surf_def_v(l)%nrs(m)
1798	IF ( ALLOCATED( surf_def_v(l)%ol ) ) &
1799	surf_v(l)%ol(mm(l)) = surf_def_v(l)%ol(m)
1800	IF ( ALLOCATED( surf_def_v(l)%rib ) ) &
1801	surf_v(l)%rib(mm(l)) = surf_def_v(l)%rib(m)
1802	IF ( ALLOCATED( surf_def_v(l)%shf ) ) &
1803	surf_v(l)%shf(mm(l)) = surf_def_v(l)%shf(m)
1804	IF ( ALLOCATED( surf_def_v(l)%qsws ) ) &
1805	surf_v(l)%qsws(mm(l)) = surf_def_v(l)%qsws(m)
1806	IF ( ALLOCATED( surf_def_v(l)%ssws ) ) &
1807	surf_v(l)%qsws(mm(l)) = surf_def_v(l)%ssws(m)
1808	IF ( ALLOCATED( surf_def_v(l)%nrsws ) ) &
1809	surf_v(l)%nrsws(mm(l)) = surf_def_v(l)%nrsws(m)
1810	IF ( ALLOCATED( surf_def_v(l)%sasws ) ) &
1811	surf_v(l)%sasws(mm(l)) = surf_def_v(l)%sasws(m)
1812	IF ( ALLOCATED( surf_def_v(l)%mom_flux_uv) ) &
1813	surf_v(l)%mom_flux_uv(mm(l)) = surf_def_v(l)%mom_flux_uv(m)
1814	IF ( ALLOCATED( surf_def_v(l)%mom_flux_w) ) &
1815	surf_v(l)%mom_flux_w(mm(l)) = surf_def_v(l)%mom_flux_w(m)
1816	IF ( ALLOCATED( surf_def_v(l)%mom_flux_tke) ) &
1817	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_def_v(l)%mom_flux_tke(0:1,m)
1818
1819	mm(l) = mm(l) + 1
1820	ENDDO
1821
1822	DO m = surf_lsm_v(l)%start_index(j,i), &
1823	surf_lsm_v(l)%end_index(j,i)
1824	IF ( ALLOCATED( surf_lsm_v(l)%us ) ) &
1825	surf_v(l)%us(mm(l)) = surf_lsm_v(l)%us(m)
1826	IF ( ALLOCATED( surf_lsm_v(l)%ts ) ) &
1827	surf_v(l)%ts(mm(l)) = surf_lsm_v(l)%ts(m)
1828	IF ( ALLOCATED( surf_lsm_v(l)%qs ) ) &
1829	surf_v(l)%qs(mm(l)) = surf_lsm_v(l)%qs(m)
1830	IF ( ALLOCATED( surf_lsm_v(l)%ss ) ) &
1831	surf_v(l)%ss(mm(l)) = surf_lsm_v(l)%ss(m)
1832	IF ( ALLOCATED( surf_lsm_v(l)%qrs ) ) &
1833	surf_v(l)%qrs(mm(l)) = surf_lsm_v(l)%qrs(m)
1834	IF ( ALLOCATED( surf_lsm_v(l)%nrs ) ) &
1835	surf_v(l)%nrs(mm(l)) = surf_lsm_v(l)%nrs(m)
1836	IF ( ALLOCATED( surf_lsm_v(l)%ol ) ) &
1837	surf_v(l)%ol(mm(l)) = surf_lsm_v(l)%ol(m)
1838	IF ( ALLOCATED( surf_lsm_v(l)%rib ) ) &
1839	surf_v(l)%rib(mm(l)) = surf_lsm_v(l)%rib(m)
1840	IF ( ALLOCATED( surf_lsm_v(l)%usws ) ) &
1841	surf_v(l)%usws(mm(l)) = surf_lsm_v(l)%usws(m)
1842	IF ( ALLOCATED( surf_lsm_v(l)%vsws ) ) &
1843	surf_v(l)%vsws(mm(l)) = surf_lsm_v(l)%vsws(m)
1844	IF ( ALLOCATED( surf_lsm_v(l)%shf ) ) &
1845	surf_v(l)%shf(mm(l)) = surf_lsm_v(l)%shf(m)
1846	IF ( ALLOCATED( surf_lsm_v(l)%qsws ) ) &
1847	surf_v(l)%qsws(mm(l)) = surf_lsm_v(l)%qsws(m)
1848	IF ( ALLOCATED( surf_lsm_v(l)%ssws ) ) &
1849	surf_v(l)%qsws(mm(l)) = surf_lsm_v(l)%ssws(m)
1850	IF ( ALLOCATED( surf_lsm_v(l)%nrsws ) ) &
1851	surf_v(l)%nrsws(mm(l)) = surf_lsm_v(l)%nrsws(m)
1852	IF ( ALLOCATED( surf_lsm_v(l)%sasws ) ) &
1853	surf_v(l)%sasws(mm(l)) = surf_lsm_v(l)%sasws(m)
1854	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_uv) ) &
1855	surf_v(l)%mom_flux_uv(mm(l)) = surf_lsm_v(l)%mom_flux_uv(m)
1856	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_w) ) &
1857	surf_v(l)%mom_flux_w(mm(l)) = surf_lsm_v(l)%mom_flux_w(m)
1858	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_tke) ) &
1859	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_lsm_v(l)%mom_flux_tke(0:1,m)
1860
1861	mm(l) = mm(l) + 1
1862	ENDDO
1863
1864	DO m = surf_usm_v(l)%start_index(j,i), &
1865	surf_usm_v(l)%end_index(j,i)
1866	IF ( ALLOCATED( surf_usm_v(l)%us ) ) &
1867	surf_v(l)%us(mm(l)) = surf_usm_v(l)%us(m)
1868	IF ( ALLOCATED( surf_usm_v(l)%ts ) ) &
1869	surf_v(l)%ts(mm(l)) = surf_usm_v(l)%ts(m)
1870	IF ( ALLOCATED( surf_usm_v(l)%qs ) ) &
1871	surf_v(l)%qs(mm(l)) = surf_usm_v(l)%qs(m)
1872	IF ( ALLOCATED( surf_usm_v(l)%ss ) ) &
1873	surf_v(l)%ss(mm(l)) = surf_usm_v(l)%ss(m)
1874	IF ( ALLOCATED( surf_usm_v(l)%qrs ) ) &
1875	surf_v(l)%qrs(mm(l)) = surf_usm_v(l)%qrs(m)
1876	IF ( ALLOCATED( surf_usm_v(l)%nrs ) ) &
1877	surf_v(l)%nrs(mm(l)) = surf_usm_v(l)%nrs(m)
1878	IF ( ALLOCATED( surf_usm_v(l)%ol ) ) &
1879	surf_v(l)%ol(mm(l)) = surf_usm_v(l)%ol(m)
1880	IF ( ALLOCATED( surf_usm_v(l)%rib ) ) &
1881	surf_v(l)%rib(mm(l)) = surf_usm_v(l)%rib(m)
1882	IF ( ALLOCATED( surf_usm_v(l)%usws ) ) &
1883	surf_v(l)%usws(mm(l)) = surf_usm_v(l)%usws(m)
1884	IF ( ALLOCATED( surf_usm_v(l)%vsws ) ) &
1885	surf_v(l)%vsws(mm(l)) = surf_usm_v(l)%vsws(m)
1886	IF ( ALLOCATED( surf_usm_v(l)%shf ) ) &
1887	surf_v(l)%shf(mm(l)) = surf_usm_v(l)%shf(m)
1888	IF ( ALLOCATED( surf_usm_v(l)%qsws ) ) &
1889	surf_v(l)%qsws(mm(l)) = surf_usm_v(l)%qsws(m)
1890	IF ( ALLOCATED( surf_usm_v(l)%ssws ) ) &
1891	surf_v(l)%qsws(mm(l)) = surf_usm_v(l)%ssws(m)
1892	IF ( ALLOCATED( surf_usm_v(l)%nrsws ) ) &
1893	surf_v(l)%nrsws(mm(l)) = surf_usm_v(l)%nrsws(m)
1894	IF ( ALLOCATED( surf_usm_v(l)%sasws ) ) &
1895	surf_v(l)%sasws(mm(l)) = surf_usm_v(l)%sasws(m)
1896	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_uv) ) &
1897	surf_v(l)%mom_flux_uv(mm(l)) = surf_usm_v(l)%mom_flux_uv(m)
1898	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_w) ) &
1899	surf_v(l)%mom_flux_w(mm(l)) = surf_usm_v(l)%mom_flux_w(m)
1900	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_tke) ) &
1901	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_usm_v(l)%mom_flux_tke(0:1,m)
1902
1903	mm(l) = mm(l) + 1
1904	ENDDO
1905
1906	ENDDO
1907	ENDDO
1908	!
1909	!-- Finally, determine start- and end-index for the respective surface
1910	surf_v(l)%start_index = MAX( surf_def_v(l)%start_index, &
1911	surf_lsm_v(l)%start_index, &
1912	surf_usm_v(l)%start_index )
1913	surf_v(l)%end_index = MAX( surf_def_v(l)%end_index, &
1914	surf_lsm_v(l)%end_index, &
1915	surf_usm_v(l)%end_index )
1916	ENDDO
1917
1918	WRITE ( 14 ) 'ns_h_on_file '
1919	WRITE ( 14 ) ns_h_on_file
1920	WRITE ( 14 ) 'ns_v_on_file '
1921	WRITE ( 14 ) ns_v_on_file
1922	!
1923	!-- Write required restart data.
1924	!-- Start with horizontal surfaces (upward-, downward-facing, and model top)
1925	DO l = 0, 2
1926	WRITE( dum, '(I1)') l
1927
1928	WRITE ( 14 ) 'surf_h(' // dum // ')%start_index '
1929	WRITE ( 14 ) surf_h(l)%start_index
1930	WRITE ( 14 ) 'surf_h(' // dum // ')%end_index '
1931	WRITE ( 14 ) surf_h(l)%end_index
1932
1933	WRITE ( 14 ) 'surf_h(' // dum // ')%us '
1934	IF ( ALLOCATED ( surf_h(l)%us ) ) THEN
1935	WRITE ( 14 ) surf_h(l)%us
1936	ENDIF
1937	WRITE ( 14 ) 'surf_h(' // dum // ')%ts '
1938	IF ( ALLOCATED ( surf_h(l)%ts ) ) THEN
1939	WRITE ( 14 ) surf_h(l)%ts
1940	ENDIF
1941	WRITE ( 14 ) 'surf_h(' // dum // ')%qs '
1942	IF ( ALLOCATED ( surf_h(l)%qs ) ) THEN
1943	WRITE ( 14 ) surf_h(l)%qs
1944	ENDIF
1945	WRITE ( 14 ) 'surf_h(' // dum // ')%ss '
1946	IF ( ALLOCATED ( surf_h(l)%ss ) ) THEN
1947	WRITE ( 14 ) surf_h(l)%ss
1948	ENDIF
1949	WRITE ( 14 ) 'surf_h(' // dum // ')%qrs '
1950	IF ( ALLOCATED ( surf_h(l)%qrs ) ) THEN
1951	WRITE ( 14 ) surf_h(l)%qrs
1952	ENDIF
1953	WRITE ( 14 ) 'surf_h(' // dum // ')%nrs '
1954	IF ( ALLOCATED ( surf_h(l)%nrs ) ) THEN
1955	WRITE ( 14 ) surf_h(l)%nrs
1956	ENDIF
1957	WRITE ( 14 ) 'surf_h(' // dum // ')%ol '
1958	IF ( ALLOCATED ( surf_h(l)%ol ) ) THEN
1959	WRITE ( 14 ) surf_h(l)%ol
1960	ENDIF
1961	WRITE ( 14 ) 'surf_h(' // dum // ')%rib '
1962	IF ( ALLOCATED ( surf_h(l)%rib ) ) THEN
1963	WRITE ( 14 ) surf_h(l)%rib
1964	ENDIF
1965	WRITE ( 14 ) 'surf_h(' // dum // ')%usws '
1966	IF ( ALLOCATED ( surf_h(l)%usws ) ) THEN
1967	WRITE ( 14 ) surf_h(l)%usws
1968	ENDIF
1969	WRITE ( 14 ) 'surf_h(' // dum // ')%vsws '
1970	IF ( ALLOCATED ( surf_h(l)%vsws ) ) THEN
1971	WRITE ( 14 ) surf_h(l)%vsws
1972	ENDIF
1973	WRITE ( 14 ) 'surf_h(' // dum // ')%shf '
1974	IF ( ALLOCATED ( surf_h(l)%shf ) ) THEN
1975	WRITE ( 14 ) surf_h(l)%shf
1976	ENDIF
1977	WRITE ( 14 ) 'surf_h(' // dum // ')%qsws '
1978	IF ( ALLOCATED ( surf_h(l)%qsws ) ) THEN
1979	WRITE ( 14 ) surf_h(l)%qsws
1980	ENDIF
1981	WRITE ( 14 ) 'surf_h(' // dum // ')%ssws '
1982	IF ( ALLOCATED ( surf_h(l)%ssws ) ) THEN
1983	WRITE ( 14 ) surf_h(l)%ssws
1984	ENDIF
1985	WRITE ( 14 ) 'surf_h(' // dum // ')%qrsws '
1986	IF ( ALLOCATED ( surf_h(l)%qrsws ) ) THEN
1987	WRITE ( 14 ) surf_h(l)%qrsws
1988	ENDIF
1989	WRITE ( 14 ) 'surf_h(' // dum // ')%nrsws '
1990	IF ( ALLOCATED ( surf_h(l)%nrsws ) ) THEN
1991	WRITE ( 14 ) surf_h(l)%nrsws
1992	ENDIF
1993	WRITE ( 14 ) 'surf_h(' // dum // ')%sasws '
1994	IF ( ALLOCATED ( surf_h(l)%sasws ) ) THEN
1995	WRITE ( 14 ) surf_h(l)%sasws
1996	ENDIF
1997	ENDDO
1998	!
1999	!-- Write vertical surfaces
2000	DO l = 0, 3
2001	WRITE( dum, '(I1)') l
2002
2003	WRITE ( 14 ) 'surf_v(' // dum // ')%start_index '
2004	WRITE ( 14 ) surf_v(l)%start_index
2005	WRITE ( 14 ) 'surf_v(' // dum // ')%end_index '
2006	WRITE ( 14 ) surf_v(l)%end_index
2007
2008	WRITE ( 14 ) 'surf_v(' // dum // ')%us '
2009	IF ( ALLOCATED ( surf_v(l)%us ) ) THEN
2010	WRITE ( 14 ) surf_v(l)%us
2011	ENDIF
2012	WRITE ( 14 ) 'surf_v(' // dum // ')%ts '
2013	IF ( ALLOCATED ( surf_v(l)%ts ) ) THEN
2014	WRITE ( 14 ) surf_v(l)%ts
2015	ENDIF
2016	WRITE ( 14 ) 'surf_v(' // dum // ')%qs '
2017	IF ( ALLOCATED ( surf_v(l)%qs ) ) THEN
2018	WRITE ( 14 ) surf_v(l)%qs
2019	ENDIF
2020	WRITE ( 14 ) 'surf_v(' // dum // ')%ss '
2021	IF ( ALLOCATED ( surf_v(l)%ss ) ) THEN
2022	WRITE ( 14 ) surf_v(l)%ss
2023	ENDIF
2024	WRITE ( 14 ) 'surf_v(' // dum // ')%qrs '
2025	IF ( ALLOCATED ( surf_v(l)%qrs ) ) THEN
2026	WRITE ( 14 ) surf_v(l)%qrs
2027	ENDIF
2028	WRITE ( 14 ) 'surf_v(' // dum // ')%nrs '
2029	IF ( ALLOCATED ( surf_v(l)%nrs ) ) THEN
2030	WRITE ( 14 ) surf_v(l)%nrs
2031	ENDIF
2032	WRITE ( 14 ) 'surf_v(' // dum // ')%ol '
2033	IF ( ALLOCATED ( surf_v(l)%ol ) ) THEN
2034	WRITE ( 14 ) surf_v(l)%ol
2035	ENDIF
2036	WRITE ( 14 ) 'surf_v(' // dum // ')%rib '
2037	IF ( ALLOCATED ( surf_v(l)%rib ) ) THEN
2038	WRITE ( 14 ) surf_v(l)%rib
2039	ENDIF
2040	WRITE ( 14 ) 'surf_v(' // dum // ')%shf '
2041	IF ( ALLOCATED ( surf_v(l)%shf ) ) THEN
2042	WRITE ( 14 ) surf_v(l)%shf
2043	ENDIF
2044	WRITE ( 14 ) 'surf_v(' // dum // ')%qsws '
2045	IF ( ALLOCATED ( surf_v(l)%qsws ) ) THEN
2046	WRITE ( 14 ) surf_v(l)%qsws
2047	ENDIF
2048	WRITE ( 14 ) 'surf_v(' // dum // ')%ssws '
2049	IF ( ALLOCATED ( surf_v(l)%ssws ) ) THEN
2050	WRITE ( 14 ) surf_v(l)%ssws
2051	ENDIF
2052	WRITE ( 14 ) 'surf_v(' // dum // ')%qrsws '
2053	IF ( ALLOCATED ( surf_v(l)%qrsws ) ) THEN
2054	WRITE ( 14 ) surf_v(l)%qrsws
2055	ENDIF
2056	WRITE ( 14 ) 'surf_v(' // dum // ')%nrsws '
2057	IF ( ALLOCATED ( surf_v(l)%nrsws ) ) THEN
2058	WRITE ( 14 ) surf_v(l)%nrsws
2059	ENDIF
2060	WRITE ( 14 ) 'surf_v(' // dum // ')%sasws '
2061	IF ( ALLOCATED ( surf_v(l)%sasws ) ) THEN
2062	WRITE ( 14 ) surf_v(l)%sasws
2063	ENDIF
2064	WRITE ( 14 ) 'surf_v(' // dum // ')%mom_uv '
2065	IF ( ALLOCATED ( surf_v(l)%mom_flux_uv ) ) THEN
2066	WRITE ( 14 ) surf_v(l)%mom_flux_uv
2067	ENDIF
2068	WRITE ( 14 ) 'surf_v(' // dum // ')%mom_w '
2069	IF ( ALLOCATED ( surf_v(l)%mom_flux_w ) ) THEN
2070	WRITE ( 14 ) surf_v(l)%mom_flux_w
2071	ENDIF
2072	WRITE ( 14 ) 'surf_v(' // dum // ')%mom_tke '
2073	IF ( ALLOCATED ( surf_v(l)%mom_flux_tke ) ) THEN
2074	WRITE ( 14 ) surf_v(l)%mom_flux_tke
2075	ENDIF
2076
2077	ENDDO
2078
2079	WRITE ( 14 ) '* end surf * '
2080
2081	END SUBROUTINE surface_write_restart_data
2082
2083
2084	!------------------------------------------------------------------------------!
2085	! Description:
2086	! ------------
2087	!> Reads surface-related restart data. Please note, restart data for a certain
2088	!> surface orientation (e.g. horizontal upward-facing) is stored in one
2089	!> array, even if surface elements may belong to different surface types
2090	!> natural or urban for example). Surface elements are redistributed into its
2091	!> respective surface types within this routine. This allows e.g. changing the
2092	!> surface type after reading the restart data, which might be required in case
2093	!> of cyclic_fill mode.
2094	!------------------------------------------------------------------------------!
2095	SUBROUTINE surface_read_restart_data( ii, &
2096	nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
2097	nynfa, nyn_on_file, nysfa, nys_on_file, &
2098	offset_xa, offset_ya, overlap_count )
2099
2100	USE pegrid, &
2101	ONLY: numprocs_previous_run
2102
2103	CHARACTER (LEN=1) :: dum !< dummy to create correct string for reading input variable
2104	CHARACTER (LEN=30) :: field_chr !< input variable
2105
2106	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
2107	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
2108	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
2109	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
2110	INTEGER(iwp) :: k !< running index along z-direction
2111	INTEGER(iwp) :: l !< index variable for surface type
2112	INTEGER(iwp) :: m !< running index for surface elements, refers to gathered array encompassing all surface types
2113	INTEGER(iwp) :: mm !< running index for surface elements, refers to individual surface types
2114
2115	INTEGER(iwp) :: ii !< running index over input files
2116	INTEGER(iwp) :: kk !< running index over previous input files covering current local domain
2117	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
2118	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
2119	INTEGER(iwp) :: nxl_on_file !< index of left boundary on former local domain
2120	INTEGER(iwp) :: nxrc !< index of right boundary on current subdomain
2121	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
2122	INTEGER(iwp) :: nxr_on_file !< index of right boundary on former local domain
2123	INTEGER(iwp) :: nync !< index of north boundary on current subdomain
2124	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
2125	INTEGER(iwp) :: nyn_on_file !< index of norht boundary on former local domain
2126	INTEGER(iwp) :: nysc !< index of south boundary on current subdomain
2127	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
2128	INTEGER(iwp) :: nys_on_file !< index of south boundary on former local domain
2129	INTEGER(iwp) :: overlap_count !< number of overlaps
2130
2131	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !<
2132	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !<
2133	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !<
2134	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !<
2135	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !<
2136	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !<
2137
2138
2139	LOGICAL :: horizontal_surface !< flag indicating horizontal surfaces
2140	LOGICAL :: surf_match_def !< flag indicating that surface element is of default type
2141	LOGICAL :: surf_match_lsm !< flag indicating that surface element is of natural type
2142	LOGICAL :: surf_match_usm !< flag indicating that surface element is of urban type
2143	LOGICAL :: vertical_surface !< flag indicating vertical surfaces
2144
2145	TYPE(surf_type), DIMENSION(0:2) :: surf_h !< horizontal surface type on file
2146	TYPE(surf_type), DIMENSION(0:3) :: surf_v !< vertical surface type on file
2147
2148	!
2149	!-- Read number of respective surface elements on file
2150	READ ( 13 ) field_chr
2151	IF ( TRIM( field_chr ) /= 'ns_h_on_file' ) THEN
2152	!
2153	!-- Add a proper error message
2154	ENDIF
2155	READ ( 13 ) ns_h_on_file
2156
2157	READ ( 13 ) field_chr
2158	IF ( TRIM( field_chr ) /= 'ns_v_on_file' ) THEN
2159	!
2160	!-- Add a proper error message
2161	ENDIF
2162	READ ( 13 ) ns_v_on_file
2163	!
2164	!-- Allocate memory for number of surface elements on file. Please note,
2165	!-- these number is not necessarily the same as the final number of surface
2166	!-- elements on local domain, which is the case if processor topology changes
2167	!-- during restart runs.
2168	!-- Horizontal upward facing
2169	surf_h(0)%ns = ns_h_on_file(0)
2170	CALL allocate_surface_attributes_h( surf_h(0), &
2171	nys_on_file, nyn_on_file, &
2172	nxl_on_file, nxr_on_file )
2173	!
2174	!-- Horizontal downward facing
2175	surf_h(1)%ns = ns_h_on_file(1)
2176	CALL allocate_surface_attributes_h( surf_h(1), &
2177	nys_on_file, nyn_on_file, &
2178	nxl_on_file, nxr_on_file )
2179	!
2180	!-- Model top
2181	surf_h(2)%ns = ns_h_on_file(2)
2182	CALL allocate_surface_attributes_h_top( surf_h(2), &
2183	nys_on_file, nyn_on_file, &
2184	nxl_on_file, nxr_on_file )
2185	!
2186	!-- Vertical surfaces
2187	DO l = 0, 3
2188	surf_v(l)%ns = ns_v_on_file(l)
2189	CALL allocate_surface_attributes_v( surf_v(l), .FALSE., &
2190	nys_on_file, nyn_on_file, &
2191	nxl_on_file, nxr_on_file )
2192	ENDDO
2193
2194	IF ( initializing_actions == 'read_restart_data' .OR. &
2195	initializing_actions == 'cyclic_fill' ) THEN
2196	!
2197	!-- Initial setting of flags for horizontal and vertical surfaces, will
2198	!-- be set after start- and end-indices are read.
2199	horizontal_surface = .FALSE.
2200	vertical_surface = .FALSE.
2201
2202	READ ( 13 ) field_chr
2203
2204	DO WHILE ( TRIM( field_chr ) /= '* end surf *' )
2205	!
2206	!-- Map data on file as often as needed (data are read only for k=1)
2207	DO kk = 1, overlap_count
2208	!
2209	!-- Get the index range of the subdomain on file which overlap with the
2210	!-- current subdomain
2211	nxlf = nxlfa(ii,kk)
2212	nxlc = nxlfa(ii,kk) + offset_xa(ii,kk)
2213	nxrf = nxrfa(ii,kk)
2214	nxrc = nxrfa(ii,kk) + offset_xa(ii,kk)
2215	nysf = nysfa(ii,kk)
2216	nysc = nysfa(ii,kk) + offset_ya(ii,kk)
2217	nynf = nynfa(ii,kk)
2218	nync = nynfa(ii,kk) + offset_ya(ii,kk)
2219
2220	SELECT CASE ( TRIM( field_chr ) )
2221
2222	CASE ( 'surf_h(0)%start_index' )
2223	IF ( kk == 1 ) &
2224	READ ( 13 ) surf_h(0)%start_index
2225	l = 0
2226	CASE ( 'surf_h(0)%end_index' )
2227	IF ( kk == 1 ) &
2228	READ ( 13 ) surf_h(0)%end_index
2229	horizontal_surface = .TRUE.
2230	vertical_surface = .FALSE.
2231	CASE ( 'surf_h(0)%us' )
2232	IF ( ALLOCATED( surf_h(0)%us ) .AND. kk == 1 ) &
2233	READ ( 13 ) surf_h(0)%us
2234	CASE ( 'surf_h(0)%ts' )
2235	IF ( ALLOCATED( surf_h(0)%ts ) .AND. kk == 1 ) &
2236	READ ( 13 ) surf_h(0)%ts
2237	CASE ( 'surf_h(0)%qs' )
2238	IF ( ALLOCATED( surf_h(0)%qs ) .AND. kk == 1 ) &
2239	READ ( 13 ) surf_h(0)%qs
2240	CASE ( 'surf_h(0)%ss' )
2241	IF ( ALLOCATED( surf_h(0)%ss ) .AND. kk == 1 ) &
2242	READ ( 13 ) surf_h(0)%ss
2243	CASE ( 'surf_h(0)%qrs' )
2244	IF ( ALLOCATED( surf_h(0)%qrs ) .AND. kk == 1 ) &
2245	READ ( 13 ) surf_h(0)%qrs
2246	CASE ( 'surf_h(0)%nrs' )
2247	IF ( ALLOCATED( surf_h(0)%nrs ) .AND. kk == 1 ) &
2248	READ ( 13 ) surf_h(0)%nrs
2249	CASE ( 'surf_h(0)%ol' )
2250	IF ( ALLOCATED( surf_h(0)%ol ) .AND. kk == 1 ) &
2251	READ ( 13 ) surf_h(0)%ol
2252	CASE ( 'surf_h(0)%rib' )
2253	IF ( ALLOCATED( surf_h(0)%rib ) .AND. kk == 1 ) &
2254	READ ( 13 ) surf_h(0)%rib
2255	CASE ( 'surf_h(0)%usws' )
2256	IF ( ALLOCATED( surf_h(0)%usws ) .AND. kk == 1 ) &
2257	READ ( 13 ) surf_h(0)%usws
2258	CASE ( 'surf_h(0)%vsws' )
2259	IF ( ALLOCATED( surf_h(0)%vsws ) .AND. kk == 1 ) &
2260	READ ( 13 ) surf_h(0)%vsws
2261	CASE ( 'surf_h(0)%shf' )
2262	IF ( ALLOCATED( surf_h(0)%shf ) .AND. kk == 1 ) &
2263	READ ( 13 ) surf_h(0)%shf
2264	CASE ( 'surf_h(0)%qsws' )
2265	IF ( ALLOCATED( surf_h(0)%qsws ) .AND. kk == 1 ) &
2266	READ ( 13 ) surf_h(0)%qsws
2267	CASE ( 'surf_h(0)%ssws' )
2268	IF ( ALLOCATED( surf_h(0)%ssws ) .AND. kk == 1 ) &
2269	READ ( 13 ) surf_h(0)%ssws
2270	CASE ( 'surf_h(0)%qrsws' )
2271	IF ( ALLOCATED( surf_h(0)%qrsws ) .AND. kk == 1 ) &
2272	READ ( 13 ) surf_h(0)%qrsws
2273	CASE ( 'surf_h(0)%nrsws' )
2274	IF ( ALLOCATED( surf_h(0)%nrsws ) .AND. kk == 1 ) &
2275	READ ( 13 ) surf_h(0)%nrsws
2276	CASE ( 'surf_h(0)%sasws' )
2277	IF ( ALLOCATED( surf_h(0)%sasws ) .AND. kk == 1 ) &
2278	READ ( 13 ) surf_h(0)%sasws
2279
2280	CASE ( 'surf_h(1)%start_index' )
2281	IF ( kk == 1 ) &
2282	READ ( 13 ) surf_h(1)%start_index
2283	l = 1
2284	CASE ( 'surf_h(1)%end_index' )
2285	IF ( kk == 1 ) &
2286	READ ( 13 ) surf_h(1)%end_index
2287	CASE ( 'surf_h(1)%us' )
2288	IF ( ALLOCATED( surf_h(1)%us ) .AND. kk == 1 ) &
2289	READ ( 13 ) surf_h(1)%us
2290	CASE ( 'surf_h(1)%ts' )
2291	IF ( ALLOCATED( surf_h(1)%ts ) .AND. kk == 1 ) &
2292	READ ( 13 ) surf_h(1)%ts
2293	CASE ( 'surf_h(1)%qs' )
2294	IF ( ALLOCATED( surf_h(1)%qs ) .AND. kk == 1 ) &
2295	READ ( 13 ) surf_h(1)%qs
2296	CASE ( 'surf_h(1)%ss' )
2297	IF ( ALLOCATED( surf_h(1)%ss ) .AND. kk == 1 ) &
2298	READ ( 13 ) surf_h(1)%ss
2299	CASE ( 'surf_h(1)%qrs' )
2300	IF ( ALLOCATED( surf_h(1)%qrs ) .AND. kk == 1 ) &
2301	READ ( 13 ) surf_h(1)%qrs
2302	CASE ( 'surf_h(1)%nrs' )
2303	IF ( ALLOCATED( surf_h(1)%nrs ) .AND. kk == 1 ) &
2304	READ ( 13 ) surf_h(1)%nrs
2305	CASE ( 'surf_h(1)%ol' )
2306	IF ( ALLOCATED( surf_h(1)%ol ) .AND. kk == 1 ) &
2307	READ ( 13 ) surf_h(1)%ol
2308	CASE ( 'surf_h(1)%rib' )
2309	IF ( ALLOCATED( surf_h(1)%rib ) .AND. kk == 1 ) &
2310	READ ( 13 ) surf_h(1)%rib
2311	CASE ( 'surf_h(1)%usws' )
2312	IF ( ALLOCATED( surf_h(1)%usws ) .AND. kk == 1 ) &
2313	READ ( 13 ) surf_h(1)%usws
2314	CASE ( 'surf_h(1)%vsws' )
2315	IF ( ALLOCATED( surf_h(1)%vsws ) .AND. kk == 1 ) &
2316	READ ( 13 ) surf_h(1)%vsws
2317	CASE ( 'surf_h(1)%shf' )
2318	IF ( ALLOCATED( surf_h(1)%shf ) .AND. kk == 1 ) &
2319	READ ( 13 ) surf_h(1)%shf
2320	CASE ( 'surf_h(1)%qsws' )
2321	IF ( ALLOCATED( surf_h(1)%qsws ) .AND. kk == 1 ) &
2322	READ ( 13 ) surf_h(1)%qsws
2323	CASE ( 'surf_h(1)%ssws' )
2324	IF ( ALLOCATED( surf_h(1)%ssws ) .AND. kk == 1 ) &
2325	READ ( 13 ) surf_h(1)%ssws
2326	CASE ( 'surf_h(1)%qrsws' )
2327	IF ( ALLOCATED( surf_h(1)%qrsws ) .AND. kk == 1 ) &
2328	READ ( 13 ) surf_h(1)%qrsws
2329	CASE ( 'surf_h(1)%nrsws' )
2330	IF ( ALLOCATED( surf_h(1)%nrsws ) .AND. kk == 1 ) &
2331	READ ( 13 ) surf_h(1)%nrsws
2332	CASE ( 'surf_h(1)%sasws' )
2333	IF ( ALLOCATED( surf_h(1)%sasws ) .AND. kk == 1 ) &
2334	READ ( 13 ) surf_h(1)%sasws
2335
2336	CASE ( 'surf_h(2)%start_index' )
2337	IF ( kk == 1 ) &
2338	READ ( 13 ) surf_h(2)%start_index
2339	l = 2
2340	CASE ( 'surf_h(2)%end_index' )
2341	IF ( kk == 1 ) &
2342	READ ( 13 ) surf_h(2)%end_index
2343	CASE ( 'surf_h(2)%us' )
2344	IF ( ALLOCATED( surf_h(2)%us ) .AND. kk == 1 ) &
2345	READ ( 13 ) surf_h(2)%us
2346	CASE ( 'surf_h(2)%ts' )
2347	IF ( ALLOCATED( surf_h(2)%ts ) .AND. kk == 1 ) &
2348	READ ( 13 ) surf_h(2)%ts
2349	CASE ( 'surf_h(2)%qs' )
2350	IF ( ALLOCATED( surf_h(2)%qs ) .AND. kk == 1 ) &
2351	READ ( 13 ) surf_h(2)%qs
2352	CASE ( 'surf_h(2)%ss' )
2353	IF ( ALLOCATED( surf_h(2)%ss ) .AND. kk == 1 ) &
2354	READ ( 13 ) surf_h(2)%ss
2355	CASE ( 'surf_h(2)%qrs' )
2356	IF ( ALLOCATED( surf_h(2)%qrs ) .AND. kk == 1 ) &
2357	READ ( 13 ) surf_h(2)%qrs
2358	CASE ( 'surf_h(2)%nrs' )
2359	IF ( ALLOCATED( surf_h(2)%nrs ) .AND. kk == 1 ) &
2360	READ ( 13 ) surf_h(2)%nrs
2361	CASE ( 'surf_h(2)%ol' )
2362	IF ( ALLOCATED( surf_h(2)%ol ) .AND. kk == 1 ) &
2363	READ ( 13 ) surf_h(2)%ol
2364	CASE ( 'surf_h(2)%rib' )
2365	IF ( ALLOCATED( surf_h(2)%rib ) .AND. kk == 1 ) &
2366	READ ( 13 ) surf_h(2)%rib
2367	CASE ( 'surf_h(2)%usws' )
2368	IF ( ALLOCATED( surf_h(2)%usws ) .AND. kk == 1 ) &
2369	READ ( 13 ) surf_h(2)%usws
2370	CASE ( 'surf_h(2)%vsws' )
2371	IF ( ALLOCATED( surf_h(2)%vsws ) .AND. kk == 1 ) &
2372	READ ( 13 ) surf_h(2)%vsws
2373	CASE ( 'surf_h(2)%shf' )
2374	IF ( ALLOCATED( surf_h(2)%shf ) .AND. kk == 1 ) &
2375	READ ( 13 ) surf_h(2)%shf
2376	CASE ( 'surf_h(2)%qsws' )
2377	IF ( ALLOCATED( surf_h(2)%qsws ) .AND. kk == 1 ) &
2378	READ ( 13 ) surf_h(2)%qsws
2379	CASE ( 'surf_h(2)%ssws' )
2380	IF ( ALLOCATED( surf_h(2)%ssws ) .AND. kk == 1 ) &
2381	READ ( 13 ) surf_h(2)%ssws
2382	CASE ( 'surf_h(2)%qrsws' )
2383	IF ( ALLOCATED( surf_h(2)%qrsws ) .AND. kk == 1 ) &
2384	READ ( 13 ) surf_h(2)%qrsws
2385	CASE ( 'surf_h(2)%nrsws' )
2386	IF ( ALLOCATED( surf_h(2)%nrsws ) .AND. kk == 1 ) &
2387	READ ( 13 ) surf_h(2)%nrsws
2388	CASE ( 'surf_h(2)%sasws' )
2389	IF ( ALLOCATED( surf_h(2)%sasws ) .AND. kk == 1 ) &
2390	READ ( 13 ) surf_h(2)%sasws
2391
2392	CASE ( 'surf_v(0)%start_index' )
2393	IF ( kk == 1 ) &
2394	READ ( 13 ) surf_v(0)%start_index
2395	l = 0
2396	horizontal_surface = .FALSE.
2397	vertical_surface = .TRUE.
2398	CASE ( 'surf_v(0)%end_index' )
2399	IF ( kk == 1 ) &
2400	READ ( 13 ) surf_v(0)%end_index
2401	CASE ( 'surf_v(0)%us' )
2402	IF ( ALLOCATED( surf_v(0)%us ) .AND. kk == 1 ) &
2403	READ ( 13 ) surf_v(0)%us
2404	CASE ( 'surf_v(0)%ts' )
2405	IF ( ALLOCATED( surf_v(0)%ts ) .AND. kk == 1 ) &
2406	READ ( 13 ) surf_v(0)%ts
2407	CASE ( 'surf_v(0)%qs' )
2408	IF ( ALLOCATED( surf_v(0)%qs ) .AND. kk == 1 ) &
2409	READ ( 13 ) surf_v(0)%qs
2410	CASE ( 'surf_v(0)%ss' )
2411	IF ( ALLOCATED( surf_v(0)%ss ) .AND. kk == 1 ) &
2412	READ ( 13 ) surf_v(0)%ss
2413	CASE ( 'surf_v(0)%qrs' )
2414	IF ( ALLOCATED( surf_v(0)%qrs ) .AND. kk == 1 ) &
2415	READ ( 13 ) surf_v(0)%qrs
2416	CASE ( 'surf_v(0)%nrs' )
2417	IF ( ALLOCATED( surf_v(0)%nrs ) .AND. kk == 1 ) &
2418	READ ( 13 ) surf_v(0)%nrs
2419	CASE ( 'surf_v(0)%ol' )
2420	IF ( ALLOCATED( surf_v(0)%ol ) .AND. kk == 1 ) &
2421	READ ( 13 ) surf_v(0)%ol
2422	CASE ( 'surf_v(0)%rib' )
2423	IF ( ALLOCATED( surf_v(0)%rib ) .AND. kk == 1 ) &
2424	READ ( 13 ) surf_v(0)%rib
2425	CASE ( 'surf_v(0)%shf' )
2426	IF ( ALLOCATED( surf_v(0)%shf ) .AND. kk == 1 ) &
2427	READ ( 13 ) surf_v(0)%shf
2428	CASE ( 'surf_v(0)%qsws' )
2429	IF ( ALLOCATED( surf_v(0)%qsws ) .AND. kk == 1 ) &
2430	READ ( 13 ) surf_v(0)%qsws
2431	CASE ( 'surf_v(0)%ssws' )
2432	IF ( ALLOCATED( surf_v(0)%ssws ) .AND. kk == 1 ) &
2433	READ ( 13 ) surf_v(0)%ssws
2434	CASE ( 'surf_v(0)%qrsws' )
2435	IF ( ALLOCATED( surf_v(0)%qrsws ) .AND. kk == 1 ) &
2436	READ ( 13 ) surf_v(0)%qrsws
2437	CASE ( 'surf_v(0)%nrsws' )
2438	IF ( ALLOCATED( surf_v(0)%nrsws ) .AND. kk == 1 ) &
2439	READ ( 13 ) surf_v(0)%nrsws
2440	CASE ( 'surf_v(0)%sasws' )
2441	IF ( ALLOCATED( surf_v(0)%sasws ) .AND. kk == 1 ) &
2442	READ ( 13 ) surf_v(0)%sasws
2443	CASE ( 'surf_v(0)%mom_uv' )
2444	IF ( ALLOCATED( surf_v(0)%mom_flux_uv ) .AND. kk == 1 )&
2445	READ ( 13 ) surf_v(0)%mom_flux_uv
2446	CASE ( 'surf_v(0)%mom_w' )
2447	IF ( ALLOCATED( surf_v(0)%mom_flux_w ) .AND. kk == 1 ) &
2448	READ ( 13 ) surf_v(0)%mom_flux_w
2449	CASE ( 'surf_v(0)%mom_tke' )
2450	IF ( ALLOCATED( surf_v(0)%mom_flux_tke ) .AND. kk == 1 )&
2451	READ ( 13 ) surf_v(0)%mom_flux_tke
2452
2453	CASE ( 'surf_v(1)%start_index' )
2454	IF ( kk == 1 ) &
2455	READ ( 13 ) surf_v(1)%start_index
2456	l = 1
2457	CASE ( 'surf_v(1)%end_index' )
2458	IF ( kk == 1 ) &
2459	READ ( 13 ) surf_v(1)%end_index
2460	CASE ( 'surf_v(1)%us' )
2461	IF ( ALLOCATED( surf_v(1)%us ) .AND. kk == 1 ) &
2462	READ ( 13 ) surf_v(1)%us
2463	CASE ( 'surf_v(1)%ts' )
2464	IF ( ALLOCATED( surf_v(1)%ts ) .AND. kk == 1 ) &
2465	READ ( 13 ) surf_v(1)%ts
2466	CASE ( 'surf_v(1)%qs' )
2467	IF ( ALLOCATED( surf_v(1)%qs ) .AND. kk == 1 ) &
2468	READ ( 13 ) surf_v(1)%qs
2469	CASE ( 'surf_v(1)%ss' )
2470	IF ( ALLOCATED( surf_v(1)%ss ) .AND. kk == 1 ) &
2471	READ ( 13 ) surf_v(1)%ss
2472	CASE ( 'surf_v(1)%qrs' )
2473	IF ( ALLOCATED( surf_v(1)%qrs ) .AND. kk == 1 ) &
2474	READ ( 13 ) surf_v(1)%qrs
2475	CASE ( 'surf_v(1)%nrs' )
2476	IF ( ALLOCATED( surf_v(1)%nrs ) .AND. kk == 1 ) &
2477	READ ( 13 ) surf_v(1)%nrs
2478	CASE ( 'surf_v(1)%ol' )
2479	IF ( ALLOCATED( surf_v(1)%ol ) .AND. kk == 1 ) &
2480	READ ( 13 ) surf_v(1)%ol
2481	CASE ( 'surf_v(1)%rib' )
2482	IF ( ALLOCATED( surf_v(1)%rib ) .AND. kk == 1 ) &
2483	READ ( 13 ) surf_v(1)%rib
2484	CASE ( 'surf_v(1)%shf' )
2485	IF ( ALLOCATED( surf_v(1)%shf ) .AND. kk == 1 ) &
2486	READ ( 13 ) surf_v(1)%shf
2487	CASE ( 'surf_v(1)%qsws' )
2488	IF ( ALLOCATED( surf_v(1)%qsws ) .AND. kk == 1 ) &
2489	READ ( 13 ) surf_v(1)%qsws
2490	CASE ( 'surf_v(1)%ssws' )
2491	IF ( ALLOCATED( surf_v(1)%ssws ) .AND. kk == 1 ) &
2492	READ ( 13 ) surf_v(1)%ssws
2493	CASE ( 'surf_v(1)%qrsws' )
2494	IF ( ALLOCATED( surf_v(1)%qrsws ) .AND. kk == 1 ) &
2495	READ ( 13 ) surf_v(1)%qrsws
2496	CASE ( 'surf_v(1)%nrsws' )
2497	IF ( ALLOCATED( surf_v(1)%nrsws ) .AND. kk == 1 ) &
2498	READ ( 13 ) surf_v(1)%nrsws
2499	CASE ( 'surf_v(1)%sasws' )
2500	IF ( ALLOCATED( surf_v(1)%sasws ) .AND. kk == 1 ) &
2501	READ ( 13 ) surf_v(1)%sasws
2502	CASE ( 'surf_v(1)%mom_uv' )
2503	IF ( ALLOCATED( surf_v(1)%mom_flux_uv ) .AND. kk == 1 )&
2504	READ ( 13 ) surf_v(1)%mom_flux_uv
2505	CASE ( 'surf_v(1)%mom_w' )
2506	IF ( ALLOCATED( surf_v(1)%mom_flux_w ) .AND. kk == 1 ) &
2507	READ ( 13 ) surf_v(1)%mom_flux_w
2508	CASE ( 'surf_v(1)%mom_tke' )
2509	IF ( ALLOCATED( surf_v(1)%mom_flux_tke ) .AND. kk == 1 )&
2510	READ ( 13 ) surf_v(1)%mom_flux_tke
2511
2512	CASE ( 'surf_v(2)%start_index' )
2513	IF ( kk == 1 ) &
2514	READ ( 13 ) surf_v(2)%start_index
2515	l = 2
2516	CASE ( 'surf_v(2)%end_index' )
2517	IF ( kk == 1 ) &
2518	READ ( 13 ) surf_v(2)%end_index
2519	CASE ( 'surf_v(2)%us' )
2520	IF ( ALLOCATED( surf_v(2)%us ) .AND. kk == 1 ) &
2521	READ ( 13 ) surf_v(2)%us
2522	CASE ( 'surf_v(2)%ts' )
2523	IF ( ALLOCATED( surf_v(2)%ts ) .AND. kk == 1 ) &
2524	READ ( 13 ) surf_v(2)%ts
2525	CASE ( 'surf_v(2)%qs' )
2526	IF ( ALLOCATED( surf_v(2)%qs ) .AND. kk == 1 ) &
2527	READ ( 13 ) surf_v(2)%qs
2528	CASE ( 'surf_v(2)%ss' )
2529	IF ( ALLOCATED( surf_v(2)%ss ) .AND. kk == 1 ) &
2530	READ ( 13 ) surf_v(2)%ss
2531	CASE ( 'surf_v(2)%qrs' )
2532	IF ( ALLOCATED( surf_v(2)%qrs ) .AND. kk == 1 ) &
2533	READ ( 13 ) surf_v(2)%qrs
2534	CASE ( 'surf_v(2)%nrs' )
2535	IF ( ALLOCATED( surf_v(2)%nrs ) .AND. kk == 1 ) &
2536	READ ( 13 ) surf_v(2)%nrs
2537	CASE ( 'surf_v(2)%ol' )
2538	IF ( ALLOCATED( surf_v(2)%ol ) .AND. kk == 1 ) &
2539	READ ( 13 ) surf_v(2)%ol
2540	CASE ( 'surf_v(2)%rib' )
2541	IF ( ALLOCATED( surf_v(2)%rib ) .AND. kk == 1 ) &
2542	READ ( 13 ) surf_v(2)%rib
2543	CASE ( 'surf_v(2)%shf' )
2544	IF ( ALLOCATED( surf_v(2)%shf ) .AND. kk == 1 ) &
2545	READ ( 13 ) surf_v(2)%shf
2546	CASE ( 'surf_v(2)%qsws' )
2547	IF ( ALLOCATED( surf_v(2)%qsws ) .AND. kk == 1 ) &
2548	READ ( 13 ) surf_v(2)%qsws
2549	CASE ( 'surf_v(2)%ssws' )
2550	IF ( ALLOCATED( surf_v(2)%ssws ) .AND. kk == 1 ) &
2551	READ ( 13 ) surf_v(2)%ssws
2552	CASE ( 'surf_v(2)%qrsws' )
2553	IF ( ALLOCATED( surf_v(2)%qrsws ) .AND. kk == 1 ) &
2554	READ ( 13 ) surf_v(2)%qrsws
2555	CASE ( 'surf_v(2)%nrsws' )
2556	IF ( ALLOCATED( surf_v(2)%nrsws ) .AND. kk == 1 ) &
2557	READ ( 13 ) surf_v(2)%nrsws
2558	CASE ( 'surf_v(2)%sasws' )
2559	IF ( ALLOCATED( surf_v(2)%sasws ) .AND. kk == 1 ) &
2560	READ ( 13 ) surf_v(2)%sasws
2561	CASE ( 'surf_v(2)%mom_uv' )
2562	IF ( ALLOCATED( surf_v(2)%mom_flux_uv ) .AND. kk == 1 )&
2563	READ ( 13 ) surf_v(2)%mom_flux_uv
2564	CASE ( 'surf_v(2)%mom_w' )
2565	IF ( ALLOCATED( surf_v(2)%mom_flux_w ) .AND. kk == 1 ) &
2566	READ ( 13 ) surf_v(2)%mom_flux_w
2567	CASE ( 'surf_v(2)%mom_tke' )
2568	IF ( ALLOCATED( surf_v(2)%mom_flux_tke ) .AND. kk == 1 )&
2569	READ ( 13 ) surf_v(2)%mom_flux_tke
2570
2571	CASE ( 'surf_v(3)%start_index' )
2572	IF ( kk == 1 ) &
2573	READ ( 13 ) surf_v(3)%start_index
2574	l = 3
2575	CASE ( 'surf_v(3)%end_index' )
2576	IF ( kk == 1 ) &
2577	READ ( 13 ) surf_v(3)%end_index
2578	CASE ( 'surf_v(3)%us' )
2579	IF ( ALLOCATED( surf_v(3)%us ) .AND. kk == 1 ) &
2580	READ ( 13 ) surf_v(3)%us
2581	CASE ( 'surf_v(3)%ts' )
2582	IF ( ALLOCATED( surf_v(3)%ts ) .AND. kk == 1 ) &
2583	READ ( 13 ) surf_v(3)%ts
2584	CASE ( 'surf_v(3)%qs' )
2585	IF ( ALLOCATED( surf_v(3)%qs ) .AND. kk == 1 ) &
2586	READ ( 13 ) surf_v(3)%qs
2587	CASE ( 'surf_v(3)%ss' )
2588	IF ( ALLOCATED( surf_v(3)%ss ) .AND. kk == 1 ) &
2589	READ ( 13 ) surf_v(3)%ss
2590	CASE ( 'surf_v(3)%qrs' )
2591	IF ( ALLOCATED( surf_v(3)%qrs ) .AND. kk == 1 ) &
2592	READ ( 13 ) surf_v(3)%qrs
2593	CASE ( 'surf_v(3)%nrs' )
2594	IF ( ALLOCATED( surf_v(3)%nrs ) .AND. kk == 1 ) &
2595	READ ( 13 ) surf_v(3)%nrs
2596	CASE ( 'surf_v(3)%ol' )
2597	IF ( ALLOCATED( surf_v(3)%ol ) .AND. kk == 1 ) &
2598	READ ( 13 ) surf_v(3)%ol
2599	CASE ( 'surf_v(3)%rib' )
2600	IF ( ALLOCATED( surf_v(3)%rib ) .AND. kk == 1 ) &
2601	READ ( 13 ) surf_v(3)%rib
2602	CASE ( 'surf_v(3)%shf' )
2603	IF ( ALLOCATED( surf_v(3)%shf ) .AND. kk == 1 ) &
2604	READ ( 13 ) surf_v(3)%shf
2605	CASE ( 'surf_v(3)%qsws' )
2606	IF ( ALLOCATED( surf_v(3)%qsws ) .AND. kk == 1 ) &
2607	READ ( 13 ) surf_v(3)%qsws
2608	CASE ( 'surf_v(3)%ssws' )
2609	IF ( ALLOCATED( surf_v(3)%ssws ) .AND. kk == 1 ) &
2610	READ ( 13 ) surf_v(3)%ssws
2611	CASE ( 'surf_v(3)%qrsws' )
2612	IF ( ALLOCATED( surf_v(3)%qrsws ) .AND. kk == 1 ) &
2613	READ ( 13 ) surf_v(3)%qrsws
2614	CASE ( 'surf_v(3)%nrsws' )
2615	IF ( ALLOCATED( surf_v(3)%nrsws ) .AND. kk == 1 ) &
2616	READ ( 13 ) surf_v(3)%nrsws
2617	CASE ( 'surf_v(3)%sasws' )
2618	IF ( ALLOCATED( surf_v(3)%sasws ) .AND. kk == 1 ) &
2619	READ ( 13 ) surf_v(3)%sasws
2620	CASE ( 'surf_v(3)%mom_uv' )
2621	IF ( ALLOCATED( surf_v(3)%mom_flux_uv ) .AND. kk == 1 )&
2622	READ ( 13 ) surf_v(3)%mom_flux_uv
2623	CASE ( 'surf_v(3)%mom_w' )
2624	IF ( ALLOCATED( surf_v(3)%mom_flux_w ) .AND. kk == 1 ) &
2625	READ ( 13 ) surf_v(3)%mom_flux_w
2626	CASE ( 'surf_v(3)%mom_tke' )
2627	IF ( ALLOCATED( surf_v(3)%mom_flux_tke ) .AND. kk == 1 )&
2628	READ ( 13 ) surf_v(3)%mom_flux_tke
2629
2630	END SELECT
2631	!
2632	!-- Redistribute surface elements on its respective type.
2633	IF ( horizontal_surface ) THEN
2634	ic = nxlc
2635	DO i = nxlf, nxrf
2636	jc = nysc
2637	DO j = nysf, nynf
2638
2639	surf_match_def = surf_def_h(l)%end_index(jc,ic) >= &
2640	surf_def_h(l)%start_index(jc,ic)
2641	surf_match_lsm = surf_lsm_h%end_index(jc,ic) >= &
2642	surf_lsm_h%start_index(jc,ic)
2643	surf_match_usm = surf_usm_h%end_index(jc,ic) >= &
2644	surf_usm_h%start_index(jc,ic)
2645
2646	IF ( surf_match_def ) THEN
2647	mm = surf_def_h(l)%start_index(jc,ic)
2648	DO m = surf_h(l)%start_index(j,i), &
2649	surf_h(l)%end_index(j,i)
2650	CALL restore_surface_elements( surf_def_h(l), &
2651	mm, surf_h(l), m )
2652	mm = mm + 1
2653	ENDDO
2654	ENDIF
2655
2656	IF ( surf_match_lsm ) THEN
2657	mm = surf_lsm_h%start_index(jc,ic)
2658	DO m = surf_h(l)%start_index(j,i), &
2659	surf_h(l)%end_index(j,i)
2660	CALL restore_surface_elements( surf_lsm_h, &
2661	mm, surf_h(l), m )
2662	mm = mm + 1
2663	ENDDO
2664	ENDIF
2665
2666	IF ( surf_match_usm ) THEN
2667	mm = surf_usm_h%start_index(jc,ic)
2668	DO m = surf_h(l)%start_index(j,i), &
2669	surf_h(l)%end_index(j,i)
2670	CALL restore_surface_elements( surf_usm_h, &
2671	mm, surf_h(l), m )
2672	mm = mm + 1
2673	ENDDO
2674	ENDIF
2675
2676	jc = jc + 1
2677	ENDDO
2678	ic = ic + 1
2679	ENDDO
2680	ELSEIF ( vertical_surface ) THEN
2681	ic = nxlc
2682	DO i = nxlf, nxrf
2683	jc = nysc
2684	DO j = nysf, nynf
2685
2686	surf_match_def = surf_def_v(l)%end_index(jc,ic) >= &
2687	surf_def_v(l)%start_index(jc,ic)
2688	surf_match_lsm = surf_lsm_v(l)%end_index(jc,ic) >= &
2689	surf_lsm_v(l)%start_index(jc,ic)
2690	surf_match_usm = surf_usm_v(l)%end_index(jc,ic) >= &
2691	surf_usm_v(l)%start_index(jc,ic)
2692
2693
2694
2695	IF ( surf_match_def ) THEN
2696	mm = surf_def_v(l)%start_index(jc,ic)
2697	DO m = surf_v(l)%start_index(j,i), &
2698	surf_v(l)%end_index(j,i)
2699	CALL restore_surface_elements( surf_def_v(l), &
2700	mm, surf_v(l), m )
2701	mm = mm + 1
2702	ENDDO
2703	ENDIF
2704
2705	IF ( surf_match_lsm ) THEN
2706	mm = surf_lsm_v(l)%start_index(jc,ic)
2707	DO m = surf_v(l)%start_index(j,i), &
2708	surf_v(l)%end_index(j,i)
2709	CALL restore_surface_elements( surf_lsm_v(l), &
2710	mm, surf_v(l), m )
2711	mm = mm + 1
2712	ENDDO
2713	ENDIF
2714
2715	IF ( surf_match_usm ) THEN
2716	mm = surf_usm_v(l)%start_index(jc,ic)
2717	DO m = surf_v(l)%start_index(j,i), &
2718	surf_v(l)%end_index(j,i)
2719	CALL restore_surface_elements( surf_usm_v(l), &
2720	mm, surf_v(l), m )
2721	mm = mm + 1
2722	ENDDO
2723	ENDIF
2724
2725	jc = jc + 1
2726	ENDDO
2727	ic = ic + 1
2728	ENDDO
2729	ENDIF
2730
2731	ENDDO
2732
2733	READ ( 13 ) field_chr
2734
2735	ENDDO
2736
2737	ENDIF
2738
2739
2740	CONTAINS
2741	!------------------------------------------------------------------------------!
2742	! Description:
2743	! ------------
2744	!> Restores surfacle elements back on its respective type.
2745	!------------------------------------------------------------------------------!
2746	SUBROUTINE restore_surface_elements( surf_target, m_target, &
2747	surf_file, m_file )
2748
2749	IMPLICIT NONE
2750
2751	INTEGER(iwp) :: m_file !< respective surface-element index of current surface array
2752	INTEGER(iwp) :: m_target !< respecitve surface-element index of surface array on file
2753
2754	TYPE( surf_type ) :: surf_target !< target surface type
2755	TYPE( surf_type ) :: surf_file !< surface type on file
2756
2757	IF ( SCAN( TRIM( field_chr ), '%us' ) /= 0 ) THEN
2758	IF ( ALLOCATED( surf_target%us ) .AND. &
2759	ALLOCATED( surf_file%us ) ) &
2760	surf_target%us(m_target) = surf_file%us(m_file)
2761	ENDIF
2762
2763	IF ( SCAN( TRIM( field_chr ), '%ol' ) /= 0 ) THEN
2764	IF ( ALLOCATED( surf_target%ol ) .AND. &
2765	ALLOCATED( surf_file%ol ) ) &
2766	surf_target%ol(m_target) = surf_file%ol(m_file)
2767	ENDIF
2768
2769	IF ( SCAN( TRIM( field_chr ), '%usws' ) /= 0 ) THEN
2770	IF ( ALLOCATED( surf_target%usws ) .AND. &
2771	ALLOCATED( surf_file%usws ) ) &
2772	surf_target%usws(m_target) = surf_file%usws(m_file)
2773	ENDIF
2774
2775	IF ( SCAN( TRIM( field_chr ), '%vsws' ) /= 0 ) THEN
2776	IF ( ALLOCATED( surf_target%vsws ) .AND. &
2777	ALLOCATED( surf_file%vsws ) ) &
2778	surf_target%vsws(m_target) = surf_file%vsws(m_file)
2779	ENDIF
2780
2781	IF ( SCAN( TRIM( field_chr ), '%ts' ) /= 0 ) THEN
2782	IF ( ALLOCATED( surf_target%ts ) .AND. &
2783	ALLOCATED( surf_file%ts ) ) &
2784	surf_target%ts(m_target) = surf_file%ts(m_file)
2785	ENDIF
2786
2787	IF ( SCAN( TRIM( field_chr ), '%shf' ) /= 0 ) THEN
2788	IF ( ALLOCATED( surf_target%shf ) .AND. &
2789	ALLOCATED( surf_file%shf ) ) &
2790	surf_target%shf(m_target) = surf_file%shf(m_file)
2791	ENDIF
2792
2793	IF ( SCAN( TRIM( field_chr ), '%qs' ) /= 0 ) THEN
2794	IF ( ALLOCATED( surf_target%qs ) .AND. &
2795	ALLOCATED( surf_file%qs ) ) &
2796	surf_target%qs(m_target) = surf_file%qs(m_file)
2797	ENDIF
2798
2799	IF ( SCAN( TRIM( field_chr ), '%qsws' ) /= 0 ) THEN
2800	IF ( ALLOCATED( surf_target%qsws ) .AND. &
2801	ALLOCATED( surf_file%qsws ) ) &
2802	surf_target%qsws(m_target) = surf_file%qsws(m_file)
2803	ENDIF
2804
2805	IF ( SCAN( TRIM( field_chr ), '%ss' ) /= 0 ) THEN
2806	IF ( ALLOCATED( surf_target%ss ) .AND. &
2807	ALLOCATED( surf_file%ss ) ) &
2808	surf_target%ss(m_target) = surf_file%ss(m_file)
2809	ENDIF
2810
2811	IF ( SCAN( TRIM( field_chr ), '%ssws' ) /= 0 ) THEN
2812	IF ( ALLOCATED( surf_target%ssws ) .AND. &
2813	ALLOCATED( surf_file%ssws ) ) &
2814	surf_target%ssws(m_target) = surf_file%ssws(m_file)
2815	ENDIF
2816
2817	IF ( SCAN( TRIM( field_chr ), '%qrs' ) /= 0 ) THEN
2818	IF ( ALLOCATED( surf_target%qrs ) .AND. &
2819	ALLOCATED( surf_file%qrs ) ) &
2820	surf_target%qrs(m_target) = surf_file%qrs(m_file)
2821	ENDIF
2822
2823	IF ( SCAN( TRIM( field_chr ), '%qrsws' ) /= 0 ) THEN
2824	IF ( ALLOCATED( surf_target%qrsws ) .AND. &
2825	ALLOCATED( surf_file%qrsws ) ) &
2826	surf_target%qrsws(m_target) = surf_file%qrsws(m_file)
2827	ENDIF
2828
2829	IF ( SCAN( TRIM( field_chr ), '%nrs' ) /= 0 ) THEN
2830	IF ( ALLOCATED( surf_target%nrs ) .AND. &
2831	ALLOCATED( surf_file%nrs ) ) &
2832	surf_target%nrs(m_target) = surf_file%nrs(m_file)
2833	ENDIF
2834
2835	IF ( SCAN( TRIM( field_chr ), '%nrsws' ) /= 0 ) THEN
2836	IF ( ALLOCATED( surf_target%nrsws ) .AND. &
2837	ALLOCATED( surf_file%nrsws ) ) &
2838	surf_target%nrsws(m_target) = surf_file%nrsws(m_file)
2839	ENDIF
2840
2841	IF ( SCAN( TRIM( field_chr ), '%sasws' ) /= 0 ) THEN
2842	IF ( ALLOCATED( surf_target%sasws ) .AND. &
2843	ALLOCATED( surf_file%sasws ) ) &
2844	surf_target%sasws(m_target) = surf_file%sasws(m_file)
2845	ENDIF
2846
2847	IF ( SCAN( TRIM( field_chr ), '%mom_uv' ) /= 0 ) THEN
2848	IF ( ALLOCATED( surf_target%mom_flux_uv ) .AND. &
2849	ALLOCATED( surf_file%mom_flux_uv ) ) &
2850	surf_target%mom_flux_uv(m_target) = &
2851	surf_file%mom_flux_uv(m_file)
2852	ENDIF
2853
2854	IF ( SCAN( TRIM( field_chr ), '%mom_w' ) /= 0 ) THEN
2855	IF ( ALLOCATED( surf_target%mom_flux_w ) .AND. &
2856	ALLOCATED( surf_file%mom_flux_w ) ) &
2857	surf_target%mom_flux_w(m_target) = &
2858	surf_file%mom_flux_w(m_file)
2859	ENDIF
2860
2861	IF ( SCAN( TRIM( field_chr ), '%mom_tke' ) /= 0 ) THEN
2862	IF ( ALLOCATED( surf_target%mom_flux_tke ) .AND. &
2863	ALLOCATED( surf_file%mom_flux_tke ) ) &
2864	surf_target%mom_flux_tke(0:1,m_target) = &
2865	surf_file%mom_flux_tke(0:1,m_file)
2866	ENDIF
2867
2868	END SUBROUTINE restore_surface_elements
2869
2870	END SUBROUTINE surface_read_restart_data
2871
2872
2873	!------------------------------------------------------------------------------!
2874	! Description:
2875	! ------------
2876	!> Counts the number of surface elements with the same facing, required for
2877	!> reading and writing restart data.
2878	!------------------------------------------------------------------------------!
2879	SUBROUTINE surface_last_actions
2880
2881	IMPLICIT NONE
2882	!
2883	!-- Horizontal surfaces
2884	ns_h_on_file(0) = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns
2885	ns_h_on_file(1) = surf_def_h(1)%ns
2886	ns_h_on_file(2) = surf_def_h(2)%ns
2887	!
2888	!-- Vertical surfaces
2889	ns_v_on_file(0) = surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns
2890	ns_v_on_file(1) = surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns
2891	ns_v_on_file(2) = surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns
2892	ns_v_on_file(3) = surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns
2893
2894	END SUBROUTINE surface_last_actions
2895
2896
2897	END MODULE surface_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |