Home

Context Navigation

init_3d_model.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 13 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

Property svn:keywords set to Id

File size: 55.3 KB

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

Download in other formats: