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

Last change on this file since 3761 was 3761, checked in by raasch, 2 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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