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

Last change on this file since 2812 was 2812, checked in by hellstea, 3 years ago

Bugfixes in computation of the interpolation loglaw-correction parameters

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