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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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