Home

Context Navigation

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

Last change on this file since 1015 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 56.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |