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

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

change default value of zeta_min; enable scalar/heat/water-fluxes at downward-facing surfaces, bugfix Makefile

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