Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |