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

Last change on this file since 4159 was 4159, checked in by suehring, 23 months ago

Revision of topography processing to have a consistent treatment of 2D and 3D buildings (init_grid, surface_mod); Bugfix in indoor model in case of non grid-resolved buildings

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