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

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

last commit documented

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