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

Last change on this file since 2693 was 2688, checked in by Giersch, 6 years ago

Error message PA0476 added. Bugfix in case of coupled runs.

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