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

Last change on this file since 2759 was 2759, checked in by suehring, 3 years ago

Major bugfix in horizontal diffusion of all scalar quantities at vertical surfaces; further bugfix: density is considered for vertical fluxes of passive scalar and salinity

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