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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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