Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |