Home

Context Navigation

surface_mod.f90 @ 3560

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