Home

Context Navigation

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

Last change on this file since 2550 was 2550, checked in by boeske, 7 years ago

enable simulations with complex terrain

Property svn:keywords set to Id
Property svn:mergeinfo set to False
/palm/branches/forwind/SOURCE/init_3d_model.f90 1564-1913

File size: 87.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |