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

Last change on this file since 4559 was 4559, checked in by raasch, 9 months ago

files re-formatted to follow the PALM coding standard

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