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

Last change on this file since 4392 was 4392, checked in by pavelkrc, 5 years ago

Merge branch resler (updated documentation, optional flux tracing, bugfix)

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