Home

Context Navigation

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

Last change on this file since 2252 was 2252, checked in by knoop, 7 years ago

air density now depending on surface_pressure even in boussinesq mode

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |