Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |