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

Last change on this file since 4156 was 4156, checked in by schwenkel, 23 months ago

Bugfix in case of cloud microphysics morrison

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