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

Last change on this file since 2977 was 2977, checked in by kanani, 3 years ago

Fixes for radiative transfer model

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