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

Last change on this file since 2735 was 2735, checked in by suehring, 3 years ago

Output of resistance also urban-type surfaces

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