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

Last change on this file since 2805 was 2805, checked in by suehring, 3 years ago

Bugfix in re-mapping surface-element data in case of restarts; bugfix in initialization of water-surface roughness; bugfix - uninitialized resistance at urban-type surfaces

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