Home

Context Navigation

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

Last change on this file since 2696 was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
Property svn:keywords set to `Id`
File size: 198.8 KB

Line
1	!> @file surface_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: surface_mod.f90 2696 2017-12-14 17:12:51Z kanani $
28	! - Missing code for css added to surf_*, handling of surface_csflux updated (FK)
29	! - Bugfixes in reading/writing restart data in case several surface types are
30	! present at the same time (MS)
31	! - Implementation of chemistry module (FK)
32	! - Allocation of pt1 and qv1 now done for all surface types (MS)
33	! - Revised classification of surface types
34	! - Introduce offset values to simplify index determination of surface elements
35	! - Dimensions of function get_topo_top_index (MS)
36	! - added variables for USM
37	! - adapted to changes in USM (RvT)
38	!
39	! 2688 2017-12-12 17:27:04Z Giersch
40	! Allocation and initialization of the latent heat flux (qsws) at the top of
41	! the ocean domain in case of coupled runs. In addtion, a double if-query has
42	! been removed.
43	!
44	! 2638 2017-11-23 12:44:23Z raasch
45	! bugfix for constant top momentumflux
46	!
47	! 2575 2017-10-24 09:57:58Z maronga
48	! Pavement parameterization revised
49	!
50	! 2547 2017-10-16 12:41:56Z schwenkel
51	! extended by cloud_droplets option
52	!
53	! 2508 2017-10-02 08:57:09Z suehring
54	! Minor formatting adjustment
55	!
56	! 2478 2017-09-18 13:37:24Z suehring
57	! Bugfixes in initializing model top
58	!
59	! 2378 2017-08-31 13:57:27Z suehring
60	! Bugfix in write restart data
61	!
62	! 2339 2017-08-07 13:55:26Z gronemeier
63	! corrected timestamp in header
64	!
65	! 2338 2017-08-07 12:15:38Z gronemeier
66	! Modularize 1D model
67	!
68	! 2318 2017-07-20 17:27:44Z suehring
69	! New function to obtain topography top index.
70	!
71	! 2317 2017-07-20 17:27:19Z suehring
72	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
73	! includes two more prognostic equations for cloud drop concentration (nc)
74	! and cloud water content (qc).
75	!
76	! 2270 2017-06-09 12:18:47Z maronga
77	! Parameters removed/added due to changes in the LSM
78	!
79	! 2269 2017-06-09 11:57:32Z suehring
80	! Formatting and description adjustments
81	!
82	! 2256 2017-06-07 13:58:08Z suehring
83	! Enable heating at downward-facing surfaces
84	!
85	! 2233 2017-05-30 18:08:54Z suehring
86	! Initial revision
87	!
88	!
89	! Description:
90	! ------------
91	!> Surface module defines derived data structures to treat surface-
92	!> bounded grid cells. Three different types of surfaces are defined:
93	!> default surfaces, natural surfaces, and urban surfaces. The module
94	!> encompasses the allocation and initialization of surface arrays, and handles
95	!> reading and writing restart data.
96	!> In addition, a further derived data structure is defined, in order to set
97	!> boundary conditions at surfaces.
98	!> @todo For the moment, downward-facing surfaces are only classified as
99	!> default type
100	!> @todo Clean up urban-surface variables (some of them are not used any more)
101	!> @todo Revise chemistry surface flux part (reduce loops?!)
102	!------------------------------------------------------------------------------!
103	MODULE surface_mod
104
105	USE arrays_3d, &
106	ONLY: zu, zw, heatflux_input_conversion, waterflux_input_conversion, &
107	momentumflux_input_conversion
108
109	#if defined( __chem )
110	USE chem_modules
111	#endif
112
113	USE control_parameters
114
115	USE indices, &
116	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, wall_flags_0
117
118	USE grid_variables, &
119	ONLY: dx, dy
120
121	USE kinds
122
123	USE model_1d_mod, &
124	ONLY: rif1d, us1d, usws1d, vsws1d
125
126
127	IMPLICIT NONE
128
129	!
130	!-- Data type used to identify grid-points where horizontal boundary conditions
131	!-- are applied
132	TYPE bc_type
133
134	INTEGER(iwp) :: ns !< number of surface elements on the PE
135
136	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid
137	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid
138	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid
139
140	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
141	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< end index within surface data type for given (j,i)
142
143	END TYPE bc_type
144	!
145	!-- Data type used to identify and treat surface-bounded grid points
146	TYPE surf_type
147
148	INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
149	INTEGER(iwp) :: joff !< offset value in y-direction, used to determine index of surface element
150	INTEGER(iwp) :: koff !< offset value in z-direction, used to determine index of surface element
151	INTEGER(iwp) :: ns !< number of surface elements on the PE
152
153	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid
154	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid
155	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid
156
157	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: facing !< Bit indicating surface orientation
158
159	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i)
160	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< End index within surface data type for given (j,i)
161
162	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_mo !< surface-layer height
163	REAL(wp), DIMENSION(:), ALLOCATABLE :: uvw_abs !< absolute surface-parallel velocity
164	REAL(wp), DIMENSION(:), ALLOCATABLE :: us !< friction velocity
165	REAL(wp), DIMENSION(:), ALLOCATABLE :: ts !< scaling parameter temerature
166	REAL(wp), DIMENSION(:), ALLOCATABLE :: qs !< scaling parameter humidity
167	REAL(wp), DIMENSION(:), ALLOCATABLE :: ss !< scaling parameter passive scalar
168	REAL(wp), DIMENSION(:), ALLOCATABLE :: qcs !< scaling parameter qc
169	REAL(wp), DIMENSION(:), ALLOCATABLE :: ncs !< scaling parameter nc
170	REAL(wp), DIMENSION(:), ALLOCATABLE :: qrs !< scaling parameter qr
171	REAL(wp), DIMENSION(:), ALLOCATABLE :: nrs !< scaling parameter nr
172
173	REAL(wp), DIMENSION(:), ALLOCATABLE :: ol !< Obukhov length
174	REAL(wp), DIMENSION(:), ALLOCATABLE :: rib !< Richardson bulk number
175
176	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0 !< roughness length for momentum
177	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0h !< roughness length for heat
178	REAL(wp), DIMENSION(:), ALLOCATABLE :: z0q !< roughness length for humidity
179
180	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt1 !< Potential temperature at first grid level
181	REAL(wp), DIMENSION(:), ALLOCATABLE :: qv1 !< Specific humidity at first grid level
182	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: css !< scaling parameter chemical species
183	!
184	!-- Define arrays for surface fluxes
185	REAL(wp), DIMENSION(:), ALLOCATABLE :: usws !< vertical momentum flux for u-component at horizontal surfaces
186	REAL(wp), DIMENSION(:), ALLOCATABLE :: vsws !< vertical momentum flux for v-component at horizontal surfaces
187
188	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf !< surface flux sensible heat
189	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws !< surface flux latent heat
190	REAL(wp), DIMENSION(:), ALLOCATABLE :: ssws !< surface flux passive scalar
191	REAL(wp), DIMENSION(:), ALLOCATABLE :: qcsws !< surface flux qc
192	REAL(wp), DIMENSION(:), ALLOCATABLE :: ncsws !< surface flux nc
193	REAL(wp), DIMENSION(:), ALLOCATABLE :: qrsws !< surface flux qr
194	REAL(wp), DIMENSION(:), ALLOCATABLE :: nrsws !< surface flux nr
195	REAL(wp), DIMENSION(:), ALLOCATABLE :: sasws !< surface flux salinity
196
197	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: cssws !< surface flux chemical species
198	!
199	!-- Required for horizontal walls in production_e
200	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_0 !< virtual velocity component (see production_e_init for further explanation)
201	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_0 !< virtual velocity component (see production_e_init for further explanation)
202
203	REAL(wp), DIMENSION(:), ALLOCATABLE :: mom_flux_uv !< momentum flux usvs and vsus at vertical surfaces (used in diffusion_u and diffusion_v)
204	REAL(wp), DIMENSION(:), ALLOCATABLE :: mom_flux_w !< momentum flux wsus and wsvs at vertical surfaces (used in diffusion_w)
205	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mom_flux_tke !< momentum flux usvs, vsus, wsus, wsvs at vertical surfaces at grid center (used in production_e)
206	!
207	!-- Variables required for LSM as well as for USM
208	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzt_pavement !< top index for pavement in soil
209	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: albedo_type !< albedo type, for each fraction (wall,green,window or vegetation,pavement water)
210
211	LOGICAL, DIMENSION(:), ALLOCATABLE :: building_surface !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
212	LOGICAL, DIMENSION(:), ALLOCATABLE :: building_covered !< flag indicating that buildings are on top of orography, only used for vertical surfaces in LSM
213	LOGICAL, DIMENSION(:), ALLOCATABLE :: pavement_surface !< flag parameter for pavements
214	LOGICAL, DIMENSION(:), ALLOCATABLE :: water_surface !< flag parameter for water surfaces
215	LOGICAL, DIMENSION(:), ALLOCATABLE :: vegetation_surface !< flag parameter for natural land surfaces
216
217	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: albedo !< broadband albedo for each surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
218	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: emissivity !< emissivity of the surface, for each fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
219	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: frac !< relative surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
220
221	REAL(wp), DIMENSION(:), ALLOCATABLE :: aldif !< albedo for longwave diffusive radiation, solar angle of 60Â°
222	REAL(wp), DIMENSION(:), ALLOCATABLE :: aldir !< albedo for longwave direct radiation, solar angle of 60Â°
223	REAL(wp), DIMENSION(:), ALLOCATABLE :: asdif !< albedo for shortwave diffusive radiation, solar angle of 60Â°
224	REAL(wp), DIMENSION(:), ALLOCATABLE :: asdir !< albedo for shortwave direct radiation, solar angle of 60Â°
225	REAL(wp), DIMENSION(:), ALLOCATABLE :: rrtm_aldif !< albedo for longwave diffusive radiation, solar angle of 60Â°
226	REAL(wp), DIMENSION(:), ALLOCATABLE :: rrtm_aldir !< albedo for longwave direct radiation, solar angle of 60Â°
227	REAL(wp), DIMENSION(:), ALLOCATABLE :: rrtm_asdif !< albedo for shortwave diffusive radiation, solar angle of 60Â°
228	REAL(wp), DIMENSION(:), ALLOCATABLE :: rrtm_asdir !< albedo for shortwave direct radiation, solar angle of 60Â°
229
230	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_surface !< skin-surface temperature
231	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net !< net radiation
232	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< net radiation, used in USM
233	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h !< heat conductivity of soil/ wall (W/m/K)
234	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h_green !< heat conductivity of green soil (W/m/K)
235	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h_window !< heat conductivity of windows (W/m/K)
236	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h_def !< default heat conductivity of soil (W/m/K)
237
238	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_lw_in !< incoming longwave radiation
239	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_lw_out !< emitted longwave radiation
240	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_sw_in !< incoming shortwave radiation
241	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_sw_out !< emitted shortwave radiation
242
243
244
245	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_liq !< liquid water coverage (of vegetated area)
246	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_veg !< vegetation coverage
247	REAL(wp), DIMENSION(:), ALLOCATABLE :: f_sw_in !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
248	REAL(wp), DIMENSION(:), ALLOCATABLE :: ghf !< ground heat flux
249	REAL(wp), DIMENSION(:), ALLOCATABLE :: g_d !< coefficient for dependence of r_canopy on water vapour pressure deficit
250	REAL(wp), DIMENSION(:), ALLOCATABLE :: lai !< leaf area index
251	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surface_u !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
252	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surface_s !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
253	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_liq !< surface flux of latent heat (liquid water portion)
254	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_soil !< surface flux of latent heat (soil portion)
255	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_veg !< surface flux of latent heat (vegetation portion)
256	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_a !< aerodynamic resistance
257	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_canopy !< canopy resistance
258	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_soil !< soil resistance
259	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_soil_min !< minimum soil resistance
260	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_s !< total surface resistance (combination of r_soil and r_canopy)
261	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_canopy_min !< minimum canopy (stomatal) resistance
262
263	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: alpha_vg !< coef. of Van Genuchten
264	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_w !< hydraulic diffusivity of soil (?)
265	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gamma_w !< hydraulic conductivity of soil (W/m/K)
266	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gamma_w_sat !< hydraulic conductivity at saturation
267	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: l_vg !< coef. of Van Genuchten
268	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_fc !< soil moisture at field capacity (m3/m3)
269	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_res !< residual soil moisture
270	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_sat !< saturation soil moisture (m3/m3)
271	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m_wilt !< soil moisture at permanent wilting point (m3/m3)
272	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_vg !< coef. Van Genuchten
273	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_total_def !< default volumetric heat capacity of the (soil) layer (J/m3/K)
274	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_total !< volumetric heat capacity of the actual soil matrix (J/m3/K)
275	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: root_fr !< root fraction within the soil layers
276	!
277	!-- Urban surface variables
278	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_types !< array of types of wall parameters
279
280	LOGICAL, DIMENSION(:), ALLOCATABLE :: isroof_surf !< flag indicating roof surfaces
281	LOGICAL, DIMENSION(:), ALLOCATABLE :: ground_level !< flag indicating ground floor level surfaces
282
283	REAL(wp), DIMENSION(:), ALLOCATABLE :: target_temp_summer !< indoor target temperature summer
284	REAL(wp), DIMENSION(:), ALLOCATABLE :: target_temp_winter !< indoor target temperature summer
285
286	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface !< heat capacity of the wall surface skin (J/m2/K)
287	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface_green !< heat capacity of the green surface skin (J/m2/K)
288	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface_window !< heat capacity of the window surface skin (J/m2/K)
289	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf !< heat conductivity between air and surface (W/m2/K)
290	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf_green !< heat conductivity between air and green surface (W/m2/K)
291	REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf_window !< heat conductivity between air and window surface (W/m2/K)
292	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_wall !< thickness of the wall, roof and soil layers
293	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_green !< thickness of the green wall, roof and soil layers
294	REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_window !< thickness of the window wall, roof and soil layers
295	REAL(wp), DIMENSION(:), ALLOCATABLE :: transmissivity !< transmissivity of windows
296
297	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< reflected shortwave radiation for local surface in i-th reflection
298	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< reflected + emitted longwave radiation for local surface in i-th reflection
299	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf !< total radiation flux incoming to minus outgoing from local surface
300
301	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K)
302	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_window_m !< window surface temperature tendency (K)
303	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_green_m !< green surface temperature tendency (K)
304	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf !< kinematic wall heat flux of sensible heat (actually no longer needed)
305	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb !< wall heat flux of sensible heat in wall normal direction
306
307	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb !< wall ground heat flux
308	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_window !< window ground heat flux
309	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_green !< green ground heat flux
310	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb !< indoor wall ground heat flux
311	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb_window !< indoor window ground heat flux
312
313	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_lw_out_change_0
314
315	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< shortwave radiation falling to local surface including radiation from reflections
316	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< total shortwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
317	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< longwave radiation falling to local surface including radiation from reflections
318	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< total longwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
319
320
321	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_wall !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)
322	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall !< wall grid spacing (center-center)
323	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall !< 1/dz_wall
324	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall_stag !< wall grid spacing (edge-edge)
325	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall_stag !< 1/dz_wall_stag
326	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array
327	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw !< wall layer depths (m)
328	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_window !< volumetric heat capacity of the window material ( J m-3 K-1 ) (= 2.19E6)
329	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_window !< window grid spacing (center-center)
330	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_window !< 1/dz_window
331	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_window_stag !< window grid spacing (edge-edge)
332	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_window_stag !< 1/dz_window_stag
333	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_window_m !< t_window prognostic array
334	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw_window !< window layer depths (m)
335	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_green !< volumetric heat capacity of the green material ( J m-3 K-1 ) (= 2.19E6)
336	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_green !< green grid spacing (center-center)
337	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_green !< 1/dz_green
338	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_green_stag !< green grid spacing (edge-edge)
339	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_green_stag !< 1/dz_green_stag
340	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_green_m !< t_green prognostic array
341	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw_green !< green layer depths (m)
342
343
344	!-- arrays for time averages
345	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l
346	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
347	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
348	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
349	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
350	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
351	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
352	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
353	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
354	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
355	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
356	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
357	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface
358	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_av !< average of wghf_eb
359	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_window_av !< average of wghf_eb window
360	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_green_av !< average of wghf_eb window
361	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb_av !< indoor average of wghf_eb
362	REAL(wp), DIMENSION(:), ALLOCATABLE :: iwghf_eb_window_av !< indoor average of wghf_eb window
363	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb_av !< average of wshf_eb
364	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_av !< average of wall surface temperature (K)
365	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_window_av !< average of window surface temperature (K)
366	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_green_av !< average of green wall surface temperature (K)
367	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_whole_av !< average of whole wall surface temperature (K)
368	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_surf_10cm_av !< average of the near surface temperature (K)
369
370	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_wall_av !< Average of t_wall
371	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_window_av !< Average of t_window
372	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_green_av !< Average of t_green
373
374	END TYPE surf_type
375
376	TYPE (bc_type), DIMENSION(0:1) :: bc_h !< boundary condition data type, horizontal upward- and downward facing surfaces
377
378	TYPE (surf_type), DIMENSION(0:2), TARGET :: surf_def_h !< horizontal default surfaces (Up, Down, and Top)
379	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_def_v !< vertical default surfaces (North, South, West, East)
380	TYPE (surf_type) , TARGET :: surf_lsm_h !< horizontal natural land surfaces, so far only upward-facing
381	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_lsm_v !< vertical land surfaces (North, South, West, East)
382	TYPE (surf_type) , TARGET :: surf_usm_h !< horizontal urban surfaces, so far only upward-facing
383	TYPE (surf_type), DIMENSION(0:3), TARGET :: surf_usm_v !< vertical urban surfaces (North, South, West, East)
384
385	INTEGER(iwp) :: ns_h_on_file(0:2) !< total number of horizontal surfaces with the same facing, required for writing restart data
386	INTEGER(iwp) :: ns_v_on_file(0:3) !< total number of vertical surfaces with the same facing, required for writing restart data
387
388
389	SAVE
390
391	PRIVATE
392
393	INTERFACE get_topography_top_index
394	MODULE PROCEDURE get_topography_top_index
395	MODULE PROCEDURE get_topography_top_index_ji
396	END INTERFACE get_topography_top_index
397
398	INTERFACE init_bc
399	MODULE PROCEDURE init_bc
400	END INTERFACE init_bc
401
402	INTERFACE init_surfaces
403	MODULE PROCEDURE init_surfaces
404	END INTERFACE init_surfaces
405
406	INTERFACE init_surface_arrays
407	MODULE PROCEDURE init_surface_arrays
408	END INTERFACE init_surface_arrays
409
410	INTERFACE surface_read_restart_data
411	MODULE PROCEDURE surface_read_restart_data
412	END INTERFACE surface_read_restart_data
413
414	INTERFACE surface_write_restart_data
415	MODULE PROCEDURE surface_write_restart_data
416	END INTERFACE surface_write_restart_data
417
418	INTERFACE surface_last_actions
419	MODULE PROCEDURE surface_last_actions
420	END INTERFACE surface_last_actions
421
422	INTERFACE surface_restore_elements
423	MODULE PROCEDURE surface_restore_elements_1d
424	MODULE PROCEDURE surface_restore_elements_2d
425	END INTERFACE surface_restore_elements
426
427	!
428	!-- Public variables
429	PUBLIC bc_h, ns_h_on_file, ns_v_on_file, surf_def_h, surf_def_v, &
430	surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type
431	!
432	!-- Public subroutines and functions
433	PUBLIC get_topography_top_index, init_bc, init_surfaces, &
434	init_surface_arrays, surface_read_restart_data, &
435	surface_restore_elements, surface_write_restart_data, &
436	surface_last_actions
437
438
439	CONTAINS
440
441	!------------------------------------------------------------------------------!
442	! Description:
443	! ------------
444	!> Initialize data type for setting boundary conditions at horizontal surfaces.
445	!------------------------------------------------------------------------------!
446	SUBROUTINE init_bc
447
448	IMPLICIT NONE
449
450	INTEGER(iwp) :: i !<
451	INTEGER(iwp) :: j !<
452	INTEGER(iwp) :: k !<
453
454	INTEGER(iwp), DIMENSION(0:1) :: num_h !<
455	INTEGER(iwp), DIMENSION(0:1) :: num_h_kji !<
456	INTEGER(iwp), DIMENSION(0:1) :: start_index_h !<
457
458	!
459	!-- First of all, count the number of upward- and downward-facing surfaces
460	num_h = 0
461	DO i = nxlg, nxrg
462	DO j = nysg, nyng
463	DO k = nzb+1, nzt
464	!
465	!-- Check if current gridpoint belongs to the atmosphere
466	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
467	!
468	!-- Upward-facing
469	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) &
470	num_h(0) = num_h(0) + 1
471	!
472	!-- Downward-facing
473	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
474	num_h(1) = num_h(1) + 1
475	ENDIF
476	ENDDO
477	ENDDO
478	ENDDO
479	!
480	!-- Save the number of surface elements
481	bc_h(0)%ns = num_h(0)
482	bc_h(1)%ns = num_h(1)
483	!
484	!-- ALLOCATE data type variables
485	!-- Upward facing
486	ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
487	ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
488	ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
489	ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
490	ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg) )
491	bc_h(0)%start_index = 1
492	bc_h(0)%end_index = 0
493	!
494	!-- Downward facing
495	ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
496	ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
497	ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
498	ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
499	ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg) )
500	bc_h(1)%start_index = 1
501	bc_h(1)%end_index = 0
502	!
503	!-- Store the respective indices on data type
504	num_h(0:1) = 1
505	start_index_h(0:1) = 1
506	DO i = nxlg, nxrg
507	DO j = nysg, nyng
508
509	num_h_kji(0:1) = 0
510	DO k = nzb+1, nzt
511	!
512	!-- Check if current gridpoint belongs to the atmosphere
513	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
514	!
515	!-- Upward-facing
516	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
517	bc_h(0)%i(num_h(0)) = i
518	bc_h(0)%j(num_h(0)) = j
519	bc_h(0)%k(num_h(0)) = k
520	num_h_kji(0) = num_h_kji(0) + 1
521	num_h(0) = num_h(0) + 1
522	ENDIF
523	!
524	!-- Downward-facing
525	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
526	bc_h(1)%i(num_h(1)) = i
527	bc_h(1)%j(num_h(1)) = j
528	bc_h(1)%k(num_h(1)) = k
529	num_h_kji(1) = num_h_kji(1) + 1
530	num_h(1) = num_h(1) + 1
531	ENDIF
532	ENDIF
533	ENDDO
534	bc_h(0)%start_index(j,i) = start_index_h(0)
535	bc_h(0)%end_index(j,i) = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
536	start_index_h(0) = bc_h(0)%end_index(j,i) + 1
537
538	bc_h(1)%start_index(j,i) = start_index_h(1)
539	bc_h(1)%end_index(j,i) = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
540	start_index_h(1) = bc_h(1)%end_index(j,i) + 1
541	ENDDO
542	ENDDO
543
544
545	END SUBROUTINE init_bc
546
547
548	!------------------------------------------------------------------------------!
549	! Description:
550	! ------------
551	!> Initialize horizontal and vertical surfaces. Counts the number of default-,
552	!> natural and urban surfaces and allocates memory, respectively.
553	!------------------------------------------------------------------------------!
554	SUBROUTINE init_surface_arrays
555
556	IMPLICIT NONE
557
558	INTEGER(iwp) :: i !< running index x-direction
559	INTEGER(iwp) :: j !< running index y-direction
560	INTEGER(iwp) :: k !< running index z-direction
561	INTEGER(iwp) :: l !< index variable for surface facing
562	INTEGER(iwp) :: num_lsm_h !< number of horizontally-aligned natural surfaces
563	INTEGER(iwp) :: num_usm_h !< number of horizontally-aligned urban surfaces
564
565	INTEGER(iwp), DIMENSION(0:2) :: num_def_h !< number of horizontally-aligned default surfaces
566	INTEGER(iwp), DIMENSION(0:3) :: num_def_v !< number of vertically-aligned default surfaces
567	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< number of vertically-aligned natural surfaces
568	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< number of vertically-aligned urban surfaces
569
570	LOGICAL :: building !< flag indicating building grid point
571	LOGICAL :: terrain !< flag indicating natural terrain grid point
572
573	num_def_h = 0
574	num_def_v = 0
575	num_lsm_h = 0
576	num_lsm_v = 0
577	num_usm_h = 0
578	num_usm_v = 0
579	!
580	!-- Surfaces are classified according to the input data read from static
581	!-- input file. If no input file is present, all surfaces are classified
582	!-- either as natural, urban, or default, depending on the setting of
583	!-- land_surface and urban_surface. To control this, use the control
584	!-- flag topo_no_distinct
585	!
586	!-- Count number of horizontal surfaces on local domain
587	DO i = nxl, nxr
588	DO j = nys, nyn
589	DO k = nzb+1, nzt
590	!
591	!-- Check if current gridpoint belongs to the atmosphere
592	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
593	!
594	!-- Check if grid point adjoins to any upward-facing horizontal
595	!-- surface, e.g. the Earth surface, plane roofs, or ceilings.
596
597	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
598	!
599	!-- Determine flags indicating terrain or building.
600	terrain = BTEST( wall_flags_0(k-1,j,i), 5 ) .OR. &
601	topo_no_distinct
602	building = BTEST( wall_flags_0(k-1,j,i), 6 ) .OR. &
603	topo_no_distinct
604	!
605	!-- Land-surface type
606	IF ( land_surface .AND. terrain ) THEN
607	num_lsm_h = num_lsm_h + 1
608	!
609	!-- Urban surface tpye
610	ELSEIF ( urban_surface .AND. building ) THEN
611	num_usm_h = num_usm_h + 1
612	!
613	!-- Default-surface type
614	ELSEIF ( .NOT. land_surface .AND. &
615	.NOT. urban_surface ) THEN
616
617	num_def_h(0) = num_def_h(0) + 1
618	!
619	!-- Unclassifified surface-grid point. Give error message.
620	ELSE
621	WRITE( message_string, * ) &
622	'Unclassified upward-facing ' // &
623	'surface element at '// &
624	'grid point (k,j,i) = ', k, j, i
625	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
626	ENDIF
627
628	ENDIF
629	!
630	!-- Check for top-fluxes
631	IF ( k == nzt .AND. use_top_fluxes ) THEN
632	num_def_h(2) = num_def_h(2) + 1
633	!
634	!-- Check for any other downward-facing surface. So far only for
635	!-- default surface type.
636	ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
637	num_def_h(1) = num_def_h(1) + 1
638	ENDIF
639
640	ENDIF
641	ENDDO
642	ENDDO
643	ENDDO
644	!
645	!-- Count number of vertical surfaces on local domain
646	DO i = nxl, nxr
647	DO j = nys, nyn
648	DO k = nzb+1, nzt
649	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
650	!
651	!-- Northward-facing
652	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) ) THEN
653	!
654	!-- Determine flags indicating terrain or building
655
656	terrain = BTEST( wall_flags_0(k,j-1,i), 5 ) .OR. &
657	topo_no_distinct
658	building = BTEST( wall_flags_0(k,j-1,i), 6 ) .OR. &
659	topo_no_distinct
660	IF ( land_surface .AND. terrain ) THEN
661	num_lsm_v(0) = num_lsm_v(0) + 1
662	ELSEIF ( urban_surface .AND. building ) THEN
663	num_usm_v(0) = num_usm_v(0) + 1
664	!
665	!-- Default-surface type
666	ELSEIF ( .NOT. land_surface .AND. &
667	.NOT. urban_surface ) THEN
668	num_def_v(0) = num_def_v(0) + 1
669	!
670	!-- Unclassifified surface-grid point. Give error message.
671	ELSE
672	WRITE( message_string, * ) &
673	'Unclassified northward-facing ' // &
674	'surface element at '// &
675	'grid point (k,j,i) = ', k, j, i
676	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
677
678	ENDIF
679	ENDIF
680	!
681	!-- Southward-facing
682	IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) THEN
683	!
684	!-- Determine flags indicating terrain or building
685	terrain = BTEST( wall_flags_0(k,j+1,i), 5 ) .OR. &
686	topo_no_distinct
687	building = BTEST( wall_flags_0(k,j+1,i), 6 ) .OR. &
688	topo_no_distinct
689	IF ( land_surface .AND. terrain ) THEN
690	num_lsm_v(1) = num_lsm_v(1) + 1
691	ELSEIF ( urban_surface .AND. building ) THEN
692	num_usm_v(1) = num_usm_v(1) + 1
693	!
694	!-- Default-surface type
695	ELSEIF ( .NOT. land_surface .AND. &
696	.NOT. urban_surface ) THEN
697	num_def_v(1) = num_def_v(1) + 1
698	!
699	!-- Unclassifified surface-grid point. Give error message.
700	ELSE
701	WRITE( message_string, * ) &
702	'Unclassified southward-facing ' // &
703	'surface element at '// &
704	'grid point (k,j,i) = ', k, j, i
705	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
706
707	ENDIF
708	ENDIF
709	!
710	!-- Eastward-facing
711	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) ) THEN
712	!
713	!-- Determine flags indicating terrain or building
714	terrain = BTEST( wall_flags_0(k,j,i-1), 5 ) .OR. &
715	topo_no_distinct
716	building = BTEST( wall_flags_0(k,j,i-1), 6 ) .OR. &
717	topo_no_distinct
718	IF ( land_surface .AND. terrain ) THEN
719	num_lsm_v(2) = num_lsm_v(2) + 1
720	ELSEIF ( urban_surface .AND. building ) THEN
721	num_usm_v(2) = num_usm_v(2) + 1
722	!
723	!-- Default-surface type
724	ELSEIF ( .NOT. land_surface .AND. &
725	.NOT. urban_surface ) THEN
726	num_def_v(2) = num_def_v(2) + 1
727	!
728	!-- Unclassifified surface-grid point. Give error message.
729	ELSE
730	WRITE( message_string, * ) &
731	'Unclassified eastward-facing ' // &
732	'surface element at '// &
733	'grid point (k,j,i) = ', k, j, i
734	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
735
736	ENDIF
737	ENDIF
738	!
739	!-- Westward-facing
740	IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) THEN
741	!
742	!-- Determine flags indicating terrain or building
743	terrain = BTEST( wall_flags_0(k,j,i+1), 5 ) .OR. &
744	topo_no_distinct
745	building = BTEST( wall_flags_0(k,j,i+1), 6 ) .OR. &
746	topo_no_distinct
747	IF ( land_surface .AND. terrain ) THEN
748	num_lsm_v(3) = num_lsm_v(3) + 1
749	ELSEIF ( urban_surface .AND. building ) THEN
750	num_usm_v(3) = num_usm_v(3) + 1
751	!
752	!-- Default-surface type
753	ELSEIF ( .NOT. land_surface .AND. &
754	.NOT. urban_surface ) THEN
755	num_def_v(3) = num_def_v(3) + 1
756	!
757	!-- Unclassifified surface-grid point. Give error message.
758	ELSE
759	WRITE( message_string, * ) &
760	'Unclassified westward-facing ' // &
761	'surface element at '// &
762	'grid point (k,j,i) = ', k, j, i
763	CALL message( 'surface_mod', 'PA0999', 1, 2, 0, 6, 0 )
764
765	ENDIF
766	ENDIF
767	ENDIF
768	ENDDO
769	ENDDO
770	ENDDO
771
772	!
773	!-- Store number of surfaces per core.
774	!-- Horizontal surface, default type, upward facing
775	surf_def_h(0)%ns = num_def_h(0)
776	!
777	!-- Horizontal surface, default type, downward facing
778	surf_def_h(1)%ns = num_def_h(1)
779	!
780	!-- Horizontal surface, default type, top downward facing
781	surf_def_h(2)%ns = num_def_h(2)
782	!
783	!-- Horizontal surface, natural type, so far only upward-facing
784	surf_lsm_h%ns = num_lsm_h
785	!
786	!-- Horizontal surface, urban type, so far only upward-facing
787	surf_usm_h%ns = num_usm_h
788	!
789	!-- Vertical surface, default type, northward facing
790	surf_def_v(0)%ns = num_def_v(0)
791	!
792	!-- Vertical surface, default type, southward facing
793	surf_def_v(1)%ns = num_def_v(1)
794	!
795	!-- Vertical surface, default type, eastward facing
796	surf_def_v(2)%ns = num_def_v(2)
797	!
798	!-- Vertical surface, default type, westward facing
799	surf_def_v(3)%ns = num_def_v(3)
800	!
801	!-- Vertical surface, natural type, northward facing
802	surf_lsm_v(0)%ns = num_lsm_v(0)
803	!
804	!-- Vertical surface, natural type, southward facing
805	surf_lsm_v(1)%ns = num_lsm_v(1)
806	!
807	!-- Vertical surface, natural type, eastward facing
808	surf_lsm_v(2)%ns = num_lsm_v(2)
809	!
810	!-- Vertical surface, natural type, westward facing
811	surf_lsm_v(3)%ns = num_lsm_v(3)
812	!
813	!-- Vertical surface, urban type, northward facing
814	surf_usm_v(0)%ns = num_usm_v(0)
815	!
816	!-- Vertical surface, urban type, southward facing
817	surf_usm_v(1)%ns = num_usm_v(1)
818	!
819	!-- Vertical surface, urban type, eastward facing
820	surf_usm_v(2)%ns = num_usm_v(2)
821	!
822	!-- Vertical surface, urban type, westward facing
823	surf_usm_v(3)%ns = num_usm_v(3)
824	!
825	!-- Allocate required attributes for horizontal surfaces - default type.
826	!-- Upward-facing (l=0) and downward-facing (l=1).
827	DO l = 0, 1
828	CALL allocate_surface_attributes_h ( surf_def_h(l), nys, nyn, nxl, nxr )
829	ENDDO
830	!
831	!-- Allocate required attributes for model top
832	CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
833	!
834	!-- Allocate required attributes for horizontal surfaces - natural type.
835	CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
836	!
837	!-- Allocate required attributes for horizontal surfaces - urban type.
838	CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
839
840	!
841	!-- Allocate required attributes for vertical surfaces.
842	!-- Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2)
843	!-- and westward-facing (l=3).
844	!-- Default type.
845	DO l = 0, 3
846	CALL allocate_surface_attributes_v ( surf_def_v(l), &
847	nys, nyn, nxl, nxr )
848	ENDDO
849	!
850	!-- Natural type
851	DO l = 0, 3
852	CALL allocate_surface_attributes_v ( surf_lsm_v(l), &
853	nys, nyn, nxl, nxr )
854	ENDDO
855	!
856	!-- Urban type
857	DO l = 0, 3
858	CALL allocate_surface_attributes_v ( surf_usm_v(l), &
859	nys, nyn, nxl, nxr )
860	ENDDO
861
862	END SUBROUTINE init_surface_arrays
863
864	!------------------------------------------------------------------------------!
865	! Description:
866	! ------------
867	!> Allocating memory for upward and downward-facing horizontal surface types,
868	!> except for top fluxes.
869	!------------------------------------------------------------------------------!
870	SUBROUTINE allocate_surface_attributes_h( surfaces, &
871	nys_l, nyn_l, nxl_l, nxr_l )
872
873	IMPLICIT NONE
874
875	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
876	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
877	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
878	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
879
880	TYPE(surf_type) :: surfaces !< respective surface type
881
882	!
883	!-- Allocate arrays for start and end index of horizontal surface type
884	!-- for each (j,i)-grid point. This is required e.g. in diffion_x, which is
885	!-- called for each (j,i). In order to find the location where the
886	!-- respective flux is store within the surface-type, start- and end-
887	!-- index are stored for each (j,i). For example, each (j,i) can have
888	!-- several entries where fluxes for horizontal surfaces might be stored,
889	!-- e.g. for overhanging structures where several upward-facing surfaces
890	!-- might exist for given (j,i).
891	!-- If no surface of respective type exist at current (j,i), set indicies
892	!-- such that loop in diffusion routines will not be entered.
893	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
894	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
895	surfaces%start_index = 0
896	surfaces%end_index = -1
897	!
898	!-- Indices to locate surface element
899	ALLOCATE ( surfaces%i(1:surfaces%ns) )
900	ALLOCATE ( surfaces%j(1:surfaces%ns) )
901	ALLOCATE ( surfaces%k(1:surfaces%ns) )
902	!
903	!-- Surface-layer height
904	ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
905	!
906	!-- Surface orientation
907	ALLOCATE ( surfaces%facing(1:surfaces%ns) )
908	!
909	!-- Surface-parallel wind velocity
910	ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
911	!
912	!-- Roughness
913	ALLOCATE ( surfaces%z0(1:surfaces%ns) )
914	ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
915	ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
916	!
917	!-- Friction velocity
918	ALLOCATE ( surfaces%us(1:surfaces%ns) )
919	!
920	!-- Stability parameter
921	ALLOCATE ( surfaces%ol(1:surfaces%ns) )
922	!
923	!-- Bulk Richardson number
924	ALLOCATE ( surfaces%rib(1:surfaces%ns) )
925	!
926	!-- Vertical momentum fluxes of u and v
927	ALLOCATE ( surfaces%usws(1:surfaces%ns) )
928	ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
929	!
930	!-- Required in production_e
931	IF ( .NOT. constant_diffusion ) THEN
932	ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
933	ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
934	ENDIF
935	!
936	!-- Characteristic temperature and surface flux of sensible heat
937	ALLOCATE ( surfaces%ts(1:surfaces%ns) )
938	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
939	!
940	!-- surface temperature
941	ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
942	!
943	!-- Characteristic humidity and surface flux of latent heat
944	IF ( humidity ) THEN
945	ALLOCATE ( surfaces%qs(1:surfaces%ns) )
946	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
947	ENDIF
948	!
949	!-- Characteristic scalar and surface flux of scalar
950	IF ( passive_scalar ) THEN
951	ALLOCATE ( surfaces%ss(1:surfaces%ns) )
952	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
953	ENDIF
954	!
955	!-- Scaling parameter (cs*) and surface flux of chemical species
956	#if defined ( __chem )
957	IF ( air_chemistry ) THEN
958	ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
959	ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
960	ENDIF
961	#endif
962	!
963	!-- Arrays for storing potential temperature and
964	!-- specific humidity at first grid level
965	ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
966	ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
967	!
968	!--
969	IF ( cloud_physics .AND. microphysics_morrison) THEN
970	ALLOCATE ( surfaces%qcs(1:surfaces%ns) )
971	ALLOCATE ( surfaces%ncs(1:surfaces%ns) )
972	ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
973	ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
974	ENDIF
975	!
976	!--
977	IF ( cloud_physics .AND. microphysics_seifert) THEN
978	ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
979	ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
980	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
981	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
982	ENDIF
983	!
984	!-- Salinity surface flux
985	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
986
987	END SUBROUTINE allocate_surface_attributes_h
988
989
990	!------------------------------------------------------------------------------!
991	! Description:
992	! ------------
993	!> Allocating memory for model-top fluxes
994	!------------------------------------------------------------------------------!
995	SUBROUTINE allocate_surface_attributes_h_top( surfaces, &
996	nys_l, nyn_l, nxl_l, nxr_l )
997
998	IMPLICIT NONE
999
1000	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1001	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1002	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1003	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1004
1005	TYPE(surf_type) :: surfaces !< respective surface type
1006
1007	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1008	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1009	surfaces%start_index = 0
1010	surfaces%end_index = -1
1011	!
1012	!-- Indices to locate surface (model-top) element
1013	ALLOCATE ( surfaces%i(1:surfaces%ns) )
1014	ALLOCATE ( surfaces%j(1:surfaces%ns) )
1015	ALLOCATE ( surfaces%k(1:surfaces%ns) )
1016
1017	IF ( .NOT. constant_diffusion ) THEN
1018	ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
1019	ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
1020	ENDIF
1021	!
1022	!-- Vertical momentum fluxes of u and v
1023	ALLOCATE ( surfaces%usws(1:surfaces%ns) )
1024	ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
1025	!
1026	!-- Sensible heat flux
1027	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
1028	!
1029	!-- Latent heat flux
1030	IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere') THEN
1031	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1032	ENDIF
1033	!
1034	!-- Scalar flux
1035	IF ( passive_scalar ) THEN
1036	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1037	ENDIF
1038	!
1039	!-- Chemical species flux
1040	#if defined( __chem )
1041	IF ( air_chemistry ) THEN
1042	ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1043	ENDIF
1044	#endif
1045	!
1046	!--
1047	IF ( cloud_physics .AND. microphysics_morrison) THEN
1048	ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1049	ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1050	ENDIF
1051	!
1052	!--
1053	IF ( cloud_physics .AND. microphysics_seifert) THEN
1054	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1055	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1056	ENDIF
1057	!
1058	!-- Salinity flux
1059	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
1060
1061	END SUBROUTINE allocate_surface_attributes_h_top
1062
1063	!------------------------------------------------------------------------------!
1064	! Description:
1065	! ------------
1066	!> Allocating memory for vertical surface types.
1067	!------------------------------------------------------------------------------!
1068	SUBROUTINE allocate_surface_attributes_v( surfaces, &
1069	nys_l, nyn_l, nxl_l, nxr_l )
1070
1071	IMPLICIT NONE
1072
1073	INTEGER(iwp) :: nyn_l !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1074	INTEGER(iwp) :: nys_l !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1075	INTEGER(iwp) :: nxl_l !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1076	INTEGER(iwp) :: nxr_l !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1077
1078	TYPE(surf_type) :: surfaces !< respective surface type
1079
1080	!
1081	!-- Allocate arrays for start and end index of vertical surface type
1082	!-- for each (j,i)-grid point. This is required in diffion_x, which is
1083	!-- called for each (j,i). In order to find the location where the
1084	!-- respective flux is store within the surface-type, start- and end-
1085	!-- index are stored for each (j,i). For example, each (j,i) can have
1086	!-- several entries where fluxes for vertical surfaces might be stored.
1087	!-- In the flat case, where no vertical walls exit, set indicies such
1088	!-- that loop in diffusion routines will not be entered.
1089	ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1090	ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1091	surfaces%start_index = 0
1092	surfaces%end_index = -1
1093	!
1094	!-- Indices to locate surface element.
1095	ALLOCATE ( surfaces%i(1:surfaces%ns) )
1096	ALLOCATE ( surfaces%j(1:surfaces%ns) )
1097	ALLOCATE ( surfaces%k(1:surfaces%ns) )
1098	!
1099	!-- Surface-layer height
1100	ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
1101	!
1102	!-- Surface orientation
1103	ALLOCATE ( surfaces%facing(1:surfaces%ns) )
1104	!
1105	!-- Surface parallel wind velocity
1106	ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
1107	!
1108	!-- Roughness
1109	ALLOCATE ( surfaces%z0(1:surfaces%ns) )
1110	ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
1111	ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
1112
1113	!
1114	!-- Friction velocity
1115	ALLOCATE ( surfaces%us(1:surfaces%ns) )
1116	!
1117	!-- Allocate Obukhov length and bulk Richardson number. Actually, at
1118	!-- vertical surfaces these are only required for natural surfaces.
1119	!-- for natural land surfaces
1120	ALLOCATE( surfaces%ol(1:surfaces%ns) )
1121	ALLOCATE( surfaces%rib(1:surfaces%ns) )
1122	!
1123	!-- Allocate arrays for surface momentum fluxes for u and v. For u at north-
1124	!-- and south-facing surfaces, for v at east- and west-facing surfaces.
1125	ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
1126	!
1127	!-- Allocate array for surface momentum flux for w - wsus and wsvs
1128	ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) )
1129	!
1130	!-- Allocate array for surface momentum flux for subgrid-scale tke wsus and
1131	!-- wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
1132	!-- on surface.
1133	ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )
1134	!
1135	!-- Characteristic temperature and surface flux of sensible heat
1136	ALLOCATE ( surfaces%ts(1:surfaces%ns) )
1137	ALLOCATE ( surfaces%shf(1:surfaces%ns) )
1138	!
1139	!-- surface temperature
1140	ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
1141	!
1142	!-- Characteristic humidity and surface flux of latent heat
1143	IF ( humidity ) THEN
1144	ALLOCATE ( surfaces%qs(1:surfaces%ns) )
1145	ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1146	ENDIF
1147	!
1148	!-- Characteristic scalar and surface flux of scalar
1149	IF ( passive_scalar ) THEN
1150	ALLOCATE ( surfaces%ss(1:surfaces%ns) )
1151	ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1152	ENDIF
1153	!
1154	!-- Scaling parameter (cs*) and surface flux of chemical species
1155	#if defined( __chem )
1156	IF ( air_chemistry ) THEN
1157	ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
1158	ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1159	ENDIF
1160	#endif
1161	!
1162	!-- Arrays for storing potential temperature and
1163	!-- specific humidity at first grid level
1164	ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
1165	ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
1166
1167	IF ( cloud_physics .AND. microphysics_seifert) THEN
1168	ALLOCATE ( surfaces%qcs(1:surfaces%ns) )
1169	ALLOCATE ( surfaces%ncs(1:surfaces%ns) )
1170	ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1171	ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1172	ENDIF
1173
1174	IF ( cloud_physics .AND. microphysics_seifert) THEN
1175	ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
1176	ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
1177	ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1178	ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1179	ENDIF
1180	!
1181	!-- Salinity surface flux
1182	IF ( ocean ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
1183
1184	END SUBROUTINE allocate_surface_attributes_v
1185
1186	!------------------------------------------------------------------------------!
1187	! Description:
1188	! ------------
1189	!> Initialize surface elements, i.e. set initial values for surface fluxes,
1190	!> friction velocity, calcuation of start/end indices, etc. .
1191	!> Please note, further initialization concerning
1192	!> special surface characteristics, e.g. soil- and vegatation type,
1193	!> building type, etc., is done in the land-surface and urban-surface module,
1194	!> respectively.
1195	!------------------------------------------------------------------------------!
1196	SUBROUTINE init_surfaces
1197
1198	IMPLICIT NONE
1199
1200	INTEGER(iwp) :: i !< running index x-direction
1201	INTEGER(iwp) :: j !< running index y-direction
1202	INTEGER(iwp) :: k !< running index z-direction
1203	INTEGER(iwp) :: l !< index variable used to distinguish surface facing
1204	INTEGER(iwp) :: m !< running index surface elements
1205
1206	INTEGER(iwp) :: start_index_lsm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal natural surfaces
1207	INTEGER(iwp) :: start_index_usm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal urban surfaces
1208
1209	INTEGER(iwp) :: num_lsm_h !< current number of horizontal surface element, natural type
1210	INTEGER(iwp) :: num_lsm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal natural surfaces
1211	INTEGER(iwp) :: num_usm_h !< current number of horizontal surface element, urban type
1212	INTEGER(iwp) :: num_usm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban surfaces
1213
1214	INTEGER(iwp), DIMENSION(0:2) :: num_def_h !< current number of horizontal surface element, default type
1215	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
1216	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
1217
1218	INTEGER(iwp), DIMENSION(0:3) :: num_def_v !< current number of vertical surface element, default type
1219	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
1220	INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< current number of vertical surface element, natural type
1221	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
1222	INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< current number of vertical surface element, urban type
1223	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
1224
1225	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
1226	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
1227	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
1228
1229	LOGICAL :: building !< flag indicating building grid point
1230	LOGICAL :: terrain !< flag indicating natural terrain grid point
1231
1232	!
1233	!-- Set offset indices, i.e. index difference between surface element and
1234	!-- surface-bounded grid point.
1235	!-- Upward facing - no horizontal offsets
1236	surf_def_h(0:2)%ioff = 0
1237	surf_def_h(0:2)%joff = 0
1238
1239	surf_lsm_h%ioff = 0
1240	surf_lsm_h%joff = 0
1241
1242	surf_usm_h%ioff = 0
1243	surf_usm_h%joff = 0
1244	!
1245	!-- Upward facing vertical offsets
1246	surf_def_h(0)%koff = -1
1247	surf_lsm_h%koff = -1
1248	surf_usm_h%koff = -1
1249	!
1250	!-- Downward facing vertical offset
1251	surf_def_h(1:2)%koff = 1
1252	!
1253	!-- Vertical surfaces - no vertical offset
1254	surf_def_v(0:3)%koff = 0
1255	surf_lsm_v(0:3)%koff = 0
1256	surf_usm_v(0:3)%koff = 0
1257	!
1258	!-- North- and southward facing - no offset in x
1259	surf_def_v(0:1)%ioff = 0
1260	surf_lsm_v(0:1)%ioff = 0
1261	surf_usm_v(0:1)%ioff = 0
1262	!
1263	!-- Northward facing offset in y
1264	surf_def_v(0)%joff = -1
1265	surf_lsm_v(0)%joff = -1
1266	surf_usm_v(0)%joff = -1
1267	!
1268	!-- Southward facing offset in y
1269	surf_def_v(1)%joff = 1
1270	surf_lsm_v(1)%joff = 1
1271	surf_usm_v(1)%joff = 1
1272
1273	!
1274	!-- East- and westward facing - no offset in y
1275	surf_def_v(2:3)%joff = 0
1276	surf_lsm_v(2:3)%joff = 0
1277	surf_usm_v(2:3)%joff = 0
1278	!
1279	!-- Eastward facing offset in x
1280	surf_def_v(2)%ioff = -1
1281	surf_lsm_v(2)%ioff = -1
1282	surf_usm_v(2)%ioff = -1
1283	!
1284	!-- Westward facing offset in y
1285	surf_def_v(3)%ioff = 1
1286	surf_lsm_v(3)%ioff = 1
1287	surf_usm_v(3)%ioff = 1
1288
1289	!
1290	!-- Initialize surface attributes, store indicies, surfaces orientation, etc.,
1291	num_def_h(0:2) = 1
1292	num_def_v(0:3) = 1
1293
1294	num_lsm_h = 1
1295	num_lsm_v(0:3) = 1
1296
1297	num_usm_h = 1
1298	num_usm_v(0:3) = 1
1299
1300	start_index_def_h(0:2) = 1
1301	start_index_def_v(0:3) = 1
1302
1303	start_index_lsm_h = 1
1304	start_index_lsm_v(0:3) = 1
1305
1306	start_index_usm_h = 1
1307	start_index_usm_v(0:3) = 1
1308
1309	DO i = nxl, nxr
1310	DO j = nys, nyn
1311
1312	num_def_h_kji = 0
1313	num_def_v_kji = 0
1314	num_lsm_h_kji = 0
1315	num_lsm_v_kji = 0
1316	num_usm_h_kji = 0
1317	num_usm_v_kji = 0
1318
1319	DO k = nzb+1, nzt
1320	!
1321	!-- Check if current gridpoint belongs to the atmosphere
1322	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1323	!
1324	!-- Upward-facing surface. Distinguish between differet surface types.
1325	!-- To do, think about method to flag natural and non-natural
1326	!-- surfaces.
1327	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN
1328	!
1329	!-- Determine flags indicating terrain or building
1330	terrain = BTEST( wall_flags_0(k-1,j,i), 5 ) .OR. &
1331	topo_no_distinct
1332	building = BTEST( wall_flags_0(k-1,j,i), 6 ) .OR. &
1333	topo_no_distinct
1334	!
1335	!-- Natural surface type
1336	IF ( land_surface .AND. terrain ) THEN
1337	CALL initialize_horizontal_surfaces( k, j, i, &
1338	surf_lsm_h, &
1339	num_lsm_h, &
1340	num_lsm_h_kji, &
1341	.TRUE., .FALSE. )
1342	!
1343	!-- Urban surface tpye
1344	ELSEIF ( urban_surface .AND. building ) THEN
1345	CALL initialize_horizontal_surfaces( k, j, i, &
1346	surf_usm_h, &
1347	num_usm_h, &
1348	num_usm_h_kji, &
1349	.TRUE., .FALSE. )
1350	!
1351	!-- Default surface type
1352	ELSE
1353	CALL initialize_horizontal_surfaces( k, j, i, &
1354	surf_def_h(0), &
1355	num_def_h(0), &
1356	num_def_h_kji(0),&
1357	.TRUE., .FALSE. )
1358	ENDIF
1359	ENDIF
1360	!
1361	!-- downward-facing surface, first, model top. Please note,
1362	!-- for the moment, downward-facing surfaces are always of
1363	!-- default type
1364	IF ( k == nzt .AND. use_top_fluxes ) THEN
1365	CALL initialize_top( k, j, i, surf_def_h(2), &
1366	num_def_h(2), num_def_h_kji(2) )
1367	!
1368	!-- Check for any other downward-facing surface. So far only for
1369	!-- default surface type.
1370	ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN
1371	CALL initialize_horizontal_surfaces( k, j, i, &
1372	surf_def_h(1), &
1373	num_def_h(1), &
1374	num_def_h_kji(1), &
1375	.FALSE., .TRUE. )
1376	ENDIF
1377	!
1378	!-- Check for vertical walls and, if required, initialize it.
1379	! Start with northward-facing surface.
1380	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) ) THEN
1381	!
1382	!-- Determine flags indicating terrain or building
1383	terrain = BTEST( wall_flags_0(k,j-1,i), 5 ) .OR. &
1384	topo_no_distinct
1385	building = BTEST( wall_flags_0(k,j-1,i), 6 ) .OR. &
1386	topo_no_distinct
1387	IF ( urban_surface .AND. building ) THEN
1388	CALL initialize_vertical_surfaces( 0, k, j, i, &
1389	surf_usm_v(0), &
1390	num_usm_v(0), &
1391	num_usm_v_kji(0), &
1392	.FALSE., .FALSE., &
1393	.FALSE., .TRUE. )
1394	ELSEIF ( land_surface .AND. terrain ) THEN
1395	CALL initialize_vertical_surfaces( 0, k, j, i, &
1396	surf_lsm_v(0), &
1397	num_lsm_v(0), &
1398	num_lsm_v_kji(0), &
1399	.FALSE., .FALSE., &
1400	.FALSE., .TRUE. )
1401	ELSE
1402	CALL initialize_vertical_surfaces( 0, k, j, i, &
1403	surf_def_v(0), &
1404	num_def_v(0), &
1405	num_def_v_kji(0), &
1406	.FALSE., .FALSE., &
1407	.FALSE., .TRUE. )
1408	ENDIF
1409	ENDIF
1410	!
1411	!-- southward-facing surface
1412	IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) THEN
1413	!
1414	!-- Determine flags indicating terrain or building
1415	terrain = BTEST( wall_flags_0(k,j+1,i), 5 ) .OR. &
1416	topo_no_distinct
1417	building = BTEST( wall_flags_0(k,j+1,i), 6 ) .OR. &
1418	topo_no_distinct
1419	IF ( urban_surface .AND. building ) THEN
1420	CALL initialize_vertical_surfaces( 1, k, j, i, &
1421	surf_usm_v(1), &
1422	num_usm_v(1), &
1423	num_usm_v_kji(1), &
1424	.FALSE., .FALSE., &
1425	.TRUE., .FALSE. )
1426	ELSEIF ( land_surface .AND. terrain ) THEN
1427	CALL initialize_vertical_surfaces( 1, k, j, i, &
1428	surf_lsm_v(1), &
1429	num_lsm_v(1), &
1430	num_lsm_v_kji(1), &
1431	.FALSE., .FALSE., &
1432	.TRUE., .FALSE. )
1433	ELSE
1434	CALL initialize_vertical_surfaces( 1, k, j, i, &
1435	surf_def_v(1), &
1436	num_def_v(1), &
1437	num_def_v_kji(1), &
1438	.FALSE., .FALSE., &
1439	.TRUE., .FALSE. )
1440	ENDIF
1441	ENDIF
1442	!
1443	!-- eastward-facing surface
1444	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) ) THEN
1445	!
1446	!-- Determine flags indicating terrain or building
1447	terrain = BTEST( wall_flags_0(k,j,i-1), 5 ) .OR. &
1448	topo_no_distinct
1449	building = BTEST( wall_flags_0(k,j,i-1), 6 ) .OR. &
1450	topo_no_distinct
1451	IF ( urban_surface .AND. building ) THEN
1452	CALL initialize_vertical_surfaces( 2, k, j, i, &
1453	surf_usm_v(2), &
1454	num_usm_v(2), &
1455	num_usm_v_kji(2), &
1456	.TRUE., .FALSE., &
1457	.FALSE., .FALSE. )
1458	ELSEIF ( land_surface .AND. terrain ) THEN
1459	CALL initialize_vertical_surfaces( 2, k, j, i, &
1460	surf_lsm_v(2), &
1461	num_lsm_v(2), &
1462	num_lsm_v_kji(2), &
1463	.TRUE., .FALSE., &
1464	.FALSE., .FALSE. )
1465	ELSE
1466	CALL initialize_vertical_surfaces( 2, k, j, i, &
1467	surf_def_v(2), &
1468	num_def_v(2), &
1469	num_def_v_kji(2), &
1470	.TRUE., .FALSE., &
1471	.FALSE., .FALSE. )
1472	ENDIF
1473	ENDIF
1474	!
1475	!-- westward-facing surface
1476	IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) THEN
1477	!
1478	!-- Determine flags indicating terrain or building
1479	terrain = BTEST( wall_flags_0(k,j,i+1), 5 ) .OR. &
1480	topo_no_distinct
1481	building = BTEST( wall_flags_0(k,j,i+1), 6 ) .OR. &
1482	topo_no_distinct
1483	IF ( urban_surface .AND. building ) THEN
1484	CALL initialize_vertical_surfaces( 3, k, j, i, &
1485	surf_usm_v(3), &
1486	num_usm_v(3), &
1487	num_usm_v_kji(3), &
1488	.FALSE., .TRUE., &
1489	.FALSE., .FALSE. )
1490	ELSEIF ( land_surface .AND. terrain ) THEN
1491	CALL initialize_vertical_surfaces( 3, k, j, i, &
1492	surf_lsm_v(3), &
1493	num_lsm_v(3), &
1494	num_lsm_v_kji(3), &
1495	.FALSE., .TRUE., &
1496	.FALSE., .FALSE. )
1497	ELSE
1498	CALL initialize_vertical_surfaces( 3, k, j, i, &
1499	surf_def_v(3), &
1500	num_def_v(3), &
1501	num_def_v_kji(3), &
1502	.FALSE., .TRUE., &
1503	.FALSE., .FALSE. )
1504	ENDIF
1505	ENDIF
1506	ENDIF
1507
1508
1509	ENDDO
1510	!
1511	!-- Determine start- and end-index at grid point (j,i). Also, for
1512	!-- horizontal surfaces more than 1 horizontal surface element can
1513	!-- exist at grid point (j,i) if overhanging structures are present.
1514	!-- Upward-facing surfaces
1515	surf_def_h(0)%start_index(j,i) = start_index_def_h(0)
1516	surf_def_h(0)%end_index(j,i) = surf_def_h(0)%start_index(j,i) + &
1517	num_def_h_kji(0) - 1
1518	start_index_def_h(0) = surf_def_h(0)%end_index(j,i) + 1
1519	!
1520	!-- Downward-facing surfaces, except model top
1521	surf_def_h(1)%start_index(j,i) = start_index_def_h(1)
1522	surf_def_h(1)%end_index(j,i) = surf_def_h(1)%start_index(j,i) + &
1523	num_def_h_kji(1) - 1
1524	start_index_def_h(1) = surf_def_h(1)%end_index(j,i) + 1
1525	!
1526	!-- Downward-facing surfaces -- model top fluxes
1527	surf_def_h(2)%start_index(j,i) = start_index_def_h(2)
1528	surf_def_h(2)%end_index(j,i) = surf_def_h(2)%start_index(j,i) + &
1529	num_def_h_kji(2) - 1
1530	start_index_def_h(2) = surf_def_h(2)%end_index(j,i) + 1
1531	!
1532	!-- Horizontal natural land surfaces
1533	surf_lsm_h%start_index(j,i) = start_index_lsm_h
1534	surf_lsm_h%end_index(j,i) = surf_lsm_h%start_index(j,i) + &
1535	num_lsm_h_kji - 1
1536	start_index_lsm_h = surf_lsm_h%end_index(j,i) + 1
1537	!
1538	!-- Horizontal urban surfaces
1539	surf_usm_h%start_index(j,i) = start_index_usm_h
1540	surf_usm_h%end_index(j,i) = surf_usm_h%start_index(j,i) + &
1541	num_usm_h_kji - 1
1542	start_index_usm_h = surf_usm_h%end_index(j,i) + 1
1543
1544	!
1545	!-- Vertical surfaces - Default type
1546	surf_def_v(0)%start_index(j,i) = start_index_def_v(0)
1547	surf_def_v(1)%start_index(j,i) = start_index_def_v(1)
1548	surf_def_v(2)%start_index(j,i) = start_index_def_v(2)
1549	surf_def_v(3)%start_index(j,i) = start_index_def_v(3)
1550	surf_def_v(0)%end_index(j,i) = start_index_def_v(0) + &
1551	num_def_v_kji(0) - 1
1552	surf_def_v(1)%end_index(j,i) = start_index_def_v(1) + &
1553	num_def_v_kji(1) - 1
1554	surf_def_v(2)%end_index(j,i) = start_index_def_v(2) + &
1555	num_def_v_kji(2) - 1
1556	surf_def_v(3)%end_index(j,i) = start_index_def_v(3) + &
1557	num_def_v_kji(3) - 1
1558	start_index_def_v(0) = surf_def_v(0)%end_index(j,i) + 1
1559	start_index_def_v(1) = surf_def_v(1)%end_index(j,i) + 1
1560	start_index_def_v(2) = surf_def_v(2)%end_index(j,i) + 1
1561	start_index_def_v(3) = surf_def_v(3)%end_index(j,i) + 1
1562	!
1563	!-- Natural type
1564	surf_lsm_v(0)%start_index(j,i) = start_index_lsm_v(0)
1565	surf_lsm_v(1)%start_index(j,i) = start_index_lsm_v(1)
1566	surf_lsm_v(2)%start_index(j,i) = start_index_lsm_v(2)
1567	surf_lsm_v(3)%start_index(j,i) = start_index_lsm_v(3)
1568	surf_lsm_v(0)%end_index(j,i) = start_index_lsm_v(0) + &
1569	num_lsm_v_kji(0) - 1
1570	surf_lsm_v(1)%end_index(j,i) = start_index_lsm_v(1) + &
1571	num_lsm_v_kji(1) - 1
1572	surf_lsm_v(2)%end_index(j,i) = start_index_lsm_v(2) + &
1573	num_lsm_v_kji(2) - 1
1574	surf_lsm_v(3)%end_index(j,i) = start_index_lsm_v(3) + &
1575	num_lsm_v_kji(3) - 1
1576	start_index_lsm_v(0) = surf_lsm_v(0)%end_index(j,i) + 1
1577	start_index_lsm_v(1) = surf_lsm_v(1)%end_index(j,i) + 1
1578	start_index_lsm_v(2) = surf_lsm_v(2)%end_index(j,i) + 1
1579	start_index_lsm_v(3) = surf_lsm_v(3)%end_index(j,i) + 1
1580	!
1581	!-- Urban type
1582	surf_usm_v(0)%start_index(j,i) = start_index_usm_v(0)
1583	surf_usm_v(1)%start_index(j,i) = start_index_usm_v(1)
1584	surf_usm_v(2)%start_index(j,i) = start_index_usm_v(2)
1585	surf_usm_v(3)%start_index(j,i) = start_index_usm_v(3)
1586	surf_usm_v(0)%end_index(j,i) = start_index_usm_v(0) + &
1587	num_usm_v_kji(0) - 1
1588	surf_usm_v(1)%end_index(j,i) = start_index_usm_v(1) + &
1589	num_usm_v_kji(1) - 1
1590	surf_usm_v(2)%end_index(j,i) = start_index_usm_v(2) + &
1591	num_usm_v_kji(2) - 1
1592	surf_usm_v(3)%end_index(j,i) = start_index_usm_v(3) + &
1593	num_usm_v_kji(3) - 1
1594	start_index_usm_v(0) = surf_usm_v(0)%end_index(j,i) + 1
1595	start_index_usm_v(1) = surf_usm_v(1)%end_index(j,i) + 1
1596	start_index_usm_v(2) = surf_usm_v(2)%end_index(j,i) + 1
1597	start_index_usm_v(3) = surf_usm_v(3)%end_index(j,i) + 1
1598
1599
1600	ENDDO
1601	ENDDO
1602
1603
1604	CONTAINS
1605
1606	!------------------------------------------------------------------------------!
1607	! Description:
1608	! ------------
1609	!> Initialize horizontal surface elements, upward- and downward-facing.
1610	!> Note, horizontal surface type alsw comprises model-top fluxes, which are,
1611	!> initialized in a different routine.
1612	!------------------------------------------------------------------------------!
1613	SUBROUTINE initialize_horizontal_surfaces( k, j, i, surf, num_h, &
1614	num_h_kji, upward_facing, &
1615	downward_facing )
1616
1617	IMPLICIT NONE
1618
1619	INTEGER(iwp) :: i !< running index x-direction
1620	INTEGER(iwp) :: j !< running index y-direction
1621	INTEGER(iwp) :: k !< running index z-direction
1622	INTEGER(iwp) :: num_h !< current number of surface element
1623	INTEGER(iwp) :: num_h_kji !< dummy increment
1624	INTEGER(iwp) :: lsp !< running index chemical species
1625	INTEGER(iwp) :: lsp_pr !< running index chemical species??
1626
1627	LOGICAL :: upward_facing !< flag indicating upward-facing surface
1628	LOGICAL :: downward_facing !< flag indicating downward-facing surface
1629
1630	TYPE( surf_type ) :: surf !< respective surface type
1631	!
1632	!-- Store indices of respective surface element
1633	surf%i(num_h) = i
1634	surf%j(num_h) = j
1635	surf%k(num_h) = k
1636	!
1637	!-- Surface orientation, bit 0 is set to 1 for upward-facing surfaces,
1638	!-- bit 1 is for downward-facing surfaces.
1639	IF ( upward_facing ) surf%facing(num_h) = IBSET( surf%facing(num_h), 0 )
1640	IF ( downward_facing ) surf%facing(num_h) = IBSET( surf%facing(num_h), 1 )
1641	!
1642	!-- Initialize surface-layer height
1643	IF ( upward_facing ) THEN
1644	surf%z_mo(num_h) = zu(k) - zw(k-1)
1645	ELSE
1646	surf%z_mo(num_h) = zw(k) - zu(k)
1647	ENDIF
1648
1649	surf%z0(num_h) = roughness_length
1650	surf%z0h(num_h) = z0h_factor * roughness_length
1651	surf%z0q(num_h) = z0h_factor * roughness_length
1652	!
1653	!-- Initialization in case of 1D pre-cursor run
1654	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )&
1655	THEN
1656	IF ( .NOT. constant_diffusion ) THEN
1657	IF ( constant_flux_layer ) THEN
1658	surf%ol(num_h) = surf%z_mo(num_h) / &
1659	( rif1d(nzb+1) + 1.0E-20_wp )
1660	surf%us(num_h) = us1d
1661	surf%usws(num_h) = usws1d
1662	surf%vsws(num_h) = vsws1d
1663	ELSE
1664	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
1665	surf%us(num_h) = 0.0_wp
1666	surf%usws(num_h) = 0.0_wp
1667	surf%vsws(num_h) = 0.0_wp
1668	ENDIF
1669	ELSE
1670	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
1671	surf%us(num_h) = 0.0_wp
1672	surf%usws(num_h) = 0.0_wp
1673	surf%vsws(num_h) = 0.0_wp
1674	ENDIF
1675	!
1676	!-- Initialization in case of constant profiles
1677	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0&
1678	.OR. INDEX(initializing_actions, 'inifor' ) /= 0 ) THEN
1679
1680	surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
1681	!
1682	!-- Very small number is required for calculation of Obukhov length
1683	!-- at first timestep
1684	surf%us(num_h) = 1E-30_wp
1685	surf%usws(num_h) = 0.0_wp
1686	surf%vsws(num_h) = 0.0_wp
1687
1688	ENDIF
1689
1690	surf%rib(num_h) = 0.0_wp
1691	surf%uvw_abs(num_h) = 0.0_wp
1692
1693	IF ( .NOT. constant_diffusion ) THEN
1694	surf%u_0(num_h) = 0.0_wp
1695	surf%v_0(num_h) = 0.0_wp
1696	ENDIF
1697
1698	surf%ts(num_h) = 0.0_wp
1699
1700	IF ( humidity ) THEN
1701	surf%qs(num_h) = 0.0_wp
1702	IF ( cloud_physics .AND. microphysics_morrison) THEN
1703	surf%qcs(num_h) = 0.0_wp
1704	surf%ncs(num_h) = 0.0_wp
1705
1706	surf%qcsws(num_h) = 0.0_wp
1707	surf%ncsws(num_h) = 0.0_wp
1708
1709	ENDIF
1710	IF ( cloud_physics .AND. microphysics_seifert) THEN
1711	surf%qrs(num_h) = 0.0_wp
1712	surf%nrs(num_h) = 0.0_wp
1713
1714	surf%qrsws(num_h) = 0.0_wp
1715	surf%nrsws(num_h) = 0.0_wp
1716
1717	surf%pt1(num_h) = 0.0_wp
1718	surf%qv1(num_h) = 0.0_wp
1719
1720	ENDIF
1721	ENDIF
1722
1723	IF ( passive_scalar ) surf%ss(num_h) = 0.0_wp
1724
1725	#if defined( __chem )
1726	DO lsp = 1, nvar
1727	IF ( air_chemistry ) surf%css(lsp,num_h) = 0.0_wp
1728	ENDDO
1729	#endif
1730	!
1731	!-- Set initial value for surface temperature
1732	surf%pt_surface(num_h) = pt_surface
1733	!
1734	!-- Inititalize surface fluxes of sensible and latent heat, as well as
1735	!-- passive scalar
1736	IF ( use_surface_fluxes ) THEN
1737
1738	IF ( upward_facing ) THEN
1739	IF ( constant_heatflux ) THEN
1740	!
1741	!-- Initialize surface heatflux. However, skip this for now if
1742	!-- if random_heatflux is set. This case, shf is initialized later.
1743	IF ( .NOT. random_heatflux ) THEN
1744	surf%shf(num_h) = surface_heatflux * &
1745	heatflux_input_conversion(nzb)
1746	!
1747	!-- Check if surface heat flux might be replaced by
1748	!-- prescribed wall heatflux
1749	IF ( k-1 /= 0 ) THEN
1750	surf%shf(num_h) = wall_heatflux(0) * &
1751	heatflux_input_conversion(k-1)
1752	ENDIF
1753	!
1754	!-- Initialize shf with data from external file LSF_DATA. Will
1755	!-- be done directly in ls_foring_surf
1756	!-- Attention: Just a workaround, need to be revised!!!
1757	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1758	! CALL ls_forcing_surf ( simulated_time )
1759	! surf%shf(num_h) = shf(j,i)
1760	ENDIF
1761	ENDIF
1762	ELSE
1763	surf%shf(num_h) = 0.0_wp
1764	ENDIF
1765	!
1766	!-- Set heat-flux at downward-facing surfaces
1767	ELSE
1768	surf%shf(num_h) = wall_heatflux(5) * &
1769	heatflux_input_conversion(k)
1770	ENDIF
1771
1772	IF ( humidity ) THEN
1773	IF ( upward_facing ) THEN
1774	IF ( constant_waterflux ) THEN
1775	surf%qsws(num_h) = surface_waterflux * &
1776	waterflux_input_conversion(nzb)
1777	IF ( k-1 /= 0 ) THEN
1778	surf%qsws(num_h) = wall_humidityflux(0) * &
1779	waterflux_input_conversion(k-1)
1780	ENDIF
1781	ELSE
1782	surf%qsws(num_h) = 0.0_wp
1783	ENDIF
1784	ELSE
1785	surf%qsws(num_h) = wall_humidityflux(5) * &
1786	heatflux_input_conversion(k)
1787	ENDIF
1788	ENDIF
1789
1790	IF ( passive_scalar ) THEN
1791	IF ( upward_facing ) THEN
1792	IF ( constant_scalarflux ) THEN
1793	surf%ssws(num_h) = surface_scalarflux
1794
1795	IF ( k-1 /= 0 ) &
1796	surf%ssws(num_h) = wall_scalarflux(0)
1797
1798	ELSE
1799	surf%ssws(num_h) = 0.0_wp
1800	ENDIF
1801	ELSE
1802	surf%ssws(num_h) = wall_scalarflux(5)
1803	ENDIF
1804	ENDIF
1805
1806	#if defined( __chem )
1807	IF ( air_chemistry ) THEN
1808	lsp_pr = 1
1809	DO WHILE ( TRIM( surface_csflux_name( lsp_pr ) ) /= 'novalue' ) !<'novalue' is the default
1810	DO lsp = 1, nvar
1811	!
1812	!-- Assign surface flux for each variable species
1813	IF ( TRIM( spc_names(lsp) ) == TRIM( surface_csflux_name(lsp_pr) ) ) THEN
1814	IF ( upward_facing ) THEN
1815	IF ( constant_csflux(lsp_pr) ) THEN
1816	surf%cssws(lsp,num_h) = surface_csflux(lsp_pr)
1817
1818	IF ( k-1 /= 0 ) &
1819	surf%cssws(lsp,num_h) = wall_csflux(lsp,0)
1820	ELSE
1821	surf%cssws(lsp,num_h) = 0.0_wp
1822	ENDIF
1823	ELSE
1824	surf%cssws(lsp,num_h) = wall_csflux(lsp,5)
1825	ENDIF
1826	ENDIF
1827	ENDDO
1828	lsp_pr = lsp_pr + 1
1829	ENDDO
1830	ENDIF
1831	#endif
1832
1833	IF ( ocean ) THEN
1834	IF ( upward_facing ) THEN
1835	surf%sasws(num_h) = bottom_salinityflux
1836	ELSE
1837	surf%sasws(num_h) = 0.0_wp
1838	ENDIF
1839	ENDIF
1840	ENDIF
1841	!
1842	!-- Increment surface indices
1843	num_h = num_h + 1
1844	num_h_kji = num_h_kji + 1
1845
1846
1847	END SUBROUTINE initialize_horizontal_surfaces
1848
1849
1850	!------------------------------------------------------------------------------!
1851	! Description:
1852	! ------------
1853	!> Initialize model-top fluxes. Currently, only the heatflux and salinity flux
1854	!> can be prescribed, latent flux is zero in this case!
1855	!------------------------------------------------------------------------------!
1856	SUBROUTINE initialize_top( k, j, i, surf, num_h, num_h_kji )
1857
1858	IMPLICIT NONE
1859
1860	INTEGER(iwp) :: i !< running index x-direction
1861	INTEGER(iwp) :: j !< running index y-direction
1862	INTEGER(iwp) :: k !< running index z-direction
1863	INTEGER(iwp) :: num_h !< current number of surface element
1864	INTEGER(iwp) :: num_h_kji !< dummy increment
1865	INTEGER(iwp) :: lsp !< running index for chemical species
1866
1867	TYPE( surf_type ) :: surf !< respective surface type
1868	!
1869	!-- Store indices of respective surface element
1870	surf%i(num_h) = i
1871	surf%j(num_h) = j
1872	surf%k(num_h) = k
1873	!
1874	!-- Initialize top heat flux
1875	IF ( constant_top_heatflux ) &
1876	surf%shf(num_h) = top_heatflux * heatflux_input_conversion(nzt+1)
1877	!
1878	!-- Initialization in case of a coupled model run
1879	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
1880	surf%shf(num_h) = 0.0_wp
1881	surf%qsws(num_h) = 0.0_wp
1882	ENDIF
1883	!
1884	!-- Prescribe latent heat flux at the top
1885	IF ( humidity ) THEN
1886	surf%qsws(num_h) = 0.0_wp
1887	IF ( cloud_physics .AND. microphysics_morrison ) THEN
1888	surf%ncsws(num_h) = 0.0_wp
1889	surf%qcsws(num_h) = 0.0_wp
1890	ENDIF
1891	IF ( cloud_physics .AND. microphysics_seifert ) THEN
1892	surf%nrsws(num_h) = 0.0_wp
1893	surf%qrsws(num_h) = 0.0_wp
1894	ENDIF
1895	ENDIF
1896	!
1897	!-- Prescribe top scalar flux
1898	IF ( passive_scalar .AND. constant_top_scalarflux ) &
1899	surf%ssws(num_h) = top_scalarflux
1900	!
1901	!-- Prescribe top chemical species' flux
1902	#if defined( __chem )
1903	DO lsp = 1, nvar
1904	IF ( air_chemistry .AND. constant_top_csflux(lsp) ) THEN
1905	surf%cssws(lsp,num_h) = top_csflux(lsp)
1906	ENDIF
1907	ENDDO
1908	#endif
1909	!
1910	!-- Prescribe top salinity flux
1911	IF ( ocean .AND. constant_top_salinityflux) &
1912	surf%sasws(num_h) = top_salinityflux
1913	!
1914	!-- Top momentum fluxes
1915	IF ( constant_top_momentumflux ) THEN
1916	surf%usws(num_h) = top_momentumflux_u * &
1917	momentumflux_input_conversion(nzt+1)
1918	surf%vsws(num_h) = top_momentumflux_v * &
1919	momentumflux_input_conversion(nzt+1)
1920	ENDIF
1921	!
1922	!-- Increment surface indices
1923	num_h = num_h + 1
1924	num_h_kji = num_h_kji + 1
1925
1926
1927	END SUBROUTINE initialize_top
1928
1929
1930	!------------------------------------------------------------------------------!
1931	! Description:
1932	! ------------
1933	!> Initialize vertical surface elements.
1934	!------------------------------------------------------------------------------!
1935	SUBROUTINE initialize_vertical_surfaces( l, k, j, i, surf, num_v, &
1936	num_v_kji, east_facing, &
1937	west_facing, south_facing, &
1938	north_facing )
1939
1940	IMPLICIT NONE
1941
1942	INTEGER(iwp) :: component !< index of wall_fluxes_ array for respective orientation
1943	INTEGER(iwp) :: i !< running index x-direction
1944	INTEGER(iwp) :: j !< running index x-direction
1945	INTEGER(iwp) :: k !< running index x-direction
1946	INTEGER(iwp) :: l !< index variable for the surface type, indicating the facing
1947	INTEGER(iwp) :: num_v !< current number of surface element
1948	INTEGER(iwp) :: num_v_kji !< current number of surface element at (j,i)
1949	INTEGER(iwp) :: lsp !< running index for chemical species
1950
1951
1952	LOGICAL :: east_facing !< flag indicating east-facing surfaces
1953	LOGICAL :: north_facing !< flag indicating north-facing surfaces
1954	LOGICAL :: south_facing !< flag indicating south-facing surfaces
1955	LOGICAL :: west_facing !< flag indicating west-facing surfaces
1956
1957	TYPE( surf_type ) :: surf !< respective surface type
1958
1959	!
1960	!-- Store indices of respective wall element
1961	surf%i(num_v) = i
1962	surf%j(num_v) = j
1963	surf%k(num_v) = k
1964	!
1965	!-- Initialize surface-layer height, or more precisely, distance to surface
1966	IF ( north_facing .OR. south_facing ) THEN
1967	surf%z_mo(num_v) = 0.5_wp * dy
1968	ELSE
1969	surf%z_mo(num_v) = 0.5_wp * dx
1970	ENDIF
1971
1972	surf%facing(num_v) = 0
1973	!
1974	!-- Surface orientation. Moreover, set component id to map wall_heatflux,
1975	!-- etc., on surface type (further below)
1976	IF ( north_facing ) THEN
1977	surf%facing(num_v) = 5 !IBSET( surf%facing(num_v), 0 )
1978	component = 4
1979	ENDIF
1980
1981	IF ( south_facing ) THEN
1982	surf%facing(num_v) = 6 !IBSET( surf%facing(num_v), 1 )
1983	component = 3
1984	ENDIF
1985
1986	IF ( east_facing ) THEN
1987	surf%facing(num_v) = 7 !IBSET( surf%facing(num_v), 2 )
1988	component = 2
1989	ENDIF
1990
1991	IF ( west_facing ) THEN
1992	surf%facing(num_v) = 8 !IBSET( surf%facing(num_v), 3 )
1993	component = 1
1994	ENDIF
1995
1996
1997	surf%z0(num_v) = roughness_length
1998	surf%z0h(num_v) = z0h_factor * roughness_length
1999	surf%z0q(num_v) = z0h_factor * roughness_length
2000
2001	surf%us(num_v) = 0.0_wp
2002	!
2003	!-- If required, initialize Obukhov length
2004	IF ( ALLOCATED( surf%ol ) ) &
2005	surf%ol(num_v) = surf%z_mo(num_v) / zeta_min
2006
2007	surf%uvw_abs(num_v) = 0.0_wp
2008
2009	surf%mom_flux_uv(num_v) = 0.0_wp
2010	surf%mom_flux_w(num_v) = 0.0_wp
2011	surf%mom_flux_tke(0:1,num_v) = 0.0_wp
2012
2013	surf%ts(num_v) = 0.0_wp
2014	surf%shf(num_v) = wall_heatflux(component)
2015	!
2016	!-- Set initial value for surface temperature
2017	surf%pt_surface(num_v) = pt_surface
2018
2019	IF ( humidity ) THEN
2020	surf%qs(num_v) = 0.0_wp
2021	surf%qsws(num_v) = wall_humidityflux(component)
2022	!
2023	!-- Following wall fluxes are assumed to be zero
2024	IF ( cloud_physics .AND. microphysics_morrison) THEN
2025	surf%qcs(num_v) = 0.0_wp
2026	surf%ncs(num_v) = 0.0_wp
2027
2028	surf%qcsws(num_v) = 0.0_wp
2029	surf%ncsws(num_v) = 0.0_wp
2030	ENDIF
2031	IF ( cloud_physics .AND. microphysics_seifert) THEN
2032	surf%qrs(num_v) = 0.0_wp
2033	surf%nrs(num_v) = 0.0_wp
2034
2035	surf%qrsws(num_v) = 0.0_wp
2036	surf%nrsws(num_v) = 0.0_wp
2037	ENDIF
2038	ENDIF
2039
2040	IF ( passive_scalar ) THEN
2041	surf%ss(num_v) = 0.0_wp
2042	surf%ssws(num_v) = wall_scalarflux(component)
2043	ENDIF
2044
2045	#if defined( __chem )
2046	IF ( air_chemistry ) THEN
2047	DO lsp = 1, nvar
2048	surf%css(lsp,num_v) = 0.0_wp
2049	surf%cssws(lsp,num_v) = wall_csflux(lsp,component)
2050	ENDDO
2051	ENDIF
2052	#endif
2053
2054	!
2055	!-- So far, salinityflux at vertical surfaces is simply zero
2056	!-- at the moment
2057	IF ( ocean ) surf%sasws(num_v) = wall_salinityflux(component)
2058	!
2059	!-- Increment wall indices
2060	num_v = num_v + 1
2061	num_v_kji = num_v_kji + 1
2062
2063	END SUBROUTINE initialize_vertical_surfaces
2064
2065	END SUBROUTINE init_surfaces
2066
2067
2068	!------------------------------------------------------------------------------!
2069	! Description:
2070	! ------------
2071	!> Determines topography-top index at given (j,i)-position.
2072	!------------------------------------------------------------------------------!
2073	FUNCTION get_topography_top_index_ji( j, i, grid )
2074
2075	IMPLICIT NONE
2076
2077	CHARACTER(LEN=*) :: grid !< flag to distinquish between staggered grids
2078	INTEGER(iwp) :: i !< grid index in x-dimension
2079	INTEGER(iwp) :: ibit !< bit position where topography information is stored on respective grid
2080	INTEGER(iwp) :: j !< grid index in y-dimension
2081	INTEGER(iwp) :: get_topography_top_index_ji !< topography top index
2082
2083	SELECT CASE ( TRIM( grid ) )
2084
2085	CASE ( 's' )
2086	ibit = 12
2087	CASE ( 'u' )
2088	ibit = 14
2089	CASE ( 'v' )
2090	ibit = 16
2091	CASE ( 'w' )
2092	ibit = 18
2093	CASE ( 's_out' )
2094	ibit = 24
2095	CASE ( 'u_out' )
2096	ibit = 26
2097	CASE ( 'v_out' )
2098	ibit = 27
2099	CASE ( 'w_out' )
2100	ibit = 28
2101	CASE DEFAULT
2102	!
2103	!-- Set default to scalar grid
2104	ibit = 12
2105
2106	END SELECT
2107
2108	get_topography_top_index_ji = MAXLOC( &
2109	MERGE( 1, 0, &
2110	BTEST( wall_flags_0(:,j,i), ibit ) &
2111	), DIM = 1 &
2112	) - 1
2113
2114	RETURN
2115
2116	END FUNCTION get_topography_top_index_ji
2117
2118	!------------------------------------------------------------------------------!
2119	! Description:
2120	! ------------
2121	!> Determines topography-top index at each (j,i)-position.
2122	!------------------------------------------------------------------------------!
2123	FUNCTION get_topography_top_index( grid )
2124
2125	IMPLICIT NONE
2126
2127	CHARACTER(LEN=*) :: grid !< flag to distinquish between staggered grids
2128	INTEGER(iwp) :: ibit !< bit position where topography information is stored on respective grid
2129	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: get_topography_top_index !< topography top index
2130
2131	SELECT CASE ( TRIM( grid ) )
2132
2133	CASE ( 's' )
2134	ibit = 12
2135	CASE ( 'u' )
2136	ibit = 14
2137	CASE ( 'v' )
2138	ibit = 16
2139	CASE ( 'w' )
2140	ibit = 18
2141	CASE ( 's_out' )
2142	ibit = 24
2143	CASE ( 'u_out' )
2144	ibit = 26
2145	CASE ( 'v_out' )
2146	ibit = 27
2147	CASE ( 'w_out' )
2148	ibit = 28
2149	CASE DEFAULT
2150	!
2151	!-- Set default to scalar grid
2152	ibit = 12
2153
2154	END SELECT
2155
2156	get_topography_top_index(nys:nyn,nxl:nxr) = MAXLOC( &
2157	MERGE( 1, 0, &
2158	BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), ibit )&
2159	), DIM = 1 &
2160	) - 1
2161
2162	RETURN
2163
2164	END FUNCTION get_topography_top_index
2165
2166	!------------------------------------------------------------------------------!
2167	! Description:
2168	! ------------
2169	!> Gathers all surface elements with the same facing (but possibly different
2170	!> type) onto a surface type, and writes binary data into restart files.
2171	!------------------------------------------------------------------------------!
2172	SUBROUTINE surface_write_restart_data
2173
2174	IMPLICIT NONE
2175
2176	CHARACTER(LEN=1) :: dum !< dummy string to create output-variable name
2177
2178	INTEGER(iwp) :: i !< running index x-direction
2179	INTEGER(iwp) :: j !< running index y-direction
2180	INTEGER(iwp) :: l !< index surface type orientation
2181	INTEGER(iwp) :: lsp !< running index chemical species
2182	INTEGER(iwp) :: m !< running index for surface elements on individual surface array
2183	INTEGER(iwp), DIMENSION(0:3) :: mm !< running index for surface elements on gathered surface array
2184
2185	TYPE(surf_type), DIMENSION(0:2) :: surf_h !< gathered horizontal surfaces, contains all surface types
2186	TYPE(surf_type), DIMENSION(0:3) :: surf_v !< gathered vertical surfaces, contains all surface types
2187
2188	!
2189	!-- Determine total number of horizontal and vertical surface elements before
2190	!-- writing var_list
2191	CALL surface_last_actions
2192	!
2193	!-- Count number of grid points with same facing and allocate attributes respectively
2194	!-- Horizontal upward facing
2195	surf_h(0)%ns = ns_h_on_file(0)
2196	CALL allocate_surface_attributes_h( surf_h(0), nys, nyn, nxl, nxr )
2197	!
2198	!-- Horizontal downward facing
2199	surf_h(1)%ns = ns_h_on_file(1)
2200	CALL allocate_surface_attributes_h( surf_h(1), nys, nyn, nxl, nxr )
2201	!
2202	!-- Model top
2203	surf_h(2)%ns = ns_h_on_file(2)
2204	CALL allocate_surface_attributes_h_top( surf_h(2), nys, nyn, nxl, nxr )
2205	!
2206	!-- Vertical surfaces
2207	DO l = 0, 3
2208	surf_v(l)%ns = ns_v_on_file(l)
2209	CALL allocate_surface_attributes_v( surf_v(l), &
2210	nys, nyn, nxl, nxr )
2211	ENDDO
2212	!
2213	!-- In the following, gather data from surfaces elements with the same
2214	!-- facing (but possibly differt type) on 1 data-type array.
2215	mm(0:2) = 1
2216	DO l = 0, 2
2217	DO i = nxl, nxr
2218	DO j = nys, nyn
2219	DO m = surf_def_h(l)%start_index(j,i), &
2220	surf_def_h(l)%end_index(j,i)
2221	IF ( ALLOCATED( surf_def_h(l)%us ) ) &
2222	surf_h(l)%us(mm(l)) = surf_def_h(l)%us(m)
2223	IF ( ALLOCATED( surf_def_h(l)%ts ) ) &
2224	surf_h(l)%ts(mm(l)) = surf_def_h(l)%ts(m)
2225	IF ( ALLOCATED( surf_def_h(l)%qs ) ) &
2226	surf_h(l)%qs(mm(l)) = surf_def_h(l)%qs(m)
2227	IF ( ALLOCATED( surf_def_h(l)%ss ) ) &
2228	surf_h(l)%ss(mm(l)) = surf_def_h(l)%ss(m)
2229	IF ( ALLOCATED( surf_def_h(l)%qcs ) ) &
2230	surf_h(l)%qcs(mm(l)) = surf_def_h(l)%qcs(m)
2231	IF ( ALLOCATED( surf_def_h(l)%ncs ) ) &
2232	surf_h(l)%ncs(mm(l)) = surf_def_h(l)%ncs(m)
2233	IF ( ALLOCATED( surf_def_h(l)%qrs ) ) &
2234	surf_h(l)%qrs(mm(l)) = surf_def_h(l)%qrs(m)
2235	IF ( ALLOCATED( surf_def_h(l)%nrs ) ) &
2236	surf_h(l)%nrs(mm(l)) = surf_def_h(l)%nrs(m)
2237	IF ( ALLOCATED( surf_def_h(l)%ol ) ) &
2238	surf_h(l)%ol(mm(l)) = surf_def_h(l)%ol(m)
2239	IF ( ALLOCATED( surf_def_h(l)%rib ) ) &
2240	surf_h(l)%rib(mm(l)) = surf_def_h(l)%rib(m)
2241	IF ( ALLOCATED( surf_def_h(l)%usws ) ) &
2242	surf_h(l)%usws(mm(l)) = surf_def_h(l)%usws(m)
2243	IF ( ALLOCATED( surf_def_h(l)%vsws ) ) &
2244	surf_h(l)%vsws(mm(l)) = surf_def_h(l)%vsws(m)
2245	IF ( ALLOCATED( surf_def_h(l)%shf ) ) &
2246	surf_h(l)%shf(mm(l)) = surf_def_h(l)%shf(m)
2247	IF ( ALLOCATED( surf_def_h(l)%qsws ) ) &
2248	surf_h(l)%qsws(mm(l)) = surf_def_h(l)%qsws(m)
2249	IF ( ALLOCATED( surf_def_h(l)%ssws ) ) &
2250	surf_h(l)%ssws(mm(l)) = surf_def_h(l)%ssws(m)
2251	#if defined( __chem )
2252	IF ( ALLOCATED( surf_def_h(l)%css ) ) THEN
2253	DO lsp = 1,nvar
2254	surf_h(l)%css(lsp,mm(l)) = surf_def_h(l)%css(lsp,m)
2255	ENDDO
2256	ENDIF
2257	IF ( ALLOCATED( surf_def_h(l)%cssws ) ) THEN
2258	DO lsp = 1,nvar
2259	surf_h(l)%cssws(lsp,mm(l)) = surf_def_h(l)%cssws(lsp,m)
2260	ENDDO
2261	ENDIF
2262	#endif
2263	IF ( ALLOCATED( surf_def_h(l)%ncsws ) ) &
2264	surf_h(l)%ncsws(mm(l)) = surf_def_h(l)%ncsws(m)
2265	IF ( ALLOCATED( surf_def_h(l)%nrsws ) ) &
2266	surf_h(l)%nrsws(mm(l)) = surf_def_h(l)%nrsws(m)
2267	IF ( ALLOCATED( surf_def_h(l)%sasws ) ) &
2268	surf_h(l)%sasws(mm(l)) = surf_def_h(l)%sasws(m)
2269
2270	mm(l) = mm(l) + 1
2271	ENDDO
2272
2273	IF ( l == 0 ) THEN
2274	DO m = surf_lsm_h%start_index(j,i), &
2275	surf_lsm_h%end_index(j,i)
2276	IF ( ALLOCATED( surf_lsm_h%us ) ) &
2277	surf_h(0)%us(mm(0)) = surf_lsm_h%us(m)
2278	IF ( ALLOCATED( surf_lsm_h%ts ) ) &
2279	surf_h(0)%ts(mm(0)) = surf_lsm_h%ts(m)
2280	IF ( ALLOCATED( surf_lsm_h%qs ) ) &
2281	surf_h(0)%qs(mm(0)) = surf_lsm_h%qs(m)
2282	IF ( ALLOCATED( surf_lsm_h%ss ) ) &
2283	surf_h(0)%ss(mm(0)) = surf_lsm_h%ss(m)
2284	IF ( ALLOCATED( surf_lsm_h%qcs ) ) &
2285	surf_h(0)%qcs(mm(0)) = surf_lsm_h%qcs(m)
2286	IF ( ALLOCATED( surf_lsm_h%ncs ) ) &
2287	surf_h(0)%ncs(mm(0)) = surf_lsm_h%ncs(m)
2288	IF ( ALLOCATED( surf_lsm_h%qrs ) ) &
2289	surf_h(0)%qrs(mm(0)) = surf_lsm_h%qrs(m)
2290	IF ( ALLOCATED( surf_lsm_h%nrs ) ) &
2291	surf_h(0)%nrs(mm(0)) = surf_lsm_h%nrs(m)
2292	IF ( ALLOCATED( surf_lsm_h%ol ) ) &
2293	surf_h(0)%ol(mm(0)) = surf_lsm_h%ol(m)
2294	IF ( ALLOCATED( surf_lsm_h%rib ) ) &
2295	surf_h(0)%rib(mm(0)) = surf_lsm_h%rib(m)
2296	IF ( ALLOCATED( surf_lsm_h%usws ) ) &
2297	surf_h(0)%usws(mm(0)) = surf_lsm_h%usws(m)
2298	IF ( ALLOCATED( surf_lsm_h%vsws ) ) &
2299	surf_h(0)%vsws(mm(0)) = surf_lsm_h%vsws(m)
2300	IF ( ALLOCATED( surf_lsm_h%shf ) ) &
2301	surf_h(0)%shf(mm(0)) = surf_lsm_h%shf(m)
2302	IF ( ALLOCATED( surf_lsm_h%qsws ) ) &
2303	surf_h(0)%qsws(mm(0)) = surf_lsm_h%qsws(m)
2304	IF ( ALLOCATED( surf_lsm_h%ssws ) ) &
2305	surf_h(0)%ssws(mm(0)) = surf_lsm_h%ssws(m)
2306	#if defined( __chem )
2307	IF ( ALLOCATED( surf_lsm_h%css ) ) THEN
2308	DO lsp = 1, nvar
2309	surf_h(0)%css(lsp,mm(0)) = surf_lsm_h%css(lsp,m)
2310	ENDDO
2311	ENDIF
2312	IF ( ALLOCATED( surf_lsm_h%cssws ) ) THEN
2313	DO lsp = 1, nvar
2314	surf_h(0)%cssws(lsp,mm(0)) = surf_lsm_h%cssws(lsp,m)
2315	ENDDO
2316	ENDIF
2317	#endif
2318	IF ( ALLOCATED( surf_lsm_h%ncsws ) ) &
2319	surf_h(0)%ncsws(mm(0)) = surf_lsm_h%ncsws(m)
2320	IF ( ALLOCATED( surf_lsm_h%nrsws ) ) &
2321	surf_h(0)%nrsws(mm(0)) = surf_lsm_h%nrsws(m)
2322	IF ( ALLOCATED( surf_lsm_h%sasws ) ) &
2323	surf_h(0)%sasws(mm(0)) = surf_lsm_h%sasws(m)
2324
2325	mm(0) = mm(0) + 1
2326
2327	ENDDO
2328
2329	DO m = surf_usm_h%start_index(j,i), &
2330	surf_usm_h%end_index(j,i)
2331	IF ( ALLOCATED( surf_usm_h%us ) ) &
2332	surf_h(0)%us(mm(0)) = surf_usm_h%us(m)
2333	IF ( ALLOCATED( surf_usm_h%ts ) ) &
2334	surf_h(0)%ts(mm(0)) = surf_usm_h%ts(m)
2335	IF ( ALLOCATED( surf_usm_h%qs ) ) &
2336	surf_h(0)%qs(mm(0)) = surf_usm_h%qs(m)
2337	IF ( ALLOCATED( surf_usm_h%ss ) ) &
2338	surf_h(0)%ss(mm(0)) = surf_usm_h%ss(m)
2339	IF ( ALLOCATED( surf_usm_h%qcs ) ) &
2340	surf_h(0)%qcs(mm(0)) = surf_usm_h%qcs(m)
2341	IF ( ALLOCATED( surf_usm_h%ncs ) ) &
2342	surf_h(0)%ncs(mm(0)) = surf_usm_h%ncs(m)
2343	IF ( ALLOCATED( surf_usm_h%qrs ) ) &
2344	surf_h(0)%qrs(mm(0)) = surf_usm_h%qrs(m)
2345	IF ( ALLOCATED( surf_usm_h%nrs ) ) &
2346	surf_h(0)%nrs(mm(0)) = surf_usm_h%nrs(m)
2347	IF ( ALLOCATED( surf_usm_h%ol ) ) &
2348	surf_h(0)%ol(mm(0)) = surf_usm_h%ol(m)
2349	IF ( ALLOCATED( surf_usm_h%rib ) ) &
2350	surf_h(0)%rib(mm(0)) = surf_usm_h%rib(m)
2351	IF ( ALLOCATED( surf_usm_h%usws ) ) &
2352	surf_h(0)%usws(mm(0)) = surf_usm_h%usws(m)
2353	IF ( ALLOCATED( surf_usm_h%vsws ) ) &
2354	surf_h(0)%vsws(mm(0)) = surf_usm_h%vsws(m)
2355	IF ( ALLOCATED( surf_usm_h%shf ) ) &
2356	surf_h(0)%shf(mm(0)) = surf_usm_h%shf(m)
2357	IF ( ALLOCATED( surf_usm_h%qsws ) ) &
2358	surf_h(0)%qsws(mm(0)) = surf_usm_h%qsws(m)
2359	IF ( ALLOCATED( surf_usm_h%ssws ) ) &
2360	surf_h(0)%ssws(mm(0)) = surf_usm_h%ssws(m)
2361	#if defined( __chem )
2362	IF ( ALLOCATED( surf_usm_h%css ) ) THEN
2363	DO lsp = 1, nvar
2364	surf_h(0)%css(lsp,mm(0)) = surf_usm_h%css(lsp,m)
2365	ENDDO
2366	ENDIF
2367	IF ( ALLOCATED( surf_usm_h%cssws ) ) THEN
2368	DO lsp = 1, nvar
2369	surf_h(0)%cssws(lsp,mm(0)) = surf_usm_h%cssws(lsp,m)
2370	ENDDO
2371	ENDIF
2372	#endif
2373	IF ( ALLOCATED( surf_usm_h%ncsws ) ) &
2374	surf_h(0)%ncsws(mm(0)) = surf_usm_h%ncsws(m)
2375	IF ( ALLOCATED( surf_usm_h%nrsws ) ) &
2376	surf_h(0)%nrsws(mm(0)) = surf_usm_h%nrsws(m)
2377	IF ( ALLOCATED( surf_usm_h%sasws ) ) &
2378	surf_h(0)%sasws(mm(0)) = surf_usm_h%sasws(m)
2379
2380	mm(0) = mm(0) + 1
2381
2382	ENDDO
2383
2384
2385	ENDIF
2386
2387	ENDDO
2388
2389	ENDDO
2390	!
2391	!-- Gather start- and end indices
2392	DO i = nxl, nxr
2393	DO j = nys, nyn
2394	IF ( surf_def_h(l)%start_index(j,i) <= &
2395	surf_def_h(l)%end_index(j,i) ) THEN
2396	surf_h(l)%start_index(j,i) = surf_def_h(l)%start_index(j,i)
2397	surf_h(l)%end_index(j,i) = surf_def_h(l)%end_index(j,i)
2398	ENDIF
2399	IF ( l == 0 ) THEN
2400	IF ( surf_lsm_h%start_index(j,i) <= &
2401	surf_lsm_h%end_index(j,i) ) THEN
2402	surf_h(l)%start_index(j,i) = surf_lsm_h%start_index(j,i)
2403	surf_h(l)%end_index(j,i) = surf_lsm_h%end_index(j,i)
2404	ENDIF
2405	IF ( surf_usm_h%start_index(j,i) <= &
2406	surf_usm_h%end_index(j,i) ) THEN
2407	surf_h(l)%start_index(j,i) = surf_usm_h%start_index(j,i)
2408	surf_h(l)%end_index(j,i) = surf_usm_h%end_index(j,i)
2409	ENDIF
2410	ENDIF
2411	ENDDO
2412	ENDDO
2413	ENDDO
2414
2415
2416	mm(0:3) = 1
2417	DO l = 0, 3
2418	DO i = nxl, nxr
2419	DO j = nys, nyn
2420	DO m = surf_def_v(l)%start_index(j,i), &
2421	surf_def_v(l)%end_index(j,i)
2422	IF ( ALLOCATED( surf_def_v(l)%us ) ) &
2423	surf_v(l)%us(mm(l)) = surf_def_v(l)%us(m)
2424	IF ( ALLOCATED( surf_def_v(l)%ts ) ) &
2425	surf_v(l)%ts(mm(l)) = surf_def_v(l)%ts(m)
2426	IF ( ALLOCATED( surf_def_v(l)%qs ) ) &
2427	surf_v(l)%qs(mm(l)) = surf_def_v(l)%qs(m)
2428	IF ( ALLOCATED( surf_def_v(l)%ss ) ) &
2429	surf_v(l)%ss(mm(l)) = surf_def_v(l)%ss(m)
2430	IF ( ALLOCATED( surf_def_v(l)%qcs ) ) &
2431	surf_v(l)%qcs(mm(l)) = surf_def_v(l)%qcs(m)
2432	IF ( ALLOCATED( surf_def_v(l)%ncs ) ) &
2433	surf_v(l)%ncs(mm(l)) = surf_def_v(l)%ncs(m)
2434	IF ( ALLOCATED( surf_def_v(l)%qrs ) ) &
2435	surf_v(l)%qrs(mm(l)) = surf_def_v(l)%qrs(m)
2436	IF ( ALLOCATED( surf_def_v(l)%nrs ) ) &
2437	surf_v(l)%nrs(mm(l)) = surf_def_v(l)%nrs(m)
2438	IF ( ALLOCATED( surf_def_v(l)%ol ) ) &
2439	surf_v(l)%ol(mm(l)) = surf_def_v(l)%ol(m)
2440	IF ( ALLOCATED( surf_def_v(l)%rib ) ) &
2441	surf_v(l)%rib(mm(l)) = surf_def_v(l)%rib(m)
2442	IF ( ALLOCATED( surf_def_v(l)%shf ) ) &
2443	surf_v(l)%shf(mm(l)) = surf_def_v(l)%shf(m)
2444	IF ( ALLOCATED( surf_def_v(l)%qsws ) ) &
2445	surf_v(l)%qsws(mm(l)) = surf_def_v(l)%qsws(m)
2446	IF ( ALLOCATED( surf_def_v(l)%ssws ) ) &
2447	surf_v(l)%ssws(mm(l)) = surf_def_v(l)%ssws(m)
2448	#if defined( __chem )
2449	IF ( ALLOCATED( surf_def_v(l)%css ) ) THEN
2450	DO lsp = 1, nvar
2451	surf_v(l)%css(lsp,mm(l)) = surf_def_v(l)%css(lsp,m)
2452	ENDDO
2453	ENDIF
2454	IF ( ALLOCATED( surf_def_v(l)%cssws ) ) THEN
2455	DO lsp = 1, nvar
2456	surf_v(l)%cssws(lsp,mm(l)) = surf_def_v(l)%cssws(lsp,m)
2457	ENDDO
2458	ENDIF
2459	#endif
2460	IF ( ALLOCATED( surf_def_v(l)%ncsws ) ) &
2461	surf_v(l)%ncsws(mm(l)) = surf_def_v(l)%ncsws(m)
2462	IF ( ALLOCATED( surf_def_v(l)%nrsws ) ) &
2463	surf_v(l)%nrsws(mm(l)) = surf_def_v(l)%nrsws(m)
2464	IF ( ALLOCATED( surf_def_v(l)%sasws ) ) &
2465	surf_v(l)%sasws(mm(l)) = surf_def_v(l)%sasws(m)
2466	IF ( ALLOCATED( surf_def_v(l)%mom_flux_uv) ) &
2467	surf_v(l)%mom_flux_uv(mm(l)) = surf_def_v(l)%mom_flux_uv(m)
2468	IF ( ALLOCATED( surf_def_v(l)%mom_flux_w) ) &
2469	surf_v(l)%mom_flux_w(mm(l)) = surf_def_v(l)%mom_flux_w(m)
2470	IF ( ALLOCATED( surf_def_v(l)%mom_flux_tke) ) &
2471	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_def_v(l)%mom_flux_tke(0:1,m)
2472
2473	mm(l) = mm(l) + 1
2474	ENDDO
2475
2476	DO m = surf_lsm_v(l)%start_index(j,i), &
2477	surf_lsm_v(l)%end_index(j,i)
2478	IF ( ALLOCATED( surf_lsm_v(l)%us ) ) &
2479	surf_v(l)%us(mm(l)) = surf_lsm_v(l)%us(m)
2480	IF ( ALLOCATED( surf_lsm_v(l)%ts ) ) &
2481	surf_v(l)%ts(mm(l)) = surf_lsm_v(l)%ts(m)
2482	IF ( ALLOCATED( surf_lsm_v(l)%qs ) ) &
2483	surf_v(l)%qs(mm(l)) = surf_lsm_v(l)%qs(m)
2484	IF ( ALLOCATED( surf_lsm_v(l)%ss ) ) &
2485	surf_v(l)%ss(mm(l)) = surf_lsm_v(l)%ss(m)
2486	IF ( ALLOCATED( surf_lsm_v(l)%qcs ) ) &
2487	surf_v(l)%qcs(mm(l)) = surf_lsm_v(l)%qcs(m)
2488	IF ( ALLOCATED( surf_lsm_v(l)%ncs ) ) &
2489	surf_v(l)%ncs(mm(l)) = surf_lsm_v(l)%ncs(m)
2490	IF ( ALLOCATED( surf_lsm_v(l)%qrs ) ) &
2491	surf_v(l)%qrs(mm(l)) = surf_lsm_v(l)%qrs(m)
2492	IF ( ALLOCATED( surf_lsm_v(l)%nrs ) ) &
2493	surf_v(l)%nrs(mm(l)) = surf_lsm_v(l)%nrs(m)
2494	IF ( ALLOCATED( surf_lsm_v(l)%ol ) ) &
2495	surf_v(l)%ol(mm(l)) = surf_lsm_v(l)%ol(m)
2496	IF ( ALLOCATED( surf_lsm_v(l)%rib ) ) &
2497	surf_v(l)%rib(mm(l)) = surf_lsm_v(l)%rib(m)
2498	IF ( ALLOCATED( surf_lsm_v(l)%usws ) ) &
2499	surf_v(l)%usws(mm(l)) = surf_lsm_v(l)%usws(m)
2500	IF ( ALLOCATED( surf_lsm_v(l)%vsws ) ) &
2501	surf_v(l)%vsws(mm(l)) = surf_lsm_v(l)%vsws(m)
2502	IF ( ALLOCATED( surf_lsm_v(l)%shf ) ) &
2503	surf_v(l)%shf(mm(l)) = surf_lsm_v(l)%shf(m)
2504	IF ( ALLOCATED( surf_lsm_v(l)%qsws ) ) &
2505	surf_v(l)%qsws(mm(l)) = surf_lsm_v(l)%qsws(m)
2506	IF ( ALLOCATED( surf_lsm_v(l)%ssws ) ) &
2507	surf_v(l)%ssws(mm(l)) = surf_lsm_v(l)%ssws(m)
2508	#if defined( __chem )
2509	IF ( ALLOCATED( surf_lsm_v(l)%css ) ) THEN
2510	DO lsp = 1, nvar
2511	surf_v(l)%css(lsp,mm(l)) = surf_lsm_v(l)%css(lsp,m)
2512	ENDDO
2513	ENDIF
2514	IF ( ALLOCATED( surf_lsm_v(l)%cssws ) ) THEN
2515	DO lsp = 1, nvar
2516	surf_v(l)%cssws(lsp,mm(l)) = surf_lsm_v(l)%cssws(lsp,m)
2517	ENDDO
2518	ENDIF
2519	#endif
2520	IF ( ALLOCATED( surf_lsm_v(l)%ncsws ) ) &
2521	surf_v(l)%ncsws(mm(l)) = surf_lsm_v(l)%ncsws(m)
2522	IF ( ALLOCATED( surf_lsm_v(l)%nrsws ) ) &
2523	surf_v(l)%nrsws(mm(l)) = surf_lsm_v(l)%nrsws(m)
2524	IF ( ALLOCATED( surf_lsm_v(l)%sasws ) ) &
2525	surf_v(l)%sasws(mm(l)) = surf_lsm_v(l)%sasws(m)
2526	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_uv) ) &
2527	surf_v(l)%mom_flux_uv(mm(l)) = surf_lsm_v(l)%mom_flux_uv(m)
2528	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_w) ) &
2529	surf_v(l)%mom_flux_w(mm(l)) = surf_lsm_v(l)%mom_flux_w(m)
2530	IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_tke) ) &
2531	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_lsm_v(l)%mom_flux_tke(0:1,m)
2532
2533	mm(l) = mm(l) + 1
2534	ENDDO
2535
2536	DO m = surf_usm_v(l)%start_index(j,i), &
2537	surf_usm_v(l)%end_index(j,i)
2538	IF ( ALLOCATED( surf_usm_v(l)%us ) ) &
2539	surf_v(l)%us(mm(l)) = surf_usm_v(l)%us(m)
2540	IF ( ALLOCATED( surf_usm_v(l)%ts ) ) &
2541	surf_v(l)%ts(mm(l)) = surf_usm_v(l)%ts(m)
2542	IF ( ALLOCATED( surf_usm_v(l)%qs ) ) &
2543	surf_v(l)%qs(mm(l)) = surf_usm_v(l)%qs(m)
2544	IF ( ALLOCATED( surf_usm_v(l)%ss ) ) &
2545	surf_v(l)%ss(mm(l)) = surf_usm_v(l)%ss(m)
2546	IF ( ALLOCATED( surf_usm_v(l)%qcs ) ) &
2547	surf_v(l)%qcs(mm(l)) = surf_usm_v(l)%qcs(m)
2548	IF ( ALLOCATED( surf_usm_v(l)%ncs ) ) &
2549	surf_v(l)%ncs(mm(l)) = surf_usm_v(l)%ncs(m)
2550	IF ( ALLOCATED( surf_usm_v(l)%qrs ) ) &
2551	surf_v(l)%qrs(mm(l)) = surf_usm_v(l)%qrs(m)
2552	IF ( ALLOCATED( surf_usm_v(l)%nrs ) ) &
2553	surf_v(l)%nrs(mm(l)) = surf_usm_v(l)%nrs(m)
2554	IF ( ALLOCATED( surf_usm_v(l)%ol ) ) &
2555	surf_v(l)%ol(mm(l)) = surf_usm_v(l)%ol(m)
2556	IF ( ALLOCATED( surf_usm_v(l)%rib ) ) &
2557	surf_v(l)%rib(mm(l)) = surf_usm_v(l)%rib(m)
2558	IF ( ALLOCATED( surf_usm_v(l)%usws ) ) &
2559	surf_v(l)%usws(mm(l)) = surf_usm_v(l)%usws(m)
2560	IF ( ALLOCATED( surf_usm_v(l)%vsws ) ) &
2561	surf_v(l)%vsws(mm(l)) = surf_usm_v(l)%vsws(m)
2562	IF ( ALLOCATED( surf_usm_v(l)%shf ) ) &
2563	surf_v(l)%shf(mm(l)) = surf_usm_v(l)%shf(m)
2564	IF ( ALLOCATED( surf_usm_v(l)%qsws ) ) &
2565	surf_v(l)%qsws(mm(l)) = surf_usm_v(l)%qsws(m)
2566	IF ( ALLOCATED( surf_usm_v(l)%ssws ) ) &
2567	surf_v(l)%ssws(mm(l)) = surf_usm_v(l)%ssws(m)
2568	#if defined( __chem )
2569	IF ( ALLOCATED( surf_usm_v(l)%css ) ) THEN
2570	DO lsp = 1, nvar
2571	surf_v(l)%css(lsp,mm(l)) = surf_usm_v(l)%css(lsp,m)
2572	ENDDO
2573	ENDIF
2574	IF ( ALLOCATED( surf_usm_v(l)%cssws ) ) THEN
2575	DO lsp = 1, nvar
2576	surf_v(l)%cssws(lsp,mm(l)) = surf_usm_v(l)%cssws(lsp,m)
2577	ENDDO
2578	ENDIF
2579	#endif
2580	IF ( ALLOCATED( surf_usm_v(l)%ncsws ) ) &
2581	surf_v(l)%ncsws(mm(l)) = surf_usm_v(l)%ncsws(m)
2582	IF ( ALLOCATED( surf_usm_v(l)%nrsws ) ) &
2583	surf_v(l)%nrsws(mm(l)) = surf_usm_v(l)%nrsws(m)
2584	IF ( ALLOCATED( surf_usm_v(l)%sasws ) ) &
2585	surf_v(l)%sasws(mm(l)) = surf_usm_v(l)%sasws(m)
2586	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_uv) ) &
2587	surf_v(l)%mom_flux_uv(mm(l)) = surf_usm_v(l)%mom_flux_uv(m)
2588	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_w) ) &
2589	surf_v(l)%mom_flux_w(mm(l)) = surf_usm_v(l)%mom_flux_w(m)
2590	IF ( ALLOCATED( surf_usm_v(l)%mom_flux_tke) ) &
2591	surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_usm_v(l)%mom_flux_tke(0:1,m)
2592
2593	mm(l) = mm(l) + 1
2594	ENDDO
2595
2596	ENDDO
2597	ENDDO
2598	!
2599	!-- Finally, determine start- and end-index for the respective surface
2600	surf_v(l)%start_index = MAX( surf_def_v(l)%start_index, &
2601	surf_lsm_v(l)%start_index, &
2602	surf_usm_v(l)%start_index )
2603	surf_v(l)%end_index = MAX( surf_def_v(l)%end_index, &
2604	surf_lsm_v(l)%end_index, &
2605	surf_usm_v(l)%end_index )
2606	ENDDO
2607
2608	WRITE ( 14 ) 'ns_h_on_file '
2609	WRITE ( 14 ) ns_h_on_file
2610	WRITE ( 14 ) 'ns_v_on_file '
2611	WRITE ( 14 ) ns_v_on_file
2612	!
2613	!-- Write required restart data.
2614	!-- Start with horizontal surfaces (upward-, downward-facing, and model top)
2615	DO l = 0, 2
2616	WRITE( dum, '(I1)') l
2617
2618	WRITE ( 14 ) 'surf_h(' // dum // ')%start_index '
2619	WRITE ( 14 ) surf_h(l)%start_index
2620	WRITE ( 14 ) 'surf_h(' // dum // ')%end_index '
2621	WRITE ( 14 ) surf_h(l)%end_index
2622
2623	WRITE ( 14 ) 'surf_h(' // dum // ')%us '
2624	IF ( ALLOCATED ( surf_h(l)%us ) ) THEN
2625	WRITE ( 14 ) surf_h(l)%us
2626	ENDIF
2627	WRITE ( 14 ) 'surf_h(' // dum // ')%ts '
2628	IF ( ALLOCATED ( surf_h(l)%ts ) ) THEN
2629	WRITE ( 14 ) surf_h(l)%ts
2630	ENDIF
2631	WRITE ( 14 ) 'surf_h(' // dum // ')%qs '
2632	IF ( ALLOCATED ( surf_h(l)%qs ) ) THEN
2633	WRITE ( 14 ) surf_h(l)%qs
2634	ENDIF
2635	WRITE ( 14 ) 'surf_h(' // dum // ')%ss '
2636	IF ( ALLOCATED ( surf_h(l)%ss ) ) THEN
2637	WRITE ( 14 ) surf_h(l)%ss
2638	ENDIF
2639	WRITE ( 14 ) 'surf_h(' // dum // ')%qcs '
2640	IF ( ALLOCATED ( surf_h(l)%qcs ) ) THEN
2641	WRITE ( 14 ) surf_h(l)%qcs
2642	ENDIF
2643	WRITE ( 14 ) 'surf_h(' // dum // ')%ncs '
2644	IF ( ALLOCATED ( surf_h(l)%ncs ) ) THEN
2645	WRITE ( 14 ) surf_h(l)%ncs
2646	ENDIF
2647	WRITE ( 14 ) 'surf_h(' // dum // ')%qrs '
2648	IF ( ALLOCATED ( surf_h(l)%qrs ) ) THEN
2649	WRITE ( 14 ) surf_h(l)%qrs
2650	ENDIF
2651	WRITE ( 14 ) 'surf_h(' // dum // ')%nrs '
2652	IF ( ALLOCATED ( surf_h(l)%nrs ) ) THEN
2653	WRITE ( 14 ) surf_h(l)%nrs
2654	ENDIF
2655	WRITE ( 14 ) 'surf_h(' // dum // ')%ol '
2656	IF ( ALLOCATED ( surf_h(l)%ol ) ) THEN
2657	WRITE ( 14 ) surf_h(l)%ol
2658	ENDIF
2659	WRITE ( 14 ) 'surf_h(' // dum // ')%rib '
2660	IF ( ALLOCATED ( surf_h(l)%rib ) ) THEN
2661	WRITE ( 14 ) surf_h(l)%rib
2662	ENDIF
2663	WRITE ( 14 ) 'surf_h(' // dum // ')%usws '
2664	IF ( ALLOCATED ( surf_h(l)%usws ) ) THEN
2665	WRITE ( 14 ) surf_h(l)%usws
2666	ENDIF
2667	WRITE ( 14 ) 'surf_h(' // dum // ')%vsws '
2668	IF ( ALLOCATED ( surf_h(l)%vsws ) ) THEN
2669	WRITE ( 14 ) surf_h(l)%vsws
2670	ENDIF
2671	WRITE ( 14 ) 'surf_h(' // dum // ')%shf '
2672	IF ( ALLOCATED ( surf_h(l)%shf ) ) THEN
2673	WRITE ( 14 ) surf_h(l)%shf
2674	ENDIF
2675	WRITE ( 14 ) 'surf_h(' // dum // ')%qsws '
2676	IF ( ALLOCATED ( surf_h(l)%qsws ) ) THEN
2677	WRITE ( 14 ) surf_h(l)%qsws
2678	ENDIF
2679	WRITE ( 14 ) 'surf_h(' // dum // ')%ssws '
2680	IF ( ALLOCATED ( surf_h(l)%ssws ) ) THEN
2681	WRITE ( 14 ) surf_h(l)%ssws
2682	ENDIF
2683	#if defined ( __chem )
2684	WRITE ( 14 ) 'surf_h(' // dum // ')%css '
2685	IF ( ALLOCATED ( surf_h(l)%css ) ) THEN
2686	WRITE ( 14 ) surf_h(l)%css
2687	ENDIF
2688	WRITE ( 14 ) 'surf_h(' // dum // ')%cssws '
2689	IF ( ALLOCATED ( surf_h(l)%cssws ) ) THEN
2690	WRITE ( 14 ) surf_h(l)%cssws
2691	ENDIF
2692	#endif
2693	WRITE ( 14 ) 'surf_h(' // dum // ')%qcsws '
2694	IF ( ALLOCATED ( surf_h(l)%qcsws ) ) THEN
2695	WRITE ( 14 ) surf_h(l)%qcsws
2696	ENDIF
2697	WRITE ( 14 ) 'surf_h(' // dum // ')%ncsws '
2698	IF ( ALLOCATED ( surf_h(l)%ncsws ) ) THEN
2699	WRITE ( 14 ) surf_h(l)%ncsws
2700	ENDIF
2701	WRITE ( 14 ) 'surf_h(' // dum // ')%qrsws '
2702	IF ( ALLOCATED ( surf_h(l)%qrsws ) ) THEN
2703	WRITE ( 14 ) surf_h(l)%qrsws
2704	ENDIF
2705	WRITE ( 14 ) 'surf_h(' // dum // ')%nrsws '
2706	IF ( ALLOCATED ( surf_h(l)%nrsws ) ) THEN
2707	WRITE ( 14 ) surf_h(l)%nrsws
2708	ENDIF
2709	WRITE ( 14 ) 'surf_h(' // dum // ')%sasws '
2710	IF ( ALLOCATED ( surf_h(l)%sasws ) ) THEN
2711	WRITE ( 14 ) surf_h(l)%sasws
2712	ENDIF
2713	ENDDO
2714	!
2715	!-- Write vertical surfaces
2716	DO l = 0, 3
2717	WRITE( dum, '(I1)') l
2718
2719	WRITE ( 14 ) 'surf_v(' // dum // ')%start_index '
2720	WRITE ( 14 ) surf_v(l)%start_index
2721	WRITE ( 14 ) 'surf_v(' // dum // ')%end_index '
2722	WRITE ( 14 ) surf_v(l)%end_index
2723
2724	WRITE ( 14 ) 'surf_v(' // dum // ')%us '
2725	IF ( ALLOCATED ( surf_v(l)%us ) ) THEN
2726	WRITE ( 14 ) surf_v(l)%us
2727	ENDIF
2728	WRITE ( 14 ) 'surf_v(' // dum // ')%ts '
2729	IF ( ALLOCATED ( surf_v(l)%ts ) ) THEN
2730	WRITE ( 14 ) surf_v(l)%ts
2731	ENDIF
2732	WRITE ( 14 ) 'surf_v(' // dum // ')%qs '
2733	IF ( ALLOCATED ( surf_v(l)%qs ) ) THEN
2734	WRITE ( 14 ) surf_v(l)%qs
2735	ENDIF
2736	WRITE ( 14 ) 'surf_v(' // dum // ')%ss '
2737	IF ( ALLOCATED ( surf_v(l)%ss ) ) THEN
2738	WRITE ( 14 ) surf_v(l)%ss
2739	ENDIF
2740	WRITE ( 14 ) 'surf_v(' // dum // ')%qcs '
2741	IF ( ALLOCATED ( surf_v(l)%qcs ) ) THEN
2742	WRITE ( 14 ) surf_v(l)%qcs
2743	ENDIF
2744	WRITE ( 14 ) 'surf_v(' // dum // ')%ncs '
2745	IF ( ALLOCATED ( surf_v(l)%ncs ) ) THEN
2746	WRITE ( 14 ) surf_v(l)%ncs
2747	ENDIF
2748	WRITE ( 14 ) 'surf_v(' // dum // ')%qrs '
2749	IF ( ALLOCATED ( surf_v(l)%qrs ) ) THEN
2750	WRITE ( 14 ) surf_v(l)%qrs
2751	ENDIF
2752	WRITE ( 14 ) 'surf_v(' // dum // ')%nrs '
2753	IF ( ALLOCATED ( surf_v(l)%nrs ) ) THEN
2754	WRITE ( 14 ) surf_v(l)%nrs
2755	ENDIF
2756	WRITE ( 14 ) 'surf_v(' // dum // ')%ol '
2757	IF ( ALLOCATED ( surf_v(l)%ol ) ) THEN
2758	WRITE ( 14 ) surf_v(l)%ol
2759	ENDIF
2760	WRITE ( 14 ) 'surf_v(' // dum // ')%rib '
2761	IF ( ALLOCATED ( surf_v(l)%rib ) ) THEN
2762	WRITE ( 14 ) surf_v(l)%rib
2763	ENDIF
2764	WRITE ( 14 ) 'surf_v(' // dum // ')%shf '
2765	IF ( ALLOCATED ( surf_v(l)%shf ) ) THEN
2766	WRITE ( 14 ) surf_v(l)%shf
2767	ENDIF
2768	WRITE ( 14 ) 'surf_v(' // dum // ')%qsws '
2769	IF ( ALLOCATED ( surf_v(l)%qsws ) ) THEN
2770	WRITE ( 14 ) surf_v(l)%qsws
2771	ENDIF
2772	WRITE ( 14 ) 'surf_v(' // dum // ')%ssws '
2773	IF ( ALLOCATED ( surf_v(l)%ssws ) ) THEN
2774	WRITE ( 14 ) surf_v(l)%ssws
2775	ENDIF
2776	#if defined( __chem )
2777	WRITE ( 14 ) 'surf_v(' // dum // ')%css '
2778	IF ( ALLOCATED ( surf_v(l)%css ) ) THEN
2779	WRITE ( 14 ) surf_v(l)%css
2780	ENDIF
2781	WRITE ( 14 ) 'surf_v(' // dum // ')%cssws '
2782	IF ( ALLOCATED ( surf_v(l)%cssws ) ) THEN
2783	WRITE ( 14 ) surf_v(l)%cssws
2784	ENDIF
2785	#endif
2786	WRITE ( 14 ) 'surf_v(' // dum // ')%qcsws '
2787	IF ( ALLOCATED ( surf_v(l)%qcsws ) ) THEN
2788	WRITE ( 14 ) surf_v(l)%qcsws
2789	ENDIF
2790	WRITE ( 14 ) 'surf_v(' // dum // ')%ncsws '
2791	IF ( ALLOCATED ( surf_v(l)%ncsws ) ) THEN
2792	WRITE ( 14 ) surf_v(l)%ncsws
2793	ENDIF
2794	WRITE ( 14 ) 'surf_v(' // dum // ')%qrsws '
2795	IF ( ALLOCATED ( surf_v(l)%qrsws ) ) THEN
2796	WRITE ( 14 ) surf_v(l)%qrsws
2797	ENDIF
2798	WRITE ( 14 ) 'surf_v(' // dum // ')%nrsws '
2799	IF ( ALLOCATED ( surf_v(l)%nrsws ) ) THEN
2800	WRITE ( 14 ) surf_v(l)%nrsws
2801	ENDIF
2802	WRITE ( 14 ) 'surf_v(' // dum // ')%sasws '
2803	IF ( ALLOCATED ( surf_v(l)%sasws ) ) THEN
2804	WRITE ( 14 ) surf_v(l)%sasws
2805	ENDIF
2806	WRITE ( 14 ) 'surf_v(' // dum // ')%mom_uv '
2807	IF ( ALLOCATED ( surf_v(l)%mom_flux_uv ) ) THEN
2808	WRITE ( 14 ) surf_v(l)%mom_flux_uv
2809	ENDIF
2810	WRITE ( 14 ) 'surf_v(' // dum // ')%mom_w '
2811	IF ( ALLOCATED ( surf_v(l)%mom_flux_w ) ) THEN
2812	WRITE ( 14 ) surf_v(l)%mom_flux_w
2813	ENDIF
2814	WRITE ( 14 ) 'surf_v(' // dum // ')%mom_tke '
2815	IF ( ALLOCATED ( surf_v(l)%mom_flux_tke ) ) THEN
2816	WRITE ( 14 ) surf_v(l)%mom_flux_tke
2817	ENDIF
2818
2819	ENDDO
2820
2821	WRITE ( 14 ) '* end surf * '
2822
2823	END SUBROUTINE surface_write_restart_data
2824
2825
2826	!------------------------------------------------------------------------------!
2827	! Description:
2828	! ------------
2829	!> Reads surface-related restart data. Please note, restart data for a certain
2830	!> surface orientation (e.g. horizontal upward-facing) is stored in one
2831	!> array, even if surface elements may belong to different surface types
2832	!> natural or urban for example). Surface elements are redistributed into its
2833	!> respective surface types within this routine. This allows e.g. changing the
2834	!> surface type after reading the restart data, which might be required in case
2835	!> of cyclic_fill mode.
2836	!------------------------------------------------------------------------------!
2837	SUBROUTINE surface_read_restart_data( ii, &
2838	nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
2839	nynfa, nyn_on_file, nysfa, nys_on_file, &
2840	offset_xa, offset_ya, overlap_count )
2841
2842	USE pegrid, &
2843	ONLY: numprocs_previous_run
2844
2845	CHARACTER (LEN=1) :: dum !< dummy to create correct string for reading input variable
2846	CHARACTER (LEN=30) :: field_chr !< input variable
2847
2848	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
2849	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
2850	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
2851	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
2852	INTEGER(iwp) :: k !< running index along z-direction
2853	INTEGER(iwp) :: l !< index variable for surface type
2854	INTEGER(iwp) :: m !< running index for surface elements, refers to gathered array encompassing all surface types
2855	INTEGER(iwp) :: mm !< running index for surface elements, refers to individual surface types
2856
2857	INTEGER(iwp) :: ii !< running index over input files
2858	INTEGER(iwp) :: kk !< running index over previous input files covering current local domain
2859	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
2860	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
2861	INTEGER(iwp) :: nxl_on_file !< index of left boundary on former local domain
2862	INTEGER(iwp) :: nxrc !< index of right boundary on current subdomain
2863	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
2864	INTEGER(iwp) :: nxr_on_file !< index of right boundary on former local domain
2865	INTEGER(iwp) :: nync !< index of north boundary on current subdomain
2866	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
2867	INTEGER(iwp) :: nyn_on_file !< index of norht boundary on former local domain
2868	INTEGER(iwp) :: nysc !< index of south boundary on current subdomain
2869	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
2870	INTEGER(iwp) :: nys_on_file !< index of south boundary on former local domain
2871	INTEGER(iwp) :: overlap_count !< number of overlaps
2872
2873	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !<
2874	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !<
2875	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !<
2876	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !<
2877	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !<
2878	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !<
2879
2880
2881	LOGICAL :: horizontal_surface !< flag indicating horizontal surfaces
2882	LOGICAL :: surf_match_def !< flag indicating that surface element is of default type
2883	LOGICAL :: surf_match_lsm !< flag indicating that surface element is of natural type
2884	LOGICAL :: surf_match_usm !< flag indicating that surface element is of urban type
2885	LOGICAL :: vertical_surface !< flag indicating vertical surfaces
2886
2887	TYPE(surf_type), DIMENSION(0:2) :: surf_h !< horizontal surface type on file
2888	TYPE(surf_type), DIMENSION(0:3) :: surf_v !< vertical surface type on file
2889
2890	!
2891	!-- Read number of respective surface elements on file
2892	READ ( 13 ) field_chr
2893	IF ( TRIM( field_chr ) /= 'ns_h_on_file' ) THEN
2894	!
2895	!-- Add a proper error message
2896	ENDIF
2897	READ ( 13 ) ns_h_on_file
2898
2899	READ ( 13 ) field_chr
2900	IF ( TRIM( field_chr ) /= 'ns_v_on_file' ) THEN
2901	!
2902	!-- Add a proper error message
2903	ENDIF
2904	READ ( 13 ) ns_v_on_file
2905	!
2906	!-- Allocate memory for number of surface elements on file. Please note,
2907	!-- these number is not necessarily the same as the final number of surface
2908	!-- elements on local domain, which is the case if processor topology changes
2909	!-- during restart runs.
2910	!-- Horizontal upward facing
2911	surf_h(0)%ns = ns_h_on_file(0)
2912	CALL allocate_surface_attributes_h( surf_h(0), &
2913	nys_on_file, nyn_on_file, &
2914	nxl_on_file, nxr_on_file )
2915	!
2916	!-- Horizontal downward facing
2917	surf_h(1)%ns = ns_h_on_file(1)
2918	CALL allocate_surface_attributes_h( surf_h(1), &
2919	nys_on_file, nyn_on_file, &
2920	nxl_on_file, nxr_on_file )
2921	!
2922	!-- Model top
2923	surf_h(2)%ns = ns_h_on_file(2)
2924	CALL allocate_surface_attributes_h_top( surf_h(2), &
2925	nys_on_file, nyn_on_file, &
2926	nxl_on_file, nxr_on_file )
2927	!
2928	!-- Vertical surfaces
2929	DO l = 0, 3
2930	surf_v(l)%ns = ns_v_on_file(l)
2931	CALL allocate_surface_attributes_v( surf_v(l), &
2932	nys_on_file, nyn_on_file, &
2933	nxl_on_file, nxr_on_file )
2934	ENDDO
2935
2936	IF ( initializing_actions == 'read_restart_data' .OR. &
2937	initializing_actions == 'cyclic_fill' ) THEN
2938	!
2939	!-- Initial setting of flags for horizontal and vertical surfaces, will
2940	!-- be set after start- and end-indices are read.
2941	horizontal_surface = .FALSE.
2942	vertical_surface = .FALSE.
2943
2944	READ ( 13 ) field_chr
2945
2946	DO WHILE ( TRIM( field_chr ) /= '* end surf *' )
2947	!
2948	!-- Map data on file as often as needed (data are read only for k=1)
2949	DO kk = 1, overlap_count
2950	!
2951	!-- Get the index range of the subdomain on file which overlap with the
2952	!-- current subdomain
2953	nxlf = nxlfa(ii,kk)
2954	nxlc = nxlfa(ii,kk) + offset_xa(ii,kk)
2955	nxrf = nxrfa(ii,kk)
2956	nxrc = nxrfa(ii,kk) + offset_xa(ii,kk)
2957	nysf = nysfa(ii,kk)
2958	nysc = nysfa(ii,kk) + offset_ya(ii,kk)
2959	nynf = nynfa(ii,kk)
2960	nync = nynfa(ii,kk) + offset_ya(ii,kk)
2961
2962	SELECT CASE ( TRIM( field_chr ) )
2963
2964	CASE ( 'surf_h(0)%start_index' )
2965	IF ( kk == 1 ) &
2966	READ ( 13 ) surf_h(0)%start_index
2967	l = 0
2968	CASE ( 'surf_h(0)%end_index' )
2969	IF ( kk == 1 ) &
2970	READ ( 13 ) surf_h(0)%end_index
2971	horizontal_surface = .TRUE.
2972	vertical_surface = .FALSE.
2973	CASE ( 'surf_h(0)%us' )
2974	IF ( ALLOCATED( surf_h(0)%us ) .AND. kk == 1 ) &
2975	READ ( 13 ) surf_h(0)%us
2976	CASE ( 'surf_h(0)%ts' )
2977	IF ( ALLOCATED( surf_h(0)%ts ) .AND. kk == 1 ) &
2978	READ ( 13 ) surf_h(0)%ts
2979	CASE ( 'surf_h(0)%qs' )
2980	IF ( ALLOCATED( surf_h(0)%qs ) .AND. kk == 1 ) &
2981	READ ( 13 ) surf_h(0)%qs
2982	CASE ( 'surf_h(0)%ss' )
2983	IF ( ALLOCATED( surf_h(0)%ss ) .AND. kk == 1 ) &
2984	READ ( 13 ) surf_h(0)%ss
2985	CASE ( 'surf_h(0)%qcs' )
2986	IF ( ALLOCATED( surf_h(0)%qcs ) .AND. kk == 1 ) &
2987	READ ( 13 ) surf_h(0)%qcs
2988	CASE ( 'surf_h(0)%ncs' )
2989	IF ( ALLOCATED( surf_h(0)%ncs ) .AND. kk == 1 ) &
2990	READ ( 13 ) surf_h(0)%ncs
2991	CASE ( 'surf_h(0)%qrs' )
2992	IF ( ALLOCATED( surf_h(0)%qrs ) .AND. kk == 1 ) &
2993	READ ( 13 ) surf_h(0)%qrs
2994	CASE ( 'surf_h(0)%nrs' )
2995	IF ( ALLOCATED( surf_h(0)%nrs ) .AND. kk == 1 ) &
2996	READ ( 13 ) surf_h(0)%nrs
2997	CASE ( 'surf_h(0)%ol' )
2998	IF ( ALLOCATED( surf_h(0)%ol ) .AND. kk == 1 ) &
2999	READ ( 13 ) surf_h(0)%ol
3000	CASE ( 'surf_h(0)%rib' )
3001	IF ( ALLOCATED( surf_h(0)%rib ) .AND. kk == 1 ) &
3002	READ ( 13 ) surf_h(0)%rib
3003	CASE ( 'surf_h(0)%usws' )
3004	IF ( ALLOCATED( surf_h(0)%usws ) .AND. kk == 1 ) &
3005	READ ( 13 ) surf_h(0)%usws
3006	CASE ( 'surf_h(0)%vsws' )
3007	IF ( ALLOCATED( surf_h(0)%vsws ) .AND. kk == 1 ) &
3008	READ ( 13 ) surf_h(0)%vsws
3009	CASE ( 'surf_h(0)%shf' )
3010	IF ( ALLOCATED( surf_h(0)%shf ) .AND. kk == 1 ) &
3011	READ ( 13 ) surf_h(0)%shf
3012	CASE ( 'surf_h(0)%qsws' )
3013	IF ( ALLOCATED( surf_h(0)%qsws ) .AND. kk == 1 ) &
3014	READ ( 13 ) surf_h(0)%qsws
3015	CASE ( 'surf_h(0)%ssws' )
3016	IF ( ALLOCATED( surf_h(0)%ssws ) .AND. kk == 1 ) &
3017	READ ( 13 ) surf_h(0)%ssws
3018	#if defined( __chem )
3019	CASE ( 'surf_h(0)%css' )
3020	IF ( ALLOCATED( surf_h(0)%css ) .AND. kk == 1 ) &
3021	READ ( 13 ) surf_h(0)%css
3022	CASE ( 'surf_h(0)%cssws' )
3023	IF ( ALLOCATED( surf_h(0)%cssws ) .AND. kk == 1 ) &
3024	READ ( 13 ) surf_h(0)%cssws
3025	#endif
3026	CASE ( 'surf_h(0)%qcsws' )
3027	IF ( ALLOCATED( surf_h(0)%qcsws ) .AND. kk == 1 ) &
3028	READ ( 13 ) surf_h(0)%qcsws
3029	CASE ( 'surf_h(0)%ncsws' )
3030	IF ( ALLOCATED( surf_h(0)%ncsws ) .AND. kk == 1 ) &
3031	READ ( 13 ) surf_h(0)%ncsws
3032	CASE ( 'surf_h(0)%qrsws' )
3033	IF ( ALLOCATED( surf_h(0)%qrsws ) .AND. kk == 1 ) &
3034	READ ( 13 ) surf_h(0)%qrsws
3035	CASE ( 'surf_h(0)%nrsws' )
3036	IF ( ALLOCATED( surf_h(0)%nrsws ) .AND. kk == 1 ) &
3037	READ ( 13 ) surf_h(0)%nrsws
3038	CASE ( 'surf_h(0)%sasws' )
3039	IF ( ALLOCATED( surf_h(0)%sasws ) .AND. kk == 1 ) &
3040	READ ( 13 ) surf_h(0)%sasws
3041
3042	CASE ( 'surf_h(1)%start_index' )
3043	IF ( kk == 1 ) &
3044	READ ( 13 ) surf_h(1)%start_index
3045	l = 1
3046	CASE ( 'surf_h(1)%end_index' )
3047	IF ( kk == 1 ) &
3048	READ ( 13 ) surf_h(1)%end_index
3049	CASE ( 'surf_h(1)%us' )
3050	IF ( ALLOCATED( surf_h(1)%us ) .AND. kk == 1 ) &
3051	READ ( 13 ) surf_h(1)%us
3052	CASE ( 'surf_h(1)%ts' )
3053	IF ( ALLOCATED( surf_h(1)%ts ) .AND. kk == 1 ) &
3054	READ ( 13 ) surf_h(1)%ts
3055	CASE ( 'surf_h(1)%qs' )
3056	IF ( ALLOCATED( surf_h(1)%qs ) .AND. kk == 1 ) &
3057	READ ( 13 ) surf_h(1)%qs
3058	CASE ( 'surf_h(1)%ss' )
3059	IF ( ALLOCATED( surf_h(1)%ss ) .AND. kk == 1 ) &
3060	READ ( 13 ) surf_h(1)%ss
3061	CASE ( 'surf_h(1)%qcs' )
3062	IF ( ALLOCATED( surf_h(1)%qcs ) .AND. kk == 1 ) &
3063	READ ( 13 ) surf_h(1)%qcs
3064	CASE ( 'surf_h(1)%ncs' )
3065	IF ( ALLOCATED( surf_h(1)%ncs ) .AND. kk == 1 ) &
3066	READ ( 13 ) surf_h(1)%ncs
3067	CASE ( 'surf_h(1)%qrs' )
3068	IF ( ALLOCATED( surf_h(1)%qrs ) .AND. kk == 1 ) &
3069	READ ( 13 ) surf_h(1)%qrs
3070	CASE ( 'surf_h(1)%nrs' )
3071	IF ( ALLOCATED( surf_h(1)%nrs ) .AND. kk == 1 ) &
3072	READ ( 13 ) surf_h(1)%nrs
3073	CASE ( 'surf_h(1)%ol' )
3074	IF ( ALLOCATED( surf_h(1)%ol ) .AND. kk == 1 ) &
3075	READ ( 13 ) surf_h(1)%ol
3076	CASE ( 'surf_h(1)%rib' )
3077	IF ( ALLOCATED( surf_h(1)%rib ) .AND. kk == 1 ) &
3078	READ ( 13 ) surf_h(1)%rib
3079	CASE ( 'surf_h(1)%usws' )
3080	IF ( ALLOCATED( surf_h(1)%usws ) .AND. kk == 1 ) &
3081	READ ( 13 ) surf_h(1)%usws
3082	CASE ( 'surf_h(1)%vsws' )
3083	IF ( ALLOCATED( surf_h(1)%vsws ) .AND. kk == 1 ) &
3084	READ ( 13 ) surf_h(1)%vsws
3085	CASE ( 'surf_h(1)%shf' )
3086	IF ( ALLOCATED( surf_h(1)%shf ) .AND. kk == 1 ) &
3087	READ ( 13 ) surf_h(1)%shf
3088	CASE ( 'surf_h(1)%qsws' )
3089	IF ( ALLOCATED( surf_h(1)%qsws ) .AND. kk == 1 ) &
3090	READ ( 13 ) surf_h(1)%qsws
3091	CASE ( 'surf_h(1)%ssws' )
3092	IF ( ALLOCATED( surf_h(1)%ssws ) .AND. kk == 1 ) &
3093	READ ( 13 ) surf_h(1)%ssws
3094	#if defined( __chem )
3095	CASE ( 'surf_h(1)%css' )
3096	IF ( ALLOCATED( surf_h(1)%css ) .AND. kk == 1 ) &
3097	READ ( 13 ) surf_h(1)%css
3098	CASE ( 'surf_h(1)%cssws' )
3099	IF ( ALLOCATED( surf_h(1)%cssws ) .AND. kk == 1 ) &
3100	READ ( 13 ) surf_h(1)%cssws
3101	#endif
3102	CASE ( 'surf_h(1)%qcsws' )
3103	IF ( ALLOCATED( surf_h(1)%qcsws ) .AND. kk == 1 ) &
3104	READ ( 13 ) surf_h(1)%qcsws
3105	CASE ( 'surf_h(1)%ncsws' )
3106	IF ( ALLOCATED( surf_h(1)%ncsws ) .AND. kk == 1 ) &
3107	READ ( 13 ) surf_h(1)%ncsws
3108	CASE ( 'surf_h(1)%qrsws' )
3109	IF ( ALLOCATED( surf_h(1)%qrsws ) .AND. kk == 1 ) &
3110	READ ( 13 ) surf_h(1)%qrsws
3111	CASE ( 'surf_h(1)%nrsws' )
3112	IF ( ALLOCATED( surf_h(1)%nrsws ) .AND. kk == 1 ) &
3113	READ ( 13 ) surf_h(1)%nrsws
3114	CASE ( 'surf_h(1)%sasws' )
3115	IF ( ALLOCATED( surf_h(1)%sasws ) .AND. kk == 1 ) &
3116	READ ( 13 ) surf_h(1)%sasws
3117
3118	CASE ( 'surf_h(2)%start_index' )
3119	IF ( kk == 1 ) &
3120	READ ( 13 ) surf_h(2)%start_index
3121	l = 2
3122	CASE ( 'surf_h(2)%end_index' )
3123	IF ( kk == 1 ) &
3124	READ ( 13 ) surf_h(2)%end_index
3125	CASE ( 'surf_h(2)%us' )
3126	IF ( ALLOCATED( surf_h(2)%us ) .AND. kk == 1 ) &
3127	READ ( 13 ) surf_h(2)%us
3128	CASE ( 'surf_h(2)%ts' )
3129	IF ( ALLOCATED( surf_h(2)%ts ) .AND. kk == 1 ) &
3130	READ ( 13 ) surf_h(2)%ts
3131	CASE ( 'surf_h(2)%qs' )
3132	IF ( ALLOCATED( surf_h(2)%qs ) .AND. kk == 1 ) &
3133	READ ( 13 ) surf_h(2)%qs
3134	CASE ( 'surf_h(2)%ss' )
3135	IF ( ALLOCATED( surf_h(2)%ss ) .AND. kk == 1 ) &
3136	READ ( 13 ) surf_h(2)%ss
3137	CASE ( 'surf_h(2)%qcs' )
3138	IF ( ALLOCATED( surf_h(2)%qcs ) .AND. kk == 1 ) &
3139	READ ( 13 ) surf_h(2)%qcs
3140	CASE ( 'surf_h(2)%ncs' )
3141	IF ( ALLOCATED( surf_h(2)%ncs ) .AND. kk == 1 ) &
3142	READ ( 13 ) surf_h(2)%ncs
3143	CASE ( 'surf_h(2)%qrs' )
3144	IF ( ALLOCATED( surf_h(2)%qrs ) .AND. kk == 1 ) &
3145	READ ( 13 ) surf_h(2)%qrs
3146	CASE ( 'surf_h(2)%nrs' )
3147	IF ( ALLOCATED( surf_h(2)%nrs ) .AND. kk == 1 ) &
3148	READ ( 13 ) surf_h(2)%nrs
3149	CASE ( 'surf_h(2)%ol' )
3150	IF ( ALLOCATED( surf_h(2)%ol ) .AND. kk == 1 ) &
3151	READ ( 13 ) surf_h(2)%ol
3152	CASE ( 'surf_h(2)%rib' )
3153	IF ( ALLOCATED( surf_h(2)%rib ) .AND. kk == 1 ) &
3154	READ ( 13 ) surf_h(2)%rib
3155	CASE ( 'surf_h(2)%usws' )
3156	IF ( ALLOCATED( surf_h(2)%usws ) .AND. kk == 1 ) &
3157	READ ( 13 ) surf_h(2)%usws
3158	CASE ( 'surf_h(2)%vsws' )
3159	IF ( ALLOCATED( surf_h(2)%vsws ) .AND. kk == 1 ) &
3160	READ ( 13 ) surf_h(2)%vsws
3161	CASE ( 'surf_h(2)%shf' )
3162	IF ( ALLOCATED( surf_h(2)%shf ) .AND. kk == 1 ) &
3163	READ ( 13 ) surf_h(2)%shf
3164	CASE ( 'surf_h(2)%qsws' )
3165	IF ( ALLOCATED( surf_h(2)%qsws ) .AND. kk == 1 ) &
3166	READ ( 13 ) surf_h(2)%qsws
3167	CASE ( 'surf_h(2)%ssws' )
3168	IF ( ALLOCATED( surf_h(2)%ssws ) .AND. kk == 1 ) &
3169	READ ( 13 ) surf_h(2)%ssws
3170	#if defined( __chem )
3171	CASE ( 'surf_h(2)%css' )
3172	IF ( ALLOCATED( surf_h(2)%css ) .AND. kk == 1 ) &
3173	READ ( 13 ) surf_h(2)%css
3174	CASE ( 'surf_h(2)%cssws' )
3175	IF ( ALLOCATED( surf_h(2)%cssws ) .AND. kk == 1 ) &
3176	READ ( 13 ) surf_h(2)%cssws
3177	#endif
3178	CASE ( 'surf_h(2)%qcsws' )
3179	IF ( ALLOCATED( surf_h(2)%qcsws ) .AND. kk == 1 ) &
3180	READ ( 13 ) surf_h(2)%qcsws
3181	CASE ( 'surf_h(2)%ncsws' )
3182	IF ( ALLOCATED( surf_h(2)%ncsws ) .AND. kk == 1 ) &
3183	READ ( 13 ) surf_h(2)%ncsws
3184	CASE ( 'surf_h(2)%qrsws' )
3185	IF ( ALLOCATED( surf_h(2)%qrsws ) .AND. kk == 1 ) &
3186	READ ( 13 ) surf_h(2)%qrsws
3187	CASE ( 'surf_h(2)%nrsws' )
3188	IF ( ALLOCATED( surf_h(2)%nrsws ) .AND. kk == 1 ) &
3189	READ ( 13 ) surf_h(2)%nrsws
3190	CASE ( 'surf_h(2)%sasws' )
3191	IF ( ALLOCATED( surf_h(2)%sasws ) .AND. kk == 1 ) &
3192	READ ( 13 ) surf_h(2)%sasws
3193
3194	CASE ( 'surf_v(0)%start_index' )
3195	IF ( kk == 1 ) &
3196	READ ( 13 ) surf_v(0)%start_index
3197	l = 0
3198	horizontal_surface = .FALSE.
3199	vertical_surface = .TRUE.
3200	CASE ( 'surf_v(0)%end_index' )
3201	IF ( kk == 1 ) &
3202	READ ( 13 ) surf_v(0)%end_index
3203	CASE ( 'surf_v(0)%us' )
3204	IF ( ALLOCATED( surf_v(0)%us ) .AND. kk == 1 ) &
3205	READ ( 13 ) surf_v(0)%us
3206	CASE ( 'surf_v(0)%ts' )
3207	IF ( ALLOCATED( surf_v(0)%ts ) .AND. kk == 1 ) &
3208	READ ( 13 ) surf_v(0)%ts
3209	CASE ( 'surf_v(0)%qs' )
3210	IF ( ALLOCATED( surf_v(0)%qs ) .AND. kk == 1 ) &
3211	READ ( 13 ) surf_v(0)%qs
3212	CASE ( 'surf_v(0)%ss' )
3213	IF ( ALLOCATED( surf_v(0)%ss ) .AND. kk == 1 ) &
3214	READ ( 13 ) surf_v(0)%ss
3215	CASE ( 'surf_v(0)%qcs' )
3216	IF ( ALLOCATED( surf_v(0)%qcs ) .AND. kk == 1 ) &
3217	READ ( 13 ) surf_v(0)%qcs
3218	CASE ( 'surf_v(0)%ncs' )
3219	IF ( ALLOCATED( surf_v(0)%ncs ) .AND. kk == 1 ) &
3220	READ ( 13 ) surf_v(0)%ncs
3221	CASE ( 'surf_v(0)%qrs' )
3222	IF ( ALLOCATED( surf_v(0)%qrs ) .AND. kk == 1 ) &
3223	READ ( 13 ) surf_v(0)%qrs
3224	CASE ( 'surf_v(0)%nrs' )
3225	IF ( ALLOCATED( surf_v(0)%nrs ) .AND. kk == 1 ) &
3226	READ ( 13 ) surf_v(0)%nrs
3227	CASE ( 'surf_v(0)%ol' )
3228	IF ( ALLOCATED( surf_v(0)%ol ) .AND. kk == 1 ) &
3229	READ ( 13 ) surf_v(0)%ol
3230	CASE ( 'surf_v(0)%rib' )
3231	IF ( ALLOCATED( surf_v(0)%rib ) .AND. kk == 1 ) &
3232	READ ( 13 ) surf_v(0)%rib
3233	CASE ( 'surf_v(0)%shf' )
3234	IF ( ALLOCATED( surf_v(0)%shf ) .AND. kk == 1 ) &
3235	READ ( 13 ) surf_v(0)%shf
3236	CASE ( 'surf_v(0)%qsws' )
3237	IF ( ALLOCATED( surf_v(0)%qsws ) .AND. kk == 1 ) &
3238	READ ( 13 ) surf_v(0)%qsws
3239	CASE ( 'surf_v(0)%ssws' )
3240	IF ( ALLOCATED( surf_v(0)%ssws ) .AND. kk == 1 ) &
3241	READ ( 13 ) surf_v(0)%ssws
3242	#if defined( __chem )
3243	CASE ( 'surf_v(0)%css' )
3244	IF ( ALLOCATED( surf_v(0)%css ) .AND. kk == 1 ) &
3245	READ ( 13 ) surf_v(0)%css
3246	CASE ( 'surf_v(0)%cssws' )
3247	IF ( ALLOCATED( surf_v(0)%cssws ) .AND. kk == 1 ) &
3248	READ ( 13 ) surf_v(0)%cssws
3249	#endif
3250	CASE ( 'surf_v(0)%qcsws' )
3251	IF ( ALLOCATED( surf_v(0)%qcsws ) .AND. kk == 1 ) &
3252	READ ( 13 ) surf_v(0)%qcsws
3253	CASE ( 'surf_v(0)%ncsws' )
3254	IF ( ALLOCATED( surf_v(0)%ncsws ) .AND. kk == 1 ) &
3255	READ ( 13 ) surf_v(0)%ncsws
3256	CASE ( 'surf_v(0)%qrsws' )
3257	IF ( ALLOCATED( surf_v(0)%qrsws ) .AND. kk == 1 ) &
3258	READ ( 13 ) surf_v(0)%qrsws
3259	CASE ( 'surf_v(0)%nrsws' )
3260	IF ( ALLOCATED( surf_v(0)%nrsws ) .AND. kk == 1 ) &
3261	READ ( 13 ) surf_v(0)%nrsws
3262	CASE ( 'surf_v(0)%sasws' )
3263	IF ( ALLOCATED( surf_v(0)%sasws ) .AND. kk == 1 ) &
3264	READ ( 13 ) surf_v(0)%sasws
3265	CASE ( 'surf_v(0)%mom_uv' )
3266	IF ( ALLOCATED( surf_v(0)%mom_flux_uv ) .AND. kk == 1 )&
3267	READ ( 13 ) surf_v(0)%mom_flux_uv
3268	CASE ( 'surf_v(0)%mom_w' )
3269	IF ( ALLOCATED( surf_v(0)%mom_flux_w ) .AND. kk == 1 ) &
3270	READ ( 13 ) surf_v(0)%mom_flux_w
3271	CASE ( 'surf_v(0)%mom_tke' )
3272	IF ( ALLOCATED( surf_v(0)%mom_flux_tke ) .AND. kk == 1 )&
3273	READ ( 13 ) surf_v(0)%mom_flux_tke
3274
3275	CASE ( 'surf_v(1)%start_index' )
3276	IF ( kk == 1 ) &
3277	READ ( 13 ) surf_v(1)%start_index
3278	l = 1
3279	CASE ( 'surf_v(1)%end_index' )
3280	IF ( kk == 1 ) &
3281	READ ( 13 ) surf_v(1)%end_index
3282	CASE ( 'surf_v(1)%us' )
3283	IF ( ALLOCATED( surf_v(1)%us ) .AND. kk == 1 ) &
3284	READ ( 13 ) surf_v(1)%us
3285	CASE ( 'surf_v(1)%ts' )
3286	IF ( ALLOCATED( surf_v(1)%ts ) .AND. kk == 1 ) &
3287	READ ( 13 ) surf_v(1)%ts
3288	CASE ( 'surf_v(1)%qs' )
3289	IF ( ALLOCATED( surf_v(1)%qs ) .AND. kk == 1 ) &
3290	READ ( 13 ) surf_v(1)%qs
3291	CASE ( 'surf_v(1)%ss' )
3292	IF ( ALLOCATED( surf_v(1)%ss ) .AND. kk == 1 ) &
3293	READ ( 13 ) surf_v(1)%ss
3294	CASE ( 'surf_v(1)%qcs' )
3295	IF ( ALLOCATED( surf_v(1)%qcs ) .AND. kk == 1 ) &
3296	READ ( 13 ) surf_v(1)%qcs
3297	CASE ( 'surf_v(1)%ncs' )
3298	IF ( ALLOCATED( surf_v(1)%ncs ) .AND. kk == 1 ) &
3299	READ ( 13 ) surf_v(1)%ncs
3300	CASE ( 'surf_v(1)%qrs' )
3301	IF ( ALLOCATED( surf_v(1)%qrs ) .AND. kk == 1 ) &
3302	READ ( 13 ) surf_v(1)%qrs
3303	CASE ( 'surf_v(1)%nrs' )
3304	IF ( ALLOCATED( surf_v(1)%nrs ) .AND. kk == 1 ) &
3305	READ ( 13 ) surf_v(1)%nrs
3306	CASE ( 'surf_v(1)%ol' )
3307	IF ( ALLOCATED( surf_v(1)%ol ) .AND. kk == 1 ) &
3308	READ ( 13 ) surf_v(1)%ol
3309	CASE ( 'surf_v(1)%rib' )
3310	IF ( ALLOCATED( surf_v(1)%rib ) .AND. kk == 1 ) &
3311	READ ( 13 ) surf_v(1)%rib
3312	CASE ( 'surf_v(1)%shf' )
3313	IF ( ALLOCATED( surf_v(1)%shf ) .AND. kk == 1 ) &
3314	READ ( 13 ) surf_v(1)%shf
3315	CASE ( 'surf_v(1)%qsws' )
3316	IF ( ALLOCATED( surf_v(1)%qsws ) .AND. kk == 1 ) &
3317	READ ( 13 ) surf_v(1)%qsws
3318	CASE ( 'surf_v(1)%ssws' )
3319	IF ( ALLOCATED( surf_v(1)%ssws ) .AND. kk == 1 ) &
3320	READ ( 13 ) surf_v(1)%ssws
3321	#if defined( __chem )
3322	CASE ( 'surf_v(1)%css' )
3323	IF ( ALLOCATED( surf_v(1)%css ) .AND. kk == 1 ) &
3324	READ ( 13 ) surf_v(1)%css
3325	CASE ( 'surf_v(1)%cssws' )
3326	IF ( ALLOCATED( surf_v(1)%cssws ) .AND. kk == 1 ) &
3327	READ ( 13 ) surf_v(1)%cssws
3328	#endif
3329	CASE ( 'surf_v(1)%qcsws' )
3330	IF ( ALLOCATED( surf_v(1)%qcsws ) .AND. kk == 1 ) &
3331	READ ( 13 ) surf_v(1)%qcsws
3332	CASE ( 'surf_v(1)%ncsws' )
3333	IF ( ALLOCATED( surf_v(1)%ncsws ) .AND. kk == 1 ) &
3334	READ ( 13 ) surf_v(1)%ncsws
3335	CASE ( 'surf_v(1)%qrsws' )
3336	IF ( ALLOCATED( surf_v(1)%qrsws ) .AND. kk == 1 ) &
3337	READ ( 13 ) surf_v(1)%qrsws
3338	CASE ( 'surf_v(1)%nrsws' )
3339	IF ( ALLOCATED( surf_v(1)%nrsws ) .AND. kk == 1 ) &
3340	READ ( 13 ) surf_v(1)%nrsws
3341	CASE ( 'surf_v(1)%sasws' )
3342	IF ( ALLOCATED( surf_v(1)%sasws ) .AND. kk == 1 ) &
3343	READ ( 13 ) surf_v(1)%sasws
3344	CASE ( 'surf_v(1)%mom_uv' )
3345	IF ( ALLOCATED( surf_v(1)%mom_flux_uv ) .AND. kk == 1 )&
3346	READ ( 13 ) surf_v(1)%mom_flux_uv
3347	CASE ( 'surf_v(1)%mom_w' )
3348	IF ( ALLOCATED( surf_v(1)%mom_flux_w ) .AND. kk == 1 ) &
3349	READ ( 13 ) surf_v(1)%mom_flux_w
3350	CASE ( 'surf_v(1)%mom_tke' )
3351	IF ( ALLOCATED( surf_v(1)%mom_flux_tke ) .AND. kk == 1 )&
3352	READ ( 13 ) surf_v(1)%mom_flux_tke
3353
3354	CASE ( 'surf_v(2)%start_index' )
3355	IF ( kk == 1 ) &
3356	READ ( 13 ) surf_v(2)%start_index
3357	l = 2
3358	CASE ( 'surf_v(2)%end_index' )
3359	IF ( kk == 1 ) &
3360	READ ( 13 ) surf_v(2)%end_index
3361	CASE ( 'surf_v(2)%us' )
3362	IF ( ALLOCATED( surf_v(2)%us ) .AND. kk == 1 ) &
3363	READ ( 13 ) surf_v(2)%us
3364	CASE ( 'surf_v(2)%ts' )
3365	IF ( ALLOCATED( surf_v(2)%ts ) .AND. kk == 1 ) &
3366	READ ( 13 ) surf_v(2)%ts
3367	CASE ( 'surf_v(2)%qs' )
3368	IF ( ALLOCATED( surf_v(2)%qs ) .AND. kk == 1 ) &
3369	READ ( 13 ) surf_v(2)%qs
3370	CASE ( 'surf_v(2)%ss' )
3371	IF ( ALLOCATED( surf_v(2)%ss ) .AND. kk == 1 ) &
3372	READ ( 13 ) surf_v(2)%ss
3373	CASE ( 'surf_v(2)%qcs' )
3374	IF ( ALLOCATED( surf_v(2)%qcs ) .AND. kk == 1 ) &
3375	READ ( 13 ) surf_v(2)%qcs
3376	CASE ( 'surf_v(2)%ncs' )
3377	IF ( ALLOCATED( surf_v(2)%ncs ) .AND. kk == 1 ) &
3378	READ ( 13 ) surf_v(2)%ncs
3379	CASE ( 'surf_v(2)%qrs' )
3380	IF ( ALLOCATED( surf_v(2)%qrs ) .AND. kk == 1 ) &
3381	READ ( 13 ) surf_v(2)%qrs
3382	CASE ( 'surf_v(2)%nrs' )
3383	IF ( ALLOCATED( surf_v(2)%nrs ) .AND. kk == 1 ) &
3384	READ ( 13 ) surf_v(2)%nrs
3385	CASE ( 'surf_v(2)%ol' )
3386	IF ( ALLOCATED( surf_v(2)%ol ) .AND. kk == 1 ) &
3387	READ ( 13 ) surf_v(2)%ol
3388	CASE ( 'surf_v(2)%rib' )
3389	IF ( ALLOCATED( surf_v(2)%rib ) .AND. kk == 1 ) &
3390	READ ( 13 ) surf_v(2)%rib
3391	CASE ( 'surf_v(2)%shf' )
3392	IF ( ALLOCATED( surf_v(2)%shf ) .AND. kk == 1 ) &
3393	READ ( 13 ) surf_v(2)%shf
3394	CASE ( 'surf_v(2)%qsws' )
3395	IF ( ALLOCATED( surf_v(2)%qsws ) .AND. kk == 1 ) &
3396	READ ( 13 ) surf_v(2)%qsws
3397	CASE ( 'surf_v(2)%ssws' )
3398	IF ( ALLOCATED( surf_v(2)%ssws ) .AND. kk == 1 ) &
3399	READ ( 13 ) surf_v(2)%ssws
3400	#if defined( __chem )
3401	CASE ( 'surf_v(2)%css' )
3402	IF ( ALLOCATED( surf_v(2)%css ) .AND. kk == 1 ) &
3403	READ ( 13 ) surf_v(2)%css
3404	CASE ( 'surf_v(2)%cssws' )
3405	IF ( ALLOCATED( surf_v(2)%cssws ) .AND. kk == 1 ) &
3406	READ ( 13 ) surf_v(2)%cssws
3407	#endif
3408	CASE ( 'surf_v(2)%qcsws' )
3409	IF ( ALLOCATED( surf_v(2)%qcsws ) .AND. kk == 1 ) &
3410	READ ( 13 ) surf_v(2)%qcsws
3411	CASE ( 'surf_v(2)%ncsws' )
3412	IF ( ALLOCATED( surf_v(2)%ncsws ) .AND. kk == 1 ) &
3413	READ ( 13 ) surf_v(2)%ncsws
3414	CASE ( 'surf_v(2)%qrsws' )
3415	IF ( ALLOCATED( surf_v(2)%qrsws ) .AND. kk == 1 ) &
3416	READ ( 13 ) surf_v(2)%qrsws
3417	CASE ( 'surf_v(2)%nrsws' )
3418	IF ( ALLOCATED( surf_v(2)%nrsws ) .AND. kk == 1 ) &
3419	READ ( 13 ) surf_v(2)%nrsws
3420	CASE ( 'surf_v(2)%sasws' )
3421	IF ( ALLOCATED( surf_v(2)%sasws ) .AND. kk == 1 ) &
3422	READ ( 13 ) surf_v(2)%sasws
3423	CASE ( 'surf_v(2)%mom_uv' )
3424	IF ( ALLOCATED( surf_v(2)%mom_flux_uv ) .AND. kk == 1 )&
3425	READ ( 13 ) surf_v(2)%mom_flux_uv
3426	CASE ( 'surf_v(2)%mom_w' )
3427	IF ( ALLOCATED( surf_v(2)%mom_flux_w ) .AND. kk == 1 ) &
3428	READ ( 13 ) surf_v(2)%mom_flux_w
3429	CASE ( 'surf_v(2)%mom_tke' )
3430	IF ( ALLOCATED( surf_v(2)%mom_flux_tke ) .AND. kk == 1 )&
3431	READ ( 13 ) surf_v(2)%mom_flux_tke
3432
3433	CASE ( 'surf_v(3)%start_index' )
3434	IF ( kk == 1 ) &
3435	READ ( 13 ) surf_v(3)%start_index
3436	l = 3
3437	CASE ( 'surf_v(3)%end_index' )
3438	IF ( kk == 1 ) &
3439	READ ( 13 ) surf_v(3)%end_index
3440	CASE ( 'surf_v(3)%us' )
3441	IF ( ALLOCATED( surf_v(3)%us ) .AND. kk == 1 ) &
3442	READ ( 13 ) surf_v(3)%us
3443	CASE ( 'surf_v(3)%ts' )
3444	IF ( ALLOCATED( surf_v(3)%ts ) .AND. kk == 1 ) &
3445	READ ( 13 ) surf_v(3)%ts
3446	CASE ( 'surf_v(3)%qs' )
3447	IF ( ALLOCATED( surf_v(3)%qs ) .AND. kk == 1 ) &
3448	READ ( 13 ) surf_v(3)%qs
3449	CASE ( 'surf_v(3)%ss' )
3450	IF ( ALLOCATED( surf_v(3)%ss ) .AND. kk == 1 ) &
3451	READ ( 13 ) surf_v(3)%ss
3452	CASE ( 'surf_v(3)%qcs' )
3453	IF ( ALLOCATED( surf_v(3)%qcs ) .AND. kk == 1 ) &
3454	READ ( 13 ) surf_v(3)%qcs
3455	CASE ( 'surf_v(3)%ncs' )
3456	IF ( ALLOCATED( surf_v(3)%ncs ) .AND. kk == 1 ) &
3457	READ ( 13 ) surf_v(3)%ncs
3458	CASE ( 'surf_v(3)%qrs' )
3459	IF ( ALLOCATED( surf_v(3)%qrs ) .AND. kk == 1 ) &
3460	READ ( 13 ) surf_v(3)%qrs
3461	CASE ( 'surf_v(3)%nrs' )
3462	IF ( ALLOCATED( surf_v(3)%nrs ) .AND. kk == 1 ) &
3463	READ ( 13 ) surf_v(3)%nrs
3464	CASE ( 'surf_v(3)%ol' )
3465	IF ( ALLOCATED( surf_v(3)%ol ) .AND. kk == 1 ) &
3466	READ ( 13 ) surf_v(3)%ol
3467	CASE ( 'surf_v(3)%rib' )
3468	IF ( ALLOCATED( surf_v(3)%rib ) .AND. kk == 1 ) &
3469	READ ( 13 ) surf_v(3)%rib
3470	CASE ( 'surf_v(3)%shf' )
3471	IF ( ALLOCATED( surf_v(3)%shf ) .AND. kk == 1 ) &
3472	READ ( 13 ) surf_v(3)%shf
3473	CASE ( 'surf_v(3)%qsws' )
3474	IF ( ALLOCATED( surf_v(3)%qsws ) .AND. kk == 1 ) &
3475	READ ( 13 ) surf_v(3)%qsws
3476	CASE ( 'surf_v(3)%ssws' )
3477	IF ( ALLOCATED( surf_v(3)%ssws ) .AND. kk == 1 ) &
3478	READ ( 13 ) surf_v(3)%ssws
3479	#if defined( __chem )
3480	CASE ( 'surf_v(3)%css' )
3481	IF ( ALLOCATED( surf_v(3)%css ) .AND. kk == 1 ) &
3482	READ ( 13 ) surf_v(3)%css
3483	CASE ( 'surf_v(3)%cssws' )
3484	IF ( ALLOCATED( surf_v(3)%cssws ) .AND. kk == 1 ) &
3485	READ ( 13 ) surf_v(3)%cssws
3486	#endif
3487	CASE ( 'surf_v(3)%qcsws' )
3488	IF ( ALLOCATED( surf_v(3)%qcsws ) .AND. kk == 1 ) &
3489	READ ( 13 ) surf_v(3)%qcsws
3490	CASE ( 'surf_v(3)%ncsws' )
3491	IF ( ALLOCATED( surf_v(3)%ncsws ) .AND. kk == 1 ) &
3492	READ ( 13 ) surf_v(3)%ncsws
3493	CASE ( 'surf_v(3)%qrsws' )
3494	IF ( ALLOCATED( surf_v(3)%qrsws ) .AND. kk == 1 ) &
3495	READ ( 13 ) surf_v(3)%qrsws
3496	CASE ( 'surf_v(3)%nrsws' )
3497	IF ( ALLOCATED( surf_v(3)%nrsws ) .AND. kk == 1 ) &
3498	READ ( 13 ) surf_v(3)%nrsws
3499	CASE ( 'surf_v(3)%sasws' )
3500	IF ( ALLOCATED( surf_v(3)%sasws ) .AND. kk == 1 ) &
3501	READ ( 13 ) surf_v(3)%sasws
3502	CASE ( 'surf_v(3)%mom_uv' )
3503	IF ( ALLOCATED( surf_v(3)%mom_flux_uv ) .AND. kk == 1 )&
3504	READ ( 13 ) surf_v(3)%mom_flux_uv
3505	CASE ( 'surf_v(3)%mom_w' )
3506	IF ( ALLOCATED( surf_v(3)%mom_flux_w ) .AND. kk == 1 ) &
3507	READ ( 13 ) surf_v(3)%mom_flux_w
3508	CASE ( 'surf_v(3)%mom_tke' )
3509	IF ( ALLOCATED( surf_v(3)%mom_flux_tke ) .AND. kk == 1 )&
3510	READ ( 13 ) surf_v(3)%mom_flux_tke
3511
3512	END SELECT
3513	!
3514	!-- Redistribute surface elements on its respective type.
3515	IF ( horizontal_surface ) THEN
3516	ic = nxlc
3517	DO i = nxlf, nxrf
3518	jc = nysc
3519	DO j = nysf, nynf
3520
3521	surf_match_def = surf_def_h(l)%end_index(jc,ic) >= &
3522	surf_def_h(l)%start_index(jc,ic)
3523	surf_match_lsm = surf_lsm_h%end_index(jc,ic) >= &
3524	surf_lsm_h%start_index(jc,ic)
3525	surf_match_usm = surf_usm_h%end_index(jc,ic) >= &
3526	surf_usm_h%start_index(jc,ic)
3527
3528	IF ( surf_match_def ) THEN
3529	mm = surf_def_h(l)%start_index(jc,ic)
3530	DO m = surf_h(l)%start_index(j,i), &
3531	surf_h(l)%end_index(j,i)
3532	CALL restore_surface_elements( surf_def_h(l), &
3533	mm, surf_h(l), m )
3534	mm = mm + 1
3535	ENDDO
3536	ENDIF
3537
3538	IF ( surf_match_lsm ) THEN
3539	mm = surf_lsm_h%start_index(jc,ic)
3540	DO m = surf_h(l)%start_index(j,i), &
3541	surf_h(l)%end_index(j,i)
3542	CALL restore_surface_elements( surf_lsm_h, &
3543	mm, surf_h(l), m )
3544	mm = mm + 1
3545	ENDDO
3546	ENDIF
3547
3548	IF ( surf_match_usm ) THEN
3549	mm = surf_usm_h%start_index(jc,ic)
3550	DO m = surf_h(l)%start_index(j,i), &
3551	surf_h(l)%end_index(j,i)
3552	CALL restore_surface_elements( surf_usm_h, &
3553	mm, surf_h(l), m )
3554	mm = mm + 1
3555	ENDDO
3556	ENDIF
3557
3558	jc = jc + 1
3559	ENDDO
3560	ic = ic + 1
3561	ENDDO
3562	ELSEIF ( vertical_surface ) THEN
3563	ic = nxlc
3564	DO i = nxlf, nxrf
3565	jc = nysc
3566	DO j = nysf, nynf
3567
3568	surf_match_def = surf_def_v(l)%end_index(jc,ic) >= &
3569	surf_def_v(l)%start_index(jc,ic)
3570	surf_match_lsm = surf_lsm_v(l)%end_index(jc,ic) >= &
3571	surf_lsm_v(l)%start_index(jc,ic)
3572	surf_match_usm = surf_usm_v(l)%end_index(jc,ic) >= &
3573	surf_usm_v(l)%start_index(jc,ic)
3574
3575
3576
3577	IF ( surf_match_def ) THEN
3578	mm = surf_def_v(l)%start_index(jc,ic)
3579	DO m = surf_v(l)%start_index(j,i), &
3580	surf_v(l)%end_index(j,i)
3581	CALL restore_surface_elements( surf_def_v(l), &
3582	mm, surf_v(l), m )
3583	mm = mm + 1
3584	ENDDO
3585	ENDIF
3586
3587	IF ( surf_match_lsm ) THEN
3588	mm = surf_lsm_v(l)%start_index(jc,ic)
3589	DO m = surf_v(l)%start_index(j,i), &
3590	surf_v(l)%end_index(j,i)
3591	CALL restore_surface_elements( surf_lsm_v(l), &
3592	mm, surf_v(l), m )
3593	mm = mm + 1
3594	ENDDO
3595	ENDIF
3596
3597	IF ( surf_match_usm ) THEN
3598	mm = surf_usm_v(l)%start_index(jc,ic)
3599	DO m = surf_v(l)%start_index(j,i), &
3600	surf_v(l)%end_index(j,i)
3601	CALL restore_surface_elements( surf_usm_v(l), &
3602	mm, surf_v(l), m )
3603	mm = mm + 1
3604	ENDDO
3605	ENDIF
3606
3607	jc = jc + 1
3608	ENDDO
3609	ic = ic + 1
3610	ENDDO
3611	ENDIF
3612
3613	ENDDO
3614
3615	READ ( 13 ) field_chr
3616
3617	ENDDO
3618
3619	ENDIF
3620
3621
3622	CONTAINS
3623	!------------------------------------------------------------------------------!
3624	! Description:
3625	! ------------
3626	!> Restores surfacle elements back on its respective type.
3627	!------------------------------------------------------------------------------!
3628	SUBROUTINE restore_surface_elements( surf_target, m_target, &
3629	surf_file, m_file )
3630
3631	IMPLICIT NONE
3632
3633	INTEGER(iwp) :: m_file !< respective surface-element index of current surface array
3634	INTEGER(iwp) :: m_target !< respecitve surface-element index of surface array on file
3635	INTEGER(iwp) :: lsp !< running index chemical species
3636
3637	TYPE( surf_type ) :: surf_target !< target surface type
3638	TYPE( surf_type ) :: surf_file !< surface type on file
3639
3640	IF ( INDEX( TRIM( field_chr ), '%us' ) /= 0 ) THEN
3641	IF ( ALLOCATED( surf_target%us ) .AND. &
3642	ALLOCATED( surf_file%us ) ) &
3643	surf_target%us(m_target) = surf_file%us(m_file)
3644	ENDIF
3645
3646	IF ( INDEX( TRIM( field_chr ), '%ol' ) /= 0 ) THEN
3647	IF ( ALLOCATED( surf_target%ol ) .AND. &
3648	ALLOCATED( surf_file%ol ) ) &
3649	surf_target%ol(m_target) = surf_file%ol(m_file)
3650	ENDIF
3651
3652	IF ( INDEX( TRIM( field_chr ), '%usws' ) /= 0 ) THEN
3653	IF ( ALLOCATED( surf_target%usws ) .AND. &
3654	ALLOCATED( surf_file%usws ) ) &
3655	surf_target%usws(m_target) = surf_file%usws(m_file)
3656	ENDIF
3657
3658	IF ( INDEX( TRIM( field_chr ), '%vsws' ) /= 0 ) THEN
3659	IF ( ALLOCATED( surf_target%vsws ) .AND. &
3660	ALLOCATED( surf_file%vsws ) ) &
3661	surf_target%vsws(m_target) = surf_file%vsws(m_file)
3662	ENDIF
3663
3664	IF ( INDEX( TRIM( field_chr ), '%ts' ) /= 0 ) THEN
3665	IF ( ALLOCATED( surf_target%ts ) .AND. &
3666	ALLOCATED( surf_file%ts ) ) &
3667	surf_target%ts(m_target) = surf_file%ts(m_file)
3668	ENDIF
3669
3670	IF ( INDEX( TRIM( field_chr ), '%shf' ) /= 0 ) THEN
3671	IF ( ALLOCATED( surf_target%shf ) .AND. &
3672	ALLOCATED( surf_file%shf ) ) &
3673	surf_target%shf(m_target) = surf_file%shf(m_file)
3674	ENDIF
3675
3676	IF ( INDEX( TRIM( field_chr ), '%qs' ) /= 0 ) THEN
3677	IF ( ALLOCATED( surf_target%qs ) .AND. &
3678	ALLOCATED( surf_file%qs ) ) &
3679	surf_target%qs(m_target) = surf_file%qs(m_file)
3680	ENDIF
3681
3682	IF ( INDEX( TRIM( field_chr ), '%qsws' ) /= 0 ) THEN
3683	IF ( ALLOCATED( surf_target%qsws ) .AND. &
3684	ALLOCATED( surf_file%qsws ) ) &
3685	surf_target%qsws(m_target) = surf_file%qsws(m_file)
3686	ENDIF
3687
3688	IF ( INDEX( TRIM( field_chr ), '%ss' ) /= 0 ) THEN
3689	IF ( ALLOCATED( surf_target%ss ) .AND. &
3690	ALLOCATED( surf_file%ss ) ) &
3691	surf_target%ss(m_target) = surf_file%ss(m_file)
3692	ENDIF
3693
3694	IF ( INDEX( TRIM( field_chr ), '%ssws' ) /= 0 ) THEN
3695	IF ( ALLOCATED( surf_target%ssws ) .AND. &
3696	ALLOCATED( surf_file%ssws ) ) &
3697	surf_target%ssws(m_target) = surf_file%ssws(m_file)
3698	ENDIF
3699
3700	#if defined( __chem )
3701	IF ( INDEX( TRIM( field_chr ), '%css' ) /= 0 ) THEN
3702	IF ( ALLOCATED( surf_target%css ) .AND. &
3703	ALLOCATED( surf_file%css ) ) THEN
3704	DO lsp = 1, nvar
3705	surf_target%css(lsp,m_target) = surf_file%css(lsp,m_file)
3706	ENDDO
3707	ENDIF
3708	ENDIF
3709	IF ( INDEX( TRIM( field_chr ), '%cssws' ) /= 0 ) THEN
3710	IF ( ALLOCATED( surf_target%cssws ) .AND. &
3711	ALLOCATED( surf_file%cssws ) ) THEN
3712	DO lsp = 1, nvar
3713	surf_target%cssws(lsp,m_target) = surf_file%cssws(lsp,m_file)
3714	ENDDO
3715	ENDIF
3716	ENDIF
3717	#endif
3718
3719	IF ( INDEX( TRIM( field_chr ), '%qcs' ) /= 0 ) THEN
3720	IF ( ALLOCATED( surf_target%qcs ) .AND. &
3721	ALLOCATED( surf_file%qcs ) ) &
3722	surf_target%qcs(m_target) = surf_file%qcs(m_file)
3723	ENDIF
3724
3725	IF ( INDEX( TRIM( field_chr ), '%qcsws' ) /= 0 ) THEN
3726	IF ( ALLOCATED( surf_target%qcsws ) .AND. &
3727	ALLOCATED( surf_file%qcsws ) ) &
3728	surf_target%qcsws(m_target) = surf_file%qcsws(m_file)
3729	ENDIF
3730
3731	IF ( INDEX( TRIM( field_chr ), '%ncs' ) /= 0 ) THEN
3732	IF ( ALLOCATED( surf_target%ncs ) .AND. &
3733	ALLOCATED( surf_file%ncs ) ) &
3734	surf_target%ncs(m_target) = surf_file%ncs(m_file)
3735	ENDIF
3736
3737	IF ( INDEX( TRIM( field_chr ), '%ncsws' ) /= 0 ) THEN
3738	IF ( ALLOCATED( surf_target%ncsws ) .AND. &
3739	ALLOCATED( surf_file%ncsws ) ) &
3740	surf_target%ncsws(m_target) = surf_file%ncsws(m_file)
3741	ENDIF
3742
3743	IF ( INDEX( TRIM( field_chr ), '%qrs' ) /= 0 ) THEN
3744	IF ( ALLOCATED( surf_target%qrs ) .AND. &
3745	ALLOCATED( surf_file%qrs ) ) &
3746	surf_target%qrs(m_target) = surf_file%qrs(m_file)
3747	ENDIF
3748
3749	IF ( INDEX( TRIM( field_chr ), '%qrsws' ) /= 0 ) THEN
3750	IF ( ALLOCATED( surf_target%qrsws ) .AND. &
3751	ALLOCATED( surf_file%qrsws ) ) &
3752	surf_target%qrsws(m_target) = surf_file%qrsws(m_file)
3753	ENDIF
3754
3755	IF ( INDEX( TRIM( field_chr ), '%nrs' ) /= 0 ) THEN
3756	IF ( ALLOCATED( surf_target%nrs ) .AND. &
3757	ALLOCATED( surf_file%nrs ) ) &
3758	surf_target%nrs(m_target) = surf_file%nrs(m_file)
3759	ENDIF
3760
3761	IF ( INDEX( TRIM( field_chr ), '%nrsws' ) /= 0 ) THEN
3762	IF ( ALLOCATED( surf_target%nrsws ) .AND. &
3763	ALLOCATED( surf_file%nrsws ) ) &
3764	surf_target%nrsws(m_target) = surf_file%nrsws(m_file)
3765	ENDIF
3766
3767	IF ( INDEX( TRIM( field_chr ), '%sasws' ) /= 0 ) THEN
3768	IF ( ALLOCATED( surf_target%sasws ) .AND. &
3769	ALLOCATED( surf_file%sasws ) ) &
3770	surf_target%sasws(m_target) = surf_file%sasws(m_file)
3771	ENDIF
3772
3773	IF ( INDEX( TRIM( field_chr ), '%mom_uv' ) /= 0 ) THEN
3774	IF ( ALLOCATED( surf_target%mom_flux_uv ) .AND. &
3775	ALLOCATED( surf_file%mom_flux_uv ) ) &
3776	surf_target%mom_flux_uv(m_target) = &
3777	surf_file%mom_flux_uv(m_file)
3778	ENDIF
3779
3780	IF ( INDEX( TRIM( field_chr ), '%mom_w' ) /= 0 ) THEN
3781	IF ( ALLOCATED( surf_target%mom_flux_w ) .AND. &
3782	ALLOCATED( surf_file%mom_flux_w ) ) &
3783	surf_target%mom_flux_w(m_target) = &
3784	surf_file%mom_flux_w(m_file)
3785	ENDIF
3786
3787	IF ( INDEX( TRIM( field_chr ), '%mom_tke' ) /= 0 ) THEN
3788	IF ( ALLOCATED( surf_target%mom_flux_tke ) .AND. &
3789	ALLOCATED( surf_file%mom_flux_tke ) ) &
3790	surf_target%mom_flux_tke(0:1,m_target) = &
3791	surf_file%mom_flux_tke(0:1,m_file)
3792	ENDIF
3793
3794	END SUBROUTINE restore_surface_elements
3795
3796	END SUBROUTINE surface_read_restart_data
3797
3798
3799	!------------------------------------------------------------------------------!
3800	! Description:
3801	! ------------
3802	!> Counts the number of surface elements with the same facing, required for
3803	!> reading and writing restart data.
3804	!------------------------------------------------------------------------------!
3805	SUBROUTINE surface_last_actions
3806
3807	IMPLICIT NONE
3808	!
3809	!-- Horizontal surfaces
3810	ns_h_on_file(0) = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns
3811	ns_h_on_file(1) = surf_def_h(1)%ns
3812	ns_h_on_file(2) = surf_def_h(2)%ns
3813	!
3814	!-- Vertical surfaces
3815	ns_v_on_file(0) = surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns
3816	ns_v_on_file(1) = surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns
3817	ns_v_on_file(2) = surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns
3818	ns_v_on_file(3) = surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns
3819
3820	END SUBROUTINE surface_last_actions
3821
3822	!------------------------------------------------------------------------------!
3823	! Description:
3824	! ------------
3825	!> Routine maps surface data read from file after restart - 1D arrays
3826	!------------------------------------------------------------------------------!
3827	SUBROUTINE surface_restore_elements_1d( surf_target, surf_file, &
3828	start_index_c, &
3829	start_index_on_file, &
3830	end_index_on_file, &
3831	nxlc, nysc, nxlf, nxrf, nysf, nynf,&
3832	nys_on_file, nyn_on_file, &
3833	nxl_on_file,nxr_on_file )
3834
3835	IMPLICIT NONE
3836
3837	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
3838	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
3839	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
3840	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
3841	INTEGER(iwp) :: m !< surface-element index on file
3842	INTEGER(iwp) :: mm !< surface-element index on current subdomain
3843	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
3844	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
3845	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
3846	INTEGER(iwp) :: nysc !< index of north boundary on current subdomain
3847	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
3848	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
3849
3850	INTEGER(iwp) :: nxl_on_file !< leftmost index on file
3851	INTEGER(iwp) :: nxr_on_file !< rightmost index on file
3852	INTEGER(iwp) :: nyn_on_file !< northmost index on file
3853	INTEGER(iwp) :: nys_on_file !< southmost index on file
3854
3855	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: start_index_c
3856	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
3857	start_index_on_file !< start index of surface elements on file
3858	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
3859	end_index_on_file !< end index of surface elements on file
3860
3861	REAL(wp), DIMENSION(:) :: surf_target !< target surface type
3862	REAL(wp), DIMENSION(:) :: surf_file !< surface type on file
3863
3864	ic = nxlc
3865	DO i = nxlf, nxrf
3866	jc = nysc
3867	DO j = nysf, nynf
3868
3869	mm = start_index_c(jc,ic)
3870	DO m = start_index_on_file(j,i), end_index_on_file(j,i)
3871	surf_target(mm) = surf_file(m)
3872	mm = mm + 1
3873	ENDDO
3874
3875	jc = jc + 1
3876	ENDDO
3877	ic = ic + 1
3878	ENDDO
3879
3880
3881	END SUBROUTINE surface_restore_elements_1d
3882
3883	!------------------------------------------------------------------------------!
3884	! Description:
3885	! ------------
3886	!> Routine maps surface data read from file after restart - 2D arrays
3887	!------------------------------------------------------------------------------!
3888	SUBROUTINE surface_restore_elements_2d( surf_target, surf_file, &
3889	start_index_c, &
3890	start_index_on_file, &
3891	end_index_on_file, &
3892	nxlc, nysc, nxlf, nxrf, nysf, nynf,&
3893	nys_on_file, nyn_on_file, &
3894	nxl_on_file,nxr_on_file )
3895
3896	IMPLICIT NONE
3897
3898	INTEGER(iwp) :: i !< running index along x-direction, refers to former domain size
3899	INTEGER(iwp) :: ic !< running index along x-direction, refers to current domain size
3900	INTEGER(iwp) :: j !< running index along y-direction, refers to former domain size
3901	INTEGER(iwp) :: jc !< running index along y-direction, refers to former domain size
3902	INTEGER(iwp) :: m !< surface-element index on file
3903	INTEGER(iwp) :: mm !< surface-element index on current subdomain
3904	INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain
3905	INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain
3906	INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain
3907	INTEGER(iwp) :: nysc !< index of north boundary on current subdomain
3908	INTEGER(iwp) :: nynf !< index of north boundary on former subdomain
3909	INTEGER(iwp) :: nysf !< index of south boundary on former subdomain
3910
3911	INTEGER(iwp) :: nxl_on_file !< leftmost index on file
3912	INTEGER(iwp) :: nxr_on_file !< rightmost index on file
3913	INTEGER(iwp) :: nyn_on_file !< northmost index on file
3914	INTEGER(iwp) :: nys_on_file !< southmost index on file
3915
3916	INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) :: start_index_c !< start index of surface type
3917	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
3918	start_index_on_file !< start index of surface elements on file
3919	INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: &
3920	end_index_on_file !< end index of surface elements on file
3921
3922	REAL(wp), DIMENSION(:,:) :: surf_target !< target surface type
3923	REAL(wp), DIMENSION(:,:) :: surf_file !< surface type on file
3924
3925	ic = nxlc
3926	DO i = nxlf, nxrf
3927	jc = nysc
3928	DO j = nysf, nynf
3929	mm = start_index_c(jc,ic)
3930	DO m = start_index_on_file(j,i), end_index_on_file(j,i)
3931	surf_target(:,mm) = surf_file(:,m)
3932	mm = mm + 1
3933	ENDDO
3934
3935	jc = jc + 1
3936	ENDDO
3937	ic = ic + 1
3938	ENDDO
3939
3940	END SUBROUTINE surface_restore_elements_2d
3941
3942
3943	END MODULE surface_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |