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

Last change on this file since 2339 was 2339, checked in by gronemeier, 7 years ago

corrected timestamp in header

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