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

Last change on this file since 3146 was 3146, checked in by maronga, 3 years ago

major bugfix in calculation of Obukhov length and subsequent handling of surface information

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