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

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

Revise and bugfix surface-element mapping and 3D soil array in mpi-io restart branch

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