Home

Context Navigation

surface_mod.f90 @ 3593

Last change on this file since 3593 was 3593, checked in by kanani, 6 years ago
Bugfix for missing array allocation (biometeorology_mod), remove degree symbol (biometeorology_mod, indoor_model_mod, multi_agent_system_mod, surface_mod, wind_turbine_model_mod)
Property svn:keywords set to `Id`
File size: 229.2 KB

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