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

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

restart mechanism for surface elements commented, some formatting

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