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

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

Enable restarts with USM with different number of PEs; some bugfixes in new surface structure in USM; formatting adjustments and descriptions in surface_mod

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