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

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

Bugfix, initialization of surface elements also in case of restart runs and cyclic fill

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