Home

Context Navigation

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

Last change on this file since 311 was 292, checked in by letzel, 16 years ago
Typographical errors (netcdf) If the inversion height calculated by the prerun is zero, inflow_damping_height must be explicitly specified (init_3d_model) Small bugfix concerning 3d 64bit netcdf output format (header) Adjustments for lcsgib and lcsgih (subjob)
Property svn:keywords set to `Id`
File size: 51.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |