Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |