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

Last change on this file since 3444 was 3444, checked in by schwenkel, 3 years ago

Bugfix in surface_wrd_local for cloud and rain water flux

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