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

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

Additional output for radiative fluxes added

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