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

Last change on this file since 4168 was 4168, checked in by suehring, 22 months ago

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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