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

Last change on this file since 3298 was 3298, checked in by kanani, 3 years ago

Merge chemistry branch at r3297 to trunk

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