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

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

Additional checks for surface_fractions and building_id; Remove redundant subroutine argument in surface_mod

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