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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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