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

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

Further adjustments for surface structure

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