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

Last change on this file since 3593 was 3593, checked in by kanani, 3 years ago

Bugfix for missing array allocation (biometeorology_mod), remove degree symbol (biometeorology_mod, indoor_model_mod, multi_agent_system_mod, surface_mod, wind_turbine_model_mod)

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