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

Last change on this file since 4366 was 4366, checked in by raasch, 4 years ago

code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length

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