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

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

Bugfix, missing deallocation of q_surface in case of restarts

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