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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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