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

Last change on this file since 2638 was 2638, checked in by raasch, 4 years ago

bugfix for cases with constant top momentumflux

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