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

Last change on this file since 2270 was 2270, checked in by maronga, 7 years ago

major revisions in land surface model

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