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

Last change on this file since 3767 was 3767, checked in by raasch, 6 years ago

unused variables removed from rrd-subroutines parameter list

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