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

Last change on this file since 3241 was 3241, checked in by raasch, 3 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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