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

Last change on this file since 3655 was 3655, checked in by knoop, 2 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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