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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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