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

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

Bugfixes concerning top fluxes and TKE production

  • Property svn:keywords set to Id
File size: 161.8 KB
Line 
1!> @file surface_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
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! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: surface_mod.f90 2478 2017-09-18 13:37:24Z suehring $
27! Bugfixes in initializing model top
28!
29! 2378 2017-08-31 13:57:27Z suehring
30! Bugfix in write restart data
31!
32! 2339 2017-08-07 13:55:26Z gronemeier
33! corrected timestamp in header
34!
35! 2338 2017-08-07 12:15:38Z gronemeier
36! Modularize 1D model
37!
38! 2318 2017-07-20 17:27:44Z suehring
39! New function to obtain topography top index.
40!
41! 2317 2017-07-20 17:27:19Z suehring
42! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
43! includes two more prognostic equations for cloud drop concentration (nc) 
44! and cloud water content (qc).
45!
46! 2270 2017-06-09 12:18:47Z maronga
47! Parameters removed/added due to changes in the LSM
48!
49! 2269 2017-06-09 11:57:32Z suehring
50! Formatting and description adjustments
51!
52! 2256 2017-06-07 13:58:08Z suehring
53! Enable heating at downward-facing surfaces
54!
55! 2233 2017-05-30 18:08:54Z suehring
56! Initial revision
57!
58!
59! Description:
60! ------------
61!> Surface module defines derived data structures to treat surface-
62!> bounded grid cells. Three different types of surfaces are defined:
63!> default surfaces, natural surfaces, and urban surfaces. The module
64!> encompasses the allocation and initialization of surface arrays, and handles
65!> reading and writing restart data.
66!> In addition, a further derived data structure is defined, in order to set
67!> boundary conditions at surfaces. 
68!------------------------------------------------------------------------------!
69 MODULE surface_mod
70
71    USE arrays_3d,                                                             &
72        ONLY:  zu, zw, heatflux_input_conversion, waterflux_input_conversion,  &
73               momentumflux_input_conversion
74
75    USE control_parameters               
76
77    USE indices,                                                               &
78        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, wall_flags_0
79
80    USE grid_variables,                                                        &
81        ONLY:  dx, dy
82
83    USE kinds
84
85    USE model_1d_mod,                                                          &
86        ONLY:  rif1d, us1d, usws1d, vsws1d
87
88
89    IMPLICIT NONE
90
91!
92!-- Data type used to identify grid-points where horizontal boundary conditions
93!-- are applied
94    TYPE bc_type
95
96       INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
97
98       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid 
99       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid   
100       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid   
101
102       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
103       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< end index within surface data type for given (j,i) 
104
105    END TYPE bc_type
106!
107!-- Data type used to identify and treat surface-bounded grid points 
108    TYPE surf_type
109
110       INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
111       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid 
112       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid   
113       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid       
114
115       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  facing  !< Bit indicating surface orientation
116     
117       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i)
118       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< End index within surface data type for given (j,i) 
119
120       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_mo      !< surface-layer height
121       REAL(wp), DIMENSION(:), ALLOCATABLE ::  uvw_abs   !< absolute surface-parallel velocity
122       REAL(wp), DIMENSION(:), ALLOCATABLE ::  us        !< friction velocity
123       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ts        !< scaling parameter temerature
124       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qs        !< scaling parameter humidity
125       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ss        !< scaling parameter passive scalar
126       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcs       !< scaling parameter qc
127       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncs       !< scaling parameter nc
128       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrs       !< scaling parameter qr
129       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrs       !< scaling parameter nr
130
131       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ol        !< Obukhov length
132       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rib       !< Richardson bulk number
133
134       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0        !< roughness length for momentum
135       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0h       !< roughness length for heat
136       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0q       !< roughness length for humidity
137
138       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< Specific humidity at first grid level (required for cloud_physics = .T.)
139       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< Potential temperature at first grid level (required for cloud_physics = .T.)
140!
141!--    Define arrays for surface fluxes
142       REAL(wp), DIMENSION(:), ALLOCATABLE ::  usws      !< vertical momentum flux for u-component at horizontal surfaces
143       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vsws      !< vertical momentum flux for v-component at horizontal surfaces
144
145       REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf       !< surface flux sensible heat
146       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat
147       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ssws      !< surface flux passive scalar
148       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcsws     !< surface flux qc
149       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncsws     !< surface flux nc
150       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrsws     !< surface flux qr
151       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws     !< surface flux nr
152       REAL(wp), DIMENSION(:), ALLOCATABLE ::  sasws     !< surface flux salinity
153!
154!--    Required for horizontal walls in production_e
155       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_0       !< virtual velocity component (see production_e_init for further explanation)
156       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_0       !< virtual velocity component (see production_e_init for further explanation)
157
158       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mom_flux_uv  !< momentum flux usvs and vsus at vertical surfaces (used in diffusion_u and diffusion_v)
159       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mom_flux_w   !< momentum flux wsus and wsvs at vertical surfaces (used in diffusion_w)
160       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  mom_flux_tke !< momentum flux usvs, vsus, wsus, wsvs at vertical surfaces at grid center (used in production_e)
161!
162!--    Variables required for LSM as well as for USM
163       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface        !< skin-surface temperature
164       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation
165       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil (W/m/K)
166       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)   
167       
168       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  building_surface    !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
169       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  pavement_surface    !< flag parameter for pavements
170       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces
171       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface     !< flag parameter for natural land surfaces
172
173       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_liq               !< liquid water coverage (of vegetated area)
174       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_veg               !< vegetation coverage
175       REAL(wp), DIMENSION(:), ALLOCATABLE ::  f_sw_in             !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
176       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ghf              !< ground heat flux
177       REAL(wp), DIMENSION(:), ALLOCATABLE ::  g_d                 !< coefficient for dependence of r_canopy on water vapour pressure deficit
178       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lai                 !< leaf area index
179       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_u    !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
180       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_s    !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
181       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq         !< surface flux of latent heat (liquid water portion)
182       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_soil        !< surface flux of latent heat (soil portion)
183       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg         !< surface flux of latent heat (vegetation portion)
184       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance
185       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy            !< canopy resistance
186       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_soil              !< soil resistance
187       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_soil_min          !< minimum soil resistance
188       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_s                 !< total surface resistance (combination of r_soil and r_canopy)
189       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy_min        !< minimum canopy (stomatal) resistance
190
191       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg          !< coef. of Van Genuchten
192       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_w          !< hydraulic diffusivity of soil (?)
193       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w           !< hydraulic conductivity of soil (W/m/K)
194       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w_sat       !< hydraulic conductivity at saturation
195       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  l_vg              !< coef. of Van Genuchten
196       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_fc              !< soil moisture at field capacity (m3/m3)
197       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_res             !< residual soil moisture
198       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_sat             !< saturation soil moisture (m3/m3)
199       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_wilt            !< soil moisture at permanent wilting point (m3/m3)
200       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg              !< coef. Van Genuchten 
201       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  rho_c_def         !< default volumetric heat capacity of the (soil) layer (J/m3/K)
202       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total       !< volumetric heat capacity of the actual soil matrix (J/m3/K)
203       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  root_fr           !< root fraction within the soil layers
204!
205!--    Urban surface variables
206       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  surface_types   !< array of types of wall parameters
207
208       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  isroof_surf         !< flag indication roof surfaces
209
210       REAL(wp), DIMENSION(:), ALLOCATABLE ::  albedo_surf         !< albedo of the surface
211       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin (J/m2/K)
212       REAL(wp), DIMENSION(:), ALLOCATABLE ::  emiss_surf          !< emissivity of the wall surface
213       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf         !< heat conductivity between air and surface (W/m2/K)
214       REAL(wp), DIMENSION(:), ALLOCATABLE ::  roughness_wall      !< roughness relative to concrete
215       REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_wall      !< thickness of the wall, roof and soil layers
216
217       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsl           !< reflected shortwave radiation for local surface in i-th reflection
218       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutll           !< reflected + emitted longwave radiation for local surface in i-th reflection
219       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf              !< total radiation flux incoming to minus outgoing from local surface
220
221       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_m        !< surface temperature tendency (K)
222       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf                !< kinematic wall heat flux of sensible heat (actually no longer needed)
223       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf_eb             !< wall heat flux of sensible heat in wall normal direction
224
225
226       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb             !< wall ground heat flux
227
228       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_in_sw           !< incoming shortwave radiation
229       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_out_sw          !< emitted shortwave radiation
230       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_in_lw           !< incoming longwave radiation
231       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_out_lw          !< emitted longwave radiation
232
233       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw            !< shortwave radiation falling to local surface including radiation from reflections
234       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsw           !< total shortwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
235       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw            !< longwave radiation falling to local surface including radiation from reflections
236       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw           !< total longwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
237
238
239       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_wall        !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)
240       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_wall           !< wall grid spacing (center-center)
241       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_wall          !< 1/dz_wall
242       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_wall_stag      !< wall grid spacing (edge-edge)
243       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_wall_stag     !< 1/dz_wall_stag
244       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m         !< t_wall prognostic array
245       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw                !< wall layer depths (m)
246
247
248!-- arrays for time averages
249       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_net_av       !< average of rad_net_l
250       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
251       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
252       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
253       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
254       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
255       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
256       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
257       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
258       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
259       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
260       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
261       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface 
262       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_av       !< average of wghf_eb
263       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf_eb_av       !< average of wshf_eb
264       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_av        !< average of wall surface temperature (K)
265
266       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_wall_av      !< Average of t_wall
267
268    END TYPE surf_type
269
270    TYPE (bc_type), DIMENSION(0:1)           ::  bc_h        !< boundary condition data type, horizontal upward- and downward facing surfaces
271
272    TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
273    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_def_v  !< vertical default surfaces (North, South, West, East)
274    TYPE (surf_type)                , TARGET ::  surf_lsm_h  !< horizontal natural land surfaces, so far only upward-facing
275    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_lsm_v  !< vertical land surfaces (North, South, West, East)
276    TYPE (surf_type)                , TARGET ::  surf_usm_h  !< horizontal urban surfaces, so far only upward-facing
277    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, West, East)
278
279    INTEGER(iwp) ::  ns_h_on_file(0:2)                       !< total number of horizontal surfaces with the same facing, required for writing restart data
280    INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data
281
282
283    SAVE
284
285    PRIVATE
286
287    INTERFACE get_topography_top_index
288       MODULE PROCEDURE get_topography_top_index
289    END  INTERFACE get_topography_top_index
290
291    INTERFACE init_bc
292       MODULE PROCEDURE init_bc
293    END INTERFACE init_bc
294
295    INTERFACE init_surfaces
296       MODULE PROCEDURE init_surfaces
297    END INTERFACE init_surfaces
298
299    INTERFACE init_surface_arrays
300       MODULE PROCEDURE init_surface_arrays
301    END INTERFACE init_surface_arrays
302
303    INTERFACE surface_read_restart_data
304       MODULE PROCEDURE surface_read_restart_data
305    END INTERFACE surface_read_restart_data
306
307    INTERFACE surface_write_restart_data
308       MODULE PROCEDURE surface_write_restart_data
309    END INTERFACE surface_write_restart_data
310
311    INTERFACE surface_last_actions
312       MODULE PROCEDURE surface_last_actions
313    END INTERFACE surface_last_actions
314
315!
316!-- Public variables
317    PUBLIC bc_h, ns_h_on_file, ns_v_on_file, surf_def_h, surf_def_v,           &
318           surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type
319!
320!-- Public subroutines and functions
321    PUBLIC get_topography_top_index, init_bc, init_surfaces,                   &
322           init_surface_arrays, surface_read_restart_data,                     &
323           surface_write_restart_data, surface_last_actions
324
325
326 CONTAINS
327
328!------------------------------------------------------------------------------!
329! Description:
330! ------------
331!> Initialize data type for setting boundary conditions at horizontal surfaces.
332!------------------------------------------------------------------------------!
333    SUBROUTINE init_bc
334
335       IMPLICIT NONE
336
337       INTEGER(iwp) ::  i         !<
338       INTEGER(iwp) ::  j         !<
339       INTEGER(iwp) ::  k         !<
340
341       INTEGER(iwp), DIMENSION(0:1) ::  num_h         !<
342       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !<
343       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !<
344
345!
346!--    First of all, count the number of upward- and downward-facing surfaces
347       num_h = 0
348       DO  i = nxlg, nxrg
349          DO  j = nysg, nyng
350             DO  k = nzb+1, nzt
351!
352!--             Check if current gridpoint belongs to the atmosphere
353                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
354!
355!--                Upward-facing
356                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )              &
357                      num_h(0) = num_h(0) + 1
358!
359!--                Downward-facing
360                   IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )              &
361                      num_h(1) = num_h(1) + 1
362                ENDIF
363             ENDDO
364          ENDDO
365       ENDDO
366!
367!--    Save the number of surface elements
368       bc_h(0)%ns = num_h(0)
369       bc_h(1)%ns = num_h(1)
370!
371!--    ALLOCATE data type variables
372!--    Upward facing
373       ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
374       ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
375       ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
376       ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
377       ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg)   )
378       bc_h(0)%start_index = 1
379       bc_h(0)%end_index   = 0
380!
381!--    Downward facing
382       ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
383       ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
384       ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
385       ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
386       ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg)   )
387       bc_h(1)%start_index = 1
388       bc_h(1)%end_index   = 0
389!
390!--    Store the respective indices on data type
391       num_h(0:1)         = 1
392       start_index_h(0:1) = 1
393       DO  i = nxlg, nxrg
394          DO  j = nysg, nyng
395
396             num_h_kji(0:1) = 0
397             DO  k = nzb+1, nzt
398!
399!--             Check if current gridpoint belongs to the atmosphere
400                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
401!
402!--                Upward-facing
403                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN
404                      bc_h(0)%i(num_h(0)) = i
405                      bc_h(0)%j(num_h(0)) = j
406                      bc_h(0)%k(num_h(0)) = k
407                      num_h_kji(0)        = num_h_kji(0) + 1
408                      num_h(0)            = num_h(0) + 1
409                   ENDIF
410!
411!--                Downward-facing
412                   IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )  THEN
413                      bc_h(1)%i(num_h(1)) = i
414                      bc_h(1)%j(num_h(1)) = j
415                      bc_h(1)%k(num_h(1)) = k
416                      num_h_kji(1)        = num_h_kji(1) + 1
417                      num_h(1)            = num_h(1) + 1
418                   ENDIF
419                ENDIF
420             ENDDO
421             bc_h(0)%start_index(j,i) = start_index_h(0)
422             bc_h(0)%end_index(j,i)   = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
423             start_index_h(0)         = bc_h(0)%end_index(j,i) + 1
424
425             bc_h(1)%start_index(j,i) = start_index_h(1)
426             bc_h(1)%end_index(j,i)   = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
427             start_index_h(1)         = bc_h(1)%end_index(j,i) + 1
428          ENDDO
429       ENDDO
430
431
432    END SUBROUTINE init_bc
433
434
435!------------------------------------------------------------------------------!
436! Description:
437! ------------
438!> Initialize horizontal and vertical surfaces. Counts the number of default-,
439!> natural and urban surfaces and allocates memory, respectively.
440!------------------------------------------------------------------------------!
441    SUBROUTINE init_surface_arrays
442
443       IMPLICIT NONE
444
445       INTEGER(iwp)                 ::  i         !< running index x-direction
446       INTEGER(iwp)                 ::  j         !< running index y-direction
447       INTEGER(iwp)                 ::  k         !< running index z-direction
448       INTEGER(iwp)                 ::  l         !< index variable for surface facing
449       INTEGER(iwp)                 ::  num_lsm_h !< number of horizontally-aligned natural surfaces
450       INTEGER(iwp)                 ::  num_usm_h !< number of horizontally-aligned urban surfaces
451
452       INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces
453       INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces
454       INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces
455       INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces
456
457
458       num_def_h = 0
459       num_def_v = 0
460       num_lsm_h = 0
461       num_lsm_v = 0
462       num_usm_h = 0
463       num_usm_v = 0
464!
465!--    Count number of horizontal surfaces on local domain
466       DO  i = nxl, nxr
467          DO  j = nys, nyn
468             DO  k = nzb+1, nzt
469!
470!--             Check if current gridpoint belongs to the atmosphere
471                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
472!
473!--                Check if grid point adjoins to any upward-facing horizontal
474!--                surface, e.g. the Earth surface, plane roofs, or ceilings.
475                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN
476!
477!--                   Land-surface type
478                      IF ( land_surface )  THEN
479                         num_lsm_h    = num_lsm_h    + 1 
480!
481!--                   Urban surface tpye
482                      ELSEIF ( urban_surface )  THEN
483                         num_usm_h    = num_usm_h    + 1 
484!
485!--                   Default-surface type
486                      ELSE
487                         num_def_h(0) = num_def_h(0) + 1 
488                      ENDIF
489
490                   ENDIF
491!
492!--                Check for top-fluxes
493                   IF ( k == nzt  .AND.  use_top_fluxes )  THEN
494                      num_def_h(2) = num_def_h(2) + 1
495!
496!--                Check for any other downward-facing surface. So far only for
497!--                default surface type.
498                   ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )  THEN
499                      num_def_h(1) = num_def_h(1) + 1
500                   ENDIF
501
502                ENDIF
503             ENDDO
504          ENDDO
505       ENDDO
506!
507!--    Count number of vertical surfaces on local domain
508       DO  i = nxl, nxr
509          DO  j = nys, nyn
510             DO  k = nzb+1, nzt
511                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
512!
513!--                Northward-facing
514                   IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) )  THEN
515                      IF ( urban_surface )  THEN
516                         num_usm_v(0) = num_usm_v(0) + 1 
517                      ELSEIF ( land_surface )  THEN
518                         num_lsm_v(0) = num_lsm_v(0) + 1 
519                      ELSE
520                         num_def_v(0) = num_def_v(0) + 1 
521                      ENDIF
522                   ENDIF
523!
524!--                Southward-facing
525                   IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )  THEN
526                      IF ( urban_surface )  THEN
527                         num_usm_v(1) = num_usm_v(1) + 1 
528                      ELSEIF ( land_surface )  THEN
529                         num_lsm_v(1) = num_lsm_v(1) + 1 
530                      ELSE
531                         num_def_v(1) = num_def_v(1) + 1 
532                      ENDIF
533                   ENDIF
534!
535!--                Eastward-facing
536                   IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) )  THEN
537                      IF ( urban_surface )  THEN
538                         num_usm_v(2) = num_usm_v(2) + 1 
539                      ELSEIF ( land_surface )  THEN
540                         num_lsm_v(2) = num_lsm_v(2) + 1 
541                      ELSE
542                         num_def_v(2) = num_def_v(2) + 1 
543                      ENDIF
544                   ENDIF
545!
546!--                Westward-facing
547                   IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )  THEN
548                      IF ( urban_surface )  THEN
549                         num_usm_v(3) = num_usm_v(3) + 1
550                      ELSEIF ( land_surface )  THEN
551                         num_lsm_v(3) = num_lsm_v(3) + 1 
552                      ELSE
553                         num_def_v(3) = num_def_v(3) + 1 
554                      ENDIF
555                   ENDIF
556                ENDIF
557             ENDDO
558          ENDDO
559       ENDDO
560
561!
562!--    Store number of surfaces per core.
563!--    Horizontal surface, default type, upward facing
564       surf_def_h(0)%ns = num_def_h(0)
565!
566!--    Horizontal surface, default type, downward facing
567       surf_def_h(1)%ns = num_def_h(1)
568!
569!--    Horizontal surface, default type, top downward facing
570       surf_def_h(2)%ns = num_def_h(2)
571!
572!--    Horizontal surface, natural type, so far only upward-facing
573       surf_lsm_h%ns    = num_lsm_h 
574!
575!--    Horizontal surface, urban type, so far only upward-facing
576       surf_usm_h%ns    = num_usm_h   
577!
578!--    Vertical surface, default type, northward facing
579       surf_def_v(0)%ns = num_def_v(0)
580!
581!--    Vertical surface, default type, southward facing
582       surf_def_v(1)%ns = num_def_v(1)
583!
584!--    Vertical surface, default type, eastward facing
585       surf_def_v(2)%ns = num_def_v(2)
586!
587!--    Vertical surface, default type, westward facing
588       surf_def_v(3)%ns = num_def_v(3)
589!
590!--    Vertical surface, natural type, northward facing
591       surf_lsm_v(0)%ns = num_lsm_v(0)
592!
593!--    Vertical surface, natural type, southward facing
594       surf_lsm_v(1)%ns = num_lsm_v(1)
595!
596!--    Vertical surface, natural type, eastward facing
597       surf_lsm_v(2)%ns = num_lsm_v(2)
598!
599!--    Vertical surface, natural type, westward facing
600       surf_lsm_v(3)%ns = num_lsm_v(3)
601!
602!--    Vertical surface, urban type, northward facing
603       surf_usm_v(0)%ns = num_usm_v(0)
604!
605!--    Vertical surface, urban type, southward facing
606       surf_usm_v(1)%ns = num_usm_v(1)
607!
608!--    Vertical surface, urban type, eastward facing
609       surf_usm_v(2)%ns = num_usm_v(2)
610!
611!--    Vertical surface, urban type, westward facing
612       surf_usm_v(3)%ns = num_usm_v(3)
613!
614!--    Allocate required attributes for horizontal surfaces - default type.
615!--    Upward-facing (l=0) and downward-facing (l=1).
616       DO  l = 0, 1
617          CALL allocate_surface_attributes_h ( surf_def_h(l), nys, nyn, nxl, nxr )
618       ENDDO
619!
620!--    Allocate required attributes for model top
621       CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
622!
623!--    Allocate required attributes for horizontal surfaces - natural type.
624       CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
625!
626!--    Allocate required attributes for horizontal surfaces - urban type.
627       CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
628
629!
630!--    Allocate required attributes for vertical surfaces.
631!--    Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2)
632!--    and westward-facing (l=3).
633!--    Default type.
634       DO  l = 0, 3
635          CALL allocate_surface_attributes_v ( surf_def_v(l), .FALSE.,         &
636                                               nys, nyn, nxl, nxr )
637       ENDDO
638!
639!--    Natural type
640       DO  l = 0, 3
641          CALL allocate_surface_attributes_v ( surf_lsm_v(l), .TRUE.,          &
642                                               nys, nyn, nxl, nxr )
643       ENDDO
644!
645!--    Urban type
646       DO  l = 0, 3
647          CALL allocate_surface_attributes_v ( surf_usm_v(l), .FALSE.,         &
648                                               nys, nyn, nxl, nxr )
649       ENDDO
650
651    END SUBROUTINE init_surface_arrays
652
653!------------------------------------------------------------------------------!
654! Description:
655! ------------
656!> Allocating memory for upward and downward-facing horizontal surface types,
657!> except for top fluxes.
658!------------------------------------------------------------------------------!
659    SUBROUTINE allocate_surface_attributes_h( surfaces,                        &
660                                              nys_l, nyn_l, nxl_l, nxr_l )
661
662       IMPLICIT NONE
663
664       INTEGER(iwp) ::  nyn_l  !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
665       INTEGER(iwp) ::  nys_l  !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
666       INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
667       INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
668
669       TYPE(surf_type) ::  surfaces  !< respective surface type
670
671!
672!--    Allocate arrays for start and end index of horizontal surface type
673!--    for each (j,i)-grid point. This is required e.g. in diffion_x, which is
674!--    called for each (j,i). In order to find the location where the
675!--    respective flux is store within the surface-type, start- and end-
676!--    index are stored for each (j,i). For example, each (j,i) can have
677!--    several entries where fluxes for horizontal surfaces might be stored,
678!--    e.g. for overhanging structures where several upward-facing surfaces
679!--    might exist for given (j,i).
680!--    If no surface of respective type exist at current (j,i), set indicies
681!--    such that loop in diffusion routines will not be entered.
682       ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
683       ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l)   )
684       surfaces%start_index = 0
685       surfaces%end_index   = -1
686!
687!--    Indices to locate surface element
688       ALLOCATE ( surfaces%i(1:surfaces%ns)  )
689       ALLOCATE ( surfaces%j(1:surfaces%ns)  )
690       ALLOCATE ( surfaces%k(1:surfaces%ns)  )
691!
692!--    Surface-layer height
693       ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
694!
695!--    Surface orientation
696       ALLOCATE ( surfaces%facing(1:surfaces%ns) )
697!
698!--    Surface-parallel wind velocity
699       ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
700!      ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
701!
702!--    Roughness
703       ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
704       ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
705       ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
706!
707!--    Friction velocity
708       ALLOCATE ( surfaces%us(1:surfaces%ns) )
709!
710!--    Stability parameter
711       ALLOCATE ( surfaces%ol(1:surfaces%ns) )
712!
713!--    Bulk Richardson number
714       ALLOCATE ( surfaces%rib(1:surfaces%ns) )
715!
716!--    Vertical momentum fluxes of u and v
717       ALLOCATE ( surfaces%usws(1:surfaces%ns) ) 
718       ALLOCATE ( surfaces%vsws(1:surfaces%ns) ) 
719!
720!--    Required in production_e
721       IF ( .NOT. constant_diffusion )  THEN   
722          ALLOCATE ( surfaces%u_0(1:surfaces%ns) ) 
723          ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
724       ENDIF
725!
726!--    Characteristic temperature and surface flux of sensible heat
727       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )   
728       ALLOCATE ( surfaces%shf(1:surfaces%ns) )   
729!
730!--    Characteristic humidity and surface flux of latent heat
731       IF ( humidity )  THEN
732          ALLOCATE ( surfaces%qs(1:surfaces%ns)   ) 
733          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )     
734       ENDIF
735!
736!--    Characteristic scalar and surface flux of scalar
737       IF ( passive_scalar )  THEN
738          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
739          ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
740       ENDIF
741!
742!--    When cloud physics is used, arrays for storing potential temperature and
743!--    specific humidity at first grid level are required.
744       IF ( cloud_physics )  THEN
745          ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
746          ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
747       ENDIF
748!
749!--       
750       IF ( cloud_physics .AND. microphysics_morrison)  THEN
751          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
752          ALLOCATE ( surfaces%ncs(1:surfaces%ns)   )
753          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
754          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
755       ENDIF
756!
757!--       
758       IF ( cloud_physics .AND. microphysics_seifert)  THEN
759          ALLOCATE ( surfaces%qrs(1:surfaces%ns)   )
760          ALLOCATE ( surfaces%nrs(1:surfaces%ns)   )
761          ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
762          ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
763       ENDIF
764!
765!--    Salinity surface flux
766       IF ( ocean )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
767
768    END SUBROUTINE allocate_surface_attributes_h
769
770
771!------------------------------------------------------------------------------!
772! Description:
773! ------------
774!> Allocating memory for model-top fluxes 
775!------------------------------------------------------------------------------!
776    SUBROUTINE allocate_surface_attributes_h_top( surfaces,                    &
777                                                  nys_l, nyn_l, nxl_l, nxr_l )
778
779       IMPLICIT NONE
780
781       INTEGER(iwp) ::  nyn_l  !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
782       INTEGER(iwp) ::  nys_l  !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
783       INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
784       INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
785
786       TYPE(surf_type) ::  surfaces !< respective surface type
787
788       ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
789       ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l)   )
790       surfaces%start_index = 0
791       surfaces%end_index   = -1
792!
793!--    Indices to locate surface (model-top) element
794       ALLOCATE ( surfaces%i(1:surfaces%ns)  )
795       ALLOCATE ( surfaces%j(1:surfaces%ns)  )
796       ALLOCATE ( surfaces%k(1:surfaces%ns)  )
797!
798!--    Vertical momentum fluxes of u and v
799       ALLOCATE ( surfaces%usws(1:surfaces%ns) ) 
800       ALLOCATE ( surfaces%vsws(1:surfaces%ns) ) 
801!
802!--    Sensible heat flux
803       ALLOCATE ( surfaces%shf(1:surfaces%ns) )   
804!
805!--    Latent heat flux
806       IF ( humidity )  THEN
807          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )     
808       ENDIF
809!
810!--    Scalar flux
811       IF ( passive_scalar )  THEN
812          ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
813       ENDIF
814!
815!--       
816       IF ( cloud_physics .AND. microphysics_morrison)  THEN
817          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
818          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
819       ENDIF
820!
821!--       
822       IF ( cloud_physics .AND. microphysics_seifert)  THEN
823          ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
824          ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
825       ENDIF
826!
827!--    Salinity flux
828       IF ( ocean )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
829
830    END SUBROUTINE allocate_surface_attributes_h_top
831
832!------------------------------------------------------------------------------!
833! Description:
834! ------------
835!> Allocating memory for vertical surface types.
836!------------------------------------------------------------------------------!
837    SUBROUTINE allocate_surface_attributes_v( surfaces, lsm,                   &
838                                              nys_l, nyn_l, nxl_l, nxr_l )
839
840       IMPLICIT NONE
841
842       INTEGER(iwp) ::  nyn_l  !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
843       INTEGER(iwp) ::  nys_l  !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
844       INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
845       INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
846
847       LOGICAL         ::  lsm      !< flag indicating data type of natural land surface
848
849       TYPE(surf_type) ::  surfaces !< respective surface type
850
851!
852!--    Allocate arrays for start and end index of vertical surface type
853!--    for each (j,i)-grid point. This is required in diffion_x, which is
854!--    called for each (j,i). In order to find the location where the
855!--    respective flux is store within the surface-type, start- and end-
856!--    index are stored for each (j,i). For example, each (j,i) can have
857!--    several entries where fluxes for vertical surfaces might be stored. 
858!--    In the flat case, where no vertical walls exit, set indicies such
859!--    that loop in diffusion routines will not be entered.
860       ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
861       ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l)   )
862       surfaces%start_index = 0
863       surfaces%end_index   = -1
864!
865!--    Indices to locate surface element.
866       ALLOCATE ( surfaces%i(1:surfaces%ns) )
867       ALLOCATE ( surfaces%j(1:surfaces%ns) )
868       ALLOCATE ( surfaces%k(1:surfaces%ns) )
869!
870!--    Surface-layer height
871       ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
872!
873!--    Surface orientation
874       ALLOCATE ( surfaces%facing(1:surfaces%ns) )
875!
876!--    Surface parallel wind velocity
877       ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
878!
879!--    Roughness
880       ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
881       ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
882       ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
883
884!
885!--    Friction velocity
886       ALLOCATE ( surfaces%us(1:surfaces%ns) )
887!
888!--    Allocate Obukhov length and bulk Richardson number. Only required
889!--    for natural land surfaces
890       IF ( lsm )  THEN
891          ALLOCATE( surfaces%ol(1:surfaces%ns)  ) 
892          ALLOCATE( surfaces%rib(1:surfaces%ns) ) 
893       ENDIF
894!
895!--    Allocate arrays for surface momentum fluxes for u and v. For u at north-
896!--    and south-facing surfaces, for v at east- and west-facing surfaces.
897       ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
898!
899!--    Allocate array for surface momentum flux for w - wsus and wsvs
900       ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) ) 
901!
902!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and
903!--    wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
904!--    on surface.
905       ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) ) 
906!
907!--    Characteristic temperature and surface flux of sensible heat
908       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )   
909       ALLOCATE ( surfaces%shf(1:surfaces%ns) )   
910!
911!--    Characteristic humidity and surface flux of latent heat
912       IF ( humidity )  THEN
913          ALLOCATE ( surfaces%qs(1:surfaces%ns)   ) 
914          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )     
915       ENDIF
916!
917!--    Characteristic scalar and surface flux of scalar
918       IF ( passive_scalar )  THEN
919          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
920          ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
921       ENDIF
922
923       IF ( cloud_physics .AND. microphysics_seifert)  THEN
924          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
925          ALLOCATE ( surfaces%ncs(1:surfaces%ns)   )
926          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
927          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
928       ENDIF
929
930       IF ( cloud_physics .AND. microphysics_seifert)  THEN
931          ALLOCATE ( surfaces%qrs(1:surfaces%ns)   )
932          ALLOCATE ( surfaces%nrs(1:surfaces%ns)   )
933          ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
934          ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
935       ENDIF
936!
937!--    Salinity surface flux
938       IF ( ocean )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
939
940    END SUBROUTINE allocate_surface_attributes_v
941
942!------------------------------------------------------------------------------!
943! Description:
944! ------------
945!> Initialize surface elements.
946!------------------------------------------------------------------------------!
947    SUBROUTINE init_surfaces
948
949       IMPLICIT NONE
950
951       INTEGER(iwp) ::  i         !< running index x-direction
952       INTEGER(iwp) ::  j         !< running index y-direction
953       INTEGER(iwp) ::  k         !< running index z-direction
954       INTEGER(iwp) ::  l         !< index variable used to distinguish surface facing
955       INTEGER(iwp) ::  m         !< running index surface elements
956
957       INTEGER(iwp)                 ::  start_index_lsm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal natural surfaces
958       INTEGER(iwp)                 ::  start_index_usm_h !< dummy to determing local start index in surface type for given (j,i), for horizontal urban surfaces
959
960       INTEGER(iwp)                 ::  num_lsm_h     !< current number of horizontal surface element, natural type
961       INTEGER(iwp)                 ::  num_lsm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal natural surfaces
962       INTEGER(iwp)                 ::  num_usm_h     !< current number of horizontal surface element, urban type
963       INTEGER(iwp)                 ::  num_usm_h_kji !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban surfaces
964
965       INTEGER(iwp), DIMENSION(0:2) ::  num_def_h     !< current number of horizontal surface element, default type
966       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
967       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
968     
969       INTEGER(iwp), DIMENSION(0:3) ::  num_def_v     !< current number of vertical surface element, default type
970       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
971       INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v     !< current number of vertical surface element, natural type
972       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
973       INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v     !< current number of vertical surface element, urban type
974       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
975
976       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
977       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
978       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
979
980
981!
982!--    Initialize surface attributes, store indicies, surfaces orientation, etc.,
983       num_def_h(0:2) = 1
984       num_def_v(0:3) = 1
985
986       num_lsm_h      = 1
987       num_lsm_v(0:3) = 1
988
989       num_usm_h      = 1
990       num_usm_v(0:3) = 1
991
992       start_index_def_h(0:2) = 1
993       start_index_def_v(0:3) = 1
994
995       start_index_lsm_h      = 1
996       start_index_lsm_v(0:3) = 1
997
998       start_index_usm_h      = 1
999       start_index_usm_v(0:3) = 1
1000
1001       DO  i = nxl, nxr
1002          DO  j = nys, nyn
1003
1004             num_def_h_kji = 0
1005             num_def_v_kji = 0
1006             num_lsm_h_kji = 0
1007             num_lsm_v_kji = 0
1008             num_usm_h_kji = 0
1009             num_usm_v_kji = 0
1010
1011             DO  k = nzb+1, nzt
1012!
1013!--             Check if current gridpoint belongs to the atmosphere
1014                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
1015!
1016!--                Upward-facing surface. Distinguish between differet surface types.
1017!--                To do, think about method to flag natural and non-natural
1018!--                surfaces. Only to ask for land_surface or urban surface
1019!--                is just a work-around.
1020                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN 
1021!
1022!--                   Natural surface type         
1023                      IF ( land_surface )  THEN
1024                         CALL initialize_horizontal_surfaces( k, j, i,         &
1025                                                              surf_lsm_h,      &
1026                                                              num_lsm_h,       &
1027                                                              num_lsm_h_kji,   &
1028                                                              .TRUE., .FALSE. ) 
1029!
1030!--                   Urban surface tpye
1031                      ELSEIF ( urban_surface )  THEN
1032                         CALL initialize_horizontal_surfaces( k, j, i,         &
1033                                                              surf_usm_h,      &
1034                                                              num_usm_h,       &
1035                                                              num_usm_h_kji,   &
1036                                                              .TRUE., .FALSE. ) 
1037!
1038!--                   Default surface type
1039                      ELSE
1040                         CALL initialize_horizontal_surfaces( k, j, i,         &
1041                                                              surf_def_h(0),   &
1042                                                              num_def_h(0),    &
1043                                                              num_def_h_kji(0),&
1044                                                              .TRUE., .FALSE. ) 
1045                      ENDIF
1046                   ENDIF 
1047!
1048!--                downward-facing surface, first, model top
1049                   IF ( k == nzt  .AND.  use_top_fluxes )  THEN
1050                      CALL initialize_top( k, j, i, surf_def_h(2),             &
1051                                           num_def_h(2), num_def_h_kji(2) )
1052!
1053!--                Check for any other downward-facing surface. So far only for
1054!--                default surface type.
1055                   ELSEIF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )  THEN
1056                      CALL initialize_horizontal_surfaces( k, j, i,            &
1057                                                           surf_def_h(1),      &
1058                                                           num_def_h(1),       &
1059                                                           num_def_h_kji(1),   &
1060                                                           .FALSE., .TRUE. )   
1061                   ENDIF
1062!
1063!--                Check for vertical walls and, if required, initialize it.
1064!                  Start with northward-facing surface.
1065                   IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) )  THEN
1066                      IF ( urban_surface )  THEN
1067                         CALL initialize_vertical_surfaces( 0, k, j, i,        &
1068                                                            surf_usm_v(0),     &
1069                                                            num_usm_v(0),      &
1070                                                            num_usm_v_kji(0),  &
1071                                                            .FALSE., .FALSE.,  &             
1072                                                            .FALSE., .TRUE. ) 
1073                      ELSEIF ( land_surface )  THEN
1074                         CALL initialize_vertical_surfaces( 0, k, j, i,        &
1075                                                            surf_lsm_v(0),     &
1076                                                            num_lsm_v(0),      &
1077                                                            num_lsm_v_kji(0),  &
1078                                                            .FALSE., .FALSE.,  &             
1079                                                            .FALSE., .TRUE. ) 
1080                      ELSE
1081                         CALL initialize_vertical_surfaces( 0, k, j, i,        &
1082                                                            surf_def_v(0),     &
1083                                                            num_def_v(0),      &
1084                                                            num_def_v_kji(0),  &
1085                                                            .FALSE., .FALSE.,  &             
1086                                                            .FALSE., .TRUE. ) 
1087                      ENDIF
1088                   ENDIF
1089!
1090!--                southward-facing surface
1091                   IF ( .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )  THEN
1092                      IF ( urban_surface )  THEN
1093                         CALL initialize_vertical_surfaces( 1, k, j, i,        &
1094                                                            surf_usm_v(1),     &
1095                                                            num_usm_v(1),      &
1096                                                            num_usm_v_kji(1),  &
1097                                                            .FALSE., .FALSE.,  &
1098                                                            .TRUE., .FALSE. )
1099                      ELSEIF ( land_surface )  THEN
1100                         CALL initialize_vertical_surfaces( 1, k, j, i,        &
1101                                                            surf_lsm_v(1),     &
1102                                                            num_lsm_v(1),      &
1103                                                            num_lsm_v_kji(1),  &
1104                                                            .FALSE., .FALSE.,  &
1105                                                            .TRUE., .FALSE. ) 
1106                      ELSE
1107                         CALL initialize_vertical_surfaces( 1, k, j, i,        &
1108                                                            surf_def_v(1),     &
1109                                                            num_def_v(1),      &
1110                                                            num_def_v_kji(1),  &
1111                                                            .FALSE., .FALSE.,  &
1112                                                            .TRUE., .FALSE. ) 
1113                      ENDIF
1114                   ENDIF
1115!
1116!--                eastward-facing surface
1117                   IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) )  THEN
1118                      IF ( urban_surface )  THEN
1119                         CALL initialize_vertical_surfaces( 2, k, j, i,        &
1120                                                            surf_usm_v(2),     &
1121                                                            num_usm_v(2),      &
1122                                                            num_usm_v_kji(2),  &
1123                                                            .TRUE., .FALSE.,   &
1124                                                            .FALSE., .FALSE. ) 
1125                      ELSEIF ( land_surface )  THEN
1126                         CALL initialize_vertical_surfaces( 2, k, j, i,        &
1127                                                            surf_lsm_v(2),     &
1128                                                            num_lsm_v(2),      &
1129                                                            num_lsm_v_kji(2),  &
1130                                                            .TRUE., .FALSE.,   &
1131                                                            .FALSE., .FALSE. ) 
1132                      ELSE
1133                         CALL initialize_vertical_surfaces( 2, k, j, i,        &
1134                                                            surf_def_v(2),     &
1135                                                            num_def_v(2),      &
1136                                                            num_def_v_kji(2),  &
1137                                                            .TRUE., .FALSE.,   &
1138                                                            .FALSE., .FALSE. ) 
1139                      ENDIF
1140                   ENDIF
1141!   
1142!--                westward-facing surface
1143                   IF ( .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )  THEN
1144                      IF ( urban_surface )  THEN
1145                         CALL initialize_vertical_surfaces( 3, k, j, i,        &
1146                                                            surf_usm_v(3),     &
1147                                                            num_usm_v(3),      &
1148                                                            num_usm_v_kji(3),  &
1149                                                           .FALSE., .TRUE.,    &
1150                                                           .FALSE., .FALSE. ) 
1151                      ELSEIF ( land_surface )  THEN
1152                         CALL initialize_vertical_surfaces( 3, k, j, i,        &
1153                                                            surf_lsm_v(3),     &
1154                                                            num_lsm_v(3),      &
1155                                                            num_lsm_v_kji(3),  &
1156                                                           .FALSE., .TRUE.,    &
1157                                                           .FALSE., .FALSE. ) 
1158                      ELSE
1159                         CALL initialize_vertical_surfaces( 3, k, j, i,        &
1160                                                            surf_def_v(3),     &
1161                                                            num_def_v(3),      &
1162                                                            num_def_v_kji(3),  &
1163                                                           .FALSE., .TRUE.,    &
1164                                                           .FALSE., .FALSE. ) 
1165                      ENDIF
1166                   ENDIF
1167                ENDIF
1168
1169 
1170             ENDDO
1171!
1172!--          Determine start- and end-index at grid point (j,i). Also, for
1173!--          horizontal surfaces more than 1 horizontal surface element can
1174!--          exist at grid point (j,i) if overhanging structures are present.
1175!--          Upward-facing surfaces
1176             surf_def_h(0)%start_index(j,i) = start_index_def_h(0)
1177             surf_def_h(0)%end_index(j,i)   = surf_def_h(0)%start_index(j,i) + &
1178                                                 num_def_h_kji(0) - 1
1179             start_index_def_h(0)           = surf_def_h(0)%end_index(j,i) + 1
1180!
1181!--          Downward-facing surfaces, except model top
1182             surf_def_h(1)%start_index(j,i) = start_index_def_h(1)                                                 
1183             surf_def_h(1)%end_index(j,i)   = surf_def_h(1)%start_index(j,i) + &
1184                                                 num_def_h_kji(1) - 1
1185             start_index_def_h(1)           = surf_def_h(1)%end_index(j,i) + 1
1186!
1187!--          Downward-facing surfaces -- model top fluxes
1188             surf_def_h(2)%start_index(j,i) = start_index_def_h(2)                                                 
1189             surf_def_h(2)%end_index(j,i)   = surf_def_h(2)%start_index(j,i) + &
1190                                                 num_def_h_kji(2) - 1
1191             start_index_def_h(2)           = surf_def_h(2)%end_index(j,i) + 1
1192!
1193!--          Horizontal natural land surfaces
1194             surf_lsm_h%start_index(j,i)    = start_index_lsm_h
1195             surf_lsm_h%end_index(j,i)      = surf_lsm_h%start_index(j,i) +    &
1196                                                 num_lsm_h_kji - 1
1197             start_index_lsm_h              = surf_lsm_h%end_index(j,i) + 1
1198!
1199!--          Horizontal urban surfaces
1200             surf_usm_h%start_index(j,i)    = start_index_usm_h
1201             surf_usm_h%end_index(j,i)      = surf_usm_h%start_index(j,i) +    &
1202                                                 num_usm_h_kji - 1
1203             start_index_usm_h              = surf_usm_h%end_index(j,i) + 1
1204
1205!
1206!--          Vertical surfaces - Default type
1207             surf_def_v(0)%start_index(j,i) = start_index_def_v(0)
1208             surf_def_v(1)%start_index(j,i) = start_index_def_v(1)
1209             surf_def_v(2)%start_index(j,i) = start_index_def_v(2)
1210             surf_def_v(3)%start_index(j,i) = start_index_def_v(3)
1211             surf_def_v(0)%end_index(j,i)   = start_index_def_v(0) +           & 
1212                                              num_def_v_kji(0) - 1
1213             surf_def_v(1)%end_index(j,i)   = start_index_def_v(1) +           &
1214                                              num_def_v_kji(1) - 1
1215             surf_def_v(2)%end_index(j,i)   = start_index_def_v(2) +           &
1216                                              num_def_v_kji(2) - 1
1217             surf_def_v(3)%end_index(j,i)   = start_index_def_v(3) +           &
1218                                              num_def_v_kji(3) - 1
1219             start_index_def_v(0)           = surf_def_v(0)%end_index(j,i) + 1
1220             start_index_def_v(1)           = surf_def_v(1)%end_index(j,i) + 1
1221             start_index_def_v(2)           = surf_def_v(2)%end_index(j,i) + 1
1222             start_index_def_v(3)           = surf_def_v(3)%end_index(j,i) + 1
1223!
1224!--          Natural type
1225             surf_lsm_v(0)%start_index(j,i) = start_index_lsm_v(0)
1226             surf_lsm_v(1)%start_index(j,i) = start_index_lsm_v(1)
1227             surf_lsm_v(2)%start_index(j,i) = start_index_lsm_v(2)
1228             surf_lsm_v(3)%start_index(j,i) = start_index_lsm_v(3)
1229             surf_lsm_v(0)%end_index(j,i)   = start_index_lsm_v(0) +           & 
1230                                              num_lsm_v_kji(0) - 1
1231             surf_lsm_v(1)%end_index(j,i)   = start_index_lsm_v(1) +           &
1232                                              num_lsm_v_kji(1) - 1
1233             surf_lsm_v(2)%end_index(j,i)   = start_index_lsm_v(2) +           &
1234                                              num_lsm_v_kji(2) - 1
1235             surf_lsm_v(3)%end_index(j,i)   = start_index_lsm_v(3) +           &
1236                                              num_lsm_v_kji(3) - 1
1237             start_index_lsm_v(0)           = surf_lsm_v(0)%end_index(j,i) + 1
1238             start_index_lsm_v(1)           = surf_lsm_v(1)%end_index(j,i) + 1
1239             start_index_lsm_v(2)           = surf_lsm_v(2)%end_index(j,i) + 1
1240             start_index_lsm_v(3)           = surf_lsm_v(3)%end_index(j,i) + 1
1241!
1242!--          Urban type
1243             surf_usm_v(0)%start_index(j,i) = start_index_usm_v(0)
1244             surf_usm_v(1)%start_index(j,i) = start_index_usm_v(1)
1245             surf_usm_v(2)%start_index(j,i) = start_index_usm_v(2)
1246             surf_usm_v(3)%start_index(j,i) = start_index_usm_v(3)
1247             surf_usm_v(0)%end_index(j,i)   = start_index_usm_v(0) +           & 
1248                                              num_usm_v_kji(0) - 1
1249             surf_usm_v(1)%end_index(j,i)   = start_index_usm_v(1) +           &
1250                                              num_usm_v_kji(1) - 1
1251             surf_usm_v(2)%end_index(j,i)   = start_index_usm_v(2) +           &
1252                                              num_usm_v_kji(2) - 1
1253             surf_usm_v(3)%end_index(j,i)   = start_index_usm_v(3) +           &
1254                                              num_usm_v_kji(3) - 1
1255             start_index_usm_v(0)           = surf_usm_v(0)%end_index(j,i) + 1
1256             start_index_usm_v(1)           = surf_usm_v(1)%end_index(j,i) + 1
1257             start_index_usm_v(2)           = surf_usm_v(2)%end_index(j,i) + 1
1258             start_index_usm_v(3)           = surf_usm_v(3)%end_index(j,i) + 1
1259
1260
1261          ENDDO
1262       ENDDO
1263
1264       CONTAINS
1265
1266!------------------------------------------------------------------------------!
1267! Description:
1268! ------------
1269!> Initialize horizontal surface elements, upward- and downward-facing.
1270!> Note, horizontal surface type alsw comprises model-top fluxes, which are,
1271!> initialized in a different routine.
1272!------------------------------------------------------------------------------!
1273          SUBROUTINE initialize_horizontal_surfaces( k, j, i, surf, num_h,     &
1274                                                     num_h_kji, upward_facing, &
1275                                                     downward_facing )       
1276
1277             IMPLICIT NONE
1278
1279             INTEGER(iwp)  ::  i                !< running index x-direction
1280             INTEGER(iwp)  ::  j                !< running index y-direction
1281             INTEGER(iwp)  ::  k                !< running index z-direction
1282             INTEGER(iwp)  ::  num_h            !< current number of surface element
1283             INTEGER(iwp)  ::  num_h_kji        !< dummy increment
1284
1285             LOGICAL       ::  upward_facing    !< flag indicating upward-facing surface
1286             LOGICAL       ::  downward_facing  !< flag indicating downward-facing surface
1287
1288             TYPE( surf_type ) :: surf          !< respective surface type
1289!
1290!--          Store indices of respective surface element
1291             surf%i(num_h) = i
1292             surf%j(num_h) = j
1293             surf%k(num_h) = k
1294!
1295!--          Surface orientation, bit 0 is set to 1 for upward-facing surfaces,
1296!--          bit 1 is for downward-facing surfaces.
1297             IF ( upward_facing   )  surf%facing(num_h) = IBSET( surf%facing(num_h), 0 )
1298             IF ( downward_facing )  surf%facing(num_h) = IBSET( surf%facing(num_h), 1 )
1299!
1300!--          Initialize surface-layer height
1301             IF ( upward_facing )  THEN
1302                surf%z_mo(num_h)  = zu(k) - zw(k-1)
1303             ELSE
1304                surf%z_mo(num_h)  = zw(k) - zu(k)
1305             ENDIF
1306 
1307             surf%z0(num_h)    = roughness_length
1308             surf%z0h(num_h)   = z0h_factor * roughness_length
1309             surf%z0q(num_h)   = z0h_factor * roughness_length         
1310!
1311!--          Initialization in case of 1D pre-cursor run
1312             IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )&
1313             THEN
1314                IF ( .NOT. constant_diffusion )  THEN
1315                   IF ( constant_flux_layer )  THEN
1316                      surf%ol(num_h)   = surf%z_mo(num_h) /                    &
1317                                            ( rif1d(nzb+1) + 1.0E-20_wp )
1318                      surf%us(num_h)   = us1d
1319                      surf%usws(num_h) = usws1d
1320                      surf%vsws(num_h) = vsws1d
1321                   ELSE
1322                      surf%ol(num_h)   = surf%z_mo(num_h) / zeta_min
1323                      surf%us(num_h)   = 0.0_wp
1324                      surf%usws(num_h) = 0.0_wp
1325                      surf%vsws(num_h) = 0.0_wp
1326                   ENDIF
1327                ELSE
1328                   surf%ol(num_h)   = surf%z_mo(num_h) / zeta_min
1329                   surf%us(num_h)   = 0.0_wp
1330                   surf%usws(num_h) = 0.0_wp
1331                   surf%vsws(num_h) = 0.0_wp
1332                ENDIF
1333!
1334!--          Initialization in case of constant profiles
1335             ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )&
1336             THEN
1337
1338                surf%ol(num_h)   = surf%z_mo(num_h) / zeta_min
1339!
1340!--             Very small number is required for calculation of Obukhov length
1341!--             at first timestep     
1342                surf%us(num_h)    = 1E-30_wp 
1343                surf%usws(num_h)  = 0.0_wp
1344                surf%vsws(num_h)  = 0.0_wp
1345       
1346             ENDIF
1347
1348             surf%rib(num_h)   = 0.0_wp 
1349             surf%uvw_abs(num_h) = 0.0_wp
1350
1351             IF ( .NOT. constant_diffusion )  THEN   
1352                surf%u_0(num_h)     = 0.0_wp 
1353                surf%v_0(num_h)     = 0.0_wp
1354             ENDIF
1355
1356             surf%ts(num_h)   = 0.0_wp
1357
1358             IF ( humidity )  THEN
1359                surf%qs(num_h)   = 0.0_wp
1360                IF ( cloud_physics .AND. microphysics_morrison)  THEN
1361                   surf%qcs(num_h) = 0.0_wp
1362                   surf%ncs(num_h) = 0.0_wp
1363   
1364                   surf%qcsws(num_h) = 0.0_wp
1365                   surf%ncsws(num_h) = 0.0_wp
1366
1367                ENDIF
1368                IF ( cloud_physics .AND. microphysics_seifert)  THEN
1369                   surf%qrs(num_h) = 0.0_wp
1370                   surf%nrs(num_h) = 0.0_wp
1371   
1372                   surf%qrsws(num_h) = 0.0_wp
1373                   surf%nrsws(num_h) = 0.0_wp
1374
1375                   surf%pt1(num_h) = 0.0_wp
1376                   surf%qv1(num_h) = 0.0_wp
1377
1378                ENDIF
1379             ENDIF
1380
1381             IF ( passive_scalar )  surf%ss(num_h) = 0.0_wp
1382!
1383!--          Inititalize surface fluxes of sensible and latent heat, as well as
1384!--          passive scalar
1385             IF ( use_surface_fluxes )  THEN
1386
1387                IF ( upward_facing )  THEN
1388                   IF ( constant_heatflux )  THEN
1389!   
1390!--                   Initialize surface heatflux. However, skip this for now if
1391!--                   if random_heatflux is set. This case, shf is initialized later.
1392                      IF ( .NOT. random_heatflux )  THEN
1393                         surf%shf(num_h) = surface_heatflux *               &
1394                                                 heatflux_input_conversion(nzb)
1395!
1396!--                      Check if surface heat flux might be replaced by
1397!--                      prescribed wall heatflux
1398                         IF ( k-1 /= 0 )  THEN
1399                            surf%shf(num_h) = wall_heatflux(0) *            &
1400                                                 heatflux_input_conversion(k-1)
1401                         ENDIF
1402!
1403!--                      Initialize shf with data from external file LSF_DATA. Will
1404!--                      be done directly in ls_foring_surf
1405!--                      Attention: Just a workaround, need to be revised!!!
1406                         IF ( large_scale_forcing .AND. lsf_surf )  THEN
1407!                             CALL ls_forcing_surf ( simulated_time )
1408!                             surf%shf(num_h) = shf(j,i) 
1409                         ENDIF
1410                      ENDIF
1411                   ELSE
1412                      surf%shf(num_h) = 0.0_wp
1413                   ENDIF
1414!
1415!--             Set heat-flux at downward-facing surfaces
1416                ELSE
1417                   surf%shf(num_h) = wall_heatflux(5) *                        &
1418                                             heatflux_input_conversion(k)
1419                ENDIF
1420
1421                IF ( humidity )  THEN
1422                   IF ( upward_facing )  THEN
1423                      IF ( constant_waterflux )  THEN
1424                         surf%qsws(num_h) = surface_waterflux *                &
1425                                                 waterflux_input_conversion(nzb)
1426                         IF ( k-1 /= 0 )  THEN
1427                            surf%qsws(num_h) = wall_humidityflux(0) *          &
1428                                                 waterflux_input_conversion(k-1)
1429                         ENDIF
1430                      ELSE
1431                         surf%qsws(num_h) = 0.0_wp
1432                      ENDIF
1433                   ELSE
1434                      surf%qsws(num_h) = wall_humidityflux(5) *                &
1435                                             heatflux_input_conversion(k)
1436                   ENDIF
1437                ENDIF
1438
1439                IF ( passive_scalar )  THEN
1440                   IF ( upward_facing )  THEN
1441                      IF ( constant_scalarflux )  THEN
1442                         surf%ssws(num_h) = surface_scalarflux
1443
1444                         IF ( k-1 /= 0 )                                       &
1445                            surf%ssws(num_h) = wall_scalarflux(0)
1446
1447                      ELSE
1448                         surf%ssws(num_h) = 0.0_wp
1449                      ENDIF
1450                   ELSE
1451                      surf%ssws(num_h) = wall_scalarflux(5)
1452                   ENDIF
1453                ENDIF
1454
1455                IF ( ocean )  THEN
1456                   IF ( upward_facing )  THEN
1457                      surf%sasws(num_h) = bottom_salinityflux
1458                   ELSE
1459                      surf%sasws(num_h) = 0.0_wp
1460                   ENDIF
1461                ENDIF
1462             ENDIF
1463!
1464!--          Increment surface indices
1465             num_h     = num_h + 1
1466             num_h_kji = num_h_kji + 1     
1467
1468
1469          END SUBROUTINE initialize_horizontal_surfaces
1470       
1471
1472!------------------------------------------------------------------------------!
1473! Description:
1474! ------------
1475!> Initialize model-top fluxes. Currently, only the heatflux and salinity flux
1476!> can be prescribed, latent flux is zero in this case!
1477!------------------------------------------------------------------------------!
1478          SUBROUTINE initialize_top( k, j, i, surf, num_h, num_h_kji )       
1479
1480             IMPLICIT NONE
1481
1482             INTEGER(iwp)  ::  i                !< running index x-direction
1483             INTEGER(iwp)  ::  j                !< running index y-direction
1484             INTEGER(iwp)  ::  k                !< running index z-direction
1485             INTEGER(iwp)  ::  num_h            !< current number of surface element
1486             INTEGER(iwp)  ::  num_h_kji        !< dummy increment
1487
1488             TYPE( surf_type ) :: surf          !< respective surface type
1489!
1490!--          Store indices of respective surface element
1491             surf%i(num_h) = i
1492             surf%j(num_h) = j
1493             surf%k(num_h) = k
1494!
1495!--          Initialize top heat flux
1496             IF ( constant_top_heatflux )                                      &
1497                surf%shf(num_h) = top_heatflux * heatflux_input_conversion(nzt+1)
1498!
1499!--          Initialization in case of a coupled model run
1500             IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1501                surf%shf(num_h) = 0.0_wp
1502             ENDIF
1503!
1504!--          Prescribe latent heat flux at the top     
1505             IF ( humidity )  THEN
1506             surf%qsws(num_h) = 0.0_wp
1507                IF ( cloud_physics  .AND.  microphysics_morrison ) THEN
1508                   surf%ncsws(num_h) = 0.0_wp
1509                   surf%qcsws(num_h) = 0.0_wp
1510                ENDIF
1511                IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
1512                   surf%nrsws(num_h) = 0.0_wp
1513                   surf%qrsws(num_h) = 0.0_wp
1514                ENDIF
1515             ENDIF
1516!
1517!--          Prescribe top scalar flux
1518             IF ( passive_scalar .AND. constant_top_scalarflux )               &
1519                surf%ssws(num_h) = top_scalarflux
1520!
1521!--          Prescribe top salinity flux
1522             IF ( ocean .AND. constant_top_salinityflux)                       &
1523                surf%sasws(num_h) = top_salinityflux
1524!
1525!--          Initialization in case of a coupled model run
1526             IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1527                surf%shf(num_h) = 0.0_wp
1528             ENDIF
1529!
1530!--          Top momentum fluxes
1531             surf%usws(num_h) = top_momentumflux_u *                           &
1532                                            momentumflux_input_conversion(nzt+1)
1533             surf%vsws(num_h) = top_momentumflux_v *                           &
1534                                            momentumflux_input_conversion(nzt+1)
1535!
1536!--          Increment surface indices
1537             num_h     = num_h + 1
1538             num_h_kji = num_h_kji + 1     
1539
1540
1541          END SUBROUTINE initialize_top
1542
1543
1544!------------------------------------------------------------------------------!
1545! Description:
1546! ------------
1547!> Initialize vertical surface elements.
1548!------------------------------------------------------------------------------!
1549          SUBROUTINE initialize_vertical_surfaces( l, k, j, i, surf, num_v,    &
1550                                                num_v_kji, east_facing,        &
1551                                                west_facing, south_facing,     &
1552                                                north_facing )       
1553
1554             IMPLICIT NONE
1555
1556             INTEGER(iwp)  ::  component !<
1557             INTEGER(iwp)  ::  i               !< running index x-direction
1558             INTEGER(iwp)  ::  j               !< running index x-direction
1559             INTEGER(iwp)  ::  k               !< running index x-direction
1560             INTEGER(iwp)  ::  l               !< index variable for the surface type, indicating the facing
1561             INTEGER(iwp)  ::  num_v           !< current number of surface element
1562             INTEGER(iwp)  ::  num_v_kji       !< current number of surface element at (j,i)
1563
1564
1565             LOGICAL       ::  east_facing     !< flag indicating east-facing surfaces
1566             LOGICAL       ::  north_facing    !< flag indicating north-facing surfaces
1567             LOGICAL       ::  south_facing    !< flag indicating south-facing surfaces
1568             LOGICAL       ::  west_facing     !< flag indicating west-facing surfaces
1569
1570             TYPE( surf_type ) :: surf         !< respective surface type
1571
1572!
1573!--          Store indices of respective wall element
1574             surf%i(num_v)   = i
1575             surf%j(num_v)   = j
1576             surf%k(num_v)   = k
1577!
1578!--          Initialize surface-layer height, or more precisely, distance to surface
1579             IF ( north_facing  .OR.  south_facing )  THEN
1580                surf%z_mo(num_v)  = 0.5_wp * dy
1581             ELSE
1582                surf%z_mo(num_v)  = 0.5_wp * dx
1583             ENDIF
1584
1585             surf%facing(num_v)  = 0
1586!
1587!--          Surface orientation. Moreover, set component id to map wall_heatflux,
1588!--          etc., on surface type (further below)
1589             IF ( north_facing )  THEN
1590                surf%facing(num_v) = IBSET( surf%facing(num_v), 0 ) 
1591                component          = 4
1592             ENDIF
1593
1594             IF ( south_facing )  THEN
1595                surf%facing(num_v) = IBSET( surf%facing(num_v), 1 ) 
1596                component          = 3
1597             ENDIF
1598
1599             IF ( east_facing )  THEN
1600                surf%facing(num_v) = IBSET( surf%facing(num_v), 2 )
1601                component          = 2
1602             ENDIF
1603
1604             IF ( west_facing )  THEN
1605                surf%facing(num_v) = IBSET( surf%facing(num_v), 3 ) 
1606                component          = 1
1607             ENDIF
1608
1609 
1610             surf%z0(num_v)  = roughness_length
1611             surf%z0h(num_v) = z0h_factor * roughness_length
1612             surf%z0q(num_v) = z0h_factor * roughness_length
1613
1614             surf%us(num_v)  = 0.0_wp
1615!
1616!--          If required, initialize Obukhov length
1617             IF ( ALLOCATED( surf%ol ) )                                       &
1618                surf%ol(num_v) = surf%z_mo(num_v) / zeta_min
1619
1620             surf%uvw_abs(num_v)   = 0.0_wp
1621
1622             surf%mom_flux_uv(num_v) = 0.0_wp
1623             surf%mom_flux_w(num_v)  = 0.0_wp
1624             surf%mom_flux_tke(0:1,num_v) = 0.0_wp
1625
1626             surf%ts(num_v)    = 0.0_wp
1627             surf%shf(num_v)   = wall_heatflux(component)
1628
1629             IF ( humidity )  THEN
1630                surf%qs(num_v)   = 0.0_wp
1631                surf%qsws(num_v) = wall_humidityflux(component)
1632!
1633!--             Following wall fluxes are assumed to be zero
1634                IF ( cloud_physics .AND. microphysics_morrison)  THEN
1635                   surf%qcs(num_v) = 0.0_wp
1636                   surf%ncs(num_v) = 0.0_wp
1637   
1638                   surf%qcsws(num_v) = 0.0_wp
1639                   surf%ncsws(num_v) = 0.0_wp
1640                ENDIF
1641                IF ( cloud_physics .AND. microphysics_seifert)  THEN
1642                   surf%qrs(num_v) = 0.0_wp
1643                   surf%nrs(num_v) = 0.0_wp
1644   
1645                   surf%qrsws(num_v) = 0.0_wp
1646                   surf%nrsws(num_v) = 0.0_wp
1647                ENDIF
1648             ENDIF
1649
1650             IF ( passive_scalar )  THEN
1651                surf%ss(num_v)   = 0.0_wp
1652                surf%ssws(num_v) = wall_scalarflux(component)
1653             ENDIF
1654!
1655!--          So far, salinityflux at vertical surfaces is simply zero
1656!--          at the moment 
1657             IF ( ocean )  surf%sasws(num_v) = wall_salinityflux(component)
1658!
1659!--          Increment wall indices
1660             num_v                 = num_v + 1
1661             num_v_kji             = num_v_kji + 1
1662
1663          END SUBROUTINE initialize_vertical_surfaces
1664
1665    END SUBROUTINE init_surfaces
1666
1667
1668!------------------------------------------------------------------------------!
1669! Description:
1670! ------------
1671!> Determines topography-top index at given (j,i)-position. 
1672!------------------------------------------------------------------------------!
1673    FUNCTION get_topography_top_index( j, i, grid )
1674
1675       USE kinds
1676
1677       IMPLICIT NONE
1678
1679       CHARACTER(LEN=*) ::  grid                      !< flag to distinquish between staggered grids
1680       INTEGER(iwp)     ::  i                         !< grid index in x-dimension
1681       INTEGER(iwp)     ::  ibit                      !< bit position where topography information is stored on respective grid
1682       INTEGER(iwp)     ::  j                         !< grid index in y-dimension
1683       INTEGER(iwp)     ::  get_topography_top_index  !< topography top index
1684
1685       SELECT CASE ( TRIM( grid ) )
1686
1687          CASE ( 's'     )
1688             ibit = 12
1689          CASE ( 'u'     )
1690             ibit = 14
1691          CASE ( 'v'     )
1692             ibit = 16
1693          CASE ( 'w'     )
1694             ibit = 18
1695          CASE ( 's_out' )
1696             ibit = 24
1697          CASE ( 'u_out' )
1698             ibit = 26
1699          CASE ( 'v_out' )
1700             ibit = 27
1701          CASE ( 'w_out' )
1702             ibit = 28
1703          CASE DEFAULT
1704!
1705!--          Set default to scalar grid
1706             ibit = 12 
1707
1708       END SELECT
1709
1710       get_topography_top_index = MAXLOC(                                      &
1711                                     MERGE( 1, 0,                              &
1712                                            BTEST( wall_flags_0(:,j,i), ibit ) &
1713                                          ), DIM = 1                           &
1714                                        ) - 1
1715
1716       RETURN
1717
1718    END FUNCTION get_topography_top_index
1719
1720!------------------------------------------------------------------------------!
1721! Description:
1722! ------------
1723!> Gathers all surface elements with the same facing (but possibly different
1724!> type) onto a surface type, and writes binary data into restart files.
1725!------------------------------------------------------------------------------!
1726    SUBROUTINE surface_write_restart_data
1727
1728       IMPLICIT NONE
1729
1730       CHARACTER(LEN=1)             ::  dum  !< dummy string to create output-variable name
1731
1732       INTEGER(iwp)                 ::  i    !< running index x-direction
1733       INTEGER(iwp)                 ::  j    !< running index y-direction
1734       INTEGER(iwp)                 ::  l    !< index surface type orientation
1735       INTEGER(iwp)                 ::  m    !< running index for surface elements on individual surface array
1736       INTEGER(iwp), DIMENSION(0:3) ::  mm   !< running index for surface elements on gathered surface array
1737
1738       TYPE(surf_type), DIMENSION(0:2) ::  surf_h !< gathered horizontal surfaces, contains all surface types
1739       TYPE(surf_type), DIMENSION(0:3) ::  surf_v !< gathered vertical surfaces, contains all surface types
1740
1741!
1742!--    Determine total number of horizontal and vertical surface elements before
1743!--    writing var_list
1744       CALL surface_last_actions
1745!
1746!--    Count number of grid points with same facing and allocate attributes respectively
1747!--    Horizontal upward facing
1748       surf_h(0)%ns = ns_h_on_file(0)
1749       CALL allocate_surface_attributes_h( surf_h(0), nys, nyn, nxl, nxr )
1750!
1751!--    Horizontal downward facing
1752       surf_h(1)%ns = ns_h_on_file(1)
1753       CALL allocate_surface_attributes_h( surf_h(1), nys, nyn, nxl, nxr )
1754!
1755!--    Model top
1756       surf_h(2)%ns = ns_h_on_file(2)
1757       CALL allocate_surface_attributes_h_top( surf_h(2), nys, nyn, nxl, nxr )
1758!
1759!--    Vertical surfaces
1760       DO  l = 0, 3
1761          surf_v(l)%ns = ns_v_on_file(l)
1762          CALL allocate_surface_attributes_v( surf_v(l), .FALSE.,              &
1763                                              nys, nyn, nxl, nxr )
1764       ENDDO
1765!
1766!--    In the following, gather data from surfaces elements with the same
1767!--    facing (but possibly differt type) on 1 data-type array.
1768       mm(0:2) = 1
1769       DO  l = 0, 2
1770          DO  i = nxl, nxr
1771             DO  j = nys, nyn
1772                DO  m = surf_def_h(l)%start_index(j,i),                        &
1773                        surf_def_h(l)%end_index(j,i)
1774                   IF ( ALLOCATED( surf_def_h(l)%us ) )                        &
1775                      surf_h(l)%us(mm(l))      = surf_def_h(l)%us(m)
1776                   IF ( ALLOCATED( surf_def_h(l)%ts ) )                        &
1777                      surf_h(l)%ts(mm(l))      = surf_def_h(l)%ts(m)
1778                   IF ( ALLOCATED( surf_def_h(l)%qs ) )                        &
1779                      surf_h(l)%qs(mm(l))      = surf_def_h(l)%qs(m)
1780                   IF ( ALLOCATED( surf_def_h(l)%ss ) )                        &
1781                      surf_h(l)%ss(mm(l))      = surf_def_h(l)%ss(m)
1782                   IF ( ALLOCATED( surf_def_h(l)%qcs ) )                       &
1783                      surf_h(l)%qcs(mm(l))     = surf_def_h(l)%qcs(m)
1784                   IF ( ALLOCATED( surf_def_h(l)%ncs ) )                       &
1785                      surf_h(l)%ncs(mm(l))     = surf_def_h(l)%ncs(m)
1786                   IF ( ALLOCATED( surf_def_h(l)%qrs ) )                       &
1787                      surf_h(l)%qrs(mm(l))     = surf_def_h(l)%qrs(m)
1788                   IF ( ALLOCATED( surf_def_h(l)%nrs ) )                       &
1789                      surf_h(l)%nrs(mm(l))     = surf_def_h(l)%nrs(m)
1790                   IF ( ALLOCATED( surf_def_h(l)%ol ) )                        &
1791                      surf_h(l)%ol(mm(l))      = surf_def_h(l)%ol(m)
1792                   IF ( ALLOCATED( surf_def_h(l)%rib ) )                       &
1793                      surf_h(l)%rib(mm(l))     = surf_def_h(l)%rib(m)
1794                   IF ( ALLOCATED( surf_def_h(l)%usws ) )                      &
1795                      surf_h(l)%usws(mm(l))    = surf_def_h(l)%usws(m)
1796                   IF ( ALLOCATED( surf_def_h(l)%vsws ) )                      &
1797                      surf_h(l)%vsws(mm(l))    = surf_def_h(l)%vsws(m)
1798                   IF ( ALLOCATED( surf_def_h(l)%shf ) )                       &
1799                      surf_h(l)%shf(mm(l))     = surf_def_h(l)%shf(m)
1800                   IF ( ALLOCATED( surf_def_h(l)%qsws ) )                      &
1801                      surf_h(l)%qsws(mm(l))    = surf_def_h(l)%qsws(m)
1802                   IF ( ALLOCATED( surf_def_h(l)%ssws ) )                      &
1803                      surf_h(l)%ssws(mm(l))    = surf_def_h(l)%ssws(m)
1804                   IF ( ALLOCATED( surf_def_h(l)%ncsws ) )                     &
1805                      surf_h(l)%ncsws(mm(l))   = surf_def_h(l)%ncsws(m)
1806                   IF ( ALLOCATED( surf_def_h(l)%nrsws ) )                     &
1807                      surf_h(l)%nrsws(mm(l))   = surf_def_h(l)%nrsws(m)
1808                   IF ( ALLOCATED( surf_def_h(l)%sasws ) )                     &
1809                      surf_h(l)%sasws(mm(l))   = surf_def_h(l)%sasws(m)
1810               
1811                   mm(l) = mm(l) + 1
1812                ENDDO
1813
1814                IF ( l == 0 )  THEN
1815                   DO  m = surf_lsm_h%start_index(j,i),                        &
1816                           surf_lsm_h%end_index(j,i)
1817                      IF ( ALLOCATED( surf_lsm_h%us ) )                        &
1818                         surf_h(0)%us(mm(0))      = surf_lsm_h%us(m)
1819                      IF ( ALLOCATED( surf_lsm_h%ts ) )                        &
1820                         surf_h(0)%ts(mm(0))      = surf_lsm_h%ts(m)
1821                      IF ( ALLOCATED( surf_lsm_h%qs ) )                        &
1822                         surf_h(0)%qs(mm(0))      = surf_lsm_h%qs(m)
1823                      IF ( ALLOCATED( surf_lsm_h%ss ) )                        &
1824                         surf_h(0)%ss(mm(0))      = surf_lsm_h%ss(m)
1825                      IF ( ALLOCATED( surf_lsm_h%qcs ) )                       &
1826                         surf_h(0)%qcs(mm(0))     = surf_lsm_h%qcs(m)
1827                      IF ( ALLOCATED( surf_lsm_h%ncs ) )                       &
1828                         surf_h(0)%ncs(mm(0))     = surf_lsm_h%ncs(m)
1829                      IF ( ALLOCATED( surf_lsm_h%qrs ) )                       &
1830                         surf_h(0)%qrs(mm(0))     = surf_lsm_h%qrs(m)
1831                      IF ( ALLOCATED( surf_lsm_h%nrs ) )                       &
1832                         surf_h(0)%nrs(mm(0))     = surf_lsm_h%nrs(m)
1833                      IF ( ALLOCATED( surf_lsm_h%ol ) )                        &
1834                         surf_h(0)%ol(mm(0))      = surf_lsm_h%ol(m)
1835                      IF ( ALLOCATED( surf_lsm_h%rib ) )                       &
1836                         surf_h(0)%rib(mm(0))     = surf_lsm_h%rib(m)
1837                      IF ( ALLOCATED( surf_lsm_h%usws ) )                      &
1838                         surf_h(0)%usws(mm(0))    = surf_lsm_h%usws(m)
1839                      IF ( ALLOCATED( surf_lsm_h%vsws ) )                      &
1840                         surf_h(0)%vsws(mm(0))    = surf_lsm_h%vsws(m)
1841                      IF ( ALLOCATED( surf_lsm_h%shf ) )                       &
1842                         surf_h(0)%shf(mm(0))     = surf_lsm_h%shf(m)
1843                      IF ( ALLOCATED( surf_lsm_h%qsws ) )                      &
1844                         surf_h(0)%qsws(mm(0))    = surf_lsm_h%qsws(m)
1845                      IF ( ALLOCATED( surf_lsm_h%ssws ) )                      &
1846                         surf_h(0)%ssws(mm(0))    = surf_lsm_h%ssws(m)
1847                      IF ( ALLOCATED( surf_lsm_h%ncsws ) )                     &
1848                         surf_h(0)%ncsws(mm(0))   = surf_lsm_h%ncsws(m)
1849                      IF ( ALLOCATED( surf_lsm_h%nrsws ) )                     &
1850                         surf_h(0)%nrsws(mm(0))   = surf_lsm_h%nrsws(m)
1851                      IF ( ALLOCATED( surf_lsm_h%sasws ) )                     &
1852                        surf_h(0)%sasws(mm(0))   = surf_lsm_h%sasws(m)
1853               
1854                      mm(0) = mm(0) + 1
1855             
1856                   ENDDO
1857
1858                   DO  m = surf_usm_h%start_index(j,i),                        &
1859                           surf_usm_h%end_index(j,i)
1860                      IF ( ALLOCATED( surf_usm_h%us ) )                        &
1861                         surf_h(0)%us(mm(0))      = surf_usm_h%us(m)
1862                      IF ( ALLOCATED( surf_usm_h%ts ) )                        &
1863                         surf_h(0)%ts(mm(0))      = surf_usm_h%ts(m)
1864                      IF ( ALLOCATED( surf_usm_h%qs ) )                        &
1865                         surf_h(0)%qs(mm(0))      = surf_usm_h%qs(m)
1866                      IF ( ALLOCATED( surf_usm_h%ss ) )                        &
1867                         surf_h(0)%ss(mm(0))      = surf_usm_h%ss(m)
1868                      IF ( ALLOCATED( surf_usm_h%qcs ) )                       &
1869                         surf_h(0)%qcs(mm(0))     = surf_usm_h%qcs(m)
1870                      IF ( ALLOCATED( surf_usm_h%ncs ) )                       &
1871                         surf_h(0)%ncs(mm(0))     = surf_usm_h%ncs(m)
1872                      IF ( ALLOCATED( surf_usm_h%qrs ) )                       &
1873                         surf_h(0)%qrs(mm(0))     = surf_usm_h%qrs(m)
1874                      IF ( ALLOCATED( surf_usm_h%nrs ) )                       &
1875                         surf_h(0)%nrs(mm(0))     = surf_usm_h%nrs(m)
1876                      IF ( ALLOCATED( surf_usm_h%ol ) )                        &
1877                         surf_h(0)%ol(mm(0))      = surf_usm_h%ol(m)
1878                      IF ( ALLOCATED( surf_usm_h%rib ) )                       &
1879                         surf_h(0)%rib(mm(0))     = surf_usm_h%rib(m)
1880                      IF ( ALLOCATED( surf_usm_h%usws ) )                      &
1881                         surf_h(0)%usws(mm(0))    = surf_usm_h%usws(m)
1882                      IF ( ALLOCATED( surf_usm_h%vsws ) )                      &
1883                         surf_h(0)%vsws(mm(0))    = surf_usm_h%vsws(m)
1884                      IF ( ALLOCATED( surf_usm_h%shf ) )                       &
1885                         surf_h(0)%shf(mm(0))     = surf_usm_h%shf(m)
1886                      IF ( ALLOCATED( surf_usm_h%qsws ) )                      &
1887                         surf_h(0)%qsws(mm(0))    = surf_usm_h%qsws(m)
1888                      IF ( ALLOCATED( surf_usm_h%ssws ) )                      &
1889                         surf_h(0)%ssws(mm(0))    = surf_usm_h%ssws(m)
1890                      IF ( ALLOCATED( surf_usm_h%ncsws ) )                     &
1891                         surf_h(0)%ncsws(mm(0))   = surf_usm_h%ncsws(m)
1892                      IF ( ALLOCATED( surf_usm_h%nrsws ) )                     &
1893                         surf_h(0)%nrsws(mm(0))   = surf_usm_h%nrsws(m)
1894                      IF ( ALLOCATED( surf_usm_h%sasws ) )                     &
1895                        surf_h(0)%sasws(mm(0))   = surf_usm_h%sasws(m)
1896               
1897                      mm(0) = mm(0) + 1
1898             
1899                   ENDDO
1900
1901
1902                ENDIF
1903
1904             ENDDO
1905
1906          ENDDO
1907          IF ( l == 0 )  THEN
1908             surf_h(l)%start_index = MAX( surf_def_h(l)%start_index,           &
1909                                          surf_lsm_h%start_index,              &
1910                                          surf_usm_h%start_index )
1911             surf_h(l)%end_index   = MAX( surf_def_h(l)%end_index,             &
1912                                          surf_lsm_h%end_index,                &
1913                                          surf_usm_h%end_index )
1914          ELSE
1915             surf_h(l)%start_index = surf_def_h(l)%start_index
1916             surf_h(l)%end_index   = surf_def_h(l)%end_index
1917          ENDIF
1918       ENDDO
1919
1920
1921       mm(0:3) = 1
1922       DO  l = 0, 3
1923          DO  i = nxl, nxr
1924             DO  j = nys, nyn
1925                DO  m = surf_def_v(l)%start_index(j,i),                        &
1926                        surf_def_v(l)%end_index(j,i)
1927                   IF ( ALLOCATED( surf_def_v(l)%us ) )                        &
1928                      surf_v(l)%us(mm(l))      = surf_def_v(l)%us(m)
1929                   IF ( ALLOCATED( surf_def_v(l)%ts ) )                        &
1930                      surf_v(l)%ts(mm(l))      = surf_def_v(l)%ts(m)
1931                   IF ( ALLOCATED( surf_def_v(l)%qs ) )                        &
1932                      surf_v(l)%qs(mm(l))      = surf_def_v(l)%qs(m)
1933                   IF ( ALLOCATED( surf_def_v(l)%ss ) )                        &
1934                      surf_v(l)%ss(mm(l))      = surf_def_v(l)%ss(m)
1935                   IF ( ALLOCATED( surf_def_v(l)%qcs ) )                       &
1936                      surf_v(l)%qcs(mm(l))     = surf_def_v(l)%qcs(m)
1937                   IF ( ALLOCATED( surf_def_v(l)%ncs ) )                       &
1938                      surf_v(l)%ncs(mm(l))     = surf_def_v(l)%ncs(m)
1939                   IF ( ALLOCATED( surf_def_v(l)%qrs ) )                       &
1940                      surf_v(l)%qrs(mm(l))     = surf_def_v(l)%qrs(m)
1941                   IF ( ALLOCATED( surf_def_v(l)%nrs ) )                       &
1942                      surf_v(l)%nrs(mm(l))     = surf_def_v(l)%nrs(m)
1943                   IF ( ALLOCATED( surf_def_v(l)%ol ) )                        &
1944                      surf_v(l)%ol(mm(l))      = surf_def_v(l)%ol(m)
1945                   IF ( ALLOCATED( surf_def_v(l)%rib ) )                       &
1946                      surf_v(l)%rib(mm(l))     = surf_def_v(l)%rib(m)
1947                   IF ( ALLOCATED( surf_def_v(l)%shf ) )                       &
1948                      surf_v(l)%shf(mm(l))     = surf_def_v(l)%shf(m)
1949                   IF ( ALLOCATED( surf_def_v(l)%qsws ) )                      &
1950                      surf_v(l)%qsws(mm(l))    = surf_def_v(l)%qsws(m)
1951                   IF ( ALLOCATED( surf_def_v(l)%ssws ) )                      &
1952                      surf_v(l)%ssws(mm(l))    = surf_def_v(l)%ssws(m)
1953                   IF ( ALLOCATED( surf_def_v(l)%ncsws ) )                     &
1954                      surf_v(l)%ncsws(mm(l))   = surf_def_v(l)%ncsws(m)
1955                   IF ( ALLOCATED( surf_def_v(l)%nrsws ) )                     &
1956                      surf_v(l)%nrsws(mm(l))   = surf_def_v(l)%nrsws(m)
1957                   IF ( ALLOCATED( surf_def_v(l)%sasws ) )                     &
1958                      surf_v(l)%sasws(mm(l))   = surf_def_v(l)%sasws(m)
1959                   IF ( ALLOCATED( surf_def_v(l)%mom_flux_uv) )                &
1960                      surf_v(l)%mom_flux_uv(mm(l))  = surf_def_v(l)%mom_flux_uv(m)
1961                   IF ( ALLOCATED( surf_def_v(l)%mom_flux_w) )                 &
1962                      surf_v(l)%mom_flux_w(mm(l))   = surf_def_v(l)%mom_flux_w(m)
1963                   IF ( ALLOCATED( surf_def_v(l)%mom_flux_tke) )               &
1964                      surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_def_v(l)%mom_flux_tke(0:1,m)
1965               
1966                   mm(l) = mm(l) + 1
1967                ENDDO
1968
1969                DO  m = surf_lsm_v(l)%start_index(j,i),                        &
1970                        surf_lsm_v(l)%end_index(j,i)
1971                   IF ( ALLOCATED( surf_lsm_v(l)%us ) )                        &
1972                      surf_v(l)%us(mm(l))      = surf_lsm_v(l)%us(m)
1973                   IF ( ALLOCATED( surf_lsm_v(l)%ts ) )                        &
1974                      surf_v(l)%ts(mm(l))      = surf_lsm_v(l)%ts(m)
1975                   IF ( ALLOCATED( surf_lsm_v(l)%qs ) )                        &
1976                      surf_v(l)%qs(mm(l))      = surf_lsm_v(l)%qs(m)
1977                   IF ( ALLOCATED( surf_lsm_v(l)%ss ) )                        &
1978                      surf_v(l)%ss(mm(l))      = surf_lsm_v(l)%ss(m)
1979                   IF ( ALLOCATED( surf_lsm_v(l)%qcs ) )                       &
1980                      surf_v(l)%qcs(mm(l))     = surf_lsm_v(l)%qcs(m)
1981                   IF ( ALLOCATED( surf_lsm_v(l)%ncs ) )                       &
1982                      surf_v(l)%ncs(mm(l))     = surf_lsm_v(l)%ncs(m)
1983                   IF ( ALLOCATED( surf_lsm_v(l)%qrs ) )                       &
1984                      surf_v(l)%qrs(mm(l))     = surf_lsm_v(l)%qrs(m)
1985                   IF ( ALLOCATED( surf_lsm_v(l)%nrs ) )                       &
1986                      surf_v(l)%nrs(mm(l))     = surf_lsm_v(l)%nrs(m)
1987                   IF ( ALLOCATED( surf_lsm_v(l)%ol ) )                        &
1988                      surf_v(l)%ol(mm(l))      = surf_lsm_v(l)%ol(m)
1989                   IF ( ALLOCATED( surf_lsm_v(l)%rib ) )                       &
1990                      surf_v(l)%rib(mm(l))     = surf_lsm_v(l)%rib(m)
1991                   IF ( ALLOCATED( surf_lsm_v(l)%usws ) )                      &
1992                      surf_v(l)%usws(mm(l))    = surf_lsm_v(l)%usws(m)
1993                   IF ( ALLOCATED( surf_lsm_v(l)%vsws ) )                      &
1994                      surf_v(l)%vsws(mm(l))    = surf_lsm_v(l)%vsws(m)
1995                   IF ( ALLOCATED( surf_lsm_v(l)%shf ) )                       &
1996                      surf_v(l)%shf(mm(l))     = surf_lsm_v(l)%shf(m)
1997                   IF ( ALLOCATED( surf_lsm_v(l)%qsws ) )                      &
1998                      surf_v(l)%qsws(mm(l))    = surf_lsm_v(l)%qsws(m)
1999                   IF ( ALLOCATED( surf_lsm_v(l)%ssws ) )                      &
2000                      surf_v(l)%ssws(mm(l))    = surf_lsm_v(l)%ssws(m)
2001                   IF ( ALLOCATED( surf_lsm_v(l)%ncsws ) )                     &
2002                      surf_v(l)%ncsws(mm(l))   = surf_lsm_v(l)%ncsws(m)
2003                   IF ( ALLOCATED( surf_lsm_v(l)%nrsws ) )                     &
2004                      surf_v(l)%nrsws(mm(l))   = surf_lsm_v(l)%nrsws(m)
2005                   IF ( ALLOCATED( surf_lsm_v(l)%sasws ) )                     &
2006                      surf_v(l)%sasws(mm(l))   = surf_lsm_v(l)%sasws(m)
2007                   IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_uv) )                &
2008                      surf_v(l)%mom_flux_uv(mm(l))  = surf_lsm_v(l)%mom_flux_uv(m)
2009                   IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_w) )                 &
2010                      surf_v(l)%mom_flux_w(mm(l))   = surf_lsm_v(l)%mom_flux_w(m)
2011                   IF ( ALLOCATED( surf_lsm_v(l)%mom_flux_tke) )               &
2012                      surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_lsm_v(l)%mom_flux_tke(0:1,m)
2013               
2014                   mm(l) = mm(l) + 1
2015                ENDDO
2016
2017                DO  m = surf_usm_v(l)%start_index(j,i),                        &
2018                        surf_usm_v(l)%end_index(j,i)
2019                   IF ( ALLOCATED( surf_usm_v(l)%us ) )                        &
2020                      surf_v(l)%us(mm(l))      = surf_usm_v(l)%us(m)
2021                   IF ( ALLOCATED( surf_usm_v(l)%ts ) )                        &
2022                      surf_v(l)%ts(mm(l))      = surf_usm_v(l)%ts(m)
2023                   IF ( ALLOCATED( surf_usm_v(l)%qs ) )                        &
2024                      surf_v(l)%qs(mm(l))      = surf_usm_v(l)%qs(m)
2025                   IF ( ALLOCATED( surf_usm_v(l)%ss ) )                        &
2026                      surf_v(l)%ss(mm(l))      = surf_usm_v(l)%ss(m)
2027                   IF ( ALLOCATED( surf_usm_v(l)%qcs ) )                       &
2028                      surf_v(l)%qcs(mm(l))     = surf_usm_v(l)%qcs(m)
2029                   IF ( ALLOCATED( surf_usm_v(l)%ncs ) )                       &
2030                      surf_v(l)%ncs(mm(l))     = surf_usm_v(l)%ncs(m)
2031                   IF ( ALLOCATED( surf_usm_v(l)%qrs ) )                       &
2032                      surf_v(l)%qrs(mm(l))     = surf_usm_v(l)%qrs(m)
2033                   IF ( ALLOCATED( surf_usm_v(l)%nrs ) )                       &
2034                      surf_v(l)%nrs(mm(l))     = surf_usm_v(l)%nrs(m)
2035                   IF ( ALLOCATED( surf_usm_v(l)%ol ) )                        &
2036                      surf_v(l)%ol(mm(l))      = surf_usm_v(l)%ol(m)
2037                   IF ( ALLOCATED( surf_usm_v(l)%rib ) )                       &
2038                      surf_v(l)%rib(mm(l))     = surf_usm_v(l)%rib(m)
2039                   IF ( ALLOCATED( surf_usm_v(l)%usws ) )                      &
2040                      surf_v(l)%usws(mm(l))    = surf_usm_v(l)%usws(m)
2041                   IF ( ALLOCATED( surf_usm_v(l)%vsws ) )                      &
2042                      surf_v(l)%vsws(mm(l))    = surf_usm_v(l)%vsws(m)
2043                   IF ( ALLOCATED( surf_usm_v(l)%shf ) )                       &
2044                      surf_v(l)%shf(mm(l))     = surf_usm_v(l)%shf(m)
2045                   IF ( ALLOCATED( surf_usm_v(l)%qsws ) )                      &
2046                      surf_v(l)%qsws(mm(l))    = surf_usm_v(l)%qsws(m)
2047                   IF ( ALLOCATED( surf_usm_v(l)%ssws ) )                      &
2048                      surf_v(l)%ssws(mm(l))    = surf_usm_v(l)%ssws(m)
2049                   IF ( ALLOCATED( surf_usm_v(l)%ncsws ) )                     &
2050                      surf_v(l)%ncsws(mm(l))   = surf_usm_v(l)%ncsws(m)
2051                   IF ( ALLOCATED( surf_usm_v(l)%nrsws ) )                     &
2052                      surf_v(l)%nrsws(mm(l))   = surf_usm_v(l)%nrsws(m)
2053                   IF ( ALLOCATED( surf_usm_v(l)%sasws ) )                     &
2054                      surf_v(l)%sasws(mm(l))   = surf_usm_v(l)%sasws(m)
2055                   IF ( ALLOCATED( surf_usm_v(l)%mom_flux_uv) )                &
2056                      surf_v(l)%mom_flux_uv(mm(l))  = surf_usm_v(l)%mom_flux_uv(m)
2057                   IF ( ALLOCATED( surf_usm_v(l)%mom_flux_w) )                 &
2058                      surf_v(l)%mom_flux_w(mm(l))   = surf_usm_v(l)%mom_flux_w(m)
2059                   IF ( ALLOCATED( surf_usm_v(l)%mom_flux_tke) )               &
2060                      surf_v(l)%mom_flux_tke(0:1,mm(l)) = surf_usm_v(l)%mom_flux_tke(0:1,m)
2061               
2062                   mm(l) = mm(l) + 1
2063                ENDDO
2064             
2065             ENDDO
2066          ENDDO
2067!
2068!--       Finally, determine start- and end-index for the respective surface
2069          surf_v(l)%start_index = MAX( surf_def_v(l)%start_index,              &
2070                                       surf_lsm_v(l)%start_index,              &
2071                                       surf_usm_v(l)%start_index )
2072          surf_v(l)%end_index   = MAX( surf_def_v(l)%end_index,                &
2073                                       surf_lsm_v(l)%end_index,                &
2074                                       surf_usm_v(l)%end_index   )
2075       ENDDO
2076
2077       WRITE ( 14 )  'ns_h_on_file                  '
2078       WRITE ( 14 )   ns_h_on_file
2079       WRITE ( 14 )  'ns_v_on_file                  '
2080       WRITE ( 14 )   ns_v_on_file
2081!
2082!--    Write required restart data.
2083!--    Start with horizontal surfaces (upward-, downward-facing, and model top)
2084       DO  l = 0, 2
2085          WRITE( dum, '(I1)')  l
2086         
2087          WRITE ( 14 )  'surf_h(' // dum // ')%start_index         ' 
2088          WRITE ( 14 )   surf_h(l)%start_index
2089          WRITE ( 14 )  'surf_h(' // dum // ')%end_index           ' 
2090          WRITE ( 14 )   surf_h(l)%end_index
2091
2092          WRITE ( 14 )  'surf_h(' // dum // ')%us                  ' 
2093          IF ( ALLOCATED ( surf_h(l)%us ) )  THEN
2094             WRITE ( 14 )  surf_h(l)%us
2095          ENDIF
2096          WRITE ( 14 )  'surf_h(' // dum // ')%ts                  ' 
2097          IF ( ALLOCATED ( surf_h(l)%ts ) )  THEN
2098             WRITE ( 14 )  surf_h(l)%ts
2099          ENDIF
2100          WRITE ( 14 )  'surf_h(' // dum // ')%qs                  ' 
2101          IF ( ALLOCATED ( surf_h(l)%qs ) )  THEN
2102             WRITE ( 14 )  surf_h(l)%qs
2103          ENDIF
2104          WRITE ( 14 )  'surf_h(' // dum // ')%ss                  ' 
2105          IF ( ALLOCATED ( surf_h(l)%ss ) )  THEN
2106             WRITE ( 14 )  surf_h(l)%ss
2107          ENDIF
2108          WRITE ( 14 )  'surf_h(' // dum // ')%qcs                 '
2109          IF ( ALLOCATED ( surf_h(l)%qcs ) )  THEN 
2110             WRITE ( 14 )  surf_h(l)%qcs
2111          ENDIF
2112          WRITE ( 14 )  'surf_h(' // dum // ')%ncs                 ' 
2113          IF ( ALLOCATED ( surf_h(l)%ncs ) )  THEN
2114             WRITE ( 14 )  surf_h(l)%ncs
2115          ENDIF
2116          WRITE ( 14 )  'surf_h(' // dum // ')%qrs                 '
2117          IF ( ALLOCATED ( surf_h(l)%qrs ) )  THEN 
2118             WRITE ( 14 )  surf_h(l)%qrs
2119          ENDIF
2120          WRITE ( 14 )  'surf_h(' // dum // ')%nrs                 ' 
2121          IF ( ALLOCATED ( surf_h(l)%nrs ) )  THEN
2122             WRITE ( 14 )  surf_h(l)%nrs
2123          ENDIF
2124          WRITE ( 14 )  'surf_h(' // dum // ')%ol                  ' 
2125          IF ( ALLOCATED ( surf_h(l)%ol ) )  THEN
2126             WRITE ( 14 )  surf_h(l)%ol
2127          ENDIF
2128          WRITE ( 14 )  'surf_h(' // dum // ')%rib                 ' 
2129          IF ( ALLOCATED ( surf_h(l)%rib ) )  THEN
2130             WRITE ( 14 )  surf_h(l)%rib
2131          ENDIF
2132          WRITE ( 14 )  'surf_h(' // dum // ')%usws                ' 
2133          IF ( ALLOCATED ( surf_h(l)%usws ) )  THEN
2134             WRITE ( 14 )  surf_h(l)%usws
2135          ENDIF
2136          WRITE ( 14 )  'surf_h(' // dum // ')%vsws                ' 
2137          IF ( ALLOCATED ( surf_h(l)%vsws ) )  THEN
2138             WRITE ( 14 )  surf_h(l)%vsws
2139          ENDIF
2140          WRITE ( 14 )  'surf_h(' // dum // ')%shf                 ' 
2141          IF ( ALLOCATED ( surf_h(l)%shf ) )  THEN
2142             WRITE ( 14 )  surf_h(l)%shf
2143          ENDIF
2144          WRITE ( 14 )  'surf_h(' // dum // ')%qsws                ' 
2145          IF ( ALLOCATED ( surf_h(l)%qsws ) )  THEN
2146             WRITE ( 14 )  surf_h(l)%qsws
2147          ENDIF
2148          WRITE ( 14 )  'surf_h(' // dum // ')%ssws                ' 
2149          IF ( ALLOCATED ( surf_h(l)%ssws ) )  THEN
2150             WRITE ( 14 )  surf_h(l)%ssws
2151          ENDIF
2152          WRITE ( 14 )  'surf_h(' // dum // ')%qcsws               ' 
2153          IF ( ALLOCATED ( surf_h(l)%qcsws ) )  THEN
2154             WRITE ( 14 )  surf_h(l)%qcsws
2155          ENDIF
2156          WRITE ( 14 )  'surf_h(' // dum // ')%ncsws               ' 
2157          IF ( ALLOCATED ( surf_h(l)%ncsws ) )  THEN
2158             WRITE ( 14 )  surf_h(l)%ncsws
2159          ENDIF
2160          WRITE ( 14 )  'surf_h(' // dum // ')%qrsws               ' 
2161          IF ( ALLOCATED ( surf_h(l)%qrsws ) )  THEN
2162             WRITE ( 14 )  surf_h(l)%qrsws
2163          ENDIF
2164          WRITE ( 14 )  'surf_h(' // dum // ')%nrsws               ' 
2165          IF ( ALLOCATED ( surf_h(l)%nrsws ) )  THEN
2166             WRITE ( 14 )  surf_h(l)%nrsws
2167          ENDIF
2168          WRITE ( 14 )  'surf_h(' // dum // ')%sasws               ' 
2169          IF ( ALLOCATED ( surf_h(l)%sasws ) )  THEN
2170             WRITE ( 14 )  surf_h(l)%sasws
2171          ENDIF
2172       ENDDO
2173!
2174!--    Write vertical surfaces
2175       DO  l = 0, 3
2176          WRITE( dum, '(I1)')  l
2177
2178          WRITE ( 14 )  'surf_v(' // dum // ')%start_index         ' 
2179          WRITE ( 14 )   surf_v(l)%start_index
2180          WRITE ( 14 )  'surf_v(' // dum // ')%end_index           ' 
2181          WRITE ( 14 )   surf_v(l)%end_index
2182
2183          WRITE ( 14 )  'surf_v(' // dum // ')%us                  ' 
2184          IF ( ALLOCATED ( surf_v(l)%us ) )  THEN
2185             WRITE ( 14 )  surf_v(l)%us
2186          ENDIF
2187          WRITE ( 14 )  'surf_v(' // dum // ')%ts                  ' 
2188          IF ( ALLOCATED ( surf_v(l)%ts ) )  THEN
2189             WRITE ( 14 )  surf_v(l)%ts
2190          ENDIF
2191          WRITE ( 14 )  'surf_v(' // dum // ')%qs                  ' 
2192          IF ( ALLOCATED ( surf_v(l)%qs ) )  THEN
2193             WRITE ( 14 )  surf_v(l)%qs
2194          ENDIF
2195          WRITE ( 14 )  'surf_v(' // dum // ')%ss                  ' 
2196          IF ( ALLOCATED ( surf_v(l)%ss ) )  THEN
2197             WRITE ( 14 )  surf_v(l)%ss
2198          ENDIF
2199          WRITE ( 14 )  'surf_v(' // dum // ')%qcs                 ' 
2200          IF ( ALLOCATED ( surf_v(l)%qcs ) )  THEN
2201             WRITE ( 14 )  surf_v(l)%qcs
2202          ENDIF
2203          WRITE ( 14 )  'surf_v(' // dum // ')%ncs                 ' 
2204          IF ( ALLOCATED ( surf_v(l)%ncs ) )  THEN
2205             WRITE ( 14 )  surf_v(l)%ncs
2206          ENDIF
2207          WRITE ( 14 )  'surf_v(' // dum // ')%qrs                 ' 
2208          IF ( ALLOCATED ( surf_v(l)%qrs ) )  THEN
2209             WRITE ( 14 )  surf_v(l)%qrs
2210          ENDIF
2211          WRITE ( 14 )  'surf_v(' // dum // ')%nrs                 ' 
2212          IF ( ALLOCATED ( surf_v(l)%nrs ) )  THEN
2213             WRITE ( 14 )  surf_v(l)%nrs
2214          ENDIF
2215          WRITE ( 14 )  'surf_v(' // dum // ')%ol                  ' 
2216          IF ( ALLOCATED ( surf_v(l)%ol ) )  THEN
2217             WRITE ( 14 )  surf_v(l)%ol
2218          ENDIF
2219          WRITE ( 14 )  'surf_v(' // dum // ')%rib                 ' 
2220          IF ( ALLOCATED ( surf_v(l)%rib ) )  THEN
2221             WRITE ( 14 )  surf_v(l)%rib
2222          ENDIF
2223          WRITE ( 14 )  'surf_v(' // dum // ')%shf                 ' 
2224          IF ( ALLOCATED ( surf_v(l)%shf ) )  THEN
2225             WRITE ( 14 )  surf_v(l)%shf
2226          ENDIF
2227          WRITE ( 14 )  'surf_v(' // dum // ')%qsws                ' 
2228          IF ( ALLOCATED ( surf_v(l)%qsws ) )  THEN
2229             WRITE ( 14 )  surf_v(l)%qsws
2230          ENDIF
2231          WRITE ( 14 )  'surf_v(' // dum // ')%ssws                ' 
2232          IF ( ALLOCATED ( surf_v(l)%ssws ) )  THEN
2233             WRITE ( 14 )  surf_v(l)%ssws
2234          ENDIF
2235          WRITE ( 14 )  'surf_v(' // dum // ')%qcsws               ' 
2236          IF ( ALLOCATED ( surf_v(l)%qcsws ) )  THEN
2237             WRITE ( 14 )  surf_v(l)%qcsws
2238          ENDIF
2239          WRITE ( 14 )  'surf_v(' // dum // ')%ncsws               ' 
2240          IF ( ALLOCATED ( surf_v(l)%ncsws ) )  THEN
2241             WRITE ( 14 )  surf_v(l)%ncsws
2242          ENDIF
2243          WRITE ( 14 )  'surf_v(' // dum // ')%qrsws               ' 
2244          IF ( ALLOCATED ( surf_v(l)%qrsws ) )  THEN
2245             WRITE ( 14 )  surf_v(l)%qrsws
2246          ENDIF
2247          WRITE ( 14 )  'surf_v(' // dum // ')%nrsws               ' 
2248          IF ( ALLOCATED ( surf_v(l)%nrsws ) )  THEN
2249             WRITE ( 14 )  surf_v(l)%nrsws
2250          ENDIF
2251          WRITE ( 14 )  'surf_v(' // dum // ')%sasws               ' 
2252          IF ( ALLOCATED ( surf_v(l)%sasws ) )  THEN
2253             WRITE ( 14 )  surf_v(l)%sasws
2254          ENDIF
2255          WRITE ( 14 )  'surf_v(' // dum // ')%mom_uv              ' 
2256          IF ( ALLOCATED ( surf_v(l)%mom_flux_uv ) )  THEN
2257             WRITE ( 14 )  surf_v(l)%mom_flux_uv
2258          ENDIF
2259          WRITE ( 14 )  'surf_v(' // dum // ')%mom_w               ' 
2260          IF ( ALLOCATED ( surf_v(l)%mom_flux_w ) )  THEN
2261             WRITE ( 14 )  surf_v(l)%mom_flux_w
2262          ENDIF
2263          WRITE ( 14 )  'surf_v(' // dum // ')%mom_tke             ' 
2264          IF ( ALLOCATED ( surf_v(l)%mom_flux_tke ) )  THEN
2265             WRITE ( 14 )  surf_v(l)%mom_flux_tke
2266          ENDIF
2267
2268       ENDDO
2269
2270       WRITE ( 14 )  '*** end surf ***              '
2271
2272    END SUBROUTINE surface_write_restart_data
2273
2274
2275!------------------------------------------------------------------------------!
2276! Description:
2277! ------------
2278!> Reads surface-related restart data. Please note, restart data for a certain
2279!> surface orientation (e.g. horizontal upward-facing) is stored in one
2280!> array, even if surface elements may belong to different surface types
2281!> natural or urban for example). Surface elements are redistributed into its
2282!> respective surface types within this routine. This allows e.g. changing the
2283!> surface type after reading the restart data, which might be required in case
2284!> of cyclic_fill mode.
2285!------------------------------------------------------------------------------!
2286    SUBROUTINE surface_read_restart_data( ii,                                  &
2287                                       nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
2288                                       nynfa, nyn_on_file, nysfa, nys_on_file, &
2289                                       offset_xa, offset_ya, overlap_count )
2290
2291       USE pegrid,                                                             &
2292           ONLY: numprocs_previous_run
2293
2294       CHARACTER (LEN=1)  ::  dum         !< dummy to create correct string for reading input variable
2295       CHARACTER (LEN=30) ::  field_chr   !< input variable
2296
2297       INTEGER(iwp)       ::  i           !< running index along x-direction, refers to former domain size
2298       INTEGER(iwp)       ::  ic          !< running index along x-direction, refers to current domain size
2299       INTEGER(iwp)       ::  j           !< running index along y-direction, refers to former domain size
2300       INTEGER(iwp)       ::  jc          !< running index along y-direction, refers to former domain size
2301       INTEGER(iwp)       ::  k           !< running index along z-direction
2302       INTEGER(iwp)       ::  l           !< index variable for surface type
2303       INTEGER(iwp)       ::  m           !< running index for surface elements, refers to gathered array encompassing all surface types
2304       INTEGER(iwp)       ::  mm          !< running index for surface elements, refers to individual surface types
2305
2306       INTEGER(iwp)       ::  ii               !< running index over input files
2307       INTEGER(iwp)       ::  kk               !< running index over previous input files covering current local domain
2308       INTEGER(iwp)       ::  nxlc             !< index of left boundary on current subdomain
2309       INTEGER(iwp)       ::  nxlf             !< index of left boundary on former subdomain
2310       INTEGER(iwp)       ::  nxl_on_file      !< index of left boundary on former local domain
2311       INTEGER(iwp)       ::  nxrc             !< index of right boundary on current subdomain
2312       INTEGER(iwp)       ::  nxrf             !< index of right boundary on former subdomain
2313       INTEGER(iwp)       ::  nxr_on_file      !< index of right boundary on former local domain 
2314       INTEGER(iwp)       ::  nync             !< index of north boundary on current subdomain
2315       INTEGER(iwp)       ::  nynf             !< index of north boundary on former subdomain
2316       INTEGER(iwp)       ::  nyn_on_file      !< index of norht boundary on former local domain 
2317       INTEGER(iwp)       ::  nysc             !< index of south boundary on current subdomain
2318       INTEGER(iwp)       ::  nysf             !< index of south boundary on former subdomain
2319       INTEGER(iwp)       ::  nys_on_file      !< index of south boundary on former local domain 
2320       INTEGER(iwp)       ::  overlap_count    !< number of overlaps
2321 
2322       INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
2323       INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !<
2324       INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !<
2325       INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !<
2326       INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<
2327       INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
2328
2329
2330       LOGICAL                         ::  horizontal_surface !< flag indicating horizontal surfaces
2331       LOGICAL                         ::  surf_match_def     !< flag indicating that surface element is of default type
2332       LOGICAL                         ::  surf_match_lsm     !< flag indicating that surface element is of natural type
2333       LOGICAL                         ::  surf_match_usm     !< flag indicating that surface element is of urban type
2334       LOGICAL                         ::  vertical_surface   !< flag indicating vertical surfaces
2335
2336       TYPE(surf_type), DIMENSION(0:2) ::  surf_h             !< horizontal surface type on file
2337       TYPE(surf_type), DIMENSION(0:3) ::  surf_v             !< vertical surface type on file
2338
2339!
2340!--    Read number of respective surface elements on file
2341       READ ( 13 )  field_chr
2342       IF ( TRIM( field_chr ) /= 'ns_h_on_file' )  THEN
2343!
2344!--       Add a proper error message
2345       ENDIF
2346       READ ( 13 ) ns_h_on_file
2347
2348       READ ( 13 )  field_chr
2349       IF ( TRIM( field_chr ) /= 'ns_v_on_file' )  THEN
2350!
2351!--       Add a proper error message
2352       ENDIF
2353       READ ( 13 ) ns_v_on_file
2354!
2355!--    Allocate memory for number of surface elements on file. Please note,
2356!--    these number is not necessarily the same as the final number of surface
2357!--    elements on local domain, which is the case if processor topology changes
2358!--    during restart runs.
2359!--    Horizontal upward facing
2360       surf_h(0)%ns = ns_h_on_file(0)
2361       CALL allocate_surface_attributes_h( surf_h(0),                          &
2362                                           nys_on_file, nyn_on_file,           &
2363                                           nxl_on_file, nxr_on_file )
2364!
2365!--    Horizontal downward facing
2366       surf_h(1)%ns = ns_h_on_file(1)
2367       CALL allocate_surface_attributes_h( surf_h(1),                          &
2368                                           nys_on_file, nyn_on_file,           &
2369                                           nxl_on_file, nxr_on_file )
2370!
2371!--    Model top
2372       surf_h(2)%ns = ns_h_on_file(2)
2373       CALL allocate_surface_attributes_h_top( surf_h(2),                      &
2374                                               nys_on_file, nyn_on_file,       &
2375                                               nxl_on_file, nxr_on_file )
2376!
2377!--    Vertical surfaces
2378       DO  l = 0, 3
2379          surf_v(l)%ns = ns_v_on_file(l)
2380          CALL allocate_surface_attributes_v( surf_v(l), .FALSE.,              &
2381                                              nys_on_file, nyn_on_file,        &
2382                                              nxl_on_file, nxr_on_file )
2383       ENDDO
2384
2385       IF ( initializing_actions == 'read_restart_data'  .OR.                  &
2386            initializing_actions == 'cyclic_fill' )  THEN
2387!
2388!--       Initial setting of flags for horizontal and vertical surfaces, will
2389!--       be set after start- and end-indices are read.
2390          horizontal_surface = .FALSE.
2391          vertical_surface   = .FALSE.
2392
2393          READ ( 13 )  field_chr
2394
2395          DO  WHILE ( TRIM( field_chr ) /= '*** end surf ***' )
2396!
2397!--          Map data on file as often as needed (data are read only for k=1)
2398             DO  kk = 1, overlap_count
2399!
2400!--             Get the index range of the subdomain on file which overlap with the
2401!--             current subdomain
2402                nxlf = nxlfa(ii,kk)
2403                nxlc = nxlfa(ii,kk) + offset_xa(ii,kk)
2404                nxrf = nxrfa(ii,kk)
2405                nxrc = nxrfa(ii,kk) + offset_xa(ii,kk)
2406                nysf = nysfa(ii,kk)
2407                nysc = nysfa(ii,kk) + offset_ya(ii,kk)
2408                nynf = nynfa(ii,kk)
2409                nync = nynfa(ii,kk) + offset_ya(ii,kk)
2410
2411                SELECT CASE ( TRIM( field_chr ) )
2412
2413                   CASE ( 'surf_h(0)%start_index' )
2414                      IF ( kk == 1 )                                           &
2415                         READ ( 13 )  surf_h(0)%start_index
2416                      l = 0
2417                   CASE ( 'surf_h(0)%end_index' )   
2418                      IF ( kk == 1 )                                           &
2419                         READ ( 13 )  surf_h(0)%end_index
2420                      horizontal_surface = .TRUE.
2421                      vertical_surface   = .FALSE.
2422                   CASE ( 'surf_h(0)%us' )         
2423                      IF ( ALLOCATED( surf_h(0)%us )  .AND.  kk == 1 )         &
2424                         READ ( 13 )  surf_h(0)%us
2425                   CASE ( 'surf_h(0)%ts' )         
2426                      IF ( ALLOCATED( surf_h(0)%ts )  .AND.  kk == 1 )         &
2427                         READ ( 13 )  surf_h(0)%ts
2428                   CASE ( 'surf_h(0)%qs' )         
2429                      IF ( ALLOCATED( surf_h(0)%qs )  .AND.  kk == 1 )         &
2430                         READ ( 13 )  surf_h(0)%qs
2431                   CASE ( 'surf_h(0)%ss' )         
2432                      IF ( ALLOCATED( surf_h(0)%ss )  .AND.  kk == 1 )         &
2433                         READ ( 13 )  surf_h(0)%ss
2434                   CASE ( 'surf_h(0)%qcs' )         
2435                      IF ( ALLOCATED( surf_h(0)%qcs )  .AND.  kk == 1 )        &
2436                         READ ( 13 )  surf_h(0)%qcs
2437                   CASE ( 'surf_h(0)%ncs' )         
2438                      IF ( ALLOCATED( surf_h(0)%ncs )  .AND.  kk == 1 )        &
2439                         READ ( 13 )  surf_h(0)%ncs
2440                   CASE ( 'surf_h(0)%qrs' )         
2441                      IF ( ALLOCATED( surf_h(0)%qrs )  .AND.  kk == 1 )        &
2442                         READ ( 13 )  surf_h(0)%qrs
2443                   CASE ( 'surf_h(0)%nrs' )         
2444                      IF ( ALLOCATED( surf_h(0)%nrs )  .AND.  kk == 1 )        &
2445                         READ ( 13 )  surf_h(0)%nrs
2446                   CASE ( 'surf_h(0)%ol' )         
2447                      IF ( ALLOCATED( surf_h(0)%ol )  .AND.  kk == 1 )         &
2448                         READ ( 13 )  surf_h(0)%ol
2449                   CASE ( 'surf_h(0)%rib' )         
2450                      IF ( ALLOCATED( surf_h(0)%rib )  .AND.  kk == 1 )        &
2451                         READ ( 13 )  surf_h(0)%rib
2452                   CASE ( 'surf_h(0)%usws' )         
2453                      IF ( ALLOCATED( surf_h(0)%usws )  .AND.  kk == 1 )       &
2454                         READ ( 13 )  surf_h(0)%usws
2455                   CASE ( 'surf_h(0)%vsws' )         
2456                      IF ( ALLOCATED( surf_h(0)%vsws )  .AND.  kk == 1 )       &
2457                         READ ( 13 )  surf_h(0)%vsws
2458                   CASE ( 'surf_h(0)%shf' )         
2459                      IF ( ALLOCATED( surf_h(0)%shf )  .AND.  kk == 1 )        &
2460                         READ ( 13 )  surf_h(0)%shf
2461                   CASE ( 'surf_h(0)%qsws' )         
2462                      IF ( ALLOCATED( surf_h(0)%qsws )  .AND.  kk == 1 )       &
2463                         READ ( 13 )  surf_h(0)%qsws
2464                   CASE ( 'surf_h(0)%ssws' )         
2465                      IF ( ALLOCATED( surf_h(0)%ssws )  .AND.  kk == 1 )       &
2466                         READ ( 13 )  surf_h(0)%ssws
2467                   CASE ( 'surf_h(0)%qcsws' )         
2468                      IF ( ALLOCATED( surf_h(0)%qcsws )  .AND.  kk == 1 )      &
2469                         READ ( 13 )  surf_h(0)%qcsws
2470                   CASE ( 'surf_h(0)%ncsws' )         
2471                      IF ( ALLOCATED( surf_h(0)%ncsws )  .AND.  kk == 1 )      &
2472                         READ ( 13 )  surf_h(0)%ncsws
2473                   CASE ( 'surf_h(0)%qrsws' )         
2474                      IF ( ALLOCATED( surf_h(0)%qrsws )  .AND.  kk == 1 )      &
2475                         READ ( 13 )  surf_h(0)%qrsws
2476                   CASE ( 'surf_h(0)%nrsws' )         
2477                      IF ( ALLOCATED( surf_h(0)%nrsws )  .AND.  kk == 1 )      &
2478                         READ ( 13 )  surf_h(0)%nrsws
2479                   CASE ( 'surf_h(0)%sasws' )         
2480                      IF ( ALLOCATED( surf_h(0)%sasws )  .AND.  kk == 1 )      &
2481                         READ ( 13 )  surf_h(0)%sasws
2482
2483                   CASE ( 'surf_h(1)%start_index' )   
2484                      IF ( kk == 1 )                                           &
2485                         READ ( 13 )  surf_h(1)%start_index
2486                      l = 1
2487                   CASE ( 'surf_h(1)%end_index' )   
2488                      IF ( kk == 1 )                                           &
2489                         READ ( 13 )  surf_h(1)%end_index
2490                   CASE ( 'surf_h(1)%us' )         
2491                      IF ( ALLOCATED( surf_h(1)%us )  .AND.  kk == 1 )         &
2492                         READ ( 13 )  surf_h(1)%us
2493                   CASE ( 'surf_h(1)%ts' )         
2494                      IF ( ALLOCATED( surf_h(1)%ts )  .AND.  kk == 1 )         &
2495                         READ ( 13 )  surf_h(1)%ts
2496                   CASE ( 'surf_h(1)%qs' )         
2497                      IF ( ALLOCATED( surf_h(1)%qs )  .AND.  kk == 1 )         &
2498                         READ ( 13 )  surf_h(1)%qs
2499                   CASE ( 'surf_h(1)%ss' )         
2500                      IF ( ALLOCATED( surf_h(1)%ss )  .AND.  kk == 1 )         &
2501                         READ ( 13 )  surf_h(1)%ss
2502                   CASE ( 'surf_h(1)%qcs' )         
2503                      IF ( ALLOCATED( surf_h(1)%qcs )  .AND.  kk == 1 )        &
2504                         READ ( 13 )  surf_h(1)%qcs
2505                   CASE ( 'surf_h(1)%ncs' )         
2506                      IF ( ALLOCATED( surf_h(1)%ncs )  .AND.  kk == 1 )        &
2507                         READ ( 13 )  surf_h(1)%ncs
2508                   CASE ( 'surf_h(1)%qrs' )         
2509                      IF ( ALLOCATED( surf_h(1)%qrs )  .AND.  kk == 1 )        &
2510                         READ ( 13 )  surf_h(1)%qrs
2511                   CASE ( 'surf_h(1)%nrs' )         
2512                      IF ( ALLOCATED( surf_h(1)%nrs )  .AND.  kk == 1 )        &
2513                         READ ( 13 )  surf_h(1)%nrs
2514                   CASE ( 'surf_h(1)%ol' )         
2515                      IF ( ALLOCATED( surf_h(1)%ol )  .AND.  kk == 1 )         &
2516                         READ ( 13 )  surf_h(1)%ol
2517                   CASE ( 'surf_h(1)%rib' )         
2518                      IF ( ALLOCATED( surf_h(1)%rib )  .AND.  kk == 1 )        &
2519                         READ ( 13 )  surf_h(1)%rib
2520                   CASE ( 'surf_h(1)%usws' )         
2521                      IF ( ALLOCATED( surf_h(1)%usws )  .AND.  kk == 1 )       &
2522                         READ ( 13 )  surf_h(1)%usws
2523                   CASE ( 'surf_h(1)%vsws' )         
2524                      IF ( ALLOCATED( surf_h(1)%vsws )  .AND.  kk == 1 )       &
2525                         READ ( 13 )  surf_h(1)%vsws
2526                   CASE ( 'surf_h(1)%shf' )         
2527                      IF ( ALLOCATED( surf_h(1)%shf )  .AND.  kk == 1 )        &
2528                         READ ( 13 )  surf_h(1)%shf
2529                   CASE ( 'surf_h(1)%qsws' )         
2530                      IF ( ALLOCATED( surf_h(1)%qsws )  .AND.  kk == 1 )       &
2531                         READ ( 13 )  surf_h(1)%qsws
2532                   CASE ( 'surf_h(1)%ssws' )         
2533                      IF ( ALLOCATED( surf_h(1)%ssws )  .AND.  kk == 1 )       &
2534                         READ ( 13 )  surf_h(1)%ssws
2535                   CASE ( 'surf_h(1)%qcsws' )         
2536                      IF ( ALLOCATED( surf_h(1)%qcsws )  .AND.  kk == 1 )      &
2537                         READ ( 13 )  surf_h(1)%qcsws
2538                   CASE ( 'surf_h(1)%ncsws' )         
2539                      IF ( ALLOCATED( surf_h(1)%ncsws )  .AND.  kk == 1 )      &
2540                         READ ( 13 )  surf_h(1)%ncsws
2541                   CASE ( 'surf_h(1)%qrsws' )         
2542                      IF ( ALLOCATED( surf_h(1)%qrsws )  .AND.  kk == 1 )      &
2543                         READ ( 13 )  surf_h(1)%qrsws
2544                   CASE ( 'surf_h(1)%nrsws' )         
2545                      IF ( ALLOCATED( surf_h(1)%nrsws )  .AND.  kk == 1 )      &
2546                         READ ( 13 )  surf_h(1)%nrsws
2547                   CASE ( 'surf_h(1)%sasws' )         
2548                      IF ( ALLOCATED( surf_h(1)%sasws )  .AND.  kk == 1 )      &
2549                         READ ( 13 )  surf_h(1)%sasws
2550
2551                   CASE ( 'surf_h(2)%start_index' )   
2552                      IF ( kk == 1 )                                           &
2553                         READ ( 13 )  surf_h(2)%start_index
2554                      l = 2
2555                   CASE ( 'surf_h(2)%end_index' )   
2556                      IF ( kk == 1 )                                           &
2557                         READ ( 13 )  surf_h(2)%end_index
2558                   CASE ( 'surf_h(2)%us' )         
2559                      IF ( ALLOCATED( surf_h(2)%us )  .AND.  kk == 1 )         &
2560                         READ ( 13 )  surf_h(2)%us
2561                   CASE ( 'surf_h(2)%ts' )         
2562                      IF ( ALLOCATED( surf_h(2)%ts )  .AND.  kk == 1 )         &
2563                         READ ( 13 )  surf_h(2)%ts
2564                   CASE ( 'surf_h(2)%qs' )       
2565                      IF ( ALLOCATED( surf_h(2)%qs )  .AND.  kk == 1 )         &
2566                         READ ( 13 )  surf_h(2)%qs
2567                   CASE ( 'surf_h(2)%ss' )         
2568                      IF ( ALLOCATED( surf_h(2)%ss )  .AND.  kk == 1 )         &
2569                         READ ( 13 )  surf_h(2)%ss
2570                   CASE ( 'surf_h(2)%qcs' )         
2571                      IF ( ALLOCATED( surf_h(2)%qcs )  .AND.  kk == 1 )        &
2572                         READ ( 13 )  surf_h(2)%qcs
2573                   CASE ( 'surf_h(2)%ncs' )         
2574                      IF ( ALLOCATED( surf_h(2)%ncs )  .AND.  kk == 1 )        &
2575                         READ ( 13 )  surf_h(2)%ncs
2576                   CASE ( 'surf_h(2)%qrs' )         
2577                      IF ( ALLOCATED( surf_h(2)%qrs )  .AND.  kk == 1 )        &
2578                         READ ( 13 )  surf_h(2)%qrs
2579                   CASE ( 'surf_h(2)%nrs' )         
2580                      IF ( ALLOCATED( surf_h(2)%nrs )  .AND.  kk == 1 )        &
2581                         READ ( 13 )  surf_h(2)%nrs
2582                   CASE ( 'surf_h(2)%ol' )         
2583                      IF ( ALLOCATED( surf_h(2)%ol )  .AND.  kk == 1 )         &
2584                         READ ( 13 )  surf_h(2)%ol
2585                   CASE ( 'surf_h(2)%rib' )         
2586                      IF ( ALLOCATED( surf_h(2)%rib )  .AND.  kk == 1 )        &
2587                         READ ( 13 )  surf_h(2)%rib
2588                   CASE ( 'surf_h(2)%usws' )         
2589                      IF ( ALLOCATED( surf_h(2)%usws )  .AND.  kk == 1 )       &
2590                         READ ( 13 )  surf_h(2)%usws
2591                   CASE ( 'surf_h(2)%vsws' )         
2592                      IF ( ALLOCATED( surf_h(2)%vsws )  .AND.  kk == 1 )       &
2593                         READ ( 13 )  surf_h(2)%vsws
2594                   CASE ( 'surf_h(2)%shf' )         
2595                      IF ( ALLOCATED( surf_h(2)%shf )  .AND.  kk == 1 )        &
2596                         READ ( 13 )  surf_h(2)%shf
2597                   CASE ( 'surf_h(2)%qsws' )         
2598                      IF ( ALLOCATED( surf_h(2)%qsws )  .AND.  kk == 1 )       &
2599                         READ ( 13 )  surf_h(2)%qsws
2600                   CASE ( 'surf_h(2)%ssws' )         
2601                      IF ( ALLOCATED( surf_h(2)%ssws )  .AND.  kk == 1 )       &
2602                         READ ( 13 )  surf_h(2)%ssws
2603                   CASE ( 'surf_h(2)%qcsws' )         
2604                      IF ( ALLOCATED( surf_h(2)%qcsws )  .AND.  kk == 1 )      &
2605                         READ ( 13 )  surf_h(2)%qcsws
2606                   CASE ( 'surf_h(2)%ncsws' )         
2607                      IF ( ALLOCATED( surf_h(2)%ncsws )  .AND.  kk == 1 )      &
2608                         READ ( 13 )  surf_h(2)%ncsws
2609                   CASE ( 'surf_h(2)%qrsws' )         
2610                      IF ( ALLOCATED( surf_h(2)%qrsws )  .AND.  kk == 1 )      &
2611                         READ ( 13 )  surf_h(2)%qrsws
2612                   CASE ( 'surf_h(2)%nrsws' )         
2613                      IF ( ALLOCATED( surf_h(2)%nrsws )  .AND.  kk == 1 )      &
2614                         READ ( 13 )  surf_h(2)%nrsws
2615                   CASE ( 'surf_h(2)%sasws' )         
2616                      IF ( ALLOCATED( surf_h(2)%sasws )  .AND.  kk == 1 )      &
2617                         READ ( 13 )  surf_h(2)%sasws
2618
2619                   CASE ( 'surf_v(0)%start_index' )   
2620                      IF ( kk == 1 )                                           &
2621                         READ ( 13 )  surf_v(0)%start_index
2622                      l = 0
2623                      horizontal_surface = .FALSE.
2624                      vertical_surface   = .TRUE.
2625                   CASE ( 'surf_v(0)%end_index' )   
2626                      IF ( kk == 1 )                                           &
2627                         READ ( 13 )  surf_v(0)%end_index
2628                   CASE ( 'surf_v(0)%us' )         
2629                      IF ( ALLOCATED( surf_v(0)%us )  .AND.  kk == 1 )         &
2630                         READ ( 13 )  surf_v(0)%us
2631                   CASE ( 'surf_v(0)%ts' )         
2632                      IF ( ALLOCATED( surf_v(0)%ts )  .AND.  kk == 1 )         &
2633                         READ ( 13 )  surf_v(0)%ts
2634                   CASE ( 'surf_v(0)%qs' )         
2635                      IF ( ALLOCATED( surf_v(0)%qs )  .AND.  kk == 1 )         &
2636                         READ ( 13 )  surf_v(0)%qs
2637                   CASE ( 'surf_v(0)%ss' )         
2638                      IF ( ALLOCATED( surf_v(0)%ss )  .AND.  kk == 1 )         &
2639                         READ ( 13 )  surf_v(0)%ss
2640                   CASE ( 'surf_v(0)%qcs' )         
2641                      IF ( ALLOCATED( surf_v(0)%qcs )  .AND.  kk == 1 )        &
2642                         READ ( 13 )  surf_v(0)%qcs
2643                   CASE ( 'surf_v(0)%ncs' )         
2644                      IF ( ALLOCATED( surf_v(0)%ncs )  .AND.  kk == 1 )        &
2645                         READ ( 13 )  surf_v(0)%ncs
2646                   CASE ( 'surf_v(0)%qrs' )         
2647                      IF ( ALLOCATED( surf_v(0)%qrs )  .AND.  kk == 1 )        &
2648                         READ ( 13 )  surf_v(0)%qrs
2649                   CASE ( 'surf_v(0)%nrs' )         
2650                      IF ( ALLOCATED( surf_v(0)%nrs )  .AND.  kk == 1 )        &
2651                         READ ( 13 )  surf_v(0)%nrs
2652                   CASE ( 'surf_v(0)%ol' )         
2653                      IF ( ALLOCATED( surf_v(0)%ol )  .AND.  kk == 1 )         &
2654                         READ ( 13 )  surf_v(0)%ol
2655                   CASE ( 'surf_v(0)%rib' )         
2656                      IF ( ALLOCATED( surf_v(0)%rib )  .AND.  kk == 1 )        &
2657                         READ ( 13 )  surf_v(0)%rib
2658                   CASE ( 'surf_v(0)%shf' )         
2659                      IF ( ALLOCATED( surf_v(0)%shf )  .AND.  kk == 1 )        &
2660                         READ ( 13 )  surf_v(0)%shf
2661                   CASE ( 'surf_v(0)%qsws' )         
2662                      IF ( ALLOCATED( surf_v(0)%qsws )  .AND.  kk == 1 )       &
2663                         READ ( 13 )  surf_v(0)%qsws
2664                   CASE ( 'surf_v(0)%ssws' )         
2665                      IF ( ALLOCATED( surf_v(0)%ssws )  .AND.  kk == 1 )       &
2666                         READ ( 13 )  surf_v(0)%ssws
2667                   CASE ( 'surf_v(0)%qcsws' )         
2668                      IF ( ALLOCATED( surf_v(0)%qcsws )  .AND.  kk == 1 )      &
2669                         READ ( 13 )  surf_v(0)%qcsws
2670                   CASE ( 'surf_v(0)%ncsws' )         
2671                      IF ( ALLOCATED( surf_v(0)%ncsws )  .AND.  kk == 1 )      &
2672                         READ ( 13 )  surf_v(0)%ncsws
2673                   CASE ( 'surf_v(0)%qrsws' )         
2674                      IF ( ALLOCATED( surf_v(0)%qrsws )  .AND.  kk == 1 )      &
2675                         READ ( 13 )  surf_v(0)%qrsws
2676                   CASE ( 'surf_v(0)%nrsws' )         
2677                      IF ( ALLOCATED( surf_v(0)%nrsws )  .AND.  kk == 1 )      &
2678                         READ ( 13 )  surf_v(0)%nrsws
2679                   CASE ( 'surf_v(0)%sasws' )         
2680                      IF ( ALLOCATED( surf_v(0)%sasws )  .AND.  kk == 1 )      &
2681                         READ ( 13 )  surf_v(0)%sasws
2682                   CASE ( 'surf_v(0)%mom_uv' )         
2683                      IF ( ALLOCATED( surf_v(0)%mom_flux_uv )  .AND.  kk == 1 )&
2684                         READ ( 13 )  surf_v(0)%mom_flux_uv
2685                   CASE ( 'surf_v(0)%mom_w' )         
2686                      IF ( ALLOCATED( surf_v(0)%mom_flux_w )  .AND.  kk == 1 ) &
2687                         READ ( 13 )  surf_v(0)%mom_flux_w
2688                   CASE ( 'surf_v(0)%mom_tke' )         
2689                      IF ( ALLOCATED( surf_v(0)%mom_flux_tke )  .AND.  kk == 1 )&
2690                         READ ( 13 )  surf_v(0)%mom_flux_tke
2691
2692                   CASE ( 'surf_v(1)%start_index' )   
2693                      IF ( kk == 1 )                                           &
2694                         READ ( 13 )  surf_v(1)%start_index
2695                      l = 1
2696                   CASE ( 'surf_v(1)%end_index' )   
2697                      IF ( kk == 1 )                                           &
2698                         READ ( 13 )  surf_v(1)%end_index
2699                   CASE ( 'surf_v(1)%us' )         
2700                      IF ( ALLOCATED( surf_v(1)%us )  .AND.  kk == 1 )         &
2701                         READ ( 13 )  surf_v(1)%us
2702                   CASE ( 'surf_v(1)%ts' )         
2703                      IF ( ALLOCATED( surf_v(1)%ts )  .AND.  kk == 1 )         &
2704                         READ ( 13 )  surf_v(1)%ts
2705                   CASE ( 'surf_v(1)%qs' )         
2706                      IF ( ALLOCATED( surf_v(1)%qs )  .AND.  kk == 1 )         &
2707                         READ ( 13 )  surf_v(1)%qs
2708                   CASE ( 'surf_v(1)%ss' )         
2709                      IF ( ALLOCATED( surf_v(1)%ss )  .AND.  kk == 1 )         &
2710                         READ ( 13 )  surf_v(1)%ss
2711                   CASE ( 'surf_v(1)%qcs' )         
2712                      IF ( ALLOCATED( surf_v(1)%qcs )  .AND.  kk == 1 )        &
2713                         READ ( 13 )  surf_v(1)%qcs
2714                   CASE ( 'surf_v(1)%ncs' )         
2715                      IF ( ALLOCATED( surf_v(1)%ncs )  .AND.  kk == 1 )        &
2716                         READ ( 13 )  surf_v(1)%ncs
2717                   CASE ( 'surf_v(1)%qrs' )         
2718                      IF ( ALLOCATED( surf_v(1)%qrs )  .AND.  kk == 1 )        &
2719                         READ ( 13 )  surf_v(1)%qrs
2720                   CASE ( 'surf_v(1)%nrs' )         
2721                      IF ( ALLOCATED( surf_v(1)%nrs )  .AND.  kk == 1 )        &
2722                         READ ( 13 )  surf_v(1)%nrs
2723                   CASE ( 'surf_v(1)%ol' )         
2724                      IF ( ALLOCATED( surf_v(1)%ol )  .AND.  kk == 1 )         &
2725                         READ ( 13 )  surf_v(1)%ol
2726                   CASE ( 'surf_v(1)%rib' )         
2727                      IF ( ALLOCATED( surf_v(1)%rib )  .AND.  kk == 1 )        &
2728                         READ ( 13 )  surf_v(1)%rib
2729                   CASE ( 'surf_v(1)%shf' )         
2730                      IF ( ALLOCATED( surf_v(1)%shf )  .AND.  kk == 1 )        &
2731                         READ ( 13 )  surf_v(1)%shf
2732                   CASE ( 'surf_v(1)%qsws' )         
2733                      IF ( ALLOCATED( surf_v(1)%qsws )  .AND.  kk == 1 )       &
2734                         READ ( 13 )  surf_v(1)%qsws
2735                   CASE ( 'surf_v(1)%ssws' )         
2736                      IF ( ALLOCATED( surf_v(1)%ssws )  .AND.  kk == 1 )       &
2737                         READ ( 13 )  surf_v(1)%ssws
2738                   CASE ( 'surf_v(1)%qcsws' )         
2739                      IF ( ALLOCATED( surf_v(1)%qcsws )  .AND.  kk == 1 )      &
2740                         READ ( 13 )  surf_v(1)%qcsws
2741                   CASE ( 'surf_v(1)%ncsws' )         
2742                      IF ( ALLOCATED( surf_v(1)%ncsws )  .AND.  kk == 1 )      &
2743                         READ ( 13 )  surf_v(1)%ncsws
2744                   CASE ( 'surf_v(1)%qrsws' )         
2745                      IF ( ALLOCATED( surf_v(1)%qrsws )  .AND.  kk == 1 )      &
2746                         READ ( 13 )  surf_v(1)%qrsws
2747                   CASE ( 'surf_v(1)%nrsws' )         
2748                      IF ( ALLOCATED( surf_v(1)%nrsws )  .AND.  kk == 1 )      &
2749                         READ ( 13 )  surf_v(1)%nrsws
2750                   CASE ( 'surf_v(1)%sasws' )         
2751                      IF ( ALLOCATED( surf_v(1)%sasws )  .AND.  kk == 1 )      &
2752                         READ ( 13 )  surf_v(1)%sasws
2753                   CASE ( 'surf_v(1)%mom_uv' )         
2754                      IF ( ALLOCATED( surf_v(1)%mom_flux_uv )  .AND.  kk == 1 )&
2755                         READ ( 13 )  surf_v(1)%mom_flux_uv
2756                   CASE ( 'surf_v(1)%mom_w' )         
2757                      IF ( ALLOCATED( surf_v(1)%mom_flux_w )  .AND.  kk == 1 ) &
2758                         READ ( 13 )  surf_v(1)%mom_flux_w
2759                   CASE ( 'surf_v(1)%mom_tke' )         
2760                      IF ( ALLOCATED( surf_v(1)%mom_flux_tke )  .AND.  kk == 1 )&
2761                         READ ( 13 )  surf_v(1)%mom_flux_tke
2762
2763                   CASE ( 'surf_v(2)%start_index' )   
2764                      IF ( kk == 1 )                                           &
2765                         READ ( 13 )  surf_v(2)%start_index
2766                      l = 2
2767                   CASE ( 'surf_v(2)%end_index' )   
2768                      IF ( kk == 1 )                                           &
2769                         READ ( 13 )  surf_v(2)%end_index
2770                   CASE ( 'surf_v(2)%us' )         
2771                      IF ( ALLOCATED( surf_v(2)%us )  .AND.  kk == 1 )         &
2772                         READ ( 13 )  surf_v(2)%us
2773                   CASE ( 'surf_v(2)%ts' )         
2774                      IF ( ALLOCATED( surf_v(2)%ts )  .AND.  kk == 1 )         &
2775                         READ ( 13 )  surf_v(2)%ts
2776                   CASE ( 'surf_v(2)%qs' )         
2777                      IF ( ALLOCATED( surf_v(2)%qs )  .AND.  kk == 1 )         &
2778                         READ ( 13 )  surf_v(2)%qs
2779                   CASE ( 'surf_v(2)%ss' )         
2780                      IF ( ALLOCATED( surf_v(2)%ss )  .AND.  kk == 1 )         &
2781                         READ ( 13 )  surf_v(2)%ss
2782                   CASE ( 'surf_v(2)%qcs' )         
2783                      IF ( ALLOCATED( surf_v(2)%qcs )  .AND.  kk == 1 )        &
2784                         READ ( 13 )  surf_v(2)%qcs
2785                   CASE ( 'surf_v(2)%ncs' )         
2786                      IF ( ALLOCATED( surf_v(2)%ncs )  .AND.  kk == 1 )        &
2787                         READ ( 13 )  surf_v(2)%ncs
2788                   CASE ( 'surf_v(2)%qrs' )         
2789                      IF ( ALLOCATED( surf_v(2)%qrs )  .AND.  kk == 1 )        &
2790                         READ ( 13 )  surf_v(2)%qrs
2791                   CASE ( 'surf_v(2)%nrs' )         
2792                      IF ( ALLOCATED( surf_v(2)%nrs )  .AND.  kk == 1 )        &
2793                         READ ( 13 )  surf_v(2)%nrs
2794                   CASE ( 'surf_v(2)%ol' )         
2795                      IF ( ALLOCATED( surf_v(2)%ol )  .AND.  kk == 1 )         &
2796                         READ ( 13 )  surf_v(2)%ol
2797                   CASE ( 'surf_v(2)%rib' )         
2798                      IF ( ALLOCATED( surf_v(2)%rib )  .AND.  kk == 1 )        &
2799                         READ ( 13 )  surf_v(2)%rib
2800                   CASE ( 'surf_v(2)%shf' )         
2801                      IF ( ALLOCATED( surf_v(2)%shf )  .AND.  kk == 1 )        &
2802                         READ ( 13 )  surf_v(2)%shf
2803                   CASE ( 'surf_v(2)%qsws' )         
2804                      IF ( ALLOCATED( surf_v(2)%qsws )  .AND.  kk == 1 )       &
2805                         READ ( 13 )  surf_v(2)%qsws
2806                   CASE ( 'surf_v(2)%ssws' )         
2807                      IF ( ALLOCATED( surf_v(2)%ssws )  .AND.  kk == 1 )       &
2808                         READ ( 13 )  surf_v(2)%ssws
2809                   CASE ( 'surf_v(2)%qcsws' )         
2810                      IF ( ALLOCATED( surf_v(2)%qcsws )  .AND.  kk == 1 )      &
2811                         READ ( 13 )  surf_v(2)%qcsws
2812                   CASE ( 'surf_v(2)%ncsws' )         
2813                      IF ( ALLOCATED( surf_v(2)%ncsws )  .AND.  kk == 1 )      &
2814                         READ ( 13 )  surf_v(2)%ncsws
2815                   CASE ( 'surf_v(2)%qrsws' )         
2816                      IF ( ALLOCATED( surf_v(2)%qrsws )  .AND.  kk == 1 )      &
2817                         READ ( 13 )  surf_v(2)%qrsws
2818                   CASE ( 'surf_v(2)%nrsws' )         
2819                      IF ( ALLOCATED( surf_v(2)%nrsws )  .AND.  kk == 1 )      &
2820                         READ ( 13 )  surf_v(2)%nrsws
2821                   CASE ( 'surf_v(2)%sasws' )         
2822                      IF ( ALLOCATED( surf_v(2)%sasws )  .AND.  kk == 1 )      &
2823                         READ ( 13 )  surf_v(2)%sasws
2824                   CASE ( 'surf_v(2)%mom_uv' )         
2825                      IF ( ALLOCATED( surf_v(2)%mom_flux_uv )  .AND.  kk == 1 )&
2826                         READ ( 13 )  surf_v(2)%mom_flux_uv
2827                   CASE ( 'surf_v(2)%mom_w' )         
2828                      IF ( ALLOCATED( surf_v(2)%mom_flux_w )  .AND.  kk == 1 ) &
2829                         READ ( 13 )  surf_v(2)%mom_flux_w
2830                   CASE ( 'surf_v(2)%mom_tke' )         
2831                      IF ( ALLOCATED( surf_v(2)%mom_flux_tke )  .AND.  kk == 1 )&
2832                         READ ( 13 )  surf_v(2)%mom_flux_tke
2833
2834                   CASE ( 'surf_v(3)%start_index' )   
2835                      IF ( kk == 1 )                                           &
2836                         READ ( 13 )  surf_v(3)%start_index
2837                      l = 3
2838                   CASE ( 'surf_v(3)%end_index' )   
2839                      IF ( kk == 1 )                                           &
2840                         READ ( 13 )  surf_v(3)%end_index
2841                   CASE ( 'surf_v(3)%us' )         
2842                      IF ( ALLOCATED( surf_v(3)%us )  .AND.  kk == 1 )         &
2843                         READ ( 13 )  surf_v(3)%us
2844                   CASE ( 'surf_v(3)%ts' )         
2845                      IF ( ALLOCATED( surf_v(3)%ts )  .AND.  kk == 1 )         &
2846                         READ ( 13 )  surf_v(3)%ts
2847                   CASE ( 'surf_v(3)%qs' )       
2848                      IF ( ALLOCATED( surf_v(3)%qs )  .AND.  kk == 1 )         &
2849                         READ ( 13 )  surf_v(3)%qs
2850                   CASE ( 'surf_v(3)%ss' )         
2851                      IF ( ALLOCATED( surf_v(3)%ss )  .AND.  kk == 1 )         &
2852                         READ ( 13 )  surf_v(3)%ss
2853                   CASE ( 'surf_v(3)%qcs' )         
2854                      IF ( ALLOCATED( surf_v(3)%qcs )  .AND.  kk == 1 )        &
2855                         READ ( 13 )  surf_v(3)%qcs
2856                   CASE ( 'surf_v(3)%ncs' )         
2857                      IF ( ALLOCATED( surf_v(3)%ncs )  .AND.  kk == 1 )        &
2858                         READ ( 13 )  surf_v(3)%ncs
2859                   CASE ( 'surf_v(3)%qrs' )         
2860                      IF ( ALLOCATED( surf_v(3)%qrs )  .AND.  kk == 1 )        &
2861                         READ ( 13 )  surf_v(3)%qrs
2862                   CASE ( 'surf_v(3)%nrs' )         
2863                      IF ( ALLOCATED( surf_v(3)%nrs )  .AND.  kk == 1 )        &
2864                         READ ( 13 )  surf_v(3)%nrs
2865                   CASE ( 'surf_v(3)%ol' )         
2866                      IF ( ALLOCATED( surf_v(3)%ol )  .AND.  kk == 1 )         &
2867                         READ ( 13 )  surf_v(3)%ol
2868                   CASE ( 'surf_v(3)%rib' )         
2869                      IF ( ALLOCATED( surf_v(3)%rib )  .AND.  kk == 1 )        &
2870                         READ ( 13 )  surf_v(3)%rib
2871                   CASE ( 'surf_v(3)%shf' )         
2872                      IF ( ALLOCATED( surf_v(3)%shf )  .AND.  kk == 1 )        &
2873                         READ ( 13 )  surf_v(3)%shf
2874                   CASE ( 'surf_v(3)%qsws' )         
2875                      IF ( ALLOCATED( surf_v(3)%qsws )  .AND.  kk == 1 )       &
2876                         READ ( 13 )  surf_v(3)%qsws
2877                   CASE ( 'surf_v(3)%ssws' )         
2878                      IF ( ALLOCATED( surf_v(3)%ssws )  .AND.  kk == 1 )       &
2879                         READ ( 13 )  surf_v(3)%ssws
2880                   CASE ( 'surf_v(3)%qcsws' )         
2881                      IF ( ALLOCATED( surf_v(3)%qcsws )  .AND.  kk == 1 )      &
2882                         READ ( 13 )  surf_v(3)%qcsws
2883                   CASE ( 'surf_v(3)%ncsws' )         
2884                      IF ( ALLOCATED( surf_v(3)%ncsws )  .AND.  kk == 1 )      &
2885                         READ ( 13 )  surf_v(3)%ncsws
2886                   CASE ( 'surf_v(3)%qrsws' )         
2887                      IF ( ALLOCATED( surf_v(3)%qrsws )  .AND.  kk == 1 )      &
2888                         READ ( 13 )  surf_v(3)%qrsws
2889                   CASE ( 'surf_v(3)%nrsws' )         
2890                      IF ( ALLOCATED( surf_v(3)%nrsws )  .AND.  kk == 1 )      &
2891                         READ ( 13 )  surf_v(3)%nrsws
2892                   CASE ( 'surf_v(3)%sasws' )         
2893                      IF ( ALLOCATED( surf_v(3)%sasws )  .AND.  kk == 1 )      &
2894                         READ ( 13 )  surf_v(3)%sasws
2895                   CASE ( 'surf_v(3)%mom_uv' )         
2896                      IF ( ALLOCATED( surf_v(3)%mom_flux_uv )  .AND.  kk == 1 )&
2897                         READ ( 13 )  surf_v(3)%mom_flux_uv
2898                   CASE ( 'surf_v(3)%mom_w' )         
2899                      IF ( ALLOCATED( surf_v(3)%mom_flux_w )  .AND.  kk == 1 ) &
2900                         READ ( 13 )  surf_v(3)%mom_flux_w
2901                   CASE ( 'surf_v(3)%mom_tke' )         
2902                      IF ( ALLOCATED( surf_v(3)%mom_flux_tke )  .AND.  kk == 1 )&
2903                         READ ( 13 )  surf_v(3)%mom_flux_tke
2904
2905                END SELECT
2906!
2907!--             Redistribute surface elements on its respective type.
2908                IF ( horizontal_surface )  THEN
2909                   ic = nxlc
2910                   DO  i = nxlf, nxrf
2911                      jc = nysc
2912                      DO  j = nysf, nynf
2913
2914                         surf_match_def  = surf_def_h(l)%end_index(jc,ic) >=   &
2915                                           surf_def_h(l)%start_index(jc,ic)
2916                         surf_match_lsm  = surf_lsm_h%end_index(jc,ic)    >=   &
2917                                           surf_lsm_h%start_index(jc,ic)
2918                         surf_match_usm  = surf_usm_h%end_index(jc,ic)    >=   &
2919                                           surf_usm_h%start_index(jc,ic)
2920
2921                         IF ( surf_match_def )  THEN
2922                            mm = surf_def_h(l)%start_index(jc,ic)
2923                            DO  m = surf_h(l)%start_index(j,i),                &
2924                                    surf_h(l)%end_index(j,i)
2925                               CALL restore_surface_elements( surf_def_h(l),   &
2926                                                              mm, surf_h(l), m )
2927                               mm = mm + 1
2928                            ENDDO
2929                         ENDIF
2930
2931                         IF ( surf_match_lsm )  THEN
2932                            mm = surf_lsm_h%start_index(jc,ic)
2933                            DO  m = surf_h(l)%start_index(j,i),                &
2934                                    surf_h(l)%end_index(j,i)
2935                               CALL restore_surface_elements( surf_lsm_h,      &
2936                                                              mm, surf_h(l), m )
2937                               mm = mm + 1
2938                            ENDDO
2939                         ENDIF
2940
2941                         IF ( surf_match_usm )  THEN
2942                            mm = surf_usm_h%start_index(jc,ic)
2943                            DO  m = surf_h(l)%start_index(j,i),                &
2944                                    surf_h(l)%end_index(j,i)
2945                               CALL restore_surface_elements( surf_usm_h,      &
2946                                                              mm, surf_h(l), m )
2947                               mm = mm + 1
2948                            ENDDO
2949                         ENDIF
2950
2951                         jc = jc + 1
2952                      ENDDO
2953                      ic = ic + 1
2954                   ENDDO
2955                ELSEIF ( vertical_surface )  THEN
2956                   ic = nxlc
2957                   DO  i = nxlf, nxrf
2958                      jc = nysc
2959                      DO  j = nysf, nynf
2960
2961                         surf_match_def  = surf_def_v(l)%end_index(jc,ic) >=   &
2962                                           surf_def_v(l)%start_index(jc,ic)
2963                         surf_match_lsm  = surf_lsm_v(l)%end_index(jc,ic) >=   &
2964                                           surf_lsm_v(l)%start_index(jc,ic)
2965                         surf_match_usm  = surf_usm_v(l)%end_index(jc,ic) >=   &
2966                                           surf_usm_v(l)%start_index(jc,ic)
2967
2968
2969
2970                         IF ( surf_match_def )  THEN
2971                            mm = surf_def_v(l)%start_index(jc,ic)
2972                            DO  m = surf_v(l)%start_index(j,i),                &
2973                                    surf_v(l)%end_index(j,i)
2974                               CALL restore_surface_elements( surf_def_v(l),   &
2975                                                              mm, surf_v(l), m )
2976                               mm = mm + 1
2977                            ENDDO
2978                         ENDIF
2979
2980                         IF ( surf_match_lsm )  THEN
2981                            mm = surf_lsm_v(l)%start_index(jc,ic)
2982                            DO  m = surf_v(l)%start_index(j,i),                &
2983                                    surf_v(l)%end_index(j,i)
2984                               CALL restore_surface_elements( surf_lsm_v(l),   &
2985                                                              mm, surf_v(l), m )
2986                               mm = mm + 1
2987                            ENDDO
2988                         ENDIF
2989   
2990                         IF ( surf_match_usm )  THEN
2991                            mm = surf_usm_v(l)%start_index(jc,ic)
2992                            DO  m = surf_v(l)%start_index(j,i),                &
2993                                    surf_v(l)%end_index(j,i)
2994                               CALL restore_surface_elements( surf_usm_v(l),   &
2995                                                              mm, surf_v(l), m )
2996                               mm = mm + 1
2997                            ENDDO
2998                         ENDIF
2999
3000                         jc = jc + 1
3001                      ENDDO
3002                      ic = ic + 1
3003                   ENDDO
3004                ENDIF
3005
3006             ENDDO
3007
3008             READ ( 13 )  field_chr
3009
3010          ENDDO
3011
3012       ENDIF
3013
3014
3015       CONTAINS
3016!------------------------------------------------------------------------------!
3017! Description:
3018! ------------
3019!> Restores surfacle elements back on its respective type.
3020!------------------------------------------------------------------------------!
3021          SUBROUTINE restore_surface_elements( surf_target, m_target,          &
3022                                               surf_file,   m_file )
3023
3024             IMPLICIT NONE
3025
3026             INTEGER(iwp)      ::  m_file      !< respective surface-element index of current surface array
3027             INTEGER(iwp)      ::  m_target    !< respecitve surface-element index of surface array on file
3028
3029             TYPE( surf_type ) ::  surf_target !< target surface type
3030             TYPE( surf_type ) ::  surf_file   !< surface type on file
3031
3032             IF ( SCAN( TRIM( field_chr ), '%us' ) /= 0 )  THEN
3033                IF ( ALLOCATED( surf_target%us )  .AND.                        &
3034                     ALLOCATED( surf_file%us   ) )                             & 
3035                   surf_target%us(m_target) = surf_file%us(m_file)
3036             ENDIF
3037
3038             IF ( SCAN( TRIM( field_chr ), '%ol' ) /= 0 )  THEN
3039                IF ( ALLOCATED( surf_target%ol )  .AND.                        &
3040                     ALLOCATED( surf_file%ol   ) )                             & 
3041                   surf_target%ol(m_target) = surf_file%ol(m_file)
3042             ENDIF
3043
3044             IF ( SCAN( TRIM( field_chr ), '%usws' ) /= 0 )  THEN
3045                IF ( ALLOCATED( surf_target%usws )  .AND.                      &
3046                     ALLOCATED( surf_file%usws   ) )                           & 
3047                   surf_target%usws(m_target) = surf_file%usws(m_file)
3048             ENDIF
3049
3050             IF ( SCAN( TRIM( field_chr ), '%vsws' ) /= 0 )  THEN
3051                IF ( ALLOCATED( surf_target%vsws )  .AND.                      &
3052                     ALLOCATED( surf_file%vsws   ) )                           & 
3053                   surf_target%vsws(m_target) = surf_file%vsws(m_file)
3054             ENDIF
3055
3056             IF ( SCAN( TRIM( field_chr ), '%ts' ) /= 0 )  THEN
3057                IF ( ALLOCATED( surf_target%ts )  .AND.                        &
3058                     ALLOCATED( surf_file%ts   ) )                             & 
3059                   surf_target%ts(m_target) = surf_file%ts(m_file)
3060             ENDIF
3061
3062             IF ( SCAN( TRIM( field_chr ), '%shf' ) /= 0 )  THEN
3063                IF ( ALLOCATED( surf_target%shf )  .AND.                       &
3064                     ALLOCATED( surf_file%shf   ) )                            & 
3065                   surf_target%shf(m_target) = surf_file%shf(m_file)
3066             ENDIF
3067
3068             IF ( SCAN( TRIM( field_chr ), '%qs' ) /= 0 )  THEN
3069                IF ( ALLOCATED( surf_target%qs )  .AND.                        &
3070                     ALLOCATED( surf_file%qs   ) )                             & 
3071                   surf_target%qs(m_target) = surf_file%qs(m_file)
3072             ENDIF
3073
3074             IF ( SCAN( TRIM( field_chr ), '%qsws' ) /= 0 )  THEN
3075                IF ( ALLOCATED( surf_target%qsws )  .AND.                      &
3076                     ALLOCATED( surf_file%qsws   ) )                           & 
3077                   surf_target%qsws(m_target) = surf_file%qsws(m_file)
3078             ENDIF
3079
3080             IF ( SCAN( TRIM( field_chr ), '%ss' ) /= 0 )  THEN
3081                IF ( ALLOCATED( surf_target%ss )  .AND.                        &
3082                     ALLOCATED( surf_file%ss   ) )                             & 
3083                   surf_target%ss(m_target) = surf_file%ss(m_file)
3084             ENDIF
3085
3086             IF ( SCAN( TRIM( field_chr ), '%ssws' ) /= 0 )  THEN
3087                IF ( ALLOCATED( surf_target%ssws )  .AND.                      &
3088                     ALLOCATED( surf_file%ssws   ) )                           & 
3089                   surf_target%ssws(m_target) = surf_file%ssws(m_file)
3090             ENDIF
3091
3092             IF ( SCAN( TRIM( field_chr ), '%qcs' ) /= 0 )  THEN
3093                IF ( ALLOCATED( surf_target%qcs )  .AND.                       &
3094                     ALLOCATED( surf_file%qcs   ) )                            & 
3095                  surf_target%qcs(m_target) = surf_file%qcs(m_file)
3096             ENDIF
3097
3098             IF ( SCAN( TRIM( field_chr ), '%qcsws' ) /= 0 )  THEN
3099                IF ( ALLOCATED( surf_target%qcsws )  .AND.                     &
3100                     ALLOCATED( surf_file%qcsws   ) )                          & 
3101                   surf_target%qcsws(m_target) = surf_file%qcsws(m_file)
3102             ENDIF
3103
3104             IF ( SCAN( TRIM( field_chr ), '%ncs' ) /= 0 )  THEN
3105                IF ( ALLOCATED( surf_target%ncs )  .AND.                       &
3106                     ALLOCATED( surf_file%ncs   ) )                            & 
3107                   surf_target%ncs(m_target) = surf_file%ncs(m_file)
3108             ENDIF
3109
3110             IF ( SCAN( TRIM( field_chr ), '%ncsws' ) /= 0 )  THEN
3111                IF ( ALLOCATED( surf_target%ncsws )  .AND.                     &
3112                     ALLOCATED( surf_file%ncsws   ) )                          & 
3113                   surf_target%ncsws(m_target) = surf_file%ncsws(m_file)
3114             ENDIF
3115
3116             IF ( SCAN( TRIM( field_chr ), '%qrs' ) /= 0 )  THEN
3117                IF ( ALLOCATED( surf_target%qrs )  .AND.                       &
3118                     ALLOCATED( surf_file%qrs   ) )                            & 
3119                  surf_target%qrs(m_target) = surf_file%qrs(m_file)
3120             ENDIF
3121
3122             IF ( SCAN( TRIM( field_chr ), '%qrsws' ) /= 0 )  THEN
3123                IF ( ALLOCATED( surf_target%qrsws )  .AND.                     &
3124                     ALLOCATED( surf_file%qrsws   ) )                          & 
3125                   surf_target%qrsws(m_target) = surf_file%qrsws(m_file)
3126             ENDIF
3127
3128             IF ( SCAN( TRIM( field_chr ), '%nrs' ) /= 0 )  THEN
3129                IF ( ALLOCATED( surf_target%nrs )  .AND.                       &
3130                     ALLOCATED( surf_file%nrs   ) )                            & 
3131                   surf_target%nrs(m_target) = surf_file%nrs(m_file)
3132             ENDIF
3133
3134             IF ( SCAN( TRIM( field_chr ), '%nrsws' ) /= 0 )  THEN
3135                IF ( ALLOCATED( surf_target%nrsws )  .AND.                     &
3136                     ALLOCATED( surf_file%nrsws   ) )                          & 
3137                   surf_target%nrsws(m_target) = surf_file%nrsws(m_file)
3138             ENDIF
3139
3140             IF ( SCAN( TRIM( field_chr ), '%sasws' ) /= 0 )  THEN
3141                IF ( ALLOCATED( surf_target%sasws )  .AND.                     &
3142                     ALLOCATED( surf_file%sasws   ) )                          & 
3143                   surf_target%sasws(m_target) = surf_file%sasws(m_file)
3144             ENDIF
3145
3146             IF ( SCAN( TRIM( field_chr ), '%mom_uv' ) /= 0 )  THEN
3147                IF ( ALLOCATED( surf_target%mom_flux_uv )  .AND.               &
3148                     ALLOCATED( surf_file%mom_flux_uv   ) )                    & 
3149                   surf_target%mom_flux_uv(m_target) =                         &
3150                                           surf_file%mom_flux_uv(m_file)
3151             ENDIF
3152
3153             IF ( SCAN( TRIM( field_chr ), '%mom_w' ) /= 0 )  THEN
3154                IF ( ALLOCATED( surf_target%mom_flux_w )  .AND.                &
3155                     ALLOCATED( surf_file%mom_flux_w   ) )                     & 
3156                   surf_target%mom_flux_w(m_target) =                          &
3157                                           surf_file%mom_flux_w(m_file)
3158             ENDIF
3159
3160             IF ( SCAN( TRIM( field_chr ), '%mom_tke' ) /= 0 )  THEN
3161                IF ( ALLOCATED( surf_target%mom_flux_tke )  .AND.              &
3162                     ALLOCATED( surf_file%mom_flux_tke   ) )                   & 
3163                   surf_target%mom_flux_tke(0:1,m_target) =                    &
3164                                           surf_file%mom_flux_tke(0:1,m_file)
3165             ENDIF
3166
3167          END SUBROUTINE restore_surface_elements
3168
3169    END SUBROUTINE surface_read_restart_data
3170
3171 
3172!------------------------------------------------------------------------------!
3173! Description:
3174! ------------
3175!> Counts the number of surface elements with the same facing, required for
3176!> reading and writing restart data.
3177!------------------------------------------------------------------------------!
3178    SUBROUTINE surface_last_actions
3179
3180       IMPLICIT NONE
3181!
3182!--    Horizontal surfaces
3183       ns_h_on_file(0) = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns
3184       ns_h_on_file(1) = surf_def_h(1)%ns
3185       ns_h_on_file(2) = surf_def_h(2)%ns
3186!
3187!--    Vertical surfaces
3188       ns_v_on_file(0) = surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns
3189       ns_v_on_file(1) = surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns
3190       ns_v_on_file(2) = surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns
3191       ns_v_on_file(3) = surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns
3192
3193    END SUBROUTINE surface_last_actions
3194
3195
3196 END MODULE surface_mod
Note: See TracBrowser for help on using the repository browser.