source: palm/tags/release-5.0/SOURCE/surface_mod.f90 @ 2704

Last change on this file since 2704 was 2698, checked in by suehring, 6 years ago

Particle reflections at downward-facing walls; revision of particle speed interpolations at walls; bugfixes in get_topography_index and in date constants

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