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

Last change on this file since 2421 was 2378, checked in by suehring, 7 years ago

Bugfix in write restart data

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