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

Last change on this file since 4354 was 4354, checked in by suehring, 3 years ago

surface_mod: bugfix in message call and define appropriate error number

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