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

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

Removal of chem directive, plus minor changes

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