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

Last change on this file since 4677 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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