Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |