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

Last change on this file since 4715 was 4680, checked in by gronemeier, 4 years ago

Add option to fix date and time of the simulation; renamed set_reference_date_time to init_date_time (palm_date_time_mod, init_3d_model, modules, parin)

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