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

Last change on this file since 4204 was 4187, checked in by suehring, 5 years ago

radiation: Take external radiation input from root domain dynamic file if no dynamic input file is provided for each nested domain; radiation: Combine MPI_ALLREDUCE calls to reduce latency; land_surface + plant_canopy: Give specific error numbers; land-surface: Adjust error messages for local checks; init_3d_model: Deallocate temporary string array for netcdf-data input since it may be re-used in other modules for different input files

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