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

Last change on this file since 2813 was 2813, checked in by suehring, 7 years ago

Further bugfixes concering restart runs

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