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

Last change on this file since 2942 was 2942, checked in by suehring, 3 years ago

Bugfix in assigning surface element data after restart.

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