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

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

Bugfix in large-scale forcing

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