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

Last change on this file since 4281 was 4245, checked in by pavelkrc, 20 months ago

Merge branch resler into trunk

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