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

Last change on this file since 2706 was 2706, checked in by suehring, 6 years ago

Bugfixes in radiation and restarts in LSM

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