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

Last change on this file since 4594 was 4594, checked in by suehring, 10 months ago

Bugfix, add acc directives for scalar-roughness length; Include k index in OMP PRIVATE statements

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