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

Last change on this file since 4598 was 4548, checked in by suehring, 11 months ago

Bugfix, move call for lsf_forcing_surf after lsf_init is called

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