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

Last change on this file since 3744 was 3744, checked in by suehring, 2 years ago

Coupling of indoor model to atmosphere; output of indoor temperatures and waste heat; enable restarts with indoor model; bugfix plant transpiration; bugfix - missing calculation of 10cm temperature

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