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

Last change on this file since 2894 was 2894, checked in by Giersch, 3 years ago

Reading/Writing? data in case of restart runs revised

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