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

Last change on this file since 3176 was 3176, checked in by suehring, 3 years ago

Major bugfix in calculation of ol and ts at building roofs; bugfix in restart data for surface elements; in 2D data output, mask latent heat flux at urban-type surfaces

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