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

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

get topograpyh top index via function call

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