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

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

Introduction of addtional surface variables indicating type and name of the surface elements; Bugfix in LSM

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