Home

Context Navigation

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

Last change on this file since 3590 was 3572, checked in by suehring, 6 years ago
Additional output for radiative fluxes added
Property svn:keywords set to `Id`
File size: 229.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |