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

Last change on this file since 3019 was 3014, checked in by maronga, 6 years ago

series of bugfixes

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