Home

Context Navigation

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

Last change on this file since 3425 was 3418, checked in by kanani, 6 years ago
Add green facades, update building data base, fix for thin walls in spinup
Property svn:keywords set to `Id`
File size: 222.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |