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

Last change on this file since 4784 was 4783, checked in by raasch, 3 years ago

bugfix for reading restart data with MPI-I/O (does not work with blockwise I/O); missed revision comments from previous commit added

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