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

Last change on this file since 4502 was 4502, checked in by schwenkel, 14 months ago

Implementation of ice microphysics

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