Home

Context Navigation

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

Last change on this file since 2273 was 2270, checked in by maronga, 7 years ago
major revisions in land surface model
Property svn:keywords set to `Id`
File size: 146.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |