source: palm/trunk/SOURCE/init_3d_model.f90

Last change on this file was 4878, checked in by suehring, 6 months ago

Rename resize_array into add_ghost_layers; remove number of passed indices; replace subroutine calls with interface name

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