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

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

document last commit

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