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

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

revised calculation of near surface air potential temperature

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