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

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

Fix accidently commented subroutine

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