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

Last change on this file since 3147 was 3147, checked in by maronga, 3 years ago

last commit documented

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