source: palm/trunk/SOURCE/surface_mod.f90

Last change on this file was 4896, checked in by raasch, 4 months ago

small re-formatting to follow the coding standard, typo in file appendix removed, more meaningful variable names assigned, redundant code removed

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