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

Last change on this file since 2920 was 2920, checked in by kanani, 3 years ago

Optimize SVF calculation, clean-up, bugfixes

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