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

Last change on this file since 3933 was 3933, checked in by kanani, 5 years ago

Bugfixes and clean-up for output quantity theta_2m*

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