Home

Context Navigation

surface_mod.f90 @ 3241

Last change on this file since 3241 was 3241, checked in by raasch, 6 years ago
various changes to avoid compiler warnings (mainly removal of unused variables)
Property svn:keywords set to `Id`
File size: 218.9 KB

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