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

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

Bugfix in reading restart data of vertical surface elements

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