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

Last change on this file since 4829 was 4829, checked in by suehring, 9 months ago

Bugfix in creating error-message string

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