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

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

variable description added some routines

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