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

Last change on this file since 2563 was 2547, checked in by schwenkel, 7 years ago

extended by cloud_droplets option

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