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

Last change on this file since 4850 was 4850, checked in by suehring, 8 months ago

Bugfix in indoor model: consider previous indoor temperature during restarts; Further bugfix in mpi-io restart mechanism for the waste-heat flux from buildings

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