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

Last change on this file since 3004 was 3004, checked in by Giersch, 3 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/forwind/SOURCE/init_3d_model.f901564-1913
    /palm/branches/palm4u/SOURCE/init_3d_model.f902540-2692
File size: 98.2 KB
Line 
1!> @file init_3d_model.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: init_3d_model.f90 3004 2018-04-27 12:33:25Z Giersch $
27! precipitation_rate removed
28!
29! 2995 2018-04-19 12:13:16Z Giersch
30! CALL radiation_control is not necessary during initialization because
31! calculation of radiative fluxes at model start is done in radiation_init
32! in any case
33!
34! 2977 2018-04-17 10:27:57Z kanani
35! Implement changes from branch radiation (r2948-2971) with minor modifications
36! (moh.hefny):
37! - set radiation_interactions according to the existence of urban/land vertical
38!   surfaces and trees to activiate RTM
39! - set average_radiation to TRUE if RTM is activiated
40!
41! 2938 2018-03-27 15:52:42Z suehring
42! - Revise Inifor initialization for geostrophic wind components
43! - Initialize synthetic turbulence generator in case of Inifor initialization 
44!
45! 2936 2018-03-27 14:49:27Z suehring
46! Synchronize parent and child models after initialization.
47! Remove obsolete masking of topography grid points for Runge-Kutta weighted
48! tendency arrays.
49!
50! 2920 2018-03-22 11:22:01Z kanani
51! Add call for precalculating apparent solar positions (moh.hefny)
52!
53! 2906 2018-03-19 08:56:40Z Giersch
54! The variables read/write_svf_on_init have been removed. Instead ENVIRONMENT
55! variables read/write_svf have been introduced. Location_message has been
56! added.
57!
58! 2894 2018-03-15 09:17:58Z Giersch
59! Renamed routines with respect to reading restart data, file 13 is closed in
60! rrd_read_parts_of_global now
61!
62! 2867 2018-03-09 09:40:23Z suehring
63! Further bugfix concerning call of user_init.
64!
65! 2864 2018-03-08 11:57:45Z suehring
66! Bugfix, move call of user_init in front of initialization of grid-point
67! arrays
68!
69! 2817 2018-02-19 16:32:21Z knoop
70! Preliminary gust module interface implemented
71!
72! 2776 2018-01-31 10:44:42Z Giersch
73! Variable use_synthetic_turbulence_generator has been abbreviated
74!
75! 2766 2018-01-22 17:17:47Z kanani
76! Removed preprocessor directive __chem
77!
78! 2758 2018-01-17 12:55:21Z suehring
79! In case of spinup of land- and urban-surface model, do not mask wind velocity
80! at first computational grid level
81!
82! 2746 2018-01-15 12:06:04Z suehring
83! Move flag plant canopy to modules
84!
85! 2718 2018-01-02 08:49:38Z maronga
86! Corrected "Former revisions" section
87!
88! 2705 2017-12-18 11:26:23Z maronga
89! Bugfix for reading initial profiles from ls/nuding file
90!
91! 2701 2017-12-15 15:40:50Z suehring
92! Changes from last commit documented
93!
94! 2700 2017-12-15 14:12:35Z suehring
95! Bugfix, missing initialization of surface attributes in case of
96! inifor-initialization branch
97!
98! 2698 2017-12-14 18:46:24Z suehring
99! Bugfix in get_topography_top_index
100!
101! 2696 2017-12-14 17:12:51Z kanani
102! Change in file header (GPL part)
103! Implementation of uv exposure model (FK)
104! Moved initialisation of diss, e, kh, km to turbulence_closure_mod (TG)
105! Added chemical emissions (FK)
106! Initialize masking arrays and number-of-grid-points arrays before initialize
107! LSM, USM and radiation module
108! Initialization with inifor (MS)
109!
110! 2618 2017-11-16 15:37:30Z suehring
111! Reorder calls of init_surfaces.
112!
113! 2564 2017-10-19 15:56:56Z Giersch
114! Variable wind_turbine was added to control_parameters.
115!
116! 2550 2017-10-16 17:12:01Z boeske
117! Modifications to cyclic fill method and turbulence recycling method in case of
118! complex terrain simulations
119!
120! 2513 2017-10-04 09:24:39Z kanani
121! Bugfix in storing initial scalar profile (wrong index)
122!
123! 2350 2017-08-15 11:48:26Z kanani
124! Bugfix in nopointer version
125!
126! 2339 2017-08-07 13:55:26Z gronemeier
127! corrected timestamp in header
128!
129! 2338 2017-08-07 12:15:38Z gronemeier
130! Modularize 1D model
131!
132! 2329 2017-08-03 14:24:56Z knoop
133! Removed temporary bugfix (r2327) as bug is properly resolved by this revision
134!
135! 2327 2017-08-02 07:40:57Z maronga
136! Temporary bugfix
137!
138! 2320 2017-07-21 12:47:43Z suehring
139! Modularize large-scale forcing and nudging
140!
141! 2292 2017-06-20 09:51:42Z schwenkel
142! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
143! includes two more prognostic equations for cloud drop concentration (nc) 
144! and cloud water content (qc).
145!
146! 2277 2017-06-12 10:47:51Z kanani
147! Removed unused variable sums_up_fraction_l
148!
149! 2270 2017-06-09 12:18:47Z maronga
150! dots_num must be increased when LSM and/or radiation is used
151!
152! 2259 2017-06-08 09:09:11Z gronemeier
153! Implemented synthetic turbulence generator
154!
155! 2252 2017-06-07 09:35:37Z knoop
156! rho_air now depending on surface_pressure even in Boussinesq mode
157!
158! 2233 2017-05-30 18:08:54Z suehring
159!
160! 2232 2017-05-30 17:47:52Z suehring
161! Adjustments to new topography and surface concept:
162!   - Modify passed parameters for disturb_field
163!   - Topography representation via flags
164!   - Remove unused arrays.
165!   - Move initialization of surface-related quantities to surface_mod
166!
167! 2172 2017-03-08 15:55:25Z knoop
168! Bugfix: moved parallel random generator initialization into its module
169!
170! 2118 2017-01-17 16:38:49Z raasch
171! OpenACC directives removed
172!
173! 2037 2016-10-26 11:15:40Z knoop
174! Anelastic approximation implemented
175!
176! 2031 2016-10-21 15:11:58Z knoop
177! renamed variable rho to rho_ocean
178!
179! 2011 2016-09-19 17:29:57Z kanani
180! Flag urban_surface is now defined in module control_parameters.
181!
182! 2007 2016-08-24 15:47:17Z kanani
183! Added support for urban surface model,
184! adjusted location_message in case of plant_canopy
185!
186! 2000 2016-08-20 18:09:15Z knoop
187! Forced header and separation lines into 80 columns
188!
189! 1992 2016-08-12 15:14:59Z suehring
190! Initializaton of scalarflux at model top
191! Bugfixes in initialization of surface and top salinity flux, top scalar and
192! humidity fluxes
193!
194! 1960 2016-07-12 16:34:24Z suehring
195! Separate humidity and passive scalar
196! Increase dimension for mean_inflow_profiles
197! Remove inadvertent write-statement
198! Bugfix, large-scale forcing is still not implemented for passive scalars
199!
200! 1957 2016-07-07 10:43:48Z suehring
201! flight module added
202!
203! 1920 2016-05-30 10:50:15Z suehring
204! Initialize us with very small number to avoid segmentation fault during
205! calculation of Obukhov length
206!
207! 1918 2016-05-27 14:35:57Z raasch
208! intermediate_timestep_count is set 0 instead 1 for first call of pres,
209! bugfix: initialization of local sum arrays are moved to the beginning of the
210!         routine because otherwise results from pres are overwritten
211!
212! 1914 2016-05-26 14:44:07Z witha
213! Added initialization of the wind turbine model
214!
215! 1878 2016-04-19 12:30:36Z hellstea
216! The zeroth element of weight_pres removed as unnecessary
217!
218! 1849 2016-04-08 11:33:18Z hoffmann
219! Adapted for modularization of microphysics.
220! precipitation_amount, precipitation_rate, prr moved to arrays_3d.
221! Initialization of nc_1d, nr_1d, pt_1d, qc_1d, qr_1d, q_1d moved to
222! microphysics_init.
223!
224! 1845 2016-04-08 08:29:13Z raasch
225! nzb_2d replaced by nzb_u|v_inner
226!
227! 1833 2016-04-07 14:23:03Z raasch
228! initialization of spectra quantities moved to spectra_mod
229!
230! 1831 2016-04-07 13:15:51Z hoffmann
231! turbulence renamed collision_turbulence
232!
233! 1826 2016-04-07 12:01:39Z maronga
234! Renamed radiation calls.
235! Renamed canopy model calls.
236!
237! 1822 2016-04-07 07:49:42Z hoffmann
238! icloud_scheme replaced by microphysics_*
239!
240! 1817 2016-04-06 15:44:20Z maronga
241! Renamed lsm calls.
242!
243! 1815 2016-04-06 13:49:59Z raasch
244! zero-settings for velocities inside topography re-activated (was deactivated
245! in r1762)
246!
247! 1788 2016-03-10 11:01:04Z maronga
248! Added z0q.
249! Syntax layout improved.
250!
251! 1783 2016-03-06 18:36:17Z raasch
252! netcdf module name changed + related changes
253!
254! 1764 2016-02-28 12:45:19Z raasch
255! bugfix: increase size of volume_flow_area_l and volume_flow_initial_l by 1
256!
257! 1762 2016-02-25 12:31:13Z hellstea
258! Introduction of nested domain feature
259!
260! 1738 2015-12-18 13:56:05Z raasch
261! calculate mean surface level height for each statistic region
262!
263! 1734 2015-12-02 12:17:12Z raasch
264! no initial disturbances in case that the disturbance energy limit has been
265! set zero
266!
267! 1707 2015-11-02 15:24:52Z maronga
268! Bugfix: transfer of Richardson number from 1D model to Obukhov length caused
269! devision by zero in neutral stratification
270!
271! 1691 2015-10-26 16:17:44Z maronga
272! Call to init_surface_layer added. rif is replaced by ol and zeta.
273!
274! 1682 2015-10-07 23:56:08Z knoop
275! Code annotations made doxygen readable
276!
277! 1615 2015-07-08 18:49:19Z suehring
278! Enable turbulent inflow for passive_scalar and humidity
279!
280! 1585 2015-04-30 07:05:52Z maronga
281! Initialization of radiation code is now done after LSM initializtion
282!
283! 1575 2015-03-27 09:56:27Z raasch
284! adjustments for psolver-queries
285!
286! 1551 2015-03-03 14:18:16Z maronga
287! Allocation of land surface arrays is now done in the subroutine lsm_init_arrays,
288! which is part of land_surface_model.
289!
290! 1507 2014-12-10 12:14:18Z suehring
291! Bugfix: set horizontal velocity components to zero inside topography
292!
293! 1496 2014-12-02 17:25:50Z maronga
294! Added initialization of the land surface and radiation schemes
295!
296! 1484 2014-10-21 10:53:05Z kanani
297! Changes due to new module structure of the plant canopy model:
298! canopy-related initialization (e.g. lad and canopy_heat_flux) moved to new
299! subroutine init_plant_canopy within the module plant_canopy_model_mod,
300! call of subroutine init_plant_canopy added.
301!
302! 1431 2014-07-15 14:47:17Z suehring
303! var_d added, in order to normalize spectra.
304!
305! 1429 2014-07-15 12:53:45Z knoop
306! Ensemble run capability added to parallel random number generator
307!
308! 1411 2014-05-16 18:01:51Z suehring
309! Initial horizontal velocity profiles were not set to zero at the first vertical
310! grid level in case of non-cyclic lateral boundary conditions.
311!
312! 1406 2014-05-16 13:47:01Z raasch
313! bugfix: setting of initial velocities at k=1 to zero not in case of a
314! no-slip boundary condition for uv
315!
316! 1402 2014-05-09 14:25:13Z raasch
317! location messages modified
318!
319! 1400 2014-05-09 14:03:54Z knoop
320! Parallel random number generator added
321!
322! 1384 2014-05-02 14:31:06Z raasch
323! location messages added
324!
325! 1361 2014-04-16 15:17:48Z hoffmann
326! tend_* removed
327! Bugfix: w_subs is not allocated anymore if it is already allocated
328!
329! 1359 2014-04-11 17:15:14Z hoffmann
330! module lpm_init_mod added to use statements, because lpm_init has become a
331! module
332!
333! 1353 2014-04-08 15:21:23Z heinze
334! REAL constants provided with KIND-attribute
335!
336! 1340 2014-03-25 19:45:13Z kanani
337! REAL constants defined as wp-kind
338!
339! 1322 2014-03-20 16:38:49Z raasch
340! REAL constants defined as wp-kind
341! module interfaces removed
342!
343! 1320 2014-03-20 08:40:49Z raasch
344! ONLY-attribute added to USE-statements,
345! kind-parameters added to all INTEGER and REAL declaration statements,
346! kinds are defined in new module kinds,
347! revision history before 2012 removed,
348! comment fields (!:) to be used for variable explanations added to
349! all variable declaration statements
350!
351! 1316 2014-03-17 07:44:59Z heinze
352! Bugfix: allocation of w_subs
353!
354! 1299 2014-03-06 13:15:21Z heinze
355! Allocate w_subs due to extension of large scale subsidence in combination
356! with large scale forcing data (LSF_DATA)
357!
358! 1241 2013-10-30 11:36:58Z heinze
359! Overwrite initial profiles in case of nudging
360! Inititialize shf and qsws in case of large_scale_forcing
361!
362! 1221 2013-09-10 08:59:13Z raasch
363! +rflags_s_inner in copyin statement, use copyin for most arrays instead of
364! copy
365!
366! 1212 2013-08-15 08:46:27Z raasch
367! array tri is allocated and included in data copy statement
368!
369! 1195 2013-07-01 12:27:57Z heinze
370! Bugfix: move allocation of ref_state to parin.f90 and read_var_list.f90
371!
372! 1179 2013-06-14 05:57:58Z raasch
373! allocate and set ref_state to be used in buoyancy terms
374!
375! 1171 2013-05-30 11:27:45Z raasch
376! diss array is allocated with full size if accelerator boards are used
377!
378! 1159 2013-05-21 11:58:22Z fricke
379! -bc_lr_dirneu, bc_lr_neudir, bc_ns_dirneu, bc_ns_neudir
380!
381! 1153 2013-05-10 14:33:08Z raasch
382! diss array is allocated with dummy elements even if it is not needed
383! (required by PGI 13.4 / CUDA 5.0)
384!
385! 1115 2013-03-26 18:16:16Z hoffmann
386! unused variables removed
387!
388! 1113 2013-03-10 02:48:14Z raasch
389! openACC directive modified
390!
391! 1111 2013-03-08 23:54:10Z raasch
392! openACC directives added for pres
393! array diss allocated only if required
394!
395! 1092 2013-02-02 11:24:22Z raasch
396! unused variables removed
397!
398! 1065 2012-11-22 17:42:36Z hoffmann
399! allocation of diss (dissipation rate) in case of turbulence = .TRUE. added
400!
401! 1053 2012-11-13 17:11:03Z hoffmann
402! allocation and initialisation of necessary data arrays for the two-moment
403! cloud physics scheme the two new prognostic equations (nr, qr):
404! +dr, lambda_r, mu_r, sed_*, xr, *s, *sws, *swst, *, *_p, t*_m, *_1, *_2, *_3,
405! +tend_*, prr
406!
407! 1036 2012-10-22 13:43:42Z raasch
408! code put under GPL (PALM 3.9)
409!
410! 1032 2012-10-21 13:03:21Z letzel
411! save memory by not allocating pt_2 in case of neutral = .T.
412!
413! 1025 2012-10-07 16:04:41Z letzel
414! bugfix: swap indices of mask for ghost boundaries
415!
416! 1015 2012-09-27 09:23:24Z raasch
417! mask is set to zero for ghost boundaries
418!
419! 1010 2012-09-20 07:59:54Z raasch
420! cpp switch __nopointer added for pointer free version
421!
422! 1003 2012-09-14 14:35:53Z raasch
423! nxra,nyna, nzta replaced ny nxr, nyn, nzt
424!
425! 1001 2012-09-13 14:08:46Z raasch
426! all actions concerning leapfrog scheme removed
427!
428! 996 2012-09-07 10:41:47Z raasch
429! little reformatting
430!
431! 978 2012-08-09 08:28:32Z fricke
432! outflow damping layer removed
433! roughness length for scalar quantites z0h added
434! damping zone for the potential temperatur in case of non-cyclic lateral
435! boundaries added
436! initialization of ptdf_x, ptdf_y
437! initialization of c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l
438!
439! 849 2012-03-15 10:35:09Z raasch
440! init_particles renamed lpm_init
441!
442! 825 2012-02-19 03:03:44Z raasch
443! wang_collision_kernel renamed wang_kernel
444!
445! Revision 1.1  1998/03/09 16:22:22  raasch
446! Initial revision
447!
448!
449! Description:
450! ------------
451!> Allocation of arrays and initialization of the 3D model via
452!> a) pre-run the 1D model
453!> or
454!> b) pre-set constant linear profiles
455!> or
456!> c) read values of a previous run
457!------------------------------------------------------------------------------!
458 SUBROUTINE init_3d_model
459 
460
461    USE advec_ws
462
463    USE arrays_3d
464
465    USE chemistry_model_mod,                                                   &
466        ONLY:  chem_emissions
467
468    USE cloud_parameters,                                                      &
469        ONLY:  cp, l_v, r_d
470
471    USE constants,                                                             &
472        ONLY:  pi
473   
474    USE control_parameters
475   
476    USE flight_mod,                                                            &
477        ONLY:  flight_init
478   
479    USE grid_variables,                                                        &
480        ONLY:  dx, dy, ddx2_mg, ddy2_mg
481
482    USE gust_mod,                                                              &
483        ONLY:  gust_init, gust_init_arrays, gust_module_enabled
484   
485    USE indices
486
487    USE lpm_init_mod,                                                          &
488        ONLY:  lpm_init
489   
490    USE kinds
491
492    USE land_surface_model_mod,                                                &
493        ONLY:  lsm_init, lsm_init_arrays
494 
495    USE lsf_nudging_mod,                                                       &
496        ONLY:  lsf_init, ls_forcing_surf, nudge_init
497
498    USE microphysics_mod,                                                      &
499        ONLY:  collision_turbulence, microphysics_init
500
501    USE model_1d_mod,                                                          &
502        ONLY:  e1d, init_1d_model, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d,  &
503               v1d, vsws1d
504
505    USE netcdf_interface,                                                      &
506        ONLY:  dots_max, dots_num, dots_unit, dots_label
507
508    USE netcdf_data_input_mod,                                                 &
509        ONLY:  init_3d, netcdf_data_input_interpolate, netcdf_data_input_init_3d
510   
511    USE particle_attributes,                                                   &
512        ONLY:  particle_advection, use_sgs_for_particles, wang_kernel
513   
514    USE pegrid
515   
516    USE plant_canopy_model_mod,                                                &
517        ONLY:  pcm_init
518
519    USE pmc_interface,                                                         &
520        ONLY:  nested_run
521
522    USE radiation_model_mod,                                                   &
523        ONLY:  average_radiation,                                              &
524               radiation_init, radiation, radiation_scheme,                    &
525               radiation_calc_svf, radiation_write_svf,                        &
526               radiation_interaction, radiation_interactions,                  &
527               radiation_interaction_init, radiation_read_svf,                 &
528               radiation_presimulate_solar_pos, radiation_interactions_on
529   
530    USE random_function_mod
531   
532    USE random_generator_parallel,                                             &
533        ONLY:  init_parallel_random_generator
534
535    USE read_restart_data_mod,                                                 &
536        ONLY:  rrd_read_parts_of_global, rrd_local                                     
537   
538    USE statistics,                                                            &
539        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask,        &
540               statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,  &
541               sums_l_l, sums_wsts_bc_l, ts_value,                             &
542               weight_pres, weight_substep
543
544    USE synthetic_turbulence_generator_mod,                                    &
545        ONLY:  stg_init, use_syn_turb_gen
546
547    USE surface_layer_fluxes_mod,                                              &
548        ONLY:  init_surface_layer_fluxes
549
550    USE surface_mod,                                                           &
551        ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
552                surf_usm_h, get_topography_top_index_ji, vertical_surfaces_exist
553   
554    USE transpose_indices
555
556    USE turbulence_closure_mod,                                                &
557        ONLY:  tcm_init_arrays, tcm_init
558
559    USE urban_surface_mod,                                                     &
560        ONLY:  usm_init_urban_surface, usm_allocate_surface
561
562    USE uv_exposure_model_mod,                                                 &
563        ONLY:  uvem_init, uvem_init_arrays
564
565    USE wind_turbine_model_mod,                                                &
566        ONLY:  wtm_init, wtm_init_arrays
567
568    IMPLICIT NONE
569
570    INTEGER(iwp) ::  i             !<
571    INTEGER(iwp) ::  ind_array(1)  !<
572    INTEGER(iwp) ::  j             !<
573    INTEGER(iwp) ::  k             !<
574    INTEGER(iwp) ::  k_surf        !< surface level index
575    INTEGER(iwp) ::  m             !< index of surface element in surface data type
576    INTEGER(iwp) ::  sr            !< index of statistic region
577
578    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !<
579
580    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !<
581    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !<
582
583    REAL(wp)     ::  t_surface !< air temperature at the surface
584
585    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
586
587    INTEGER(iwp) ::  l       !< loop variable
588    INTEGER(iwp) ::  nzt_l   !< index of top PE boundary for multigrid level
589    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
590    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
591
592    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !<
593    REAL(wp), DIMENSION(1:3) ::  volume_flow_initial_l  !<
594
595    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l    !<
596    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l    !<
597    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !<
598
599    INTEGER(iwp) ::  nz_u_shift   !<
600    INTEGER(iwp) ::  nz_v_shift   !<
601    INTEGER(iwp) ::  nz_w_shift   !<
602    INTEGER(iwp) ::  nz_s_shift   !<
603    INTEGER(iwp) ::  nz_u_shift_l !<
604    INTEGER(iwp) ::  nz_v_shift_l !<
605    INTEGER(iwp) ::  nz_w_shift_l !<
606    INTEGER(iwp) ::  nz_s_shift_l !<
607
608    CALL location_message( 'allocating arrays', .FALSE. )
609!
610!-- Allocate arrays
611    ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
612              mean_surface_level_height_l(0:statistic_regions),                &
613              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
614              ngp_3d(0:statistic_regions),                                     &
615              ngp_3d_inner(0:statistic_regions),                               &
616              ngp_3d_inner_l(0:statistic_regions),                             &
617              ngp_3d_inner_tmp(0:statistic_regions),                           &
618              sums_divnew_l(0:statistic_regions),                              &
619              sums_divold_l(0:statistic_regions) )
620    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
621    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
622              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
623              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
624              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
625              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
626              sums(nzb:nzt+1,pr_palm+max_pr_user),                             &
627              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),      &
628              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
629              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                   &
630              ts_value(dots_max,0:statistic_regions) )
631    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
632
633    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
634              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
635              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
636
637#if defined( __nopointer )
638    ALLOCATE( pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
639              pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
640              u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
641              u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
642              v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
643              v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
644              w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
645              w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
646              tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
647              tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
648              tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
649              tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
650#else
651    ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
652              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
653              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
654              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
655              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
656              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
657              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
658              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
659              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
660              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
661              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
662    IF (  .NOT.  neutral )  THEN
663       ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
664    ENDIF
665#endif
666
667!
668!-- Following array is required for perturbation pressure within the iterative
669!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
670!-- the weighted average of the substeps and cannot be used in the Poisson
671!-- solver.
672    IF ( psolver == 'sor' )  THEN
673       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
674    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
675!
676!--    For performance reasons, multigrid is using one ghost layer only
677       ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
678    ENDIF
679
680!
681!-- Array for storing constant coeffficients of the tridiagonal solver
682    IF ( psolver == 'poisfft' )  THEN
683       ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
684       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
685    ENDIF
686
687    IF ( humidity )  THEN
688!
689!--    3D-humidity
690#if defined( __nopointer )
691       ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
692                 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
693                 tq_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
694#else
695       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
696                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
697                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
698#endif
699
700!
701!--    3D-arrays needed for humidity
702       IF ( humidity )  THEN
703#if defined( __nopointer )
704          ALLOCATE( vpt(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
705#else
706          ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
707#endif
708
709          IF ( cloud_physics )  THEN
710!
711!--          Liquid water content
712#if defined( __nopointer )
713             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
714#else
715             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
716#endif
717
718!
719!--          3D-cloud water content
720             IF ( .NOT. microphysics_morrison )  THEN
721#if defined( __nopointer )
722                ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
723#else
724                ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
725#endif
726             ENDIF
727!
728!--          Precipitation amount and rate (only needed if output is switched)
729             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) )
730
731!
732!--          3d-precipitation rate
733             ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
734
735             IF ( microphysics_morrison )  THEN
736!
737!--             3D-cloud drop water content, cloud drop concentration arrays
738#if defined( __nopointer )
739                ALLOCATE( nc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
740                          nc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
741                          qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
742                          qc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
743                          tnc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             & 
744                          tqc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
745#else
746                ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
747                          nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
748                          nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
749                          qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
750                          qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
751                          qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
752#endif
753             ENDIF
754
755             IF ( microphysics_seifert )  THEN
756!
757!--             3D-rain water content, rain drop concentration arrays
758#if defined( __nopointer )
759                ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
760                          nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
761                          qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
762                          qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
763                          tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
764                          tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
765#else
766                ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
767                          nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
768                          nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
769                          qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
770                          qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
771                          qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
772#endif
773             ENDIF
774
775          ENDIF
776
777          IF ( cloud_droplets )  THEN
778!
779!--          Liquid water content, change in liquid water content
780#if defined( __nopointer )
781             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
782                        ql_c(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
783#else
784             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
785                        ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
786#endif
787!
788!--          Real volume of particles (with weighting), volume of particles
789             ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
790                        ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
791          ENDIF
792
793       ENDIF
794
795    ENDIF
796   
797   
798    IF ( passive_scalar )  THEN
799
800!
801!--    3D scalar arrays
802#if defined( __nopointer )
803       ALLOCATE( s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
804                 s_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
805                 ts_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
806#else
807       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
808                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
809                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
810#endif
811    ENDIF
812
813    IF ( ocean )  THEN
814#if defined( __nopointer )
815       ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
816                 rho_ocean(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
817                 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
818                 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
819                 tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
820#else
821       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
822                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
823                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
824                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
825                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
826       prho => prho_1
827       rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
828                      ! density to be apointer
829#endif
830    ENDIF
831
832!
833!-- Allocation of anelastic and Boussinesq approximation specific arrays
834    ALLOCATE( p_hydrostatic(nzb:nzt+1) )
835    ALLOCATE( rho_air(nzb:nzt+1) )
836    ALLOCATE( rho_air_zw(nzb:nzt+1) )
837    ALLOCATE( drho_air(nzb:nzt+1) )
838    ALLOCATE( drho_air_zw(nzb:nzt+1) )
839
840!
841!-- Density profile calculation for anelastic approximation
842    t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp )
843    IF ( TRIM( approximation ) == 'anelastic' ) THEN
844       DO  k = nzb, nzt+1
845          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
846                                ( 1 - ( g * zu(k) ) / ( cp * t_surface )       &
847                                )**( cp / r_d )
848          rho_air(k)          = ( p_hydrostatic(k) *                           &
849                                  ( 100000.0_wp / p_hydrostatic(k)             &
850                                  )**( r_d / cp )                              &
851                                ) / ( r_d * pt_init(k) )
852       ENDDO
853       DO  k = nzb, nzt
854          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
855       ENDDO
856       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
857                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
858    ELSE
859       DO  k = nzb, nzt+1
860          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
861                                ( 1 - ( g * zu(nzb) ) / ( cp * t_surface )       &
862                                )**( cp / r_d )
863          rho_air(k)          = ( p_hydrostatic(k) *                           &
864                                  ( 100000.0_wp / p_hydrostatic(k)             &
865                                  )**( r_d / cp )                              &
866                                ) / ( r_d * pt_init(nzb) )
867       ENDDO
868       DO  k = nzb, nzt
869          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
870       ENDDO
871       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
872                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
873    ENDIF
874!
875!-- compute the inverse density array in order to avoid expencive divisions
876    drho_air    = 1.0_wp / rho_air
877    drho_air_zw = 1.0_wp / rho_air_zw
878
879!
880!-- Allocation of flux conversion arrays
881    ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
882    ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
883    ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
884    ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
885    ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
886    ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
887
888!
889!-- calculate flux conversion factors according to approximation and in-/output mode
890    DO  k = nzb, nzt+1
891
892        IF ( TRIM( flux_input_mode ) == 'kinematic' )  THEN
893            heatflux_input_conversion(k)      = rho_air_zw(k)
894            waterflux_input_conversion(k)     = rho_air_zw(k)
895            momentumflux_input_conversion(k)  = rho_air_zw(k)
896        ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
897            heatflux_input_conversion(k)      = 1.0_wp / cp
898            waterflux_input_conversion(k)     = 1.0_wp / l_v
899            momentumflux_input_conversion(k)  = 1.0_wp
900        ENDIF
901
902        IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
903            heatflux_output_conversion(k)     = drho_air_zw(k)
904            waterflux_output_conversion(k)    = drho_air_zw(k)
905            momentumflux_output_conversion(k) = drho_air_zw(k)
906        ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
907            heatflux_output_conversion(k)     = cp
908            waterflux_output_conversion(k)    = l_v
909            momentumflux_output_conversion(k) = 1.0_wp
910        ENDIF
911
912        IF ( .NOT. humidity ) THEN
913            waterflux_input_conversion(k)  = 1.0_wp
914            waterflux_output_conversion(k) = 1.0_wp
915        ENDIF
916
917    ENDDO
918
919!
920!-- In case of multigrid method, compute grid lengths and grid factors for the
921!-- grid levels with respective density on each grid
922    IF ( psolver(1:9) == 'multigrid' )  THEN
923
924       ALLOCATE( ddx2_mg(maximum_grid_level) )
925       ALLOCATE( ddy2_mg(maximum_grid_level) )
926       ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
927       ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
928       ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
929       ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
930       ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
931       ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
932       ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
933
934       dzu_mg(:,maximum_grid_level) = dzu
935       rho_air_mg(:,maximum_grid_level) = rho_air
936!       
937!--    Next line to ensure an equally spaced grid.
938       dzu_mg(1,maximum_grid_level) = dzu(2)
939       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) +                     &
940                                             (rho_air(nzb) - rho_air(nzb+1))
941
942       dzw_mg(:,maximum_grid_level) = dzw
943       rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
944       nzt_l = nzt
945       DO  l = maximum_grid_level-1, 1, -1
946           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
947           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
948           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
949           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))
950           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
951           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
952           nzt_l = nzt_l / 2
953           DO  k = 2, nzt_l+1
954              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
955              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
956              rho_air_mg(k,l)    = rho_air_mg(2*k-1,l+1)
957              rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
958           ENDDO
959       ENDDO
960
961       nzt_l = nzt
962       dx_l  = dx
963       dy_l  = dy
964       DO  l = maximum_grid_level, 1, -1
965          ddx2_mg(l) = 1.0_wp / dx_l**2
966          ddy2_mg(l) = 1.0_wp / dy_l**2
967          DO  k = nzb+1, nzt_l
968             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
969             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
970             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
971                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
972          ENDDO
973          nzt_l = nzt_l / 2
974          dx_l  = dx_l * 2.0_wp
975          dy_l  = dy_l * 2.0_wp
976       ENDDO
977
978    ENDIF
979
980!
981!-- 1D-array for large scale subsidence velocity
982    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
983       ALLOCATE ( w_subs(nzb:nzt+1) )
984       w_subs = 0.0_wp
985    ENDIF
986
987!
988!-- Arrays to store velocity data from t-dt and the phase speeds which
989!-- are needed for radiation boundary conditions
990    IF ( outflow_l )  THEN
991       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
992                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
993                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
994    ENDIF
995    IF ( outflow_r )  THEN
996       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
997                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
998                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
999    ENDIF
1000    IF ( outflow_l  .OR.  outflow_r )  THEN
1001       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
1002                 c_w(nzb:nzt+1,nysg:nyng) )
1003    ENDIF
1004    IF ( outflow_s )  THEN
1005       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
1006                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
1007                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
1008    ENDIF
1009    IF ( outflow_n )  THEN
1010       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
1011                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
1012                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
1013    ENDIF
1014    IF ( outflow_s  .OR.  outflow_n )  THEN
1015       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
1016                 c_w(nzb:nzt+1,nxlg:nxrg) )
1017    ENDIF
1018    IF ( outflow_l  .OR.  outflow_r  .OR.  outflow_s  .OR.  outflow_n )  THEN
1019       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
1020       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
1021    ENDIF
1022
1023
1024#if ! defined( __nopointer )
1025!
1026!-- Initial assignment of the pointers
1027    IF ( .NOT. neutral )  THEN
1028       pt => pt_1;  pt_p => pt_2;  tpt_m => pt_3
1029    ELSE
1030       pt => pt_1;  pt_p => pt_1;  tpt_m => pt_3
1031    ENDIF
1032    u  => u_1;   u_p  => u_2;   tu_m  => u_3
1033    v  => v_1;   v_p  => v_2;   tv_m  => v_3
1034    w  => w_1;   w_p  => w_2;   tw_m  => w_3
1035
1036    IF ( humidity )  THEN
1037       q => q_1;  q_p => q_2;  tq_m => q_3
1038       IF ( humidity )  THEN
1039          vpt  => vpt_1   
1040          IF ( cloud_physics )  THEN
1041             ql => ql_1
1042             IF ( .NOT. microphysics_morrison )  THEN
1043                qc => qc_1
1044             ENDIF
1045             IF ( microphysics_morrison )  THEN
1046                qc => qc_1;  qc_p  => qc_2;  tqc_m  => qc_3
1047                nc => nc_1;  nc_p  => nc_2;  tnc_m  => nc_3
1048             ENDIF
1049             IF ( microphysics_seifert )  THEN
1050                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
1051                nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
1052             ENDIF
1053          ENDIF
1054       ENDIF
1055       IF ( cloud_droplets )  THEN
1056          ql   => ql_1
1057          ql_c => ql_2
1058       ENDIF
1059    ENDIF
1060   
1061    IF ( passive_scalar )  THEN
1062       s => s_1;  s_p => s_2;  ts_m => s_3
1063    ENDIF   
1064
1065    IF ( ocean )  THEN
1066       sa => sa_1;  sa_p => sa_2;  tsa_m => sa_3
1067    ENDIF
1068#endif
1069!
1070!-- Initialize arrays for turbulence closure
1071    CALL tcm_init_arrays
1072!
1073!-- Initialize surface arrays
1074    CALL init_surface_arrays
1075!
1076!-- Allocate land surface model arrays
1077    IF ( land_surface )  THEN
1078       CALL lsm_init_arrays
1079    ENDIF
1080
1081!
1082!-- Allocate wind turbine model arrays
1083    IF ( wind_turbine )  THEN
1084       CALL wtm_init_arrays
1085    ENDIF
1086!
1087!-- Allocate gust module arrays
1088    IF ( gust_module_enabled )  THEN
1089       CALL gust_init_arrays
1090    ENDIF
1091
1092!
1093!-- Initialize virtual flight measurements
1094    IF ( virtual_flight )  THEN
1095       CALL flight_init
1096    ENDIF
1097
1098!
1099!-- Read uv exposure input data
1100    IF ( uv_exposure )  THEN
1101       CALL uvem_init
1102    ENDIF
1103!
1104!-- Allocate uv exposure arrays
1105    IF ( uv_exposure )  THEN
1106       CALL uvem_init_arrays
1107    ENDIF
1108
1109!
1110!-- Initialize nudging if required
1111    IF ( nudging )  THEN
1112       CALL nudge_init
1113    ENDIF
1114
1115!
1116!-- Initialize reading of large scale forcing from external file - if required
1117    IF ( large_scale_forcing  .OR.  forcing )  THEN
1118       CALL lsf_init
1119    ENDIF
1120
1121!
1122!-- Allocate arrays containing the RK coefficient for calculation of
1123!-- perturbation pressure and turbulent fluxes. At this point values are
1124!-- set for pressure calculation during initialization (where no timestep
1125!-- is done). Further below the values needed within the timestep scheme
1126!-- will be set.
1127    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
1128              weight_pres(1:intermediate_timestep_count_max) )
1129    weight_substep = 1.0_wp
1130    weight_pres    = 1.0_wp
1131    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
1132       
1133    CALL location_message( 'finished', .TRUE. )
1134
1135!
1136!-- Initialize local summation arrays for routine flow_statistics.
1137!-- This is necessary because they may not yet have been initialized when they
1138!-- are called from flow_statistics (or - depending on the chosen model run -
1139!-- are never initialized)
1140    sums_divnew_l      = 0.0_wp
1141    sums_divold_l      = 0.0_wp
1142    sums_l_l           = 0.0_wp
1143    sums_wsts_bc_l     = 0.0_wp
1144
1145
1146
1147!
1148!-- Initialize model variables
1149    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1150         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1151!
1152!--    Initialization with provided input data derived from larger-scale model
1153       IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
1154          CALL location_message( 'initializing with INIFOR', .FALSE. )
1155!
1156!--       Read initial 1D profiles from NetCDF file if available.
1157!--       At the moment, only u, v, w, pt and q are provided.
1158          CALL netcdf_data_input_init_3d
1159!
1160!--       Please note, at the moment INIFOR assumes only an equidistant vertical
1161!--       grid. In case of vertical grid stretching, input of inital data
1162!--       need to be inter- and/or extrapolated.
1163!--       Therefore, check if zu grid on file is identical to numeric zw grid.
1164!--       Please note 
1165          IF ( ANY( zu(1:nzt+1) /= init_3d%zu_atmos(1:init_3d%nzu) ) )  THEN
1166
1167             CALL netcdf_data_input_interpolate( init_3d%u_init(nzb+1:nzt+1),  &
1168                                                 zu(nzb+1:nzt+1),              &
1169                                                 init_3d%zu_atmos )
1170             CALL netcdf_data_input_interpolate( init_3d%v_init(nzb+1:nzt+1),  &
1171                                                 zu(nzb+1:nzt+1),              &
1172                                                 init_3d%zu_atmos )
1173!              CALL netcdf_data_input_interpolate( init_3d%w_init(nzb+1:nzt),    &
1174!                                                  zw(nzb+1:nzt),                &
1175!                                                  init_3d%zw_atmos )
1176             IF ( .NOT. neutral )                                              &
1177                CALL netcdf_data_input_interpolate(                            &
1178                                             init_3d%pt_init(nzb+1:nzt+1),     &
1179                                             zu(nzb+1:nzt+1),                  &
1180                                             init_3d%zu_atmos )
1181             IF ( humidity )                                                   &
1182                CALL netcdf_data_input_interpolate(                            &
1183                                             init_3d%q_init(nzb+1:nzt+1),      &
1184                                             zu(nzb+1:nzt+1),                  &
1185                                             init_3d%zu_atmos )
1186          ENDIF
1187
1188          u_init = init_3d%u_init
1189          v_init = init_3d%v_init   
1190          IF( .NOT. neutral )  pt_init = init_3d%pt_init
1191          IF( humidity      )  q_init  = init_3d%q_init
1192
1193!
1194!--       Please note, Inifor provides data from nzb+1 to nzt+1.
1195!--       Initialize pt and q with Neumann condition at nzb.
1196          IF( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
1197          IF( humidity      )  q_init(nzb)  = q_init(nzb+1)
1198          DO  i = nxlg, nxrg
1199             DO  j = nysg, nyng
1200                u(:,j,i) = u_init(:)
1201                v(:,j,i) = v_init(:)
1202                IF( .NOT. neutral )  pt(:,j,i) = pt_init(:)
1203                IF( humidity      )  q(:,j,i)  = q_init(:)
1204             ENDDO
1205          ENDDO
1206!
1207!--       MS: What about the geostrophic wind profiles? Actually these
1208!--           are not identical to the initial wind profiles in this case.
1209!--           This need to be further revised.
1210          IF ( init_3d%from_file_ug )  THEN
1211             ug(:) = init_3d%ug_init(:)
1212          ENDIF
1213          IF ( init_3d%from_file_vg )  THEN
1214             vg(:) = init_3d%vg_init(:)
1215          ENDIF
1216
1217          ug(nzt+1) = ug(nzt)
1218          vg(nzt+1) = vg(nzt)
1219
1220!
1221!--       Set inital w to 0
1222          w = 0.0_wp
1223!
1224!--       Initialize the remaining quantities
1225          IF ( humidity )  THEN
1226             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1227                DO  i = nxlg, nxrg
1228                   DO  j = nysg, nyng
1229                      qc(:,j,i) = 0.0_wp
1230                      nc(:,j,i) = 0.0_wp
1231                   ENDDO
1232                ENDDO
1233             ENDIF
1234
1235             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1236                DO  i = nxlg, nxrg
1237                   DO  j = nysg, nyng
1238                      qr(:,j,i) = 0.0_wp
1239                      nr(:,j,i) = 0.0_wp
1240                   ENDDO
1241                ENDDO
1242             ENDIF
1243
1244          ENDIF
1245
1246          IF ( passive_scalar )  THEN
1247             DO  i = nxlg, nxrg
1248                DO  j = nysg, nyng
1249                   s(:,j,i) = s_init
1250                ENDDO
1251             ENDDO
1252          ENDIF
1253
1254          IF ( ocean )  THEN
1255             DO  i = nxlg, nxrg
1256                DO  j = nysg, nyng
1257                   sa(:,j,i) = sa_init
1258                ENDDO
1259             ENDDO
1260          ENDIF
1261
1262!
1263!--       Set velocity components at non-atmospheric / oceanic grid points to
1264!--       zero.
1265          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1266          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1267          w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
1268!
1269!--       Initialize surface variables, e.g. friction velocity, momentum
1270!--       fluxes, etc.
1271          CALL init_surfaces
1272!
1273!--       Initialize turbulence generator
1274          IF( use_syn_turb_gen )  CALL stg_init
1275
1276          CALL location_message( 'finished', .TRUE. )
1277!
1278!--    Initialization via computed 1D-model profiles
1279       ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1280
1281          CALL location_message( 'initializing with 1D model profiles', .FALSE. )
1282!
1283!--       Use solutions of the 1D model as initial profiles,
1284!--       start 1D model
1285          CALL init_1d_model
1286!
1287!--       Transfer initial profiles to the arrays of the 3D model
1288          DO  i = nxlg, nxrg
1289             DO  j = nysg, nyng
1290                pt(:,j,i) = pt_init
1291                u(:,j,i)  = u1d
1292                v(:,j,i)  = v1d
1293             ENDDO
1294          ENDDO
1295
1296          IF ( humidity )  THEN
1297             DO  i = nxlg, nxrg
1298                DO  j = nysg, nyng
1299                   q(:,j,i) = q_init
1300                ENDDO
1301             ENDDO
1302             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1303                DO  i = nxlg, nxrg
1304                   DO  j = nysg, nyng
1305                      qc(:,j,i) = 0.0_wp
1306                      nc(:,j,i) = 0.0_wp
1307                   ENDDO
1308                ENDDO
1309             ENDIF
1310             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1311                DO  i = nxlg, nxrg
1312                   DO  j = nysg, nyng
1313                      qr(:,j,i) = 0.0_wp
1314                      nr(:,j,i) = 0.0_wp
1315                   ENDDO
1316                ENDDO
1317             ENDIF
1318          ENDIF
1319
1320          IF ( passive_scalar )  THEN
1321             DO  i = nxlg, nxrg
1322                DO  j = nysg, nyng
1323                   s(:,j,i) = s_init
1324                ENDDO
1325             ENDDO   
1326          ENDIF
1327!
1328!--          Store initial profiles for output purposes etc.
1329          IF ( .NOT. constant_diffusion )  THEN
1330             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
1331          ENDIF
1332!
1333!--       Set velocities back to zero
1334          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1335          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )         
1336!
1337!--       WARNING: The extra boundary conditions set after running the
1338!--       -------  1D model impose an error on the divergence one layer
1339!--                below the topography; need to correct later
1340!--       ATTENTION: Provisional correction for Piacsek & Williams
1341!--       ---------  advection scheme: keep u and v zero one layer below
1342!--                  the topography.
1343          IF ( ibc_uv_b == 1 )  THEN
1344!
1345!--          Neumann condition
1346             DO  i = nxl-1, nxr+1
1347                DO  j = nys-1, nyn+1
1348                   u(nzb,j,i) = u(nzb+1,j,i)
1349                   v(nzb,j,i) = v(nzb+1,j,i)
1350                ENDDO
1351             ENDDO
1352
1353          ENDIF
1354!
1355!--       Initialize surface variables, e.g. friction velocity, momentum
1356!--       fluxes, etc.
1357          CALL init_surfaces
1358
1359          CALL location_message( 'finished', .TRUE. )
1360
1361       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
1362       THEN
1363
1364          CALL location_message( 'initializing with constant profiles', .FALSE. )
1365!
1366!--       Overwrite initial profiles in case of synthetic turbulence generator
1367          IF( use_syn_turb_gen )  CALL stg_init
1368
1369!
1370!--       Use constructed initial profiles (velocity constant with height,
1371!--       temperature profile with constant gradient)
1372          DO  i = nxlg, nxrg
1373             DO  j = nysg, nyng
1374                pt(:,j,i) = pt_init
1375                u(:,j,i)  = u_init
1376                v(:,j,i)  = v_init
1377             ENDDO
1378          ENDDO
1379!
1380!--       Mask topography
1381          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1382          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1383!
1384!--       Set initial horizontal velocities at the lowest computational grid
1385!--       levels to zero in order to avoid too small time steps caused by the
1386!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
1387!--       in the limiting formula!).
1388!--       Please note, in case land- or urban-surface model is used and a
1389!--       spinup is applied, masking the lowest computational level is not
1390!--       possible as MOST as well as energy-balance parametrizations will not
1391!--       work with zero wind velocity.
1392          IF ( ibc_uv_b /= 1  .AND.  .NOT.  spinup )  THEN
1393             DO  i = nxlg, nxrg
1394                DO  j = nysg, nyng
1395                   DO  k = nzb, nzt
1396                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp,                      &
1397                                        BTEST( wall_flags_0(k,j,i), 20 ) )
1398                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp,                      &
1399                                        BTEST( wall_flags_0(k,j,i), 21 ) )
1400                   ENDDO
1401                ENDDO
1402             ENDDO
1403          ENDIF
1404
1405          IF ( humidity )  THEN
1406             DO  i = nxlg, nxrg
1407                DO  j = nysg, nyng
1408                   q(:,j,i) = q_init
1409                ENDDO
1410             ENDDO
1411             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1412                DO  i = nxlg, nxrg
1413                   DO  j = nysg, nyng
1414                      qc(:,j,i) = 0.0_wp
1415                      nc(:,j,i) = 0.0_wp
1416                   ENDDO
1417                ENDDO
1418             ENDIF
1419
1420             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1421                DO  i = nxlg, nxrg
1422                   DO  j = nysg, nyng
1423                      qr(:,j,i) = 0.0_wp
1424                      nr(:,j,i) = 0.0_wp
1425                   ENDDO
1426                ENDDO
1427             ENDIF
1428
1429          ENDIF
1430         
1431          IF ( passive_scalar )  THEN
1432             DO  i = nxlg, nxrg
1433                DO  j = nysg, nyng
1434                   s(:,j,i) = s_init
1435                ENDDO
1436             ENDDO
1437          ENDIF
1438
1439          IF ( ocean )  THEN
1440             DO  i = nxlg, nxrg
1441                DO  j = nysg, nyng
1442                   sa(:,j,i) = sa_init
1443                ENDDO
1444             ENDDO
1445          ENDIF
1446!
1447!--       Compute initial temperature field and other constants used in case
1448!--       of a sloping surface
1449          IF ( sloping_surface )  CALL init_slope
1450!
1451!--       Initialize surface variables, e.g. friction velocity, momentum
1452!--       fluxes, etc.
1453          CALL init_surfaces
1454
1455          CALL location_message( 'finished', .TRUE. )
1456
1457       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
1458       THEN
1459
1460          CALL location_message( 'initializing by user', .FALSE. )
1461!
1462!--       Pre-initialize surface variables, i.e. setting start- and end-indices
1463!--       at each (j,i)-location. Please note, this does not supersede
1464!--       user-defined initialization of surface quantities.
1465          CALL init_surfaces
1466!
1467!--       Initialization will completely be done by the user
1468          CALL user_init_3d_model
1469
1470          CALL location_message( 'finished', .TRUE. )
1471
1472       ENDIF
1473
1474       CALL location_message( 'initializing statistics, boundary conditions, etc.', &
1475                              .FALSE. )
1476
1477!
1478!--    Bottom boundary
1479       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
1480          u(nzb,:,:) = 0.0_wp
1481          v(nzb,:,:) = 0.0_wp
1482       ENDIF
1483
1484!
1485!--    Apply channel flow boundary condition
1486       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
1487          u(nzt+1,:,:) = 0.0_wp
1488          v(nzt+1,:,:) = 0.0_wp
1489       ENDIF
1490
1491!
1492!--    Calculate virtual potential temperature
1493       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
1494
1495!
1496!--    Store initial profiles for output purposes etc.. Please note, in case of
1497!--    initialization of u, v, w, pt, and q via output data derived from larger
1498!--    scale models, data will not be horizontally homogeneous. Actually, a mean
1499!--    profile should be calculated before.   
1500       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
1501       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
1502       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
1503          hom(nzb,1,5,:) = 0.0_wp
1504          hom(nzb,1,6,:) = 0.0_wp
1505       ENDIF
1506       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
1507
1508
1509!
1510!--    Store initial salinity profile
1511       IF ( ocean )  THEN
1512          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
1513       ENDIF
1514
1515       IF ( humidity )  THEN
1516!
1517!--       Store initial profile of total water content, virtual potential
1518!--       temperature
1519          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
1520          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
1521!
1522!--       Store initial profile of specific humidity and potential
1523!--       temperature
1524          IF ( cloud_physics  .OR.  cloud_droplets ) THEN
1525             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
1526             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
1527          ENDIF
1528       ENDIF
1529
1530!
1531!--    Store initial scalar profile
1532       IF ( passive_scalar )  THEN
1533          hom(:,1,121,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
1534       ENDIF
1535
1536!
1537!--    Initialize the random number generators (from numerical recipes)
1538       CALL random_function_ini
1539       
1540       IF ( random_generator == 'random-parallel' )  THEN
1541          CALL init_parallel_random_generator(nx, ny, nys, nyn, nxl, nxr)
1542       ENDIF
1543!
1544!--    Set the reference state to be used in the buoyancy terms (for ocean runs
1545!--    the reference state will be set (overwritten) in init_ocean)
1546       IF ( use_single_reference_value )  THEN
1547          IF (  .NOT.  humidity )  THEN
1548             ref_state(:) = pt_reference
1549          ELSE
1550             ref_state(:) = vpt_reference
1551          ENDIF
1552       ELSE
1553          IF (  .NOT.  humidity )  THEN
1554             ref_state(:) = pt_init(:)
1555          ELSE
1556             ref_state(:) = vpt(:,nys,nxl)
1557          ENDIF
1558       ENDIF
1559
1560!
1561!--    For the moment, vertical velocity is zero
1562       w = 0.0_wp
1563
1564!
1565!--    Initialize array sums (must be defined in first call of pres)
1566       sums = 0.0_wp
1567
1568!
1569!--    In case of iterative solvers, p must get an initial value
1570       IF ( psolver(1:9) == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
1571
1572!
1573!--    Treating cloud physics, liquid water content and precipitation amount
1574!--    are zero at beginning of the simulation
1575       IF ( cloud_physics )  THEN
1576          ql = 0.0_wp
1577          qc = 0.0_wp
1578
1579          precipitation_amount = 0.0_wp
1580       ENDIF
1581!
1582!--    Impose vortex with vertical axis on the initial velocity profile
1583       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
1584          CALL init_rankine
1585       ENDIF
1586
1587!
1588!--    Impose temperature anomaly (advection test only)
1589       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
1590          CALL init_pt_anomaly
1591       ENDIF
1592
1593!
1594!--    If required, change the surface temperature at the start of the 3D run
1595       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
1596          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
1597       ENDIF
1598
1599!
1600!--    If required, change the surface humidity/scalar at the start of the 3D
1601!--    run
1602       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
1603          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
1604         
1605       IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
1606          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
1607       
1608
1609!
1610!--    Initialize old and new time levels.
1611       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1612       pt_p = pt; u_p = u; v_p = v; w_p = w
1613
1614       IF ( humidity  )  THEN
1615          tq_m = 0.0_wp
1616          q_p = q
1617          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1618             tqc_m = 0.0_wp
1619             qc_p  = qc
1620             tnc_m = 0.0_wp
1621             nc_p  = nc
1622          ENDIF
1623          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1624             tqr_m = 0.0_wp
1625             qr_p  = qr
1626             tnr_m = 0.0_wp
1627             nr_p  = nr
1628          ENDIF
1629       ENDIF
1630       
1631       IF ( passive_scalar )  THEN
1632          ts_m = 0.0_wp
1633          s_p  = s
1634       ENDIF       
1635
1636       IF ( ocean )  THEN
1637          tsa_m = 0.0_wp
1638          sa_p  = sa
1639       ENDIF
1640       
1641       CALL location_message( 'finished', .TRUE. )
1642
1643    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1644             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
1645    THEN
1646
1647       CALL location_message( 'initializing in case of restart / cyclic_fill', &
1648                              .FALSE. )
1649!
1650!--    Initialize surface elements and its attributes, e.g. heat- and
1651!--    momentumfluxes, roughness, scaling parameters. As number of surface
1652!--    elements might be different between runs, e.g. in case of cyclic fill,
1653!--    and not all surface elements are read, surface elements need to be
1654!--    initialized before.     
1655       CALL init_surfaces
1656!
1657!--    When reading data for cyclic fill of 3D prerun data files, read
1658!--    some of the global variables from the restart file which are required
1659!--    for initializing the inflow
1660       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1661
1662          DO  i = 0, io_blocks-1
1663             IF ( i == io_group )  THEN
1664                CALL rrd_read_parts_of_global
1665             ENDIF
1666#if defined( __parallel )
1667             CALL MPI_BARRIER( comm2d, ierr )
1668#endif
1669          ENDDO
1670
1671       ENDIF
1672
1673!
1674!--    Read processor specific binary data from restart file
1675       DO  i = 0, io_blocks-1
1676          IF ( i == io_group )  THEN
1677             CALL rrd_local
1678          ENDIF
1679#if defined( __parallel )
1680          CALL MPI_BARRIER( comm2d, ierr )
1681#endif
1682       ENDDO
1683
1684!
1685!--    In case of complex terrain and cyclic fill method as initialization,
1686!--    shift initial data in the vertical direction for each point in the
1687!--    x-y-plane depending on local surface height
1688       IF ( complex_terrain  .AND.                                             &
1689            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1690          DO  i = nxlg, nxrg
1691             DO  j = nysg, nyng
1692                nz_u_shift = get_topography_top_index_ji( j, i, 'u' )
1693                nz_v_shift = get_topography_top_index_ji( j, i, 'v' )
1694                nz_w_shift = get_topography_top_index_ji( j, i, 'w' )
1695                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
1696
1697                u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)               
1698
1699                v(nz_v_shift:nzt+1,j,i)  = v(0:nzt+1-nz_v_shift,j,i)
1700
1701                w(nz_w_shift:nzt+1,j,i)  = w(0:nzt+1-nz_w_shift,j,i)
1702
1703                p(nz_s_shift:nzt+1,j,i)  =  p(0:nzt+1-nz_s_shift,j,i)
1704                pt(nz_s_shift:nzt+1,j,i) = pt(0:nzt+1-nz_s_shift,j,i)
1705             ENDDO
1706          ENDDO
1707       ENDIF
1708
1709!
1710!--    Initialization of the turbulence recycling method
1711       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
1712            turbulent_inflow )  THEN
1713!
1714!--       First store the profiles to be used at the inflow.
1715!--       These profiles are the (temporally) and horizontally averaged vertical
1716!--       profiles from the prerun. Alternatively, prescribed profiles
1717!--       for u,v-components can be used.
1718          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,7) )
1719
1720          IF ( use_prescribed_profile_data )  THEN
1721             mean_inflow_profiles(:,1) = u_init            ! u
1722             mean_inflow_profiles(:,2) = v_init            ! v
1723          ELSE
1724             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
1725             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
1726          ENDIF
1727          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
1728          IF ( humidity )                                                      &
1729             mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
1730          IF ( passive_scalar )                                                &
1731             mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
1732!
1733!--       In case of complex terrain, determine vertical displacement at inflow
1734!--       boundary and adjust mean inflow profiles
1735          IF ( complex_terrain )  THEN
1736             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
1737                nz_u_shift_l = get_topography_top_index_ji( 0, 0, 'u' )
1738                nz_v_shift_l = get_topography_top_index_ji( 0, 0, 'v' )
1739                nz_w_shift_l = get_topography_top_index_ji( 0, 0, 'w' )
1740                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
1741             ELSE
1742                nz_u_shift_l = 0
1743                nz_v_shift_l = 0
1744                nz_w_shift_l = 0
1745                nz_s_shift_l = 0
1746             ENDIF
1747
1748#if defined( __parallel )
1749             CALL MPI_ALLREDUCE(nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER,      &
1750                                MPI_MAX, comm2d, ierr)
1751             CALL MPI_ALLREDUCE(nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER,      &
1752                                MPI_MAX, comm2d, ierr)
1753             CALL MPI_ALLREDUCE(nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER,      & 
1754                                MPI_MAX, comm2d, ierr)
1755             CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
1756                                MPI_MAX, comm2d, ierr)
1757#else
1758             nz_u_shift = nz_u_shift_l
1759             nz_v_shift = nz_v_shift_l
1760             nz_w_shift = nz_w_shift_l
1761             nz_s_shift = nz_s_shift_l
1762#endif
1763
1764             mean_inflow_profiles(:,1) = 0.0_wp
1765             mean_inflow_profiles(nz_u_shift:nzt+1,1) = hom_sum(0:nzt+1-nz_u_shift,1,0)  ! u
1766
1767             mean_inflow_profiles(:,2) = 0.0_wp
1768             mean_inflow_profiles(nz_v_shift:nzt+1,2) = hom_sum(0:nzt+1-nz_v_shift,2,0)  ! v
1769
1770             mean_inflow_profiles(nz_s_shift:nzt+1,4) = hom_sum(0:nzt+1-nz_s_shift,4,0)  ! pt
1771
1772          ENDIF
1773
1774!
1775!--       If necessary, adjust the horizontal flow field to the prescribed
1776!--       profiles
1777          IF ( use_prescribed_profile_data )  THEN
1778             DO  i = nxlg, nxrg
1779                DO  j = nysg, nyng
1780                   DO  k = nzb, nzt+1
1781                      u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
1782                      v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
1783                   ENDDO
1784                ENDDO
1785             ENDDO
1786          ENDIF
1787
1788!
1789!--       Use these mean profiles at the inflow (provided that Dirichlet
1790!--       conditions are used)
1791          IF ( inflow_l )  THEN
1792             DO  j = nysg, nyng
1793                DO  k = nzb, nzt+1
1794                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
1795                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
1796                   w(k,j,nxlg:-1)  = 0.0_wp
1797                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
1798                   IF ( humidity )                                             &
1799                      q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
1800                   IF ( passive_scalar )                                       &
1801                      s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                     
1802                ENDDO
1803             ENDDO
1804          ENDIF
1805
1806!
1807!--       Calculate the damping factors to be used at the inflow. For a
1808!--       turbulent inflow the turbulent fluctuations have to be limited
1809!--       vertically because otherwise the turbulent inflow layer will grow
1810!--       in time.
1811          IF ( inflow_damping_height == 9999999.9_wp )  THEN
1812!
1813!--          Default: use the inversion height calculated by the prerun; if
1814!--          this is zero, inflow_damping_height must be explicitly
1815!--          specified.
1816             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
1817                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
1818             ELSE
1819                WRITE( message_string, * ) 'inflow_damping_height must be ',   &
1820                     'explicitly specified because&the inversion height ',     &
1821                     'calculated by the prerun is zero.'
1822                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1823             ENDIF
1824
1825          ENDIF
1826
1827          IF ( inflow_damping_width == 9999999.9_wp )  THEN
1828!
1829!--          Default for the transition range: one tenth of the undamped
1830!--          layer
1831             inflow_damping_width = 0.1_wp * inflow_damping_height
1832
1833          ENDIF
1834
1835          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1836
1837          DO  k = nzb, nzt+1
1838
1839             IF ( zu(k) <= inflow_damping_height )  THEN
1840                inflow_damping_factor(k) = 1.0_wp
1841             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
1842                inflow_damping_factor(k) = 1.0_wp -                            &
1843                                           ( zu(k) - inflow_damping_height ) / &
1844                                           inflow_damping_width
1845             ELSE
1846                inflow_damping_factor(k) = 0.0_wp
1847             ENDIF
1848
1849          ENDDO
1850
1851       ENDIF
1852
1853!
1854!--    Inside buildings set velocities back to zero
1855       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
1856            topography /= 'flat' )  THEN
1857!
1858!--       Inside buildings set velocities back to zero.
1859!--       Other scalars (pt, q, s, p, sa, ...) are ignored at present,
1860!--       maybe revise later.
1861          DO  i = nxlg, nxrg
1862             DO  j = nysg, nyng
1863                DO  k = nzb, nzt
1864                   u(k,j,i)     = MERGE( u(k,j,i), 0.0_wp,                     &
1865                                         BTEST( wall_flags_0(k,j,i), 1 ) )
1866                   v(k,j,i)     = MERGE( v(k,j,i), 0.0_wp,                     &
1867                                         BTEST( wall_flags_0(k,j,i), 2 ) )
1868                   w(k,j,i)     = MERGE( w(k,j,i), 0.0_wp,                     &
1869                                         BTEST( wall_flags_0(k,j,i), 3 ) )
1870                ENDDO
1871             ENDDO
1872          ENDDO
1873
1874       ENDIF
1875
1876!
1877!--    Calculate initial temperature field and other constants used in case
1878!--    of a sloping surface
1879       IF ( sloping_surface )  CALL init_slope
1880
1881!
1882!--    Initialize new time levels (only done in order to set boundary values
1883!--    including ghost points)
1884       pt_p = pt; u_p = u; v_p = v; w_p = w
1885       IF ( humidity )  THEN
1886          q_p = q
1887          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1888             qc_p = qc
1889             nc_p = nc
1890          ENDIF
1891          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1892             qr_p = qr
1893             nr_p = nr
1894          ENDIF
1895       ENDIF
1896       IF ( passive_scalar )  s_p  = s
1897       IF ( ocean          )  sa_p = sa
1898
1899!
1900!--    Allthough tendency arrays are set in prognostic_equations, they have
1901!--    have to be predefined here because they are used (but multiplied with 0)
1902!--    there before they are set.
1903       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1904       IF ( humidity )  THEN
1905          tq_m = 0.0_wp
1906          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1907             tqc_m = 0.0_wp
1908             tnc_m = 0.0_wp
1909          ENDIF
1910          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1911             tqr_m = 0.0_wp
1912             tnr_m = 0.0_wp
1913          ENDIF
1914       ENDIF
1915       IF ( passive_scalar )  ts_m  = 0.0_wp
1916       IF ( ocean          )  tsa_m = 0.0_wp
1917!
1918!--    Initialize synthetic turbulence generator in case of restart.
1919       IF ( TRIM( initializing_actions ) == 'read_restart_data'  .AND.         &
1920            use_syn_turb_gen )  CALL stg_init
1921
1922       CALL location_message( 'finished', .TRUE. )
1923
1924    ELSE
1925!
1926!--    Actually this part of the programm should not be reached
1927       message_string = 'unknown initializing problem'
1928       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1929    ENDIF
1930
1931!
1932!-- Initialize TKE, Kh and Km
1933    CALL tcm_init
1934
1935
1936    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1937!
1938!--    Initialize old timelevels needed for radiation boundary conditions
1939       IF ( outflow_l )  THEN
1940          u_m_l(:,:,:) = u(:,:,1:2)
1941          v_m_l(:,:,:) = v(:,:,0:1)
1942          w_m_l(:,:,:) = w(:,:,0:1)
1943       ENDIF
1944       IF ( outflow_r )  THEN
1945          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1946          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1947          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1948       ENDIF
1949       IF ( outflow_s )  THEN
1950          u_m_s(:,:,:) = u(:,0:1,:)
1951          v_m_s(:,:,:) = v(:,1:2,:)
1952          w_m_s(:,:,:) = w(:,0:1,:)
1953       ENDIF
1954       IF ( outflow_n )  THEN
1955          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1956          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1957          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1958       ENDIF
1959       
1960    ENDIF
1961
1962!
1963!-- Calculate the initial volume flow at the right and north boundary
1964    IF ( conserve_volume_flow )  THEN
1965
1966       IF ( use_prescribed_profile_data )  THEN
1967
1968          volume_flow_initial_l = 0.0_wp
1969          volume_flow_area_l    = 0.0_wp
1970
1971          IF ( nxr == nx )  THEN
1972             DO  j = nys, nyn
1973                DO  k = nzb+1, nzt
1974                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1975                                              u_init(k) * dzw(k)               &
1976                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1977                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
1978                                            )
1979
1980                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1981                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1982                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
1983                                            )
1984                ENDDO
1985             ENDDO
1986          ENDIF
1987         
1988          IF ( nyn == ny )  THEN
1989             DO  i = nxl, nxr
1990                DO  k = nzb+1, nzt
1991                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1992                                              v_init(k) * dzw(k)               &       
1993                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1994                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
1995                                            )
1996                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1997                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1998                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
1999                                            )
2000                ENDDO
2001             ENDDO
2002          ENDIF
2003
2004#if defined( __parallel )
2005          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
2006                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
2007          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
2008                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
2009
2010#else
2011          volume_flow_initial = volume_flow_initial_l
2012          volume_flow_area    = volume_flow_area_l
2013#endif 
2014
2015       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
2016
2017          volume_flow_initial_l = 0.0_wp
2018          volume_flow_area_l    = 0.0_wp
2019
2020          IF ( nxr == nx )  THEN
2021             DO  j = nys, nyn
2022                DO  k = nzb+1, nzt
2023                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
2024                                              hom_sum(k,1,0) * dzw(k)          &
2025                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2026                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
2027                                            )
2028                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
2029                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2030                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
2031                                            )
2032                ENDDO
2033             ENDDO
2034          ENDIF
2035         
2036          IF ( nyn == ny )  THEN
2037             DO  i = nxl, nxr
2038                DO  k = nzb+1, nzt
2039                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
2040                                              hom_sum(k,2,0) * dzw(k)          &       
2041                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2042                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
2043                                            )
2044                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
2045                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2046                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
2047                                            )
2048                ENDDO
2049             ENDDO
2050          ENDIF
2051
2052#if defined( __parallel )
2053          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
2054                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
2055          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
2056                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
2057
2058#else
2059          volume_flow_initial = volume_flow_initial_l
2060          volume_flow_area    = volume_flow_area_l
2061#endif 
2062
2063       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2064
2065          volume_flow_initial_l = 0.0_wp
2066          volume_flow_area_l    = 0.0_wp
2067
2068          IF ( nxr == nx )  THEN
2069             DO  j = nys, nyn
2070                DO  k = nzb+1, nzt
2071                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
2072                                              u(k,j,nx) * dzw(k)               &
2073                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2074                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
2075                                            )
2076                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
2077                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2078                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
2079                                            )
2080                ENDDO
2081             ENDDO
2082          ENDIF
2083         
2084          IF ( nyn == ny )  THEN
2085             DO  i = nxl, nxr
2086                DO  k = nzb+1, nzt
2087                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
2088                                              v(k,ny,i) * dzw(k)               &       
2089                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2090                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
2091                                            )
2092                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
2093                                     * MERGE( 1.0_wp, 0.0_wp,                  &
2094                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
2095                                            )
2096                ENDDO
2097             ENDDO
2098          ENDIF
2099
2100#if defined( __parallel )
2101          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
2102                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
2103          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
2104                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
2105
2106#else
2107          volume_flow_initial = volume_flow_initial_l
2108          volume_flow_area    = volume_flow_area_l
2109#endif 
2110
2111       ENDIF
2112
2113!
2114!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
2115!--    from u|v_bulk instead
2116       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
2117          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
2118          volume_flow_initial(2) = v_bulk * volume_flow_area(2)
2119       ENDIF
2120
2121    ENDIF
2122!
2123!-- Finally, if random_heatflux is set, disturb shf at horizontal
2124!-- surfaces. Actually, this should be done in surface_mod, where all other
2125!-- initializations of surface quantities are done. However, this
2126!-- would create a ring dependency, hence, it is done here. Maybe delete
2127!-- disturb_heatflux and tranfer the respective code directly into the
2128!-- initialization in surface_mod.         
2129    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
2130         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2131 
2132       IF ( use_surface_fluxes  .AND.  constant_heatflux  .AND.                &
2133            random_heatflux )  THEN
2134          IF ( surf_def_h(0)%ns >= 1 )  CALL disturb_heatflux( surf_def_h(0) )
2135          IF ( surf_lsm_h%ns    >= 1 )  CALL disturb_heatflux( surf_lsm_h    )
2136          IF ( surf_usm_h%ns    >= 1 )  CALL disturb_heatflux( surf_usm_h    )
2137       ENDIF
2138    ENDIF
2139
2140!
2141!-- Before initializing further modules, compute total sum of active mask
2142!-- grid points and the mean surface level height for each statistic region.
2143!-- ngp_2dh: number of grid points of a horizontal cross section through the
2144!--          total domain
2145!-- ngp_3d:  number of grid points of the total domain
2146    ngp_2dh_outer_l   = 0
2147    ngp_2dh_outer     = 0
2148    ngp_2dh_s_inner_l = 0
2149    ngp_2dh_s_inner   = 0
2150    ngp_2dh_l         = 0
2151    ngp_2dh           = 0
2152    ngp_3d_inner_l    = 0.0_wp
2153    ngp_3d_inner      = 0
2154    ngp_3d            = 0
2155    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
2156
2157    mean_surface_level_height   = 0.0_wp
2158    mean_surface_level_height_l = 0.0_wp
2159!
2160!-- Pre-set masks for regional statistics. Default is the total model domain.
2161!-- Ghost points are excluded because counting values at the ghost boundaries
2162!-- would bias the statistics
2163    rmask = 1.0_wp
2164    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
2165    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
2166!
2167!-- User-defined initializing actions
2168    CALL user_init
2169!
2170!-- To do: New concept for these non-topography grid points!
2171    DO  sr = 0, statistic_regions
2172       DO  i = nxl, nxr
2173          DO  j = nys, nyn
2174             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
2175!
2176!--             All xy-grid points
2177                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
2178!
2179!--             Determine mean surface-level height. In case of downward-
2180!--             facing walls are present, more than one surface level exist.
2181!--             In this case, use the lowest surface-level height.
2182                IF ( surf_def_h(0)%start_index(j,i) <=                         &
2183                     surf_def_h(0)%end_index(j,i) )  THEN
2184                   m = surf_def_h(0)%start_index(j,i)
2185                   k = surf_def_h(0)%k(m)
2186                   mean_surface_level_height_l(sr) =                           &
2187                                       mean_surface_level_height_l(sr) + zw(k-1)
2188                ENDIF
2189                IF ( surf_lsm_h%start_index(j,i) <=                            &
2190                     surf_lsm_h%end_index(j,i) )  THEN
2191                   m = surf_lsm_h%start_index(j,i)
2192                   k = surf_lsm_h%k(m)
2193                   mean_surface_level_height_l(sr) =                           &
2194                                       mean_surface_level_height_l(sr) + zw(k-1)
2195                ENDIF
2196                IF ( surf_usm_h%start_index(j,i) <=                            &
2197                     surf_usm_h%end_index(j,i) )  THEN
2198                   m = surf_usm_h%start_index(j,i)
2199                   k = surf_usm_h%k(m)
2200                   mean_surface_level_height_l(sr) =                           &
2201                                       mean_surface_level_height_l(sr) + zw(k-1)
2202                ENDIF
2203
2204                k_surf = k - 1
2205
2206                DO  k = nzb, nzt+1
2207!
2208!--                xy-grid points above topography
2209                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr)     +         &
2210                                  MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 24 ) )
2211
2212                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) +         &
2213                                  MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 22 ) )
2214
2215                ENDDO
2216!
2217!--             All grid points of the total domain above topography
2218                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + ( nz - k_surf + 2 )
2219
2220
2221
2222             ENDIF
2223          ENDDO
2224       ENDDO
2225    ENDDO
2226!
2227!-- Initialize arrays encompassing number of grid-points in inner and outer
2228!-- domains, statistic regions, etc. Mainly used for horizontal averaging
2229!-- of turbulence statistics. Please note, user_init must be called before
2230!-- doing this.   
2231    sr = statistic_regions + 1
2232#if defined( __parallel )
2233    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2234    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
2235                        comm2d, ierr )
2236    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2237    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
2238                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
2239    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2240    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
2241                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
2242    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2243    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
2244                        MPI_SUM, comm2d, ierr )
2245    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
2246    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2247    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
2248                        mean_surface_level_height(0), sr, MPI_REAL,            &
2249                        MPI_SUM, comm2d, ierr )
2250    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
2251#else
2252    ngp_2dh         = ngp_2dh_l
2253    ngp_2dh_outer   = ngp_2dh_outer_l
2254    ngp_2dh_s_inner = ngp_2dh_s_inner_l
2255    ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
2256    mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
2257#endif
2258
2259    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
2260             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
2261
2262!
2263!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
2264!-- buoyancy, etc. A zero value will occur for cases where all grid points of
2265!-- the respective subdomain lie below the surface topography
2266    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
2267    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
2268                           ngp_3d_inner(:) )
2269    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
2270
2271    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
2272                ngp_3d_inner_l, ngp_3d_inner_tmp )
2273
2274!
2275!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
2276!-- initialize heat-fluxes, etc. via datatype. Revise it later!
2277    IF ( large_scale_forcing .AND. lsf_surf )  THEN
2278       IF ( use_surface_fluxes  .AND.  constant_heatflux )  THEN
2279          CALL ls_forcing_surf ( simulated_time )
2280       ENDIF
2281    ENDIF
2282!
2283!-- Initialize quantities for special advections schemes
2284    CALL init_advec
2285
2286!
2287!-- Impose random perturbation on the horizontal velocity field and then
2288!-- remove the divergences from the velocity field at the initial stage
2289    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
2290         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
2291         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2292
2293       CALL location_message( 'creating initial disturbances', .FALSE. )
2294       CALL disturb_field( 'u', tend, u )
2295       CALL disturb_field( 'v', tend, v )
2296       CALL location_message( 'finished', .TRUE. )
2297
2298       CALL location_message( 'calling pressure solver', .FALSE. )
2299       n_sor = nsor_ini
2300       CALL pres
2301       n_sor = nsor
2302       CALL location_message( 'finished', .TRUE. )
2303
2304    ENDIF
2305
2306!
2307!-- If required, initialize quantities needed for the plant canopy model
2308    IF ( plant_canopy )  THEN
2309       CALL location_message( 'initializing plant canopy model', .FALSE. )   
2310       CALL pcm_init
2311       CALL location_message( 'finished', .TRUE. )
2312    ENDIF
2313
2314!
2315!-- If required, initialize dvrp-software
2316    IF ( dt_dvrp /= 9999999.9_wp )  CALL init_dvrp
2317
2318    IF ( ocean )  THEN
2319!
2320!--    Initialize quantities needed for the ocean model
2321       CALL init_ocean
2322
2323    ELSE
2324!
2325!--    Initialize quantities for handling cloud physics
2326!--    This routine must be called before lpm_init, because
2327!--    otherwise, array pt_d_t, needed in data_output_dvrp (called by
2328!--    lpm_init) is not defined.
2329       CALL init_cloud_physics
2330!
2331!--    Initialize bulk cloud microphysics
2332       CALL microphysics_init
2333    ENDIF
2334
2335!
2336!-- If required, initialize particles
2337    IF ( particle_advection )  CALL lpm_init
2338
2339!
2340!-- If required, initialize quantities needed for the LSM
2341    IF ( land_surface )  THEN
2342       CALL location_message( 'initializing land surface model', .FALSE. )
2343       CALL lsm_init
2344       CALL location_message( 'finished', .TRUE. )
2345    ENDIF
2346
2347!
2348!-- If required, allocate USM and LSM surfaces
2349    IF ( urban_surface )  THEN
2350       CALL location_message( 'initializing and allocating urban surfaces', .FALSE. )
2351       CALL usm_allocate_surface
2352       CALL location_message( 'finished', .TRUE. )
2353    ENDIF
2354!
2355!-- If required, initialize urban surface model
2356    IF ( urban_surface )  THEN
2357       CALL location_message( 'initializing urban surface model', .FALSE. )
2358       CALL usm_init_urban_surface
2359       CALL location_message( 'finished', .TRUE. )
2360    ENDIF
2361
2362!
2363!-- Initialize surface layer (done after LSM as roughness length are required
2364!-- for initialization
2365    IF ( constant_flux_layer )  THEN
2366       CALL location_message( 'initializing surface layer', .FALSE. )
2367       CALL init_surface_layer_fluxes
2368       CALL location_message( 'finished', .TRUE. )
2369    ENDIF
2370
2371!
2372!-- If required, set chemical emissions
2373!-- (todo(FK): This should later on be CALLed time-dependently in init_3d_model)
2374    IF ( air_chemistry )  THEN
2375       CALL chem_emissions
2376    ENDIF
2377
2378!
2379!-- Initialize radiation processes
2380    IF ( radiation )  THEN
2381!
2382!--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees.
2383!--    The namelist parameter radiation_interactions_on can override this behavior.
2384!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
2385!--    init_surface_arrays.)
2386       IF ( radiation_interactions_on )  THEN
2387          IF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
2388             radiation_interactions    = .TRUE.
2389             average_radiation         = .TRUE.
2390          ELSE
2391             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions
2392                                                   !< calculations necessary in case of flat surface
2393          ENDIF
2394       ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
2395          message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
2396                           'vertical surfaces and/or trees exist. The model will run ' // &
2397                           'without RTM (no shadows, no radiation reflections)'
2398          CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
2399       ENDIF
2400!
2401!--    If required, initialize radiation interactions between surfaces
2402!--    via sky-view factors. This must be done before radiation is initialized.
2403       IF ( radiation_interactions )  CALL radiation_interaction_init
2404
2405!
2406!--    Initialize radiation model
2407       CALL location_message( 'initializing radiation model', .FALSE. )
2408       CALL radiation_init
2409       CALL location_message( 'finished', .TRUE. )
2410
2411!
2412!--    Find all discretized apparent solar positions for radiation interaction.
2413!--    This must be done after radiation_init.
2414       IF ( radiation_interactions )  CALL radiation_presimulate_solar_pos
2415
2416!
2417!--    If required, read or calculate and write out the SVF
2418       IF ( radiation_interactions .AND. read_svf)  THEN
2419!
2420!--       Read sky-view factors and further required data from file
2421          CALL location_message( '    Start reading SVF from file', .FALSE. )
2422          CALL radiation_read_svf()
2423          CALL location_message( '    Reading SVF from file has finished', .TRUE. )
2424
2425       ELSEIF ( radiation_interactions .AND. .NOT. read_svf)  THEN
2426!
2427!--       calculate SFV and CSF
2428          CALL location_message( '    Start calculation of SVF', .FALSE. )
2429          CALL radiation_calc_svf()
2430          CALL location_message( '    Calculation of SVF has finished', .TRUE. )
2431       ENDIF
2432
2433       IF ( radiation_interactions .AND. write_svf)  THEN
2434!
2435!--       Write svf, csf svfsurf and csfsurf data to file
2436          CALL location_message( '    Start writing SVF in file', .FALSE. )
2437          CALL radiation_write_svf()
2438          CALL location_message( '    Writing SVF in file has finished', .TRUE. )
2439       ENDIF
2440
2441!
2442!--    Adjust radiative fluxes. In case of urban and land surfaces, also
2443!--    call an initial interaction.
2444       IF ( radiation_interactions )  THEN
2445          CALL radiation_interaction
2446       ENDIF
2447    ENDIF
2448
2449!
2450!-- Temporary solution to add LSM and radiation time series to the default
2451!-- output
2452    IF ( land_surface  .OR.  radiation )  THEN
2453       IF ( TRIM( radiation_scheme ) == 'rrtmg' )  THEN
2454          dots_num = dots_num + 15
2455       ELSE
2456          dots_num = dots_num + 11
2457       ENDIF
2458    ENDIF
2459   
2460
2461
2462!
2463!-- If required, initialize quantities needed for the wind turbine model
2464    IF ( wind_turbine )  THEN
2465       CALL location_message( 'initializing wind turbine model', .FALSE. )
2466       CALL wtm_init
2467       CALL location_message( 'finished', .TRUE. )
2468    ENDIF
2469
2470!
2471!-- If required, initialize quantities needed for the gust module
2472    IF ( gust_module_enabled )  THEN
2473       CALL gust_init( dots_label, dots_unit, dots_num, dots_max )
2474    ENDIF
2475
2476!
2477!-- Initialize the ws-scheme.   
2478    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init       
2479
2480!
2481!-- Setting weighting factors for calculation of perturbation pressure
2482!-- and turbulent quantities from the RK substeps
2483    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
2484
2485       weight_substep(1) = 1._wp/6._wp
2486       weight_substep(2) = 3._wp/10._wp
2487       weight_substep(3) = 8._wp/15._wp
2488
2489       weight_pres(1)    = 1._wp/3._wp
2490       weight_pres(2)    = 5._wp/12._wp
2491       weight_pres(3)    = 1._wp/4._wp
2492
2493    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
2494
2495       weight_substep(1) = 1._wp/2._wp
2496       weight_substep(2) = 1._wp/2._wp
2497         
2498       weight_pres(1)    = 1._wp/2._wp
2499       weight_pres(2)    = 1._wp/2._wp       
2500
2501    ELSE                                     ! for Euler-method
2502
2503       weight_substep(1) = 1.0_wp     
2504       weight_pres(1)    = 1.0_wp                   
2505
2506    ENDIF
2507
2508!
2509!-- Initialize Rayleigh damping factors
2510    rdf    = 0.0_wp
2511    rdf_sc = 0.0_wp
2512    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
2513       IF (  .NOT.  ocean )  THEN
2514          DO  k = nzb+1, nzt
2515             IF ( zu(k) >= rayleigh_damping_height )  THEN
2516                rdf(k) = rayleigh_damping_factor *                             &
2517                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
2518                             / ( zu(nzt) - rayleigh_damping_height ) )         &
2519                      )**2
2520             ENDIF
2521          ENDDO
2522       ELSE
2523          DO  k = nzt, nzb+1, -1
2524             IF ( zu(k) <= rayleigh_damping_height )  THEN
2525                rdf(k) = rayleigh_damping_factor *                             &
2526                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
2527                             / ( rayleigh_damping_height - zu(nzb+1) ) )       &
2528                      )**2
2529             ENDIF
2530          ENDDO
2531       ENDIF
2532    ENDIF
2533    IF ( scalar_rayleigh_damping )  rdf_sc = rdf
2534
2535!
2536!-- Initialize the starting level and the vertical smoothing factor used for
2537!-- the external pressure gradient
2538    dp_smooth_factor = 1.0_wp
2539    IF ( dp_external )  THEN
2540!
2541!--    Set the starting level dp_level_ind_b only if it has not been set before
2542!--    (e.g. in init_grid).
2543       IF ( dp_level_ind_b == 0 )  THEN
2544          ind_array = MINLOC( ABS( dp_level_b - zu ) )
2545          dp_level_ind_b = ind_array(1) - 1 + nzb
2546                                        ! MINLOC uses lower array bound 1
2547       ENDIF
2548       IF ( dp_smooth )  THEN
2549          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
2550          DO  k = dp_level_ind_b+1, nzt
2551             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *               &
2552                        ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
2553                          REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
2554          ENDDO
2555       ENDIF
2556    ENDIF
2557
2558!
2559!-- Initialize damping zone for the potential temperature in case of
2560!-- non-cyclic lateral boundaries. The damping zone has the maximum value
2561!-- at the inflow boundary and decreases to zero at pt_damping_width.
2562    ptdf_x = 0.0_wp
2563    ptdf_y = 0.0_wp
2564    IF ( bc_lr_dirrad )  THEN
2565       DO  i = nxl, nxr
2566          IF ( ( i * dx ) < pt_damping_width )  THEN
2567             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
2568                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
2569                            REAL( pt_damping_width, KIND=wp ) ) ) )**2 
2570          ENDIF
2571       ENDDO
2572    ELSEIF ( bc_lr_raddir )  THEN
2573       DO  i = nxl, nxr
2574          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
2575             ptdf_x(i) = pt_damping_factor *                                   &
2576                         SIN( pi * 0.5_wp *                                    &
2577                                 ( ( i - nx ) * dx + pt_damping_width ) /      &
2578                                 REAL( pt_damping_width, KIND=wp ) )**2
2579          ENDIF
2580       ENDDO
2581    ELSEIF ( bc_ns_dirrad )  THEN
2582       DO  j = nys, nyn
2583          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
2584             ptdf_y(j) = pt_damping_factor *                                   &
2585                         SIN( pi * 0.5_wp *                                    &
2586                                 ( ( j - ny ) * dy + pt_damping_width ) /      &
2587                                 REAL( pt_damping_width, KIND=wp ) )**2
2588          ENDIF
2589       ENDDO
2590    ELSEIF ( bc_ns_raddir )  THEN
2591       DO  j = nys, nyn
2592          IF ( ( j * dy ) < pt_damping_width )  THEN
2593             ptdf_y(j) = pt_damping_factor *                                   &
2594                         SIN( pi * 0.5_wp *                                    &
2595                                ( pt_damping_width - j * dy ) /                &
2596                                REAL( pt_damping_width, KIND=wp ) )**2
2597          ENDIF
2598       ENDDO
2599    ENDIF
2600!
2601!-- Check if maximum number of allowed timeseries is exceeded
2602    IF ( dots_num > dots_max )  THEN
2603       WRITE( message_string, * ) 'number of time series quantities exceeds',  &
2604                                  ' its maximum of dots_max = ', dots_max,     &
2605                                  ' &Please increase dots_max in modules.f90.'
2606       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )   
2607    ENDIF
2608
2609!
2610!-- Input binary data file is not needed anymore. This line must be placed
2611!-- after call of user_init!
2612    CALL close_file( 13 )
2613!
2614!-- In case of nesting, put an barrier to assure that all parent and child
2615!-- domains finished initialization.
2616#if defined( __parallel )
2617    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
2618#endif
2619
2620
2621    CALL location_message( 'leaving init_3d_model', .TRUE. )
2622
2623 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.