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

Last change on this file since 3560 was 3560, checked in by raasch, 3 years ago

bugfix to guarantee correct particle releases in case that the release interval is smaller than the model timestep

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