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

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

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

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