source: palm/trunk/SOURCE/init_3d_model.f90 @ 4652

Last change on this file since 4652 was 4648, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to (toggle deleted branches)
    /palm/branches/chemistry/SOURCE/init_3d_model.f902047-3190,​3218-3297
    /palm/branches/rans/SOURCE/init_3d_model.f902078-3128
    /palm/branches/resler/SOURCE/init_3d_model.f902023-4492
    /palm/branches/salsa/SOURCE/init_3d_model.f902503-3581
    /palm/branches/forwind/SOURCE/init_3d_model.f901564-1913
    /palm/branches/mosaik_M2/init_3d_model.f902360-3471
    /palm/branches/palm4u/SOURCE/init_3d_model.f902540-2692
File size: 84.7 KB
Line 
1!> @file init_3d_model.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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: init_3d_model.f90 4648 2020-08-25 07:52:08Z raasch $
26! file re-formatted to follow the PALM coding standard
27!
28! 4548 2020-05-28 19:36:45Z suehring
29! Bugfix, move call for lsf_forcing_surf after lsf_init is called
30!
31! 4514 2020-04-30 16:29:59Z suehring
32! Add possibility to initialize surface sensible and latent heat fluxes via a static driver.
33!
34! 4493 2020-04-10 09:49:43Z pavelkrc
35! Overwrite u_init, v_init, pt_init, q_init and s_init with hom for all cyclic_fill-cases, not only
36! for turbulent_inflow = .TRUE.
37!
38! 4360 2020-01-07 11:25:50Z suehring
39! Introduction of wall_flags_total_0, which currently sets bits based on static topography
40! information used in wall_flags_static_0
41!
42! 4329 2019-12-10 15:46:36Z motisi
43! Renamed wall_flags_0 to wall_flags_static_0
44!
45! 4286 2019-10-30 16:01:14Z resler
46! implement new palm_date_time_mod
47!
48! 4223 2019-09-10 09:20:47Z gronemeier
49! Deallocate temporary string array since it may be re-used to read different input data in other
50! modules
51!
52! 4186 2019-08-23 16:06:14Z suehring
53! Design change, use variables defined in netcdf_data_input_mod to read netcd variables rather than
54! define local ones.
55!
56! 4185 2019-08-23 13:49:38Z oliver.maas
57! For initializing_actions = ' cyclic_fill':
58! Overwrite u_init, v_init, pt_init, q_init and s_init with the (temporally) and horizontally
59! averaged vertical profiles from the end of the prerun, because these profiles shall be used as the
60! basic state for the rayleigh damping and the pt_damping.
61!
62! 4182 2019-08-22 15:20:23Z scharf
63! Corrected "Former revisions" section
64!
65! 4168 2019-08-16 13:50:17Z suehring
66! Replace function get_topography_top_index by topo_top_ind
67!
68! 4151 2019-08-09 08:24:30Z suehring
69! Add netcdf directive around input calls (fix for last commit)
70!
71! 4150 2019-08-08 20:00:47Z suehring
72! Input of additional surface variables independent on land- or urban-surface model
73!
74! 4131 2019-08-02 11:06:18Z monakurppa
75! Allocate sums and sums_l to allow profile output for salsa variables.
76!
77! 4130 2019-08-01 13:04:13Z suehring
78! Effectively reduce 3D initialization to 1D initial profiles. This is because 3D initialization
79! produces structures in the w-component that are correlated with the processor grid for some
80! unknown reason
81!
82! 4090 2019-07-11 15:06:47Z Giersch
83! Unused variables removed
84!
85! 4088 2019-07-11 13:57:56Z Giersch
86! Pressure and density profile calculations revised using basic functions
87!
88! 4048 2019-06-21 21:00:21Z knoop
89! Further modularization of particle code components
90!
91! 4017 2019-06-06 12:16:46Z schwenkel
92! Convert most location messages to debug messages to reduce output in job logfile to a minimum
93!
94! unused variable removed
95!
96! 3937 2019-04-29 15:09:07Z suehring
97! Move initialization of synthetic turbulence generator behind initialization of offline nesting.
98! Remove call for stg_adjust, as this is now already done in stg_init.
99!
100! 3900 2019-04-16 15:17:43Z suehring
101! Fix problem with LOD = 2 initialization
102!
103! 3885 2019-04-11 11:29:34Z kanani
104! Changes related to global restructuring of location messages and introduction of additional debug
105! messages
106!
107! 3849 2019-04-01 16:35:16Z knoop
108! Move initialization of rmask before initializing user_init_arrays
109!
110! 3711 2019-01-31 13:44:26Z knoop
111! Introduced module_interface_init_checks for post-init checks in modules
112!
113! 3700 2019-01-26 17:03:42Z knoop
114! Some interface calls moved to module_interface + cleanup
115!
116! 3648 2019-01-02 16:35:46Z suehring
117! Rename subroutines for surface-data output
118!
119! Revision 1.1  1998/03/09 16:22:22  raasch
120! Initial revision
121!
122!
123! Description:
124! ------------
125!> Allocation of arrays and initialization of the 3D model via
126!> a) pre-run the 1D model
127!> or
128!> b) pre-set constant linear profiles
129!> or
130!> c) read values of a previous run
131!--------------------------------------------------------------------------------------------------!
132 SUBROUTINE init_3d_model
133
134
135    USE advec_ws
136
137    USE arrays_3d
138
139    USE basic_constants_and_equations_mod,                                                         &
140        ONLY:  barometric_formula, c_p, exner_function, exner_function_invers, g,                  &
141               ideal_gas_law_rho, ideal_gas_law_rho_pt, l_v, pi
142
143    USE bulk_cloud_model_mod,                                                                      &
144        ONLY:  bulk_cloud_model
145
146    USE chem_modules,                                                                              &
147        ONLY:  max_pr_cs ! ToDo: this dependency needs to be removed cause it is ugly #new_dom
148
149    USE control_parameters
150
151    USE grid_variables,                                                                            &
152        ONLY:  dx, dy, ddx2_mg, ddy2_mg
153
154    USE indices
155
156    USE kinds
157
158    USE lsf_nudging_mod,                                                                           &
159        ONLY:  ls_forcing_surf
160
161    USE model_1d_mod,                                                                              &
162        ONLY:  init_1d_model, l1d, u1d, v1d
163
164    USE module_interface,                                                                          &
165        ONLY:  module_interface_init_arrays,                                                       &
166               module_interface_init,                                                              &
167               module_interface_init_checks
168
169    USE multi_agent_system_mod,                                                                    &
170        ONLY:  agents_active, mas_init
171
172    USE netcdf_interface,                                                                          &
173        ONLY:  dots_max
174
175    USE netcdf_data_input_mod,                                                                     &
176        ONLY:  char_fill,                                                                          &
177               check_existence,                                                                    &
178               close_input_file,                                                                   &
179               get_attribute,                                                                      &
180               get_variable,                                                                       &
181               init_3d,                                                                            &
182               input_pids_static,                                                                  &
183               inquire_num_variables,                                                              &
184               inquire_variable_names,                                                             &
185               input_file_static,                                                                  &
186               netcdf_data_input_init_3d,                                                          &
187               num_var_pids,                                                                       &
188               open_read_file,                                                                     &
189               pids_id,                                                                            &
190               real_2d,                                                                            &
191               vars_pids
192
193    USE nesting_offl_mod,                                                                          &
194        ONLY:  nesting_offl_init
195
196    USE palm_date_time_mod,                                                                        &
197        ONLY:  set_reference_date_time
198
199    USE pegrid
200
201#if defined( __parallel )
202    USE pmc_interface,                                                                             &
203        ONLY:  nested_run
204#endif
205
206    USE random_function_mod
207
208    USE random_generator_parallel,                                                                 &
209        ONLY:  init_parallel_random_generator
210
211    USE read_restart_data_mod,                                                                     &
212        ONLY:  rrd_local, rrd_read_parts_of_global
213
214    USE statistics,                                                                                &
215        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask, statistic_regions, sums,   &
216               sums_divnew_l, sums_divold_l, sums_l, sums_l_l, sums_wsts_bc_l, ts_value,           &
217               weight_pres, weight_substep
218
219    USE synthetic_turbulence_generator_mod,                                                        &
220        ONLY:  stg_init, use_syn_turb_gen
221
222    USE surface_layer_fluxes_mod,                                                                  &
223        ONLY:  init_surface_layer_fluxes
224
225    USE surface_mod,                                                                               &
226        ONLY :  init_single_surface_properties,                                                    &
227                init_surface_arrays,                                                               &
228                init_surfaces,                                                                     &
229                surf_def_h,                                                                        &
230                surf_def_v,                                                                        &
231                surf_lsm_h,                                                                        &
232                surf_usm_h
233
234#if defined( _OPENACC )
235    USE surface_mod,                                                                               &
236        ONLY :  bc_h
237#endif
238
239    USE surface_data_output_mod,                                                                   &
240        ONLY:  surface_data_output_init
241
242    USE transpose_indices
243
244
245    IMPLICIT NONE
246
247    INTEGER(iwp) ::  i                    !< grid index in x direction
248    INTEGER(iwp) ::  ind_array(1)         !< dummy used to determine start index for external pressure forcing
249    INTEGER(iwp) ::  j                    !< grid index in y direction
250    INTEGER(iwp) ::  k                    !< grid index in z direction
251    INTEGER(iwp) ::  k_surf               !< surface level index
252    INTEGER(iwp) ::  l                    !< running index over surface orientation
253    INTEGER(iwp) ::  m                    !< index of surface element in surface data type
254    INTEGER(iwp) ::  nz_u_shift           !< topography-top index on u-grid, used to vertically shift initial profiles
255    INTEGER(iwp) ::  nz_v_shift           !< topography-top index on v-grid, used to vertically shift initial profiles
256    INTEGER(iwp) ::  nz_w_shift           !< topography-top index on w-grid, used to vertically shift initial profiles
257    INTEGER(iwp) ::  nz_s_shift           !< topography-top index on scalar-grid, used to vertically shift initial profiles
258    INTEGER(iwp) ::  nz_u_shift_l         !< topography-top index on u-grid, used to vertically shift initial profiles
259    INTEGER(iwp) ::  nz_v_shift_l         !< topography-top index on v-grid, used to vertically shift initial profiles
260    INTEGER(iwp) ::  nz_w_shift_l         !< topography-top index on w-grid, used to vertically shift initial profiles
261    INTEGER(iwp) ::  nz_s_shift_l         !< topography-top index on scalar-grid, used to vertically shift initial profiles
262    INTEGER(iwp) ::  nzt_l                !< index of top PE boundary for multigrid level
263    INTEGER(iwp) ::  sr                   !< index of statistic region
264
265    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !< toal number of horizontal grid points in statistical region on
266                                                             !< subdomain
267
268    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !< number of horizontal non-wall bounded grid points on
269                                                                     !< subdomain
270    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !< number of horizontal non-topography grid points on
271                                                                     !< subdomain
272
273    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
274    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
275
276    REAL(wp), DIMENSION(:), ALLOCATABLE ::  init_l                       !< dummy array used for averaging 3D data to obtain
277                                                                         !< inital profiles
278    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l  !< mean surface level height on subdomain
279    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l               !< total number of non-topography grid points on subdomain
280    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp             !< total number of non-topography grid points
281    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic                !< hydrostatic pressure
282
283    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !< area of lateral and top model domain surface on local subdomain
284    REAL(wp), DIMENSION(1:3) ::  volume_flow_initial_l  !< initial volume flow into model domain
285
286    TYPE(real_2d) ::  tmp_2d !< temporary variable to input additional surface-data from static file
287
288    CALL location_message( 'model initialization', 'start' )
289!
290!-- Set reference date-time
291    CALL set_reference_date_time( date_time_str=origin_date_time )
292
293    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'start' )
294!
295!-- Allocate arrays
296    ALLOCATE( mean_surface_level_height(0:statistic_regions),                                      &
297              mean_surface_level_height_l(0:statistic_regions),                                    &
298              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),                        &
299              ngp_3d(0:statistic_regions),                                                         &
300              ngp_3d_inner(0:statistic_regions),                                                   &
301              ngp_3d_inner_l(0:statistic_regions),                                                 &
302              ngp_3d_inner_tmp(0:statistic_regions),                                               &
303              sums_divnew_l(0:statistic_regions),                                                  &
304              sums_divold_l(0:statistic_regions) )
305    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
306    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                                        &
307              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                                      &
308              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                                      &
309              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                                    &
310              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                                      &
311              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa),                          &
312              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:threads_per_task-1),   &
313              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),                        &
314              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
315    ALLOCATE( ts_value(dots_max,0:statistic_regions) )
316    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
317
318    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                                        &
319              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                    &
320              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
321
322    ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                 &
323              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                 &
324              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
325              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
326              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
327              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
328              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
329              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
330              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
331              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
332              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
333    IF (  .NOT.  neutral )  THEN
334       ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
335    ENDIF
336!
337!-- Pre-set masks for regional statistics. Default is the total model domain.
338!-- Ghost points are excluded because counting values at the ghost boundaries would bias the
339!-- statistics.
340    rmask = 1.0_wp
341    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
342    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
343!
344!-- Following array is required for perturbation pressure within the iterative pressure solvers. For
345!-- the multistep schemes (Runge-Kutta), array p holds the weighted average of the substeps and
346!-- cannot be used in the Poisson solver.
347    IF ( psolver == 'sor' )  THEN
348       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
349    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
350!
351!--    For performance reasons, multigrid is using one ghost layer only
352       ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
353    ENDIF
354
355!
356!-- Array for storing constant coeffficients of the tridiagonal solver
357    IF ( psolver == 'poisfft' )  THEN
358       ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
359       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
360    ENDIF
361
362    IF ( humidity )  THEN
363!
364!--    3D-humidity
365       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
366                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
367                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
368                 vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
369    ENDIF
370
371    IF ( passive_scalar )  THEN
372
373!
374!--    3D scalar arrays
375       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
376                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
377                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
378
379    ENDIF
380
381!
382!-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to non-zero values later
383!-- in ocean_init.
384    ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1),                                      &
385              v_stokes_zu(nzb:nzt+1), v_stokes_zw(nzb:nzt+1) )
386    u_stokes_zu(:) = 0.0_wp
387    u_stokes_zw(:) = 0.0_wp
388    v_stokes_zu(:) = 0.0_wp
389    v_stokes_zw(:) = 0.0_wp
390
391!
392!-- Allocation of anelastic and Boussinesq approximation specific arrays
393    ALLOCATE( p_hydrostatic(nzb:nzt+1) )
394    ALLOCATE( rho_air(nzb:nzt+1) )
395    ALLOCATE( rho_air_zw(nzb:nzt+1) )
396    ALLOCATE( drho_air(nzb:nzt+1) )
397    ALLOCATE( drho_air_zw(nzb:nzt+1) )
398!
399!-- Density profile calculation for anelastic and Boussinesq approximation.
400!-- In case of a Boussinesq approximation, a constant density is calculated mainly for output
401!-- purposes. This density does not need to be considered in the model's system of equations.
402    IF ( TRIM( approximation ) == 'anelastic' )  THEN
403       DO  k = nzb, nzt+1
404          p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *                                &
405                                                exner_function(surface_pressure * 100.0_wp),       &
406                                                surface_pressure * 100.0_wp)
407
408          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
409       ENDDO
410
411       DO  k = nzb, nzt
412          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
413       ENDDO
414
415       rho_air_zw(nzt+1)  = rho_air_zw(nzt) + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
416
417    ELSE
418       DO  k = nzb, nzt+1
419          p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *                              &
420                                                exner_function(surface_pressure * 100.0_wp),       &
421                                                surface_pressure * 100.0_wp)
422
423          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(nzb))
424       ENDDO
425
426       DO  k = nzb, nzt
427          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
428       ENDDO
429
430       rho_air_zw(nzt+1)  = rho_air_zw(nzt) + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
431
432    ENDIF
433!
434!-- Compute the inverse density array in order to avoid expencive divisions
435    drho_air    = 1.0_wp / rho_air
436    drho_air_zw = 1.0_wp / rho_air_zw
437
438!
439!-- Allocation of flux conversion arrays
440    ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
441    ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
442    ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
443    ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
444    ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
445    ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
446
447!
448!-- Calculate flux conversion factors according to approximation and in-/output mode
449    DO  k = nzb, nzt+1
450
451        IF ( TRIM( flux_input_mode ) == 'kinematic' )  THEN
452            heatflux_input_conversion(k)      = rho_air_zw(k)
453            waterflux_input_conversion(k)     = rho_air_zw(k)
454            momentumflux_input_conversion(k)  = rho_air_zw(k)
455        ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
456            heatflux_input_conversion(k)      = 1.0_wp / c_p
457            waterflux_input_conversion(k)     = 1.0_wp / l_v
458            momentumflux_input_conversion(k)  = 1.0_wp
459        ENDIF
460
461        IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
462            heatflux_output_conversion(k)     = drho_air_zw(k)
463            waterflux_output_conversion(k)    = drho_air_zw(k)
464            momentumflux_output_conversion(k) = drho_air_zw(k)
465        ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
466            heatflux_output_conversion(k)     = c_p
467            waterflux_output_conversion(k)    = l_v
468            momentumflux_output_conversion(k) = 1.0_wp
469        ENDIF
470
471        IF ( .NOT. humidity ) THEN
472            waterflux_input_conversion(k)  = 1.0_wp
473            waterflux_output_conversion(k) = 1.0_wp
474        ENDIF
475
476    ENDDO
477
478!
479!-- In case of multigrid method, compute grid lengths and grid factors for the grid levels with
480!-- respective density on each grid.
481    IF ( psolver(1:9) == 'multigrid' )  THEN
482
483       ALLOCATE( ddx2_mg(maximum_grid_level) )
484       ALLOCATE( ddy2_mg(maximum_grid_level) )
485       ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
486       ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
487       ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
488       ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
489       ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
490       ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
491       ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
492
493       dzu_mg(:,maximum_grid_level) = dzu
494       rho_air_mg(:,maximum_grid_level) = rho_air
495!
496!--    Next line to ensure an equally spaced grid.
497       dzu_mg(1,maximum_grid_level) = dzu(2)
498       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) + (rho_air(nzb) - rho_air(nzb+1))
499
500       dzw_mg(:,maximum_grid_level) = dzw
501       rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
502       nzt_l = nzt
503       DO  l = maximum_grid_level-1, 1, -1
504           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
505           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
506           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1)    + ( rho_air_mg(nzb,l+1)    -              &
507                                                             rho_air_mg(nzb+1,l+1)    )
508           rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + ( rho_air_zw_mg(nzb,l+1) -              &
509                                                             rho_air_zw_mg(nzb+1,l+1) )
510           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
511           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
512           nzt_l = nzt_l / 2
513           DO  k = 2, nzt_l+1
514              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
515              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
516              rho_air_mg(k,l)    = rho_air_mg(2*k-1,l+1)
517              rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
518           ENDDO
519       ENDDO
520
521       nzt_l = nzt
522       dx_l  = dx
523       dy_l  = dy
524       DO  l = maximum_grid_level, 1, -1
525          ddx2_mg(l) = 1.0_wp / dx_l**2
526          ddy2_mg(l) = 1.0_wp / dy_l**2
527          DO  k = nzb+1, nzt_l
528             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
529             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
530             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) )                                     &
531                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
532          ENDDO
533          nzt_l = nzt_l / 2
534          dx_l  = dx_l * 2.0_wp
535          dy_l  = dy_l * 2.0_wp
536       ENDDO
537
538    ENDIF
539
540!
541!-- 1D-array for large scale subsidence velocity
542    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
543       ALLOCATE ( w_subs(nzb:nzt+1) )
544       w_subs = 0.0_wp
545    ENDIF
546
547!
548!-- Arrays to store velocity data from t-dt and the phase speeds which are needed for radiation
549!-- boundary conditions.
550    IF ( bc_radiation_l )  THEN
551       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                                                   &
552                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                                                   &
553                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
554    ENDIF
555    IF ( bc_radiation_r )  THEN
556       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                                               &
557                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                                               &
558                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
559    ENDIF
560    IF ( bc_radiation_l  .OR.  bc_radiation_r )  THEN
561       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), c_w(nzb:nzt+1,nysg:nyng) )
562    ENDIF
563    IF ( bc_radiation_s )  THEN
564       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                                                   &
565                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                                                   &
566                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
567    ENDIF
568    IF ( bc_radiation_n )  THEN
569       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                                               &
570                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                                               &
571                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
572    ENDIF
573    IF ( bc_radiation_s  .OR.  bc_radiation_n )  THEN
574       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), c_w(nzb:nzt+1,nxlg:nxrg) )
575    ENDIF
576    IF ( bc_radiation_l  .OR.  bc_radiation_r  .OR.  bc_radiation_s  .OR.  bc_radiation_n )  THEN
577       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )
578       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
579    ENDIF
580
581!
582!-- Initial assignment of the pointers
583    IF ( .NOT. neutral )  THEN
584       pt => pt_1;  pt_p => pt_2;  tpt_m => pt_3
585    ELSE
586       pt => pt_1;  pt_p => pt_1;  tpt_m => pt_3
587    ENDIF
588    u  => u_1;   u_p  => u_2;   tu_m  => u_3
589    v  => v_1;   v_p  => v_2;   tv_m  => v_3
590    w  => w_1;   w_p  => w_2;   tw_m  => w_3
591
592    IF ( humidity )  THEN
593       q => q_1;  q_p => q_2;  tq_m => q_3
594       vpt  => vpt_1
595    ENDIF
596
597    IF ( passive_scalar )  THEN
598       s => s_1;  s_p => s_2;  ts_m => s_3
599    ENDIF
600
601!
602!-- Initialize surface arrays
603    CALL init_surface_arrays
604!
605!-- Allocate arrays for other modules
606    CALL module_interface_init_arrays
607
608
609!
610!-- Allocate arrays containing the RK coefficient for calculation of perturbation pressure and
611!-- turbulent fluxes. At this point values are set for pressure calculation during initialization
612!-- (where no timestep is done). Further below the values needed within the timestep scheme will be
613!-- set.
614    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),                                   &
615              weight_pres(1:intermediate_timestep_count_max) )
616    weight_substep = 1.0_wp
617    weight_pres    = 1.0_wp
618    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
619
620    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'end' )
621
622!
623!-- Initialize time series
624    ts_value = 0.0_wp
625
626!
627!-- Initialize local summation arrays for routine flow_statistics.
628!-- This is necessary because they may not yet have been initialized when they are called from
629!-- flow_statistics (or - depending on the chosen model run - are never initialized)
630    sums_divnew_l  = 0.0_wp
631    sums_divold_l  = 0.0_wp
632    sums_l_l       = 0.0_wp
633    sums_wsts_bc_l = 0.0_wp
634
635!
636!-- Initialize model variables
637    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
638         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
639!
640!--    Initialization with provided input data derived from larger-scale model
641       IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
642          IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'start' )
643!
644!--       Read initial 1D profiles or 3D data from NetCDF file, depending on the provided
645!--       level-of-detail.
646!--       At the moment, only u, v, w, pt and q are provided.
647          CALL netcdf_data_input_init_3d
648!
649!--       Please note, Inifor provides data from nzb+1 to nzt.
650!--       Bottom and top boundary conditions for Inifor profiles are already set (just after
651!--       reading), so that this is not necessary here.
652!--       Depending on the provided level-of-detail, initial Inifor data is either stored on data
653!--       type (lod=1), or directly on 3D arrays (lod=2).
654!--       In order to obtain also initial profiles in case of lod=2 (which is required for e.g.
655!--       damping), average over 3D data.
656          IF( init_3d%lod_u == 1 )  THEN
657             u_init = init_3d%u_init
658          ELSEIF( init_3d%lod_u == 2 )  THEN
659             ALLOCATE( init_l(nzb:nzt+1) )
660             DO  k = nzb, nzt+1
661                init_l(k) = SUM( u(k,nys:nyn,nxl:nxr) )
662             ENDDO
663             init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
664
665#if defined( __parallel )
666             CALL MPI_ALLREDUCE( init_l, u_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
667#else
668             u_init = init_l
669#endif
670             DEALLOCATE( init_l )
671
672          ENDIF
673
674          IF( init_3d%lod_v == 1 )  THEN
675             v_init = init_3d%v_init
676          ELSEIF( init_3d%lod_v == 2 )  THEN
677             ALLOCATE( init_l(nzb:nzt+1) )
678             DO  k = nzb, nzt+1
679                init_l(k) = SUM( v(k,nys:nyn,nxl:nxr) )
680             ENDDO
681             init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
682
683#if defined( __parallel )
684             CALL MPI_ALLREDUCE( init_l, v_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
685#else
686             v_init = init_l
687#endif
688             DEALLOCATE( init_l )
689          ENDIF
690          IF( .NOT. neutral )  THEN
691             IF( init_3d%lod_pt == 1 )  THEN
692                pt_init = init_3d%pt_init
693             ELSEIF( init_3d%lod_pt == 2 )  THEN
694                ALLOCATE( init_l(nzb:nzt+1) )
695                DO  k = nzb, nzt+1
696                   init_l(k) = SUM( pt(k,nys:nyn,nxl:nxr) )
697                ENDDO
698                init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
699
700#if defined( __parallel )
701                CALL MPI_ALLREDUCE( init_l, pt_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
702#else
703                pt_init = init_l
704#endif
705                DEALLOCATE( init_l )
706             ENDIF
707          ENDIF
708
709
710          IF( humidity )  THEN
711             IF( init_3d%lod_q == 1 )  THEN
712                q_init = init_3d%q_init
713             ELSEIF( init_3d%lod_q == 2 )  THEN
714                ALLOCATE( init_l(nzb:nzt+1) )
715                DO  k = nzb, nzt+1
716                   init_l(k) = SUM( q(k,nys:nyn,nxl:nxr) )
717                ENDDO
718                init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
719
720#if defined( __parallel )
721                CALL MPI_ALLREDUCE( init_l, q_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
722#else
723                q_init = init_l
724#endif
725                DEALLOCATE( init_l )
726             ENDIF
727          ENDIF
728
729!
730!--       Write initial profiles onto 3D arrays.
731!--       Work-around, 3D initialization of u,v,w creates artificial structures which correlate with
732!--       the processor grid. The reason for this is still unknown. To work-around this, 3D
733!--       initialization will be effectively reduce to a 1D initialization where no such artificial
734!--       structures appear.
735          DO  i = nxlg, nxrg
736             DO  j = nysg, nyng
737                IF( init_3d%lod_u == 1  .OR.  init_3d%lod_u == 2 )  u(:,j,i) = u_init(:)
738                IF( init_3d%lod_v == 1  .OR.  init_3d%lod_u == 2 )  v(:,j,i) = v_init(:)
739                IF( .NOT. neutral  .AND.  ( init_3d%lod_pt == 1  .OR.  init_3d%lod_pt == 2 ) )     &
740                   pt(:,j,i) = pt_init(:)
741                IF( humidity  .AND.  ( init_3d%lod_q == 1  .OR.  init_3d%lod_q == 2 ) )            &
742                   q(:,j,i) = q_init(:)
743             ENDDO
744          ENDDO
745!
746!--       Set geostrophic wind components.
747          IF ( init_3d%from_file_ug )  THEN
748             ug(:) = init_3d%ug_init(:)
749          ENDIF
750          IF ( init_3d%from_file_vg )  THEN
751             vg(:) = init_3d%vg_init(:)
752          ENDIF
753!
754!--       Set bottom and top boundary condition for geostrophic wind
755          ug(nzt+1) = ug(nzt)
756          vg(nzt+1) = vg(nzt)
757          ug(nzb)   = ug(nzb+1)
758          vg(nzb)   = vg(nzb+1)
759!
760!--       Set inital w to 0
761          w = 0.0_wp
762
763          IF ( passive_scalar )  THEN
764             DO  i = nxlg, nxrg
765                DO  j = nysg, nyng
766                   s(:,j,i) = s_init
767                ENDDO
768             ENDDO
769          ENDIF
770
771!
772!--       Set velocity components at non-atmospheric / oceanic grid points to zero.
773          u = MERGE( u, 0.0_wp, BTEST( wall_flags_total_0, 1 ) )
774          v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )
775          w = MERGE( w, 0.0_wp, BTEST( wall_flags_total_0, 3 ) )
776!
777!--       Initialize surface variables, e.g. friction velocity, momentum fluxes, etc.
778          CALL  init_surfaces
779
780          IF ( debug_output )  CALL  debug_message( 'initializing with INIFOR', 'end' )
781!
782!--    Initialization via computed 1D-model profiles
783       ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
784
785          IF ( debug_output )  CALL  debug_message( 'initializing with 1D model profiles', 'start' )
786!
787!--       Use solutions of the 1D model as initial profiles,
788!--       start 1D model
789          CALL init_1d_model
790!
791!--       Transfer initial profiles to the arrays of the 3D model
792          DO  i = nxlg, nxrg
793             DO  j = nysg, nyng
794                pt(:,j,i) = pt_init
795                u(:,j,i)  = u1d
796                v(:,j,i)  = v1d
797             ENDDO
798          ENDDO
799
800          IF ( humidity )  THEN
801             DO  i = nxlg, nxrg
802                DO  j = nysg, nyng
803                   q(:,j,i) = q_init
804                ENDDO
805             ENDDO
806          ENDIF
807
808          IF ( passive_scalar )  THEN
809             DO  i = nxlg, nxrg
810                DO  j = nysg, nyng
811                   s(:,j,i) = s_init
812                ENDDO
813             ENDDO
814          ENDIF
815!
816!--          Store initial profiles for output purposes etc.
817          IF ( .NOT. constant_diffusion )  THEN
818             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
819          ENDIF
820!
821!--       Set velocities back to zero
822          u = MERGE( u, 0.0_wp, BTEST( wall_flags_total_0, 1 ) )
823          v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )
824!
825!--       WARNING: The extra boundary conditions set after running the 1D model impose an error on
826!--       -------- the divergence one layer below the topography; need to correct later
827!--       ATTENTION: Provisional correction for Piacsek & Williams advection scheme: keep u and v
828!--       ---------- zero one layer below the topography.
829          IF ( ibc_uv_b == 1 )  THEN
830!
831!--          Neumann condition
832             DO  i = nxl-1, nxr+1
833                DO  j = nys-1, nyn+1
834                   u(nzb,j,i) = u(nzb+1,j,i)
835                   v(nzb,j,i) = v(nzb+1,j,i)
836                ENDDO
837             ENDDO
838
839          ENDIF
840!
841!--       Initialize surface variables, e.g. friction velocity, momentum fluxes, etc.
842          CALL init_surfaces
843
844          IF ( debug_output )  CALL  debug_message( 'initializing with 1D model profiles', 'end' )
845
846       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
847
848          IF ( debug_output )  CALL  debug_message( 'initializing with constant profiles', 'start' )
849
850!
851!--       Use constructed initial profiles (velocity constant with height, temperature profile with
852!--       constant gradient)
853          DO  i = nxlg, nxrg
854             DO  j = nysg, nyng
855                pt(:,j,i) = pt_init
856                u(:,j,i)  = u_init
857                v(:,j,i)  = v_init
858             ENDDO
859          ENDDO
860!
861!--       Mask topography
862          u = MERGE( u, 0.0_wp, BTEST( wall_flags_total_0, 1 ) )
863          v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )
864!
865!--       Set initial horizontal velocities at the lowest computational grid levels to zero in order
866!--       to avoid too small time steps caused by the diffusion limit in the initial phase of a run
867!--       (at k=1, dz/2 occurs in the limiting formula!).
868!--       Please note, in case land- or urban-surface model is used and a spinup is applied, masking
869!--       the lowest computational level is not possible as MOST as well as energy-balance
870!--       parametrizations will not work with zero wind velocity.
871          IF ( ibc_uv_b /= 1  .AND.  .NOT. spinup )  THEN
872             DO  i = nxlg, nxrg
873                DO  j = nysg, nyng
874                   DO  k = nzb, nzt
875                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 20 ) )
876                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 21 ) )
877                   ENDDO
878                ENDDO
879             ENDDO
880          ENDIF
881
882          IF ( humidity )  THEN
883             DO  i = nxlg, nxrg
884                DO  j = nysg, nyng
885                   q(:,j,i) = q_init
886                ENDDO
887             ENDDO
888          ENDIF
889
890          IF ( passive_scalar )  THEN
891             DO  i = nxlg, nxrg
892                DO  j = nysg, nyng
893                   s(:,j,i) = s_init
894                ENDDO
895             ENDDO
896          ENDIF
897
898!
899!--       Compute initial temperature field and other constants used in case of a sloping surface.
900          IF ( sloping_surface )  CALL init_slope
901!
902!--       Initialize surface variables, e.g. friction velocity, momentum fluxes, etc.
903          CALL init_surfaces
904
905          IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'end' )
906
907       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )  THEN
908
909          IF ( debug_output )  CALL debug_message( 'initializing by user', 'start' )
910!
911!--       Pre-initialize surface variables, i.e. setting start- and end-indices at each
912!--       (j,i)-location. Please note, this does not supersede user-defined initialization of
913!--       surface quantities.
914          CALL init_surfaces
915!
916!--       Initialization will completely be done by the user
917          CALL user_init_3d_model
918
919          IF ( debug_output )  CALL debug_message( 'initializing by user', 'end' )
920
921       ENDIF
922
923       IF ( debug_output )  THEN
924          CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
925       ENDIF
926
927!
928!--    Bottom boundary
929       IF ( ibc_uv_b == 0  .OR.  ibc_uv_b == 2  )  THEN
930          u(nzb,:,:) = 0.0_wp
931          v(nzb,:,:) = 0.0_wp
932       ENDIF
933
934!
935!--    Apply channel flow boundary condition
936       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
937          u(nzt+1,:,:) = 0.0_wp
938          v(nzt+1,:,:) = 0.0_wp
939       ENDIF
940
941!
942!--    Calculate virtual potential temperature
943       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
944
945!
946!--    Store initial profiles for output purposes etc.. Please note, in case of initialization of u,
947!--    v, w, pt, and q via output data derived from larger scale models, data will not be
948!--    horizontally homogeneous. Actually, a mean profile should be calculated before.
949       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
950       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
951       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
952          hom(nzb,1,5,:) = 0.0_wp
953          hom(nzb,1,6,:) = 0.0_wp
954       ENDIF
955       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
956
957       IF ( humidity )  THEN
958!
959!--       Store initial profile of total water content, virtual potential temperature
960          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
961          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
962!
963!--       Store initial profile of mixing ratio and potential temperature
964          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
965             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
966             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
967          ENDIF
968       ENDIF
969
970!
971!--    Store initial scalar profile
972       IF ( passive_scalar )  THEN
973          hom(:,1,121,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
974       ENDIF
975
976!
977!--    Initialize the random number generators (from numerical recipes)
978       CALL random_function_ini
979
980       IF ( random_generator == 'random-parallel' )  THEN
981          CALL init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
982       ENDIF
983!
984!--    Set the reference state to be used in the buoyancy terms (for ocean runs the reference state
985!--    will be set (overwritten) in init_ocean).
986       IF ( use_single_reference_value )  THEN
987          IF ( .NOT. humidity )  THEN
988             ref_state(:) = pt_reference
989          ELSE
990             ref_state(:) = vpt_reference
991          ENDIF
992       ELSE
993          IF ( .NOT. humidity )  THEN
994             ref_state(:) = pt_init(:)
995          ELSE
996             ref_state(:) = vpt(:,nys,nxl)
997          ENDIF
998       ENDIF
999
1000!
1001!--    For the moment, vertical velocity is zero
1002       w = 0.0_wp
1003
1004!
1005!--    Initialize array sums (must be defined in first call of pres)
1006       sums = 0.0_wp
1007
1008!
1009!--    In case of iterative solvers, p must get an initial value
1010       IF ( psolver(1:9) == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
1011!
1012!--    Impose vortex with vertical axis on the initial velocity profile
1013       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
1014          CALL init_rankine
1015       ENDIF
1016
1017!
1018!--    Impose temperature anomaly (advection test only) or warm air bubble close to surface.
1019       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0  .OR.                          &
1020            INDEX( initializing_actions, 'initialize_bubble' ) /= 0  )  THEN
1021          CALL init_pt_anomaly
1022       ENDIF
1023
1024!
1025!--    If required, change the surface temperature at the start of the 3D run
1026       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
1027          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
1028       ENDIF
1029
1030!
1031!--    If required, change the surface humidity/scalar at the start of the 3D
1032!--    run
1033       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )                                  &
1034          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
1035
1036       IF ( passive_scalar  .AND.  s_surface_initial_change /= 0.0_wp )                            &
1037          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
1038
1039
1040!
1041!--    Initialize old and new time levels.
1042       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1043       pt_p = pt; u_p = u; v_p = v; w_p = w
1044
1045       IF ( humidity  )  THEN
1046          tq_m = 0.0_wp
1047          q_p = q
1048       ENDIF
1049
1050       IF ( passive_scalar )  THEN
1051          ts_m = 0.0_wp
1052          s_p  = s
1053       ENDIF
1054
1055       IF ( debug_output )  THEN
1056          CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
1057       ENDIF
1058
1059    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1060             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
1061    THEN
1062
1063       IF ( debug_output )  THEN
1064          CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
1065       ENDIF
1066!
1067!--    Initialize surface elements and its attributes, e.g. heat- and momentumfluxes, roughness,
1068!--    scaling parameters. As number of surface elements might be different between runs, e.g. in
1069!--    case of cyclic fill, and not all surface elements are read, surface elements need to be
1070!--    initialized before.
1071!--    Please note, in case of cyclic fill, surfaces should be initialized after restart data is
1072!--    read, else, individual settings of surface parameters will be overwritten from data of
1073!--    precursor run, hence, init_surfaces is called a second time after reading the restart data.
1074       CALL init_surfaces
1075!
1076!--    When reading data for cyclic fill of 3D prerun data files, read some of the global variables
1077!--    from the restart file which are required for initializing the inflow
1078       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1079
1080          DO  i = 0, io_blocks-1
1081             IF ( i == io_group )  THEN
1082                CALL rrd_read_parts_of_global
1083             ENDIF
1084#if defined( __parallel )
1085             CALL MPI_BARRIER( comm2d, ierr )
1086#endif
1087          ENDDO
1088
1089       ENDIF
1090
1091!
1092!--    Read processor specific binary data from restart file
1093       DO  i = 0, io_blocks-1
1094          IF ( i == io_group )  THEN
1095             CALL rrd_local
1096          ENDIF
1097#if defined( __parallel )
1098          CALL MPI_BARRIER( comm2d, ierr )
1099#endif
1100       ENDDO
1101
1102
1103       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1104
1105!
1106!--       In case of cyclic fill, call init_surfaces a second time, so that surface properties such
1107!--       as heat fluxes are initialized as prescribed.
1108          CALL init_surfaces
1109
1110!
1111!--       Overwrite u_init, v_init, pt_init, q_init and s_init with the horizontally mean (hom)
1112!--       vertical profiles from the end of the prerun, because these profiles shall be used as the
1113!--       reference state for the rayleigh damping and the pt_damping. This is especially important
1114!--       for the use of large_scale_subsidence, because the reference temperature in the free
1115!--       atmosphere changes in time.
1116          u_init(:) = hom_sum(:,1,0)
1117          v_init(:) = hom_sum(:,2,0)
1118          pt_init(:) = hom_sum(:,4,0)
1119          IF ( humidity )  q_init(:) = hom_sum(:,41,0)
1120          IF ( passive_scalar )  s_init(:) = hom_sum(:,115,0)
1121       ENDIF
1122!
1123!--    In case of complex terrain and cyclic fill method as initialization, shift initial data in
1124!--    the vertical direction for each point in the x-y-plane depending on local surface height.
1125       IF ( complex_terrain  .AND.  TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1126          DO  i = nxlg, nxrg
1127             DO  j = nysg, nyng
1128                nz_u_shift = topo_top_ind(j,i,1)
1129                nz_v_shift = topo_top_ind(j,i,2)
1130                nz_w_shift = topo_top_ind(j,i,3)
1131                nz_s_shift = topo_top_ind(j,i,0)
1132
1133                u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)
1134
1135                v(nz_v_shift:nzt+1,j,i)  = v(0:nzt+1-nz_v_shift,j,i)
1136
1137                w(nz_w_shift:nzt+1,j,i)  = w(0:nzt+1-nz_w_shift,j,i)
1138
1139                p(nz_s_shift:nzt+1,j,i)  =  p(0:nzt+1-nz_s_shift,j,i)
1140                pt(nz_s_shift:nzt+1,j,i) = pt(0:nzt+1-nz_s_shift,j,i)
1141             ENDDO
1142          ENDDO
1143       ENDIF
1144!
1145!--    Initialization of the turbulence recycling method
1146       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  turbulent_inflow )  THEN
1147!
1148!--       First store the profiles to be used at the inflow.
1149!--       These profiles are the (temporally) and horizontally averaged vertical profiles from the
1150!--       prerun. Alternatively, prescribed profiles for u,v-components can be used.
1151          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
1152
1153          IF ( use_prescribed_profile_data )  THEN
1154             mean_inflow_profiles(:,1) = u_init            ! u
1155             mean_inflow_profiles(:,2) = v_init            ! v
1156          ELSE
1157             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
1158             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
1159          ENDIF
1160          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
1161          IF ( humidity )  mean_inflow_profiles(:,6) = hom_sum(:,41,0)          ! q
1162          IF ( passive_scalar )  mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
1163
1164!
1165!--       In case of complex terrain, determine vertical displacement at inflow boundary and adjust
1166!--       mean inflow profiles
1167          IF ( complex_terrain )  THEN
1168             IF ( nxlg <= 0  .AND.  nxrg >= 0  .AND.  nysg <= 0  .AND.  nyng >= 0 )  THEN
1169                nz_u_shift_l = topo_top_ind(j,i,1)
1170                nz_v_shift_l = topo_top_ind(j,i,2)
1171                nz_w_shift_l = topo_top_ind(j,i,3)
1172                nz_s_shift_l = topo_top_ind(j,i,0)
1173             ELSE
1174                nz_u_shift_l = 0
1175                nz_v_shift_l = 0
1176                nz_w_shift_l = 0
1177                nz_s_shift_l = 0
1178             ENDIF
1179
1180#if defined( __parallel )
1181             CALL MPI_ALLREDUCE( nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
1182             CALL MPI_ALLREDUCE( nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
1183             CALL MPI_ALLREDUCE( nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
1184             CALL MPI_ALLREDUCE( nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
1185#else
1186             nz_u_shift = nz_u_shift_l
1187             nz_v_shift = nz_v_shift_l
1188             nz_w_shift = nz_w_shift_l
1189             nz_s_shift = nz_s_shift_l
1190#endif
1191
1192             mean_inflow_profiles(:,1) = 0.0_wp
1193             mean_inflow_profiles(nz_u_shift:nzt+1,1) = hom_sum(0:nzt+1-nz_u_shift,1,0)  ! u
1194
1195             mean_inflow_profiles(:,2) = 0.0_wp
1196             mean_inflow_profiles(nz_v_shift:nzt+1,2) = hom_sum(0:nzt+1-nz_v_shift,2,0)  ! v
1197
1198             mean_inflow_profiles(nz_s_shift:nzt+1,4) = hom_sum(0:nzt+1-nz_s_shift,4,0)  ! pt
1199
1200          ENDIF
1201
1202!
1203!--       If necessary, adjust the horizontal flow field to the prescribed profiles
1204          IF ( use_prescribed_profile_data )  THEN
1205             DO  i = nxlg, nxrg
1206                DO  j = nysg, nyng
1207                   DO  k = nzb, nzt+1
1208                      u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
1209                      v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
1210                   ENDDO
1211                ENDDO
1212             ENDDO
1213          ENDIF
1214
1215!
1216!--       Use these mean profiles at the inflow (provided that Dirichlet conditions are used)
1217          IF ( bc_dirichlet_l )  THEN
1218             DO  j = nysg, nyng
1219                DO  k = nzb, nzt+1
1220                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
1221                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
1222                   w(k,j,nxlg:-1)  = 0.0_wp
1223                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
1224                   IF ( humidity )  q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
1225                   IF ( passive_scalar )  s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)
1226                ENDDO
1227             ENDDO
1228          ENDIF
1229
1230!
1231!--       Calculate the damping factors to be used at the inflow. For a turbulent inflow the
1232!--       turbulent fluctuations have to be limited vertically because otherwise the turbulent
1233!--       inflow layer will grow in time.
1234          IF ( inflow_damping_height == 9999999.9_wp )  THEN
1235!
1236!--          Default: use the inversion height calculated by the prerun; if this is zero,
1237!--          inflow_damping_height must be explicitly specified.
1238             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
1239                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
1240             ELSE
1241                WRITE( message_string, * ) 'inflow_damping_height must be ',                       &
1242                                           'explicitly specified because&the inversion height ',   &
1243                                           'calculated by the prerun is zero.'
1244                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1245             ENDIF
1246
1247          ENDIF
1248
1249          IF ( inflow_damping_width == 9999999.9_wp )  THEN
1250!
1251!--          Default for the transition range: one tenth of the undamped layer
1252             inflow_damping_width = 0.1_wp * inflow_damping_height
1253
1254          ENDIF
1255
1256          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1257
1258          DO  k = nzb, nzt+1
1259
1260             IF ( zu(k) <= inflow_damping_height )  THEN
1261                inflow_damping_factor(k) = 1.0_wp
1262             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
1263                inflow_damping_factor(k) = 1.0_wp -                                                &
1264                                           ( zu(k) - inflow_damping_height ) / inflow_damping_width
1265             ELSE
1266                inflow_damping_factor(k) = 0.0_wp
1267             ENDIF
1268
1269          ENDDO
1270
1271       ENDIF
1272
1273!
1274!--    Inside buildings set velocities back to zero
1275       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  topography /= 'flat' )  THEN
1276!
1277!--       Inside buildings set velocities back to zero.
1278!--       Other scalars (pt, q, s, p, sa, ...) are ignored at present,
1279!--       maybe revise later.
1280          DO  i = nxlg, nxrg
1281             DO  j = nysg, nyng
1282                DO  k = nzb, nzt
1283                   u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
1284                   v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
1285                   w(k,j,i) = MERGE( w(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
1286                ENDDO
1287             ENDDO
1288          ENDDO
1289
1290       ENDIF
1291
1292!
1293!--    Calculate initial temperature field and other constants used in case of a sloping surface
1294       IF ( sloping_surface )  CALL init_slope
1295
1296!
1297!--    Initialize new time levels (only done in order to set boundary values including ghost points)
1298       pt_p = pt; u_p = u; v_p = v; w_p = w
1299       IF ( humidity )  THEN
1300          q_p = q
1301       ENDIF
1302       IF ( passive_scalar )  s_p  = s
1303!
1304!--    Allthough tendency arrays are set in prognostic_equations, they have have to be predefined
1305!--    here because they are used (but multiplied with 0) there before they are set.
1306       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1307       IF ( humidity )  THEN
1308          tq_m = 0.0_wp
1309       ENDIF
1310       IF ( passive_scalar )  ts_m  = 0.0_wp
1311
1312       IF ( debug_output )  THEN
1313          CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
1314       ENDIF
1315
1316    ELSE
1317!
1318!--    Actually this part of the programm should not be reached
1319       message_string = 'unknown initializing problem'
1320       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1321    ENDIF
1322
1323
1324    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1325!
1326!--    Initialize old timelevels needed for radiation boundary conditions
1327       IF ( bc_radiation_l )  THEN
1328          u_m_l(:,:,:) = u(:,:,1:2)
1329          v_m_l(:,:,:) = v(:,:,0:1)
1330          w_m_l(:,:,:) = w(:,:,0:1)
1331       ENDIF
1332       IF ( bc_radiation_r )  THEN
1333          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1334          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1335          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1336       ENDIF
1337       IF ( bc_radiation_s )  THEN
1338          u_m_s(:,:,:) = u(:,0:1,:)
1339          v_m_s(:,:,:) = v(:,1:2,:)
1340          w_m_s(:,:,:) = w(:,0:1,:)
1341       ENDIF
1342       IF ( bc_radiation_n )  THEN
1343          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1344          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1345          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1346       ENDIF
1347
1348    ENDIF
1349
1350!
1351!-- Calculate the initial volume flow at the right and north boundary
1352    IF ( conserve_volume_flow )  THEN
1353
1354       IF ( use_prescribed_profile_data )  THEN
1355
1356          volume_flow_initial_l = 0.0_wp
1357          volume_flow_area_l    = 0.0_wp
1358
1359          IF ( nxr == nx )  THEN
1360             DO  j = nys, nyn
1361                DO  k = nzb+1, nzt
1362                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +                           &
1363                                              u_init(k) * dzw(k)                                   &
1364                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1365                                                       BTEST( wall_flags_total_0(k,j,nxr), 1 )     &
1366                                                     )
1367
1368                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)                       &
1369                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1370                                                       BTEST( wall_flags_total_0(k,j,nxr), 1 )     &
1371                                                     )
1372                ENDDO
1373             ENDDO
1374          ENDIF
1375
1376          IF ( nyn == ny )  THEN
1377             DO  i = nxl, nxr
1378                DO  k = nzb+1, nzt
1379                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +                           &
1380                                              v_init(k) * dzw(k)                                   &
1381                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1382                                                       BTEST( wall_flags_total_0(k,nyn,i), 2 )     &
1383                                                     )
1384                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)                       &
1385                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1386                                                       BTEST( wall_flags_total_0(k,nyn,i), 2 )     &
1387                                                     )
1388                ENDDO
1389             ENDDO
1390          ENDIF
1391
1392#if defined( __parallel )
1393          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), 2, MPI_REAL,       &
1394                              MPI_SUM, comm2d, ierr )
1395          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), 2, MPI_REAL, MPI_SUM,    &
1396                              comm2d, ierr )
1397
1398#else
1399          volume_flow_initial = volume_flow_initial_l
1400          volume_flow_area    = volume_flow_area_l
1401#endif
1402
1403       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1404
1405          volume_flow_initial_l = 0.0_wp
1406          volume_flow_area_l    = 0.0_wp
1407
1408          IF ( nxr == nx )  THEN
1409             DO  j = nys, nyn
1410                DO  k = nzb+1, nzt
1411                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +                           &
1412                                              hom_sum(k,1,0) * dzw(k)                              &
1413                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1414                                                       BTEST( wall_flags_total_0(k,j,nx), 1 )      &
1415                                                     )
1416                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)                       &
1417                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1418                                                       BTEST( wall_flags_total_0(k,j,nx), 1 )      &
1419                                                     )
1420                ENDDO
1421             ENDDO
1422          ENDIF
1423
1424          IF ( nyn == ny )  THEN
1425             DO  i = nxl, nxr
1426                DO  k = nzb+1, nzt
1427                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +                           &
1428                                              hom_sum(k,2,0) * dzw(k)                              &
1429                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1430                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
1431                                                     )
1432                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)                       &
1433                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1434                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
1435                                                     )
1436                ENDDO
1437             ENDDO
1438          ENDIF
1439
1440#if defined( __parallel )
1441          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), 2, MPI_REAL,       &
1442                              MPI_SUM, comm2d, ierr )
1443          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), 2, MPI_REAL, MPI_SUM,    &
1444                              comm2d, ierr )
1445
1446#else
1447          volume_flow_initial = volume_flow_initial_l
1448          volume_flow_area    = volume_flow_area_l
1449#endif
1450
1451       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1452
1453          volume_flow_initial_l = 0.0_wp
1454          volume_flow_area_l    = 0.0_wp
1455
1456          IF ( nxr == nx )  THEN
1457             DO  j = nys, nyn
1458                DO  k = nzb+1, nzt
1459                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +                           &
1460                                              u(k,j,nx) * dzw(k)                                   &
1461                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1462                                                       BTEST( wall_flags_total_0(k,j,nx), 1 )      &
1463                                                     )
1464                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)                       &
1465                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1466                                                   BTEST( wall_flags_total_0(k,j,nx), 1 )          &
1467                                                     )
1468                ENDDO
1469             ENDDO
1470          ENDIF
1471
1472          IF ( nyn == ny )  THEN
1473             DO  i = nxl, nxr
1474                DO  k = nzb+1, nzt
1475                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +                           &
1476                                              v(k,ny,i) * dzw(k)                                   &
1477                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1478                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
1479                                                     )
1480                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)                       &
1481                                              * MERGE( 1.0_wp, 0.0_wp,                             &
1482                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
1483                                                     )
1484                ENDDO
1485             ENDDO
1486          ENDIF
1487
1488#if defined( __parallel )
1489          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), 2, MPI_REAL,       &
1490                              MPI_SUM, comm2d, ierr )
1491          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), 2, MPI_REAL, MPI_SUM,    &
1492                              comm2d, ierr )
1493
1494#else
1495          volume_flow_initial = volume_flow_initial_l
1496          volume_flow_area    = volume_flow_area_l
1497#endif
1498
1499       ENDIF
1500
1501!
1502!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated from u|v_bulk instead
1503       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
1504          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
1505          volume_flow_initial(2) = v_bulk * volume_flow_area(2)
1506       ENDIF
1507
1508    ENDIF
1509!
1510!-- In the following, surface properties can be further initialized with input from static driver
1511!-- file.
1512!-- At the moment this affects only default surfaces. For example, roughness length or sensible /
1513!-- latent heat fluxes can be initialized heterogeneously for default surfaces. Therefore, a generic
1514!-- routine from netcdf_data_input_mod is called to read a 2D array.
1515    IF ( input_pids_static )  THEN
1516!
1517!--    Allocate memory for possible static input
1518       ALLOCATE( tmp_2d%var(nys:nyn,nxl:nxr) )
1519       tmp_2d%var = 0.0_wp
1520!
1521!--    Open the static input file
1522#if defined( __netcdf )
1523       CALL open_read_file( TRIM( input_file_static ) //                                           &
1524                            TRIM( coupling_char ), pids_id )
1525
1526       CALL inquire_num_variables( pids_id, num_var_pids )
1527!
1528!--    Allocate memory to store variable names and read them
1529       ALLOCATE( vars_pids(1:num_var_pids) )
1530       CALL inquire_variable_names( pids_id, vars_pids )
1531!
1532!--    Input roughness length.
1533       IF ( check_existence( vars_pids, 'z0' ) )  THEN
1534!
1535!--       Read _FillValue attribute
1536          CALL get_attribute( pids_id, char_fill, tmp_2d%fill, .FALSE., 'z0' )
1537!
1538!--       Read variable
1539          CALL get_variable( pids_id, 'z0', tmp_2d%var, nxl, nxr, nys, nyn )
1540!
1541!--       Initialize roughness length. Note, z0 will be only initialized at default-type surfaces.
1542!--       At natural or urban z0 is implicitly initialized by the respective parameter lists.
1543!--       Initialize horizontal surface elements.
1544          CALL init_single_surface_properties( surf_def_h(0)%z0, tmp_2d%var, surf_def_h(0)%ns,     &
1545                                               tmp_2d%fill, surf_def_h(0)%i, surf_def_h(0)%j )
1546!
1547!--       Initialize roughness also at vertical surface elements.
1548!--       Note, the actual 2D input arrays are only defined on the subdomain. Therefore, pass the
1549!--       index arrays with their respective offset values.
1550          DO  l = 0, 3
1551             CALL init_single_surface_properties( surf_def_v(l)%z0, tmp_2d%var, surf_def_v(l)%ns,  &
1552                                                  tmp_2d%fill, surf_def_v(l)%i+surf_def_v(l)%ioff, &
1553                                                  surf_def_v(l)%j+surf_def_v(l)%joff )
1554          ENDDO
1555
1556       ENDIF
1557!
1558!--    Input surface sensible heat flux.
1559       IF ( check_existence( vars_pids, 'shf' ) )  THEN
1560!
1561!--       Read _FillValue attribute
1562          CALL get_attribute( pids_id, char_fill, tmp_2d%fill, .FALSE., 'shf' )
1563!
1564!--       Read variable
1565          CALL get_variable( pids_id, 'shf', tmp_2d%var, nxl, nxr, nys, nyn )
1566!
1567!--       Initialize heat flux. Note, shf will be only initialized at default-type surfaces. At
1568!--       natural or urban shf is implicitly initialized by the respective parameter lists.
1569!--       Initialize horizontal surface elements.
1570          CALL init_single_surface_properties( surf_def_h(0)%shf, tmp_2d%var, surf_def_h(0)%ns,    &
1571                                               tmp_2d%fill, surf_def_h(0)%i, surf_def_h(0)%j )
1572!
1573!--       Initialize heat flux also at vertical surface elements.
1574!--       Note, the actual 2D input arrays are only defined on the subdomain. Therefore, pass the
1575!--       index arrays with their respective offset values.
1576          DO  l = 0, 3
1577             CALL init_single_surface_properties( surf_def_v(l)%shf, tmp_2d%var, surf_def_v(l)%ns, &
1578                                                  tmp_2d%fill, surf_def_v(l)%i+surf_def_v(l)%ioff, &
1579                                                  surf_def_v(l)%j+surf_def_v(l)%joff )
1580          ENDDO
1581
1582       ENDIF
1583!
1584!--    Input surface sensible heat flux.
1585       IF ( check_existence( vars_pids, 'qsws' ) )  THEN
1586!
1587!--       Read _FillValue attribute
1588          CALL get_attribute( pids_id, char_fill, tmp_2d%fill,                 &
1589                              .FALSE., 'qsws' )
1590!
1591!--       Read variable
1592          CALL get_variable( pids_id, 'qsws', tmp_2d%var,                      &
1593                             nxl, nxr, nys, nyn )
1594!
1595!--       Initialize latent heat flux. Note, qsws will be only initialized at default-type surfaces.
1596!--       At natural or urban qsws is implicitly initialized by the respective parameter lists.
1597!--       Initialize horizontal surface elements.
1598          CALL init_single_surface_properties( surf_def_h(0)%qsws, tmp_2d%var, surf_def_h(0)%ns,   &
1599                                               tmp_2d%fill, surf_def_h(0)%i, surf_def_h(0)%j )
1600!
1601!--       Initialize latent heat flux also at vertical surface elements.
1602!--       Note, the actual 2D input arrays are only defined on the subdomain. Therefore, pass the
1603!--       index arrays with their respective offset values.
1604          DO  l = 0, 3
1605             CALL init_single_surface_properties( surf_def_v(l)%qsws, tmp_2d%var, surf_def_v(l)%ns,&
1606                                                  tmp_2d%fill, surf_def_v(l)%i+surf_def_v(l)%ioff, &
1607                                                  surf_def_v(l)%j+surf_def_v(l)%joff )
1608          ENDDO
1609
1610       ENDIF
1611!
1612!--    Additional variables, can be initialized the
1613!--    same way.
1614
1615!
1616!--    Finally, close the input file and deallocate temporary arrays
1617       DEALLOCATE( vars_pids )
1618
1619       CALL close_input_file( pids_id )
1620#endif
1621       DEALLOCATE( tmp_2d%var )
1622    ENDIF
1623!
1624!-- Finally, if random_heatflux is set, disturb shf at horizontal surfaces. Actually, this should be
1625!-- done in surface_mod, where all other initializations of surface quantities are done. However,
1626!-- this would create a ring dependency, hence, it is done here. Maybe delete disturb_heatflux and
1627!-- tranfer the respective code directly into the initialization in surface_mod.
1628    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
1629         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1630
1631       IF ( use_surface_fluxes  .AND.  constant_heatflux  .AND.  random_heatflux )  THEN
1632          IF ( surf_def_h(0)%ns >= 1 )  CALL disturb_heatflux( surf_def_h(0) )
1633          IF ( surf_lsm_h%ns    >= 1 )  CALL disturb_heatflux( surf_lsm_h    )
1634          IF ( surf_usm_h%ns    >= 1 )  CALL disturb_heatflux( surf_usm_h    )
1635       ENDIF
1636    ENDIF
1637
1638!
1639!-- Compute total sum of grid points and the mean surface level height for each statistic region.
1640!-- These are mainly used for horizontal averaging of turbulence statistics.
1641!-- ngp_2dh: number of grid points of a horizontal cross section through the respective statistic
1642!--          region
1643!-- ngp_3d:  number of grid points of the respective statistic region
1644    ngp_2dh_outer_l   = 0
1645    ngp_2dh_outer     = 0
1646    ngp_2dh_s_inner_l = 0
1647    ngp_2dh_s_inner   = 0
1648    ngp_2dh_l         = 0
1649    ngp_2dh           = 0
1650    ngp_3d_inner_l    = 0.0_wp
1651    ngp_3d_inner      = 0
1652    ngp_3d            = 0
1653    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
1654
1655    mean_surface_level_height   = 0.0_wp
1656    mean_surface_level_height_l = 0.0_wp
1657!
1658!-- To do: New concept for these non-topography grid points!
1659    DO  sr = 0, statistic_regions
1660       DO  i = nxl, nxr
1661          DO  j = nys, nyn
1662             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
1663!
1664!--             All xy-grid points
1665                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
1666!
1667!--             Determine mean surface-level height. In case of downward-facing walls are present,
1668!--             more than one surface level exist.
1669!--             In this case, use the lowest surface-level height.
1670                IF ( surf_def_h(0)%start_index(j,i) <= surf_def_h(0)%end_index(j,i) )  THEN
1671                   m = surf_def_h(0)%start_index(j,i)
1672                   k = surf_def_h(0)%k(m)
1673                   mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + zw(k-1)
1674                ENDIF
1675                IF ( surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) )  THEN
1676                   m = surf_lsm_h%start_index(j,i)
1677                   k = surf_lsm_h%k(m)
1678                   mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + zw(k-1)
1679                ENDIF
1680                IF ( surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) )  THEN
1681                   m = surf_usm_h%start_index(j,i)
1682                   k = surf_usm_h%k(m)
1683                   mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + zw(k-1)
1684                ENDIF
1685
1686                k_surf = k - 1
1687
1688                DO  k = nzb, nzt+1
1689!
1690!--                xy-grid points above topography
1691                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr)     +                             &
1692                                           MERGE( 1, 0, BTEST( wall_flags_total_0(k,j,i), 24 ) )
1693
1694                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) +                             &
1695                                             MERGE( 1, 0, BTEST( wall_flags_total_0(k,j,i), 22 ) )
1696
1697                ENDDO
1698!
1699!--             All grid points of the total domain above topography
1700                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + ( nz - k_surf + 2 )
1701
1702
1703
1704             ENDIF
1705          ENDDO
1706       ENDDO
1707    ENDDO
1708
1709    sr = statistic_regions + 1
1710#if defined( __parallel )
1711    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1712    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1713    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1714    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, MPI_INTEGER, MPI_SUM, &
1715                        comm2d, ierr )
1716    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1717    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), (nz+2)*sr, MPI_INTEGER,      &
1718                        MPI_SUM, comm2d, ierr )
1719    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1720    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, MPI_SUM, comm2d,     &
1721                        ierr )
1722    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
1723    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1724    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0), mean_surface_level_height(0), sr, MPI_REAL,&
1725                        MPI_SUM, comm2d, ierr )
1726    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
1727#else
1728    ngp_2dh         = ngp_2dh_l
1729    ngp_2dh_outer   = ngp_2dh_outer_l
1730    ngp_2dh_s_inner = ngp_2dh_s_inner_l
1731    ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
1732    mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
1733#endif
1734
1735    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) *                                     &
1736             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
1737
1738!
1739!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics, buoyancy, etc. A
1740!-- zero value will occur for cases where all grid points of the respective subdomain lie below the
1741!-- surface topography
1742    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   )
1743    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), ngp_3d_inner(:) )
1744    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
1745
1746    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l,           &
1747                ngp_3d_inner_tmp )
1748
1749!
1750!-- Initializae 3D offline nesting in COSMO model and read data from
1751!-- external NetCDF file.
1752    IF ( nesting_offline )  CALL nesting_offl_init
1753!
1754!-- Initialize quantities for special advections schemes
1755    CALL init_advec
1756
1757!
1758!-- Impose random perturbation on the horizontal velocity field and then
1759!-- remove the divergences from the velocity field at the initial stage
1760    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND.                     &
1761         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
1762         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1763
1764       IF ( debug_output )  THEN
1765          CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
1766       ENDIF
1767!
1768!--    Needed for both disturb_field and pres
1769!$ACC DATA &
1770!$ACC CREATE(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1771!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1772!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
1773
1774       CALL disturb_field( 'u', tend, u )
1775       CALL disturb_field( 'v', tend, v )
1776
1777!$ACC DATA &
1778!$ACC CREATE(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
1779!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1780!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1781!$ACC COPYIN(rho_air(nzb:nzt+1), rho_air_zw(nzb:nzt+1)) &
1782!$ACC COPYIN(ddzu(1:nzt+1), ddzw(1:nzt+1)) &
1783!$ACC COPYIN(wall_flags_total_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1784!$ACC COPYIN(bc_h(0:1)) &
1785!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
1786!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
1787!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
1788!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
1789!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
1790!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
1791
1792       n_sor = nsor_ini
1793       CALL pres
1794       n_sor = nsor
1795
1796!$ACC END DATA
1797!$ACC END DATA
1798
1799       IF ( debug_output )  THEN
1800          CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
1801       ENDIF
1802
1803    ENDIF
1804
1805    IF ( .NOT. ocean_mode )  THEN
1806
1807       ALLOCATE( hyp(nzb:nzt+1) )
1808       ALLOCATE( d_exner(nzb:nzt+1) )
1809       ALLOCATE( exner(nzb:nzt+1) )
1810       ALLOCATE( hyrho(nzb:nzt+1) )
1811!
1812!--    Check temperature in case of too large domain height
1813       DO  k = nzb, nzt+1
1814          IF ( ( pt_surface * exner_function( surface_pressure * 100.0_wp ) - g/c_p * zu(k) )      &
1815                 < 0.0_wp )  THEN
1816             WRITE( message_string, * )  'absolute temperature < 0.0 at zu(', k, ') = ', zu(k)
1817             CALL message( 'init_3d_model', 'PA0142', 1, 2, 0, 6, 0 )
1818          ENDIF
1819       ENDDO
1820
1821!
1822!--    Calculate vertical profile of the hydrostatic pressure (hyp)
1823       hyp    = barometric_formula( zu, pt_surface * exner_function( surface_pressure * 100.0_wp ),&
1824                                    surface_pressure * 100.0_wp )
1825       d_exner = exner_function_invers( hyp )
1826       exner = 1.0_wp / exner_function_invers( hyp )
1827       hyrho  = ideal_gas_law_rho_pt( hyp, pt_init )
1828!
1829!--    Compute reference density
1830       rho_surface = ideal_gas_law_rho( surface_pressure * 100.0_wp,                               &
1831                                        pt_surface * exner_function( surface_pressure * 100.0_wp ) )
1832
1833    ENDIF
1834
1835!
1836!-- If required, initialize particles
1837    IF ( agents_active )  CALL mas_init
1838!
1839!-- Initialization of synthetic turbulence generator
1840    IF ( use_syn_turb_gen )  CALL stg_init
1841!
1842!-- Initializing actions for all other modules
1843    CALL module_interface_init
1844!
1845!-- Initialize surface layer (done after LSM as roughness length are required for initialization
1846    IF ( constant_flux_layer )  CALL init_surface_layer_fluxes
1847!
1848!-- Initialize surface data output
1849    IF ( surface_output )  CALL surface_data_output_init
1850!
1851!-- Initialize the ws-scheme.
1852    IF ( ws_scheme_sca  .OR.  ws_scheme_mom )  CALL ws_init
1853!
1854!-- Perform post-initializing checks for all other modules
1855    CALL module_interface_init_checks
1856
1857!
1858!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
1859!-- initialize heat-fluxes, etc. via datatype. Revise it later!
1860    IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1861       IF ( use_surface_fluxes  .AND.  constant_heatflux )  THEN
1862          CALL ls_forcing_surf( simulated_time )
1863       ENDIF
1864    ENDIF
1865!
1866!-- Setting weighting factors for calculation of perturbation pressure and turbulent quantities from
1867!-- the RK substeps.
1868    IF ( TRIM( timestep_scheme ) == 'runge-kutta-3' )  THEN      ! for RK3-method
1869
1870       weight_substep(1) = 1._wp/6._wp
1871       weight_substep(2) = 3._wp/10._wp
1872       weight_substep(3) = 8._wp/15._wp
1873
1874       weight_pres(1)    = 1._wp/3._wp
1875       weight_pres(2)    = 5._wp/12._wp
1876       weight_pres(3)    = 1._wp/4._wp
1877
1878    ELSEIF ( TRIM( timestep_scheme ) == 'runge-kutta-2' )  THEN  ! for RK2-method
1879
1880       weight_substep(1) = 1._wp/2._wp
1881       weight_substep(2) = 1._wp/2._wp
1882
1883       weight_pres(1)    = 1._wp/2._wp
1884       weight_pres(2)    = 1._wp/2._wp
1885
1886    ELSE                                     ! for Euler-method
1887
1888       weight_substep(1) = 1.0_wp
1889       weight_pres(1)    = 1.0_wp
1890
1891    ENDIF
1892
1893!
1894!-- Initialize Rayleigh damping factors
1895    rdf    = 0.0_wp
1896    rdf_sc = 0.0_wp
1897    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
1898
1899       IF (  .NOT.  ocean_mode )  THEN
1900          DO  k = nzb+1, nzt
1901             IF ( zu(k) >= rayleigh_damping_height )  THEN
1902                rdf(k) = rayleigh_damping_factor *                                                 &
1903                         ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height )                  &
1904                                / ( zu(nzt) - rayleigh_damping_height ) )                          &
1905                         )**2
1906             ENDIF
1907          ENDDO
1908       ELSE
1909!
1910!--       In ocean mode, rayleigh damping is applied in the lower part of the model domain
1911          DO  k = nzt, nzb+1, -1
1912             IF ( zu(k) <= rayleigh_damping_height )  THEN
1913                rdf(k) = rayleigh_damping_factor *                                                 &
1914                         ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) )                  &
1915                                / ( rayleigh_damping_height - zu(nzb+1) ) )                        &
1916                         )**2
1917             ENDIF
1918          ENDDO
1919       ENDIF
1920
1921    ENDIF
1922    IF ( scalar_rayleigh_damping )  rdf_sc = rdf
1923
1924!
1925!-- Initialize the starting level and the vertical smoothing factor used for the external pressure
1926!-- gradient
1927    dp_smooth_factor = 1.0_wp
1928    IF ( dp_external )  THEN
1929!
1930!--    Set the starting level dp_level_ind_b only if it has not been set before (e.g. in init_grid).
1931       IF ( dp_level_ind_b == 0 )  THEN
1932          ind_array = MINLOC( ABS( dp_level_b - zu ) )
1933          dp_level_ind_b = ind_array(1) - 1 + nzb
1934                                        ! MINLOC uses lower array bound 1
1935       ENDIF
1936       IF ( dp_smooth )  THEN
1937          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
1938          DO  k = dp_level_ind_b+1, nzt
1939             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *                                   &
1940                                             ( REAL( k - dp_level_ind_b, KIND=wp ) /               &
1941                                               REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
1942          ENDDO
1943       ENDIF
1944    ENDIF
1945
1946!
1947!-- Initialize damping zone for the potential temperature in case of non-cyclic lateral boundaries.
1948!-- The damping zone has the maximum value at the inflow boundary and decreases to zero at
1949!-- pt_damping_width.
1950    ptdf_x = 0.0_wp
1951    ptdf_y = 0.0_wp
1952    IF ( bc_lr_dirrad )  THEN
1953       DO  i = nxl, nxr
1954          IF ( ( i * dx ) < pt_damping_width )  THEN
1955             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *                                  &
1956                                                    REAL( pt_damping_width - i * dx, KIND=wp ) /   &
1957                                                    REAL( pt_damping_width, KIND=wp ) ) )**2
1958                                  ENDIF
1959       ENDDO
1960    ELSEIF ( bc_lr_raddir )  THEN
1961       DO  i = nxl, nxr
1962          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
1963             ptdf_x(i) = pt_damping_factor * SIN( pi * 0.5_wp *                                    &
1964                                                  ( ( i - nx ) * dx + pt_damping_width ) /         &
1965                                                  REAL( pt_damping_width, KIND=wp ) )**2
1966          ENDIF
1967       ENDDO
1968    ELSEIF ( bc_ns_dirrad )  THEN
1969       DO  j = nys, nyn
1970          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
1971             ptdf_y(j) = pt_damping_factor * SIN( pi * 0.5_wp *                                    &
1972                                                  ( ( j - ny ) * dy + pt_damping_width ) /         &
1973                                                  REAL( pt_damping_width, KIND=wp ) )**2
1974          ENDIF
1975       ENDDO
1976    ELSEIF ( bc_ns_raddir )  THEN
1977       DO  j = nys, nyn
1978          IF ( ( j * dy ) < pt_damping_width )  THEN
1979             ptdf_y(j) = pt_damping_factor * SIN( pi * 0.5_wp *                                    &
1980                                                  ( pt_damping_width - j * dy ) /                  &
1981                                                  REAL( pt_damping_width, KIND=wp ) )**2
1982          ENDIF
1983       ENDDO
1984    ENDIF
1985
1986!
1987!-- Input binary data file is not needed anymore. This line must be placed after call of user_init!
1988    CALL close_file( 13 )
1989!
1990!-- In case of nesting, put an barrier to assure that all parent and child domains finished
1991!-- initialization.
1992#if defined( __parallel )
1993    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
1994#endif
1995
1996
1997    CALL location_message( 'model initialization', 'finished' )
1998
1999 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.