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

Last change on this file since 4262 was 4245, checked in by pavelkrc, 4 years ago

Merge branch resler into trunk

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