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

Last change on this file since 2995 was 2995, checked in by Giersch, 6 years ago

Bugfixes related to spinup and radiation model

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