Home

Context Navigation

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

Last change on this file since 16 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 28.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |