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

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago

restart data handling with MPI-IO added, first part

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