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

Last change on this file since 3943 was 3943, checked in by maronga, 2 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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