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

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

Assure that emissions are only at natural (pavement) surfaces in case of parametrized emission mode, not on urban surfaces

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