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

Last change on this file since 3833 was 3833, checked in by forkel, 2 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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