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

Last change on this file since 4104 was 4104, checked in by suehring, 2 years ago

initialization of index space for boundary data structure accidantly run over ghost points, causing a segmentation fault

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