Home

Context Navigation

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

Last change on this file since 106 was 106, checked in by raasch, 18 years ago
preliminary update of bugfixes and extensions for non-cyclic BCs
Property svn:keywords set to `Id`
File size: 36.2 KB

Line
1	#if defined( __ibmy_special )
2	@PROCESS NOOPTimize
3	#endif
4	SUBROUTINE init_3d_model
5
6	!------------------------------------------------------------------------------!
7	! Actual revisions:
8	! -----------------
9	! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
10	! +arrays for phase speed c_u, c_v, c_w, indices for u\|v\|w_m_l\|r changed
11	!
12	! Former revisions:
13	! -----------------
14	! $Id: init_3d_model.f90 106 2007-08-16 14:30:26Z raasch $
15	!
16	! 97 2007-06-21 08:23:15Z raasch
17	! Initialization of salinity, call of init_ocean
18	!
19	! 87 2007-05-22 15:46:47Z raasch
20	! var_hom and var_sum renamed pr_palm
21	!
22	! 75 2007-03-22 09:54:05Z raasch
23	! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
24	! bugfix for cases with the outflow damping layer extending over more than one
25	! subdomain, moisture renamed humidity,
26	! new initializing action "by_user" calls user_init_3d_model,
27	! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
28	! initial velocities at nzb+1 are regarded for volume
29	! flow control in case they have been set zero before (to avoid small timesteps)
30	! -uvmean_outflow, uxrp, vynp eliminated
31	!
32	! 19 2007-02-23 04:53:48Z raasch
33	! +handling of top fluxes
34	!
35	! RCS Log replace by Id keyword, revision history cleaned up
36	!
37	! Revision 1.49 2006/08/22 15:59:07 raasch
38	! No optimization of this file on the ibmy (Yonsei Univ.)
39	!
40	! Revision 1.1 1998/03/09 16:22:22 raasch
41	! Initial revision
42	!
43	!
44	! Description:
45	! ------------
46	! Allocation of arrays and initialization of the 3D model via
47	! a) pre-run the 1D model
48	! or
49	! b) pre-set constant linear profiles
50	! or
51	! c) read values of a previous run
52	!------------------------------------------------------------------------------!
53
54	USE arrays_3d
55	USE averaging
56	USE cloud_parameters
57	USE constants
58	USE control_parameters
59	USE cpulog
60	USE indices
61	USE interfaces
62	USE model_1d
63	USE netcdf_control
64	USE particle_attributes
65	USE pegrid
66	USE profil_parameter
67	USE random_function_mod
68	USE statistics
69
70	IMPLICIT NONE
71
72	INTEGER :: i, j, k, sr
73
74	INTEGER, DIMENSION(:), ALLOCATABLE :: ngp_2dh_l, ngp_3d_inner_l
75
76	INTEGER, DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_outer_l
77
78	REAL, DIMENSION(1:2) :: volume_flow_area_l, volume_flow_initial_l
79
80
81	!
82	!-- Allocate arrays
83	ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
84	ngp_3d(0:statistic_regions), &
85	ngp_3d_inner(0:statistic_regions), &
86	ngp_3d_inner_l(0:statistic_regions), &
87	sums_divnew_l(0:statistic_regions), &
88	sums_divold_l(0:statistic_regions) )
89	ALLOCATE( rdf(nzb+1:nzt) )
90	ALLOCATE( hom_sum(nzb:nzt+1,pr_palm+max_pr_user,0:statistic_regions), &
91	ngp_2dh_outer(nzb:nzt+1,0:statistic_regions), &
92	ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions), &
93	rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions), &
94	sums(nzb:nzt+1,pr_palm+max_pr_user), &
95	sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1), &
96	sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
97	sums_up_fraction_l(10,3,0:statistic_regions), &
98	sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), &
99	ts_value(var_ts,0:statistic_regions) )
100	ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
101
102	ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
103	ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1), &
104	us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1), &
105	uswst_1(nys-1:nyn+1,nxl-1:nxr+1), &
106	vsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
107	vswst_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
108
109	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
110	!
111	!-- Leapfrog scheme needs two timelevels of diffusion quantities
112	ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1), &
113	shf_2(nys-1:nyn+1,nxl-1:nxr+1), &
114	tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
115	usws_2(nys-1:nyn+1,nxl-1:nxr+1), &
116	uswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
117	vswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
118	vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
119	ENDIF
120
121	ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra), &
122	e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
123	e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
124	e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
125	kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
126	km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
127	p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
128	pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
129	pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
130	pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
131	tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
132	u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
133	u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
134	u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
135	v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
136	v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
137	v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
138	w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
139	w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
140	w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
141
142	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
143	ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
144	km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
145	ENDIF
146
147	IF ( humidity .OR. passive_scalar ) THEN
148	!
149	!-- 2D-humidity/scalar arrays
150	ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1), &
151	qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
152	qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
153
154	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
155	ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
156	qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
157	ENDIF
158	!
159	!-- 3D-humidity/scalar arrays
160	ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
161	q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
162	q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
163
164	!
165	!-- 3D-arrays needed for humidity only
166	IF ( humidity ) THEN
167	ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
168
169	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
170	ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
171	ENDIF
172
173	IF ( cloud_physics ) THEN
174	!
175	!-- Liquid water content
176	ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
177	!
178	!-- Precipitation amount and rate (only needed if output is switched)
179	ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
180	precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
181	ENDIF
182
183	IF ( cloud_droplets ) THEN
184	!
185	!-- Liquid water content, change in liquid water content,
186	!-- real volume of particles (with weighting), volume of particles
187	ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
188	ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
189	ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
190	ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
191	ENDIF
192
193	ENDIF
194
195	ENDIF
196
197	IF ( ocean ) THEN
198	ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
199	saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
200	ALLOCATE( rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
201	sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
202	sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
203	sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
204	rho => rho_1 ! routine calc_mean_profile requires density to be a
205	! pointer
206	ENDIF
207
208	!
209	!-- 3D-array for storing the dissipation, needed for calculating the sgs
210	!-- particle velocities
211	IF ( use_sgs_for_particles ) THEN
212	ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
213	ENDIF
214
215	IF ( dt_dosp /= 9999999.9 ) THEN
216	ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
217	spectrum_y( 1:ny/2, 1:10, 1:10 ) )
218	ENDIF
219
220	!
221	!-- 4D-array for storing the Rif-values at vertical walls
222	IF ( topography /= 'flat' ) THEN
223	ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
224	rif_wall = 0.0
225	ENDIF
226
227	!
228	!-- Velocities at nzb+1 needed for volume flow control
229	IF ( conserve_volume_flow ) THEN
230	ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
231	u_nzb_p1_for_vfc = 0.0
232	v_nzb_p1_for_vfc = 0.0
233	ENDIF
234
235	!
236	!-- Arrays to store velocity data from t-dt and the phase speeds which
237	!-- are needed for radiation boundary conditions
238	IF ( outflow_l ) THEN
239	ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
240	v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
241	w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
242	ENDIF
243	IF ( outflow_r ) THEN
244	ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
245	v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
246	w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
247	ENDIF
248	IF ( outflow_l .OR. outflow_r ) THEN
249	ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
250	c_w(nzb:nzt+1,nys-1:nyn+1) )
251	ENDIF
252	IF ( outflow_s ) THEN
253	ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
254	v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
255	w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
256	ENDIF
257	IF ( outflow_n ) THEN
258	ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
259	v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
260	w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
261	ENDIF
262	IF ( outflow_s .OR. outflow_n ) THEN
263	ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
264	c_w(nzb:nzt+1,nxl-1:nxr+1) )
265	ENDIF
266
267	!
268	!-- Initial assignment of the pointers
269	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
270
271	rif_m => rif_1; rif => rif_2
272	shf_m => shf_1; shf => shf_2
273	tswst_m => tswst_1; tswst => tswst_2
274	usws_m => usws_1; usws => usws_2
275	uswst_m => uswst_1; uswst => uswst_2
276	vsws_m => vsws_1; vsws => vsws_2
277	vswst_m => vswst_1; vswst => vswst_2
278	e_m => e_1; e => e_2; e_p => e_3; te_m => e_3
279	kh_m => kh_1; kh => kh_2
280	km_m => km_1; km => km_2
281	pt_m => pt_1; pt => pt_2; pt_p => pt_3; tpt_m => pt_3
282	u_m => u_1; u => u_2; u_p => u_3; tu_m => u_3
283	v_m => v_1; v => v_2; v_p => v_3; tv_m => v_3
284	w_m => w_1; w => w_2; w_p => w_3; tw_m => w_3
285
286	IF ( humidity .OR. passive_scalar ) THEN
287	qsws_m => qsws_1; qsws => qsws_2
288	qswst_m => qswst_1; qswst => qswst_2
289	q_m => q_1; q => q_2; q_p => q_3; tq_m => q_3
290	IF ( humidity ) vpt_m => vpt_1; vpt => vpt_2
291	IF ( cloud_physics ) ql => ql_1
292	IF ( cloud_droplets ) THEN
293	ql => ql_1
294	ql_c => ql_2
295	ENDIF
296	ENDIF
297
298	ELSE
299
300	rif => rif_1
301	shf => shf_1
302	tswst => tswst_1
303	usws => usws_1
304	uswst => uswst_1
305	vsws => vsws_1
306	vswst => vswst_1
307	e => e_1; e_p => e_2; te_m => e_3; e_m => e_3
308	kh => kh_1
309	km => km_1
310	pt => pt_1; pt_p => pt_2; tpt_m => pt_3; pt_m => pt_3
311	u => u_1; u_p => u_2; tu_m => u_3; u_m => u_3
312	v => v_1; v_p => v_2; tv_m => v_3; v_m => v_3
313	w => w_1; w_p => w_2; tw_m => w_3; w_m => w_3
314
315	IF ( humidity .OR. passive_scalar ) THEN
316	qsws => qsws_1
317	qswst => qswst_1
318	q => q_1; q_p => q_2; tq_m => q_3; q_m => q_3
319	IF ( humidity ) vpt => vpt_1
320	IF ( cloud_physics ) ql => ql_1
321	IF ( cloud_droplets ) THEN
322	ql => ql_1
323	ql_c => ql_2
324	ENDIF
325	ENDIF
326
327	IF ( ocean ) THEN
328	saswsb => saswsb_1
329	saswst => saswst_1
330	sa => sa_1; sa_p => sa_2; tsa_m => sa_3
331	ENDIF
332
333	ENDIF
334
335	!
336	!-- Initialize model variables
337	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
338	!
339	!-- First model run of a possible job queue.
340	!-- Initial profiles of the variables must be computes.
341	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
342	!
343	!-- Use solutions of the 1D model as initial profiles,
344	!-- start 1D model
345	CALL init_1d_model
346	!
347	!-- Transfer initial profiles to the arrays of the 3D model
348	DO i = nxl-1, nxr+1
349	DO j = nys-1, nyn+1
350	e(:,j,i) = e1d
351	kh(:,j,i) = kh1d
352	km(:,j,i) = km1d
353	pt(:,j,i) = pt_init
354	u(:,j,i) = u1d
355	v(:,j,i) = v1d
356	ENDDO
357	ENDDO
358
359	IF ( humidity .OR. passive_scalar ) THEN
360	DO i = nxl-1, nxr+1
361	DO j = nys-1, nyn+1
362	q(:,j,i) = q_init
363	ENDDO
364	ENDDO
365	ENDIF
366
367	IF ( .NOT. constant_diffusion ) THEN
368	DO i = nxl-1, nxr+1
369	DO j = nys-1, nyn+1
370	e(:,j,i) = e1d
371	ENDDO
372	ENDDO
373	!
374	!-- Store initial profiles for output purposes etc.
375	hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
376
377	IF ( prandtl_layer ) THEN
378	rif = rif1d(nzb+1)
379	ts = 0.0 ! could actually be computed more accurately in the
380	! 1D model. Update when opportunity arises.
381	us = us1d
382	usws = usws1d
383	vsws = vsws1d
384	ELSE
385	ts = 0.0 ! must be set, because used in
386	rif = 0.0 ! flowste
387	us = 0.0
388	usws = 0.0
389	vsws = 0.0
390	ENDIF
391
392	ELSE
393	e = 0.0 ! must be set, because used in
394	rif = 0.0 ! flowste
395	ts = 0.0
396	us = 0.0
397	usws = 0.0
398	vsws = 0.0
399	ENDIF
400	uswst = top_momentumflux_u
401	vswst = top_momentumflux_v
402
403	!
404	!-- In every case qs = 0.0 (see also pt)
405	!-- This could actually be computed more accurately in the 1D model.
406	!-- Update when opportunity arises!
407	IF ( humidity .OR. passive_scalar ) qs = 0.0
408
409	!
410	!-- inside buildings set velocities back to zero
411	IF ( topography /= 'flat' ) THEN
412	DO i = nxl-1, nxr+1
413	DO j = nys-1, nyn+1
414	u(nzb:nzb_u_inner(j,i),j,i) = 0.0
415	v(nzb:nzb_v_inner(j,i),j,i) = 0.0
416	ENDDO
417	ENDDO
418	!
419	!-- WARNING: The extra boundary conditions set after running the
420	!-- ------- 1D model impose an error on the divergence one layer
421	!-- below the topography; need to correct later
422	!-- ATTENTION: Provisional correction for Piacsek & Williams
423	!-- --------- advection scheme: keep u and v zero one layer below
424	!-- the topography.
425	IF ( ibc_uv_b == 0 ) THEN
426	!
427	!-- Satisfying the Dirichlet condition with an extra layer below
428	!-- the surface where the u and v component change their sign.
429	DO i = nxl-1, nxr+1
430	DO j = nys-1, nyn+1
431	IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
432	IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
433	ENDDO
434	ENDDO
435
436	ELSE
437	!
438	!-- Neumann condition
439	DO i = nxl-1, nxr+1
440	DO j = nys-1, nyn+1
441	IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
442	IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
443	ENDDO
444	ENDDO
445
446	ENDIF
447
448	ENDIF
449
450	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
451	THEN
452	!
453	!-- Use constructed initial profiles (velocity constant with height,
454	!-- temperature profile with constant gradient)
455	DO i = nxl-1, nxr+1
456	DO j = nys-1, nyn+1
457	pt(:,j,i) = pt_init
458	u(:,j,i) = u_init
459	v(:,j,i) = v_init
460	ENDDO
461	ENDDO
462
463	!
464	!-- Set initial horizontal velocities at the lowest computational grid levels
465	!-- to zero in order to avoid too small time steps caused by the diffusion
466	!-- limit in the initial phase of a run (at k=1, dz/2 occurs in the
467	!-- limiting formula!). The original values are stored to be later used for
468	!-- volume flow control.
469	DO i = nxl-1, nxr+1
470	DO j = nys-1, nyn+1
471	u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
472	v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
473	ENDDO
474	ENDDO
475	IF ( conserve_volume_flow ) THEN
476	IF ( nxr == nx ) THEN
477	DO j = nys, nyn
478	k = nzb_u_inner(j,nx) + 1
479	u_nzb_p1_for_vfc(j) = u_init(k) * dzu(k)
480	ENDDO
481	ENDIF
482	IF ( nyn == ny ) THEN
483	DO i = nxl, nxr
484	k = nzb_v_inner(ny,i) + 1
485	v_nzb_p1_for_vfc(i) = v_init(k) * dzu(k)
486	ENDDO
487	ENDIF
488	ENDIF
489
490	IF ( humidity .OR. passive_scalar ) THEN
491	DO i = nxl-1, nxr+1
492	DO j = nys-1, nyn+1
493	q(:,j,i) = q_init
494	ENDDO
495	ENDDO
496	ENDIF
497
498	IF ( ocean ) THEN
499	DO i = nxl-1, nxr+1
500	DO j = nys-1, nyn+1
501	sa(:,j,i) = sa_init
502	ENDDO
503	ENDDO
504	ENDIF
505
506	IF ( constant_diffusion ) THEN
507	km = km_constant
508	kh = km / prandtl_number
509	ELSE
510	kh = 0.01 ! there must exist an initial diffusion, because
511	km = 0.01 ! otherwise no TKE would be produced by the
512	! production terms, as long as not yet
513	! e = (u/cm)*2 at k=nzb+1
514	ENDIF
515	e = 0.0
516	rif = 0.0
517	ts = 0.0
518	us = 0.0
519	usws = 0.0
520	uswst = top_momentumflux_u
521	vsws = 0.0
522	vswst = top_momentumflux_v
523	IF ( humidity .OR. passive_scalar ) qs = 0.0
524
525	!
526	!-- Compute initial temperature field and other constants used in case
527	!-- of a sloping surface
528	IF ( sloping_surface ) CALL init_slope
529
530	ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
531	THEN
532	!
533	!-- Initialization will completely be done by the user
534	CALL user_init_3d_model
535
536	ENDIF
537
538	!
539	!-- Calculate virtual potential temperature
540	IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
541
542	!
543	!-- Store initial profiles for output purposes etc.
544	hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
545	hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
546	IF ( ibc_uv_b == 0 ) THEN
547	hom(nzb,1,5,:) = -hom(nzb+1,1,5,:) ! due to satisfying the Dirichlet
548	hom(nzb,1,6,:) = -hom(nzb+1,1,6,:) ! condition with an extra layer
549	! below the surface where the u and v component change their sign
550	ENDIF
551	hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
552	hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
553	hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
554
555	IF ( ocean ) THEN
556	!
557	!-- Store initial salinity profile
558	hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
559	ENDIF
560
561	IF ( humidity ) THEN
562	!
563	!-- Store initial profile of total water content, virtual potential
564	!-- temperature
565	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
566	hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
567	IF ( cloud_physics .OR. cloud_droplets ) THEN
568	!
569	!-- Store initial profile of specific humidity and potential
570	!-- temperature
571	hom(:,1,27,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
572	hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
573	ENDIF
574	ENDIF
575
576	IF ( passive_scalar ) THEN
577	!
578	!-- Store initial scalar profile
579	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
580	ENDIF
581
582	!
583	!-- Initialize fluxes at bottom surface
584	IF ( use_surface_fluxes ) THEN
585
586	IF ( constant_heatflux ) THEN
587	!
588	!-- Heat flux is prescribed
589	IF ( random_heatflux ) THEN
590	CALL disturb_heatflux
591	ELSE
592	shf = surface_heatflux
593	!
594	!-- Over topography surface_heatflux is replaced by wall_heatflux(0)
595	IF ( TRIM( topography ) /= 'flat' ) THEN
596	DO i = nxl-1, nxr+1
597	DO j = nys-1, nyn+1
598	IF ( nzb_s_inner(j,i) /= 0 ) THEN
599	shf(j,i) = wall_heatflux(0)
600	ENDIF
601	ENDDO
602	ENDDO
603	ENDIF
604	ENDIF
605	IF ( ASSOCIATED( shf_m ) ) shf_m = shf
606	ENDIF
607
608	!
609	!-- Determine the near-surface water flux
610	IF ( humidity .OR. passive_scalar ) THEN
611	IF ( constant_waterflux ) THEN
612	qsws = surface_waterflux
613	IF ( ASSOCIATED( qsws_m ) ) qsws_m = qsws
614	ENDIF
615	ENDIF
616
617	ENDIF
618
619	!
620	!-- Initialize fluxes at top surface
621	!-- Currently, only the heatflux and salinity flux can be prescribed.
622	!-- The latent flux is zero in this case!
623	IF ( use_top_fluxes ) THEN
624
625	IF ( constant_top_heatflux ) THEN
626	!
627	!-- Heat flux is prescribed
628	tswst = top_heatflux
629	IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst
630
631	IF ( humidity .OR. passive_scalar ) THEN
632	qswst = 0.0
633	IF ( ASSOCIATED( qswst_m ) ) qswst_m = qswst
634	ENDIF
635
636	IF ( ocean ) THEN
637	saswsb = bottom_salinityflux
638	saswst = top_salinityflux
639	ENDIF
640	ENDIF
641
642	!
643	!-- Initialization in case of a coupled model run
644	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
645	tswst = 0.0
646	IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst
647	ENDIF
648
649	ENDIF
650
651	!
652	!-- Initialize Prandtl layer quantities
653	IF ( prandtl_layer ) THEN
654
655	z0 = roughness_length
656
657	IF ( .NOT. constant_heatflux ) THEN
658	!
659	!-- Surface temperature is prescribed. Here the heat flux cannot be
660	!-- simply estimated, because therefore rif, u* and theta* would have
661	!-- to be computed by iteration. This is why the heat flux is assumed
662	!-- to be zero before the first time step. It approaches its correct
663	!-- value in the course of the first few time steps.
664	shf = 0.0
665	IF ( ASSOCIATED( shf_m ) ) shf_m = 0.0
666	ENDIF
667
668	IF ( humidity .OR. passive_scalar ) THEN
669	IF ( .NOT. constant_waterflux ) THEN
670	qsws = 0.0
671	IF ( ASSOCIATED( qsws_m ) ) qsws_m = 0.0
672	ENDIF
673	ENDIF
674
675	ENDIF
676
677	!
678	!-- Calculate the initial volume flow at the right and north boundary
679	IF ( conserve_volume_flow ) THEN
680
681	volume_flow_initial_l = 0.0
682	volume_flow_area_l = 0.0
683
684	IF ( nxr == nx ) THEN
685	DO j = nys, nyn
686	DO k = nzb_2d(j,nx) + 1, nzt
687	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
688	u(k,j,nx) * dzu(k)
689	volume_flow_area_l(1) = volume_flow_area_l(1) + dzu(k)
690	ENDDO
691	!
692	!-- Correction if velocity at nzb+1 has been set zero further above
693	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
694	u_nzb_p1_for_vfc(j)
695	ENDDO
696	ENDIF
697
698	IF ( nyn == ny ) THEN
699	DO i = nxl, nxr
700	DO k = nzb_2d(ny,i) + 1, nzt
701	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
702	v(k,ny,i) * dzu(k)
703	volume_flow_area_l(2) = volume_flow_area_l(2) + dzu(k)
704	ENDDO
705	!
706	!-- Correction if velocity at nzb+1 has been set zero further above
707	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
708	v_nzb_p1_for_vfc(i)
709	ENDDO
710	ENDIF
711
712	#if defined( __parallel )
713	CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
714	2, MPI_REAL, MPI_SUM, comm2d, ierr )
715	CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &
716	2, MPI_REAL, MPI_SUM, comm2d, ierr )
717	#else
718	volume_flow_initial = volume_flow_initial_l
719	volume_flow_area = volume_flow_area_l
720	#endif
721	ENDIF
722
723	!
724	!-- For the moment, perturbation pressure and vertical velocity are zero
725	p = 0.0; w = 0.0
726
727	!
728	!-- Initialize array sums (must be defined in first call of pres)
729	sums = 0.0
730
731	!
732	!-- Treating cloud physics, liquid water content and precipitation amount
733	!-- are zero at beginning of the simulation
734	IF ( cloud_physics ) THEN
735	ql = 0.0
736	IF ( precipitation ) precipitation_amount = 0.0
737	ENDIF
738
739	!
740	!-- Initialize spectra
741	IF ( dt_dosp /= 9999999.9 ) THEN
742	spectrum_x = 0.0
743	spectrum_y = 0.0
744	ENDIF
745
746	!
747	!-- Impose vortex with vertical axis on the initial velocity profile
748	IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 ) THEN
749	CALL init_rankine
750	ENDIF
751
752	!
753	!-- Impose temperature anomaly (advection test only)
754	IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 ) THEN
755	CALL init_pt_anomaly
756	ENDIF
757
758	!
759	!-- If required, change the surface temperature at the start of the 3D run
760	IF ( pt_surface_initial_change /= 0.0 ) THEN
761	pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
762	ENDIF
763
764	!
765	!-- If required, change the surface humidity/scalar at the start of the 3D
766	!-- run
767	IF ( ( humidity .OR. passive_scalar ) .AND. &
768	q_surface_initial_change /= 0.0 ) THEN
769	q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
770	ENDIF
771
772	!
773	!-- Initialize the random number generator (from numerical recipes)
774	CALL random_function_ini
775
776	!
777	!-- Impose random perturbation on the horizontal velocity field and then
778	!-- remove the divergences from the velocity field
779	IF ( create_disturbances ) THEN
780	CALL disturb_field( nzb_u_inner, tend, u )
781	CALL disturb_field( nzb_v_inner, tend, v )
782	n_sor = nsor_ini
783	CALL pres
784	n_sor = nsor
785	ENDIF
786
787	!
788	!-- Once again set the perturbation pressure explicitly to zero in order to
789	!-- assure that it does not generate any divergences in the first time step.
790	!-- At t=0 the velocity field is free of divergence (as constructed above).
791	!-- Divergences being created during a time step are not yet known and thus
792	!-- cannot be corrected during the time step yet.
793	p = 0.0
794
795	!
796	!-- Initialize old and new time levels.
797	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
798	e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
799	ELSE
800	te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
801	ENDIF
802	e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
803
804	IF ( humidity .OR. passive_scalar ) THEN
805	IF ( ASSOCIATED( q_m ) ) q_m = q
806	IF ( timestep_scheme(1:5) == 'runge' ) tq_m = 0.0
807	q_p = q
808	IF ( humidity .AND. ASSOCIATED( vpt_m ) ) vpt_m = vpt
809	ENDIF
810
811	IF ( ocean ) THEN
812	tsa_m = 0.0
813	sa_p = sa
814	ENDIF
815
816	!
817	!-- Initialize old timelevels needed for radiation boundary conditions
818	IF ( outflow_l ) THEN
819	u_m_l(:,:,:) = u(:,:,1:2)
820	v_m_l(:,:,:) = v(:,:,0:1)
821	w_m_l(:,:,:) = w(:,:,0:1)
822	ENDIF
823	IF ( outflow_r ) THEN
824	u_m_r(:,:,:) = u(:,:,nx-1:nx)
825	v_m_r(:,:,:) = v(:,:,nx-1:nx)
826	w_m_r(:,:,:) = w(:,:,nx-1:nx)
827	ENDIF
828	IF ( outflow_s ) THEN
829	u_m_s(:,:,:) = u(:,0:1,:)
830	v_m_s(:,:,:) = v(:,1:2,:)
831	w_m_s(:,:,:) = w(:,0:1,:)
832	ENDIF
833	IF ( outflow_n ) THEN
834	u_m_n(:,:,:) = u(:,ny-1:ny,:)
835	v_m_n(:,:,:) = v(:,ny-1:ny,:)
836	w_m_n(:,:,:) = w(:,ny-1:ny,:)
837	ENDIF
838
839	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
840	THEN
841	!
842	!-- Read binary data from restart file
843	CALL read_3d_binary
844
845	!
846	!-- Calculate initial temperature field and other constants used in case
847	!-- of a sloping surface
848	IF ( sloping_surface ) CALL init_slope
849
850	!
851	!-- Initialize new time levels (only done in order to set boundary values
852	!-- including ghost points)
853	e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
854	IF ( humidity .OR. passive_scalar ) q_p = q
855	IF ( ocean ) sa_p = sa
856
857	ELSE
858	!
859	!-- Actually this part of the programm should not be reached
860	IF ( myid == 0 ) PRINT*,'+++ init_3d_model: unknown initializing ', &
861	'problem'
862	CALL local_stop
863	ENDIF
864
865	!
866	!-- If required, initialize dvrp-software
867	IF ( dt_dvrp /= 9999999.9 ) CALL init_dvrp
868
869	IF ( ocean ) THEN
870	!
871	!-- Initialize quantities needed for the ocean model
872	CALL init_ocean
873	ELSE
874	!
875	!-- Initialize quantities for handling cloud physics
876	!-- This routine must be called before init_particles, because
877	!-- otherwise, array pt_d_t, needed in data_output_dvrp (called by
878	!-- init_particles) is not defined.
879	CALL init_cloud_physics
880	ENDIF
881
882	!
883	!-- If required, initialize particles
884	IF ( particle_advection ) CALL init_particles
885
886	!
887	!-- Initialize quantities for special advections schemes
888	CALL init_advec
889
890	!
891	!-- Initialize Rayleigh damping factors
892	rdf = 0.0
893	IF ( rayleigh_damping_factor /= 0.0 ) THEN
894	DO k = nzb+1, nzt
895	IF ( zu(k) >= rayleigh_damping_height ) THEN
896	rdf(k) = rayleigh_damping_factor * &
897	( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height ) &
898	/ ( zu(nzt) - rayleigh_damping_height ) )&
899	)**2
900	ENDIF
901	ENDDO
902	ENDIF
903
904	!
905	!-- Initialize diffusivities used within the outflow damping layer in case of
906	!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
907	!-- half of the width of the damping layer
908	IF ( bc_lr == 'dirichlet/radiation' ) THEN
909
910	DO i = nxl-1, nxr+1
911	IF ( i >= nx - outflow_damping_width ) THEN
912	km_damp_x(i) = km_damp_max * MIN( 1.0, &
913	( i - ( nx - outflow_damping_width ) ) / &
914	REAL( outflow_damping_width/2 ) &
915	)
916	ELSE
917	km_damp_x(i) = 0.0
918	ENDIF
919	ENDDO
920
921	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
922
923	DO i = nxl-1, nxr+1
924	IF ( i <= outflow_damping_width ) THEN
925	km_damp_x(i) = km_damp_max * MIN( 1.0, &
926	( outflow_damping_width - i ) / &
927	REAL( outflow_damping_width/2 ) &
928	)
929	ELSE
930	km_damp_x(i) = 0.0
931	ENDIF
932	ENDDO
933
934	ENDIF
935
936	IF ( bc_ns == 'dirichlet/radiation' ) THEN
937
938	DO j = nys-1, nyn+1
939	IF ( j >= ny - outflow_damping_width ) THEN
940	km_damp_y(j) = km_damp_max * MIN( 1.0, &
941	( j - ( ny - outflow_damping_width ) ) / &
942	REAL( outflow_damping_width/2 ) &
943	)
944	ELSE
945	km_damp_y(j) = 0.0
946	ENDIF
947	ENDDO
948
949	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
950
951	DO j = nys-1, nyn+1
952	IF ( j <= outflow_damping_width ) THEN
953	km_damp_y(j) = km_damp_max * MIN( 1.0, &
954	( outflow_damping_width - j ) / &
955	REAL( outflow_damping_width/2 ) &
956	)
957	ELSE
958	km_damp_y(j) = 0.0
959	ENDIF
960	ENDDO
961
962	ENDIF
963
964	!
965	!-- Initialize local summation arrays for UP flow_statistics. This is necessary
966	!-- because they may not yet have been initialized when they are called from
967	!-- flow_statistics (or - depending on the chosen model run - are never
968	!-- initialized)
969	sums_divnew_l = 0.0
970	sums_divold_l = 0.0
971	sums_l_l = 0.0
972	sums_up_fraction_l = 0.0
973	sums_wsts_bc_l = 0.0
974
975	!
976	!-- Pre-set masks for regional statistics. Default is the total model domain.
977	rmask = 1.0
978
979	!
980	!-- User-defined initializing actions. Check afterwards, if maximum number
981	!-- of allowed timeseries is not exceeded
982	CALL user_init
983
984	IF ( dots_num > dots_max ) THEN
985	IF ( myid == 0 ) THEN
986	PRINT*, '+++ user_init: number of time series quantities exceeds', &
987	' its maximum of dots_max = ', dots_max
988	PRINT*, ' Please increase dots_max in modules.f90.'
989	ENDIF
990	CALL local_stop
991	ENDIF
992
993	!
994	!-- Input binary data file is not needed anymore. This line must be placed
995	!-- after call of user_init!
996	CALL close_file( 13 )
997
998	!
999	!-- Compute total sum of active mask grid points
1000	!-- ngp_2dh: number of grid points of a horizontal cross section through the
1001	!-- total domain
1002	!-- ngp_3d: number of grid points of the total domain
1003	!-- Note: The lower vertical index nzb_s_outer imposes a small error on the 2D
1004	!-- ---- averages of staggered variables such as u and v due to the topography
1005	!-- arrangement on the staggered grid. Maybe revise later.
1006	ngp_2dh_outer_l = 0
1007	ngp_2dh_outer = 0
1008	ngp_2dh_l = 0
1009	ngp_2dh = 0
1010	ngp_3d_inner_l = 0
1011	ngp_3d_inner = 0
1012	ngp_3d = 0
1013	ngp_sums = ( nz + 2 ) * ( pr_palm + max_pr_user )
1014
1015	DO sr = 0, statistic_regions
1016	DO i = nxl, nxr
1017	DO j = nys, nyn
1018	IF ( rmask(j,i,sr) == 1.0 ) THEN
1019	!
1020	!-- All xy-grid points
1021	ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
1022	!
1023	!-- xy-grid points above topography
1024	DO k = nzb_s_outer(j,i), nz + 1
1025	ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
1026	ENDDO
1027	!
1028	!-- All grid points of the total domain above topography
1029	ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
1030	( nz - nzb_s_inner(j,i) + 2 )
1031	ENDIF
1032	ENDDO
1033	ENDDO
1034	ENDDO
1035
1036	sr = statistic_regions + 1
1037	#if defined( __parallel )
1038	CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, &
1039	comm2d, ierr )
1040	CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
1041	MPI_INTEGER, MPI_SUM, comm2d, ierr )
1042	CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner(0), sr, MPI_INTEGER, &
1043	MPI_SUM, comm2d, ierr )
1044	#else
1045	ngp_2dh = ngp_2dh_l
1046	ngp_2dh_outer = ngp_2dh_outer_l
1047	ngp_3d_inner = ngp_3d_inner_l
1048	#endif
1049
1050	ngp_3d = ngp_2dh * ( nz + 2 )
1051
1052	!
1053	!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
1054	!-- buoyancy, etc. A zero value will occur for cases where all grid points of
1055	!-- the respective subdomain lie below the surface topography
1056	ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) )
1057	ngp_3d_inner = MAX( 1, ngp_3d_inner(:) )
1058
1059	DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l )
1060
1061
1062	END SUBROUTINE init_3d_model

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |