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

Last change on this file since 4150 was 4150, checked in by suehring, 23 months ago

Framework to input single surface variables independent on land/urban-surface model via static input file provided. This way, input and initialization of heterogeneous roughness was already added.

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