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

Last change on this file since 3049 was 3026, checked in by schwenkel, 3 years ago

Changed the name specific humidity to mixing ratio

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