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

Last change on this file since 4819 was 4791, checked in by suehring, 4 years ago

Add possibility to initialize surface flux of passive scalar via static driver

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