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

Last change on this file since 4186 was 4186, checked in by suehring, 2 years ago

Enable limitation of Obukhov length so that it does not become zero; to read input data from netcdf in init_3d_model use provided module variables instead of defining local ones; tests updated because changes in Obukhov lengths causes small differences during the initial phase of the run

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