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

Last change on this file since 4331 was 4331, checked in by suehring, 4 years ago

New diagnostic output for 10-m wind speed; Diagnostic output of 2-m potential temperature moved to diagnostic output

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