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

Last change on this file since 3421 was 3418, checked in by kanani, 3 years ago

Add green facades, update building data base, fix for thin walls in spinup

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