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

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

Do not overwrite values of albedo in radiation_init in case albedo has been already initialized in the urban-surface model via ASCII input

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