Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |