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

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

Branch salsa @3446 re-integrated into trunk

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