Home

Context Navigation

source: palm/tags/release-6.0/SOURCE/surface_mod.f90 @ 4103

Last change on this file since 4103 was 3467, checked in by suehring, 6 years ago
Branch salsa @3446 re-integrated into trunk
Property svn:keywords set to `Id`
File size: 225.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |