Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |