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

Last change on this file since 2707 was 2707, checked in by suehring, 4 years ago

changes documented

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