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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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