Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |