Home

Context Navigation

init_3d_model.f90 @ 1179

Last change on this file since 1179 was 1179, checked in by raasch, 11 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:

small bugfix for background communication (time_integration)

Property svn:keywords set to Id

File size: 64.9 KB

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

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

Context Navigation

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

Download in other formats: