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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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