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

Last change on this file since 2963 was 2963, checked in by suehring, 7 years ago

Minor revision of static input file checks, bugfix in initialization of surface-fractions in LSM; minor bugfix in initialization of albedo at window-surfaces; for clearer access of albedo and emissivity introduce index for vegetation/wall, pavement/green-wall and water/window surfaces

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