Home

Context Navigation

source: palm/trunk/SOURCE/time_integration.f90 @ 73

Last change on this file since 73 was 73, checked in by raasch, 17 years ago
preliminary changes for radiation conditions
Property svn:keywords set to `Id`
File size: 16.3 KB

Line
1	SUBROUTINE time_integration
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! Move call of user_actions( 'after_integration' ) below increment of times
7	! and counters,
8	! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
9	! .._vector, these calls are now controlled by switch loop_optimization
10	!
11	! Former revisions:
12	! -----------------
13	! $Id: time_integration.f90 73 2007-03-20 08:33:14Z raasch $
14	! RCS Log replace by Id keyword, revision history cleaned up
15	!
16	! Revision 1.8 2006/08/22 14:16:05 raasch
17	! Disturbances are imposed only for the last Runge-Kutta-substep
18	!
19	! Revision 1.2 2004/04/30 13:03:40 raasch
20	! decalpha-specific warning removed, routine name changed to time_integration,
21	! particle advection is carried out only once during the intermediate steps,
22	! impulse_advec renamed momentum_advec
23	!
24	! Revision 1.1 1997/08/11 06:19:04 raasch
25	! Initial revision
26	!
27	!
28	! Description:
29	! ------------
30	! Integration in time of the model equations, statistical analysis and graphic
31	! output
32	!------------------------------------------------------------------------------!
33
34	USE arrays_3d
35	USE averaging
36	USE control_parameters
37	USE cpulog
38	#if defined( __dvrp_graphics )
39	USE DVRP
40	#endif
41	USE grid_variables
42	USE indices
43	USE interaction_droplets_ptq_mod
44	USE interfaces
45	USE particle_attributes
46	USE pegrid
47	USE prognostic_equations_mod
48	USE statistics
49	USE user_actions_mod
50
51	IMPLICIT NONE
52
53	CHARACTER (LEN=9) :: time_to_string
54	INTEGER :: i, j, k
55
56	!
57	!-- At the beginning of a simulation determine the time step as well as
58	!-- determine and print out the run control parameters
59	IF ( simulated_time == 0.0 ) CALL timestep
60	CALL run_control
61
62	#if defined( __dvrp_graphics )
63	!
64	!-- Time measurement with dvrp software
65	CALL DVRP_LOG_EVENT( 2, current_timestep_number )
66	#endif
67
68	!
69	!-- Start of the time loop
70	DO WHILE ( simulated_time < end_time .AND. .NOT. stop_dt .AND. &
71	.NOT. terminate_run )
72
73	CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
74
75	!
76	!-- Determine size of next time step
77	IF ( simulated_time /= 0.0 ) CALL timestep
78
79	!
80	!-- Execute the user-defined actions
81	CALL user_actions( 'before_timestep' )
82
83	!
84	!-- Start of intermediate step loop
85	intermediate_timestep_count = 0
86	DO WHILE ( intermediate_timestep_count < &
87	intermediate_timestep_count_max )
88
89	intermediate_timestep_count = intermediate_timestep_count + 1
90
91	!
92	!-- Set the steering factors for the prognostic equations which depend
93	!-- on the timestep scheme
94	CALL timestep_scheme_steering
95
96	!
97	!-- Solve the prognostic equations. A fast cache optimized version with
98	!-- only one single loop is used in case of Piascek-Williams advection
99	!-- scheme. NEC vector machines use a different version, because
100	!-- in the other versions a good vectorization is prohibited due to
101	!-- inlining problems.
102	IF ( loop_optimization == 'vector' ) THEN
103	CALL prognostic_equations_vector
104	ELSE
105	IF ( momentum_advec == 'ups-scheme' .OR. &
106	scalar_advec == 'ups-scheme' .OR. &
107	scalar_advec == 'bc-scheme' ) &
108	THEN
109	CALL prognostic_equations_noopt
110	ELSE
111	CALL prognostic_equations_cache
112	ENDIF
113	ENDIF
114
115	!
116	!-- Particle advection (only once during intermediate steps, because
117	!-- it uses an Euler-step)
118	IF ( particle_advection .AND. &
119	simulated_time >= particle_advection_start .AND. &
120	intermediate_timestep_count == 1 ) THEN
121	CALL advec_particles
122	first_call_advec_particles = .FALSE.
123	ENDIF
124
125	!
126	!-- Interaction of droplets with temperature and specific humidity.
127	!-- Droplet condensation and evaporation is calculated within
128	!-- advec_particles.
129	IF ( cloud_droplets .AND. &
130	intermediate_timestep_count == intermediate_timestep_count_max )&
131	THEN
132	CALL interaction_droplets_ptq
133	ENDIF
134
135	!
136	!-- Exchange of ghost points (lateral boundary conditions)
137	CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
138	CALL exchange_horiz( u_p, uxrp, 0 )
139	CALL exchange_horiz( v_p, 0, vynp )
140	CALL exchange_horiz( w_p, 0, 0 )
141	CALL exchange_horiz( pt_p, 0, 0 )
142	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e_p, 0, 0)
143	IF ( moisture .OR. passive_scalar ) CALL exchange_horiz( q_p, 0, 0)
144	IF ( cloud_droplets ) THEN
145	CALL exchange_horiz( ql, 0, 0 )
146	CALL exchange_horiz( ql_c, 0, 0 )
147	CALL exchange_horiz( ql_v, 0, 0 )
148	CALL exchange_horiz( ql_vp, 0, 0 )
149	ENDIF
150
151	CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
152
153	!
154	!-- Apply time filter in case of leap-frog timestep
155	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
156	CALL asselin_filter
157	ENDIF
158
159	!
160	!-- Boundary conditions for the prognostic quantities (except of the
161	!-- velocities at the outflow in case of a non-cyclic lateral wall)
162	CALL boundary_conds( 'main' )
163
164	!
165	!-- Swap the time levels in preparation for the next time step.
166	CALL swap_timelevel
167
168	!
169	!-- Temperature offset must be imposed at cyclic boundaries in x-direction
170	!-- when a sloping surface is used
171	IF ( sloping_surface ) THEN
172	IF ( nxl == 0 ) pt(:,:,nxl-1) = pt(:,:,nxl-1) - pt_slope_offset
173	IF ( nxr == nx ) pt(:,:,nxr+1) = pt(:,:,nxr+1) + pt_slope_offset
174	ENDIF
175
176	!
177	!-- Impose a random perturbation on the horizontal velocity field
178	IF ( create_disturbances .AND. &
179	intermediate_timestep_count == intermediate_timestep_count_max )&
180	THEN
181	time_disturb = time_disturb + dt_3d
182	IF ( time_disturb >= dt_disturb ) THEN
183	IF ( hom(nzb+5,1,var_hom,0) < disturbance_energy_limit ) THEN
184	CALL disturb_field( nzb_u_inner, tend, u, uxrp, 0 )
185	CALL disturb_field( nzb_v_inner, tend, v , 0, vynp )
186	ELSEIF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN
187	!
188	!-- Runs with a non-cyclic lateral wall need perturbations
189	!-- near the inflow throughout the whole simulation
190	dist_range = 1
191	CALL disturb_field( nzb_u_inner, tend, u, uxrp, 0 )
192	CALL disturb_field( nzb_v_inner, tend, v , 0, vynp )
193	dist_range = 0
194	ENDIF
195	time_disturb = time_disturb - dt_disturb
196	ENDIF
197	ENDIF
198
199	!
200	!-- Reduce the velocity divergence via the equation for perturbation
201	!-- pressure.
202	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
203	.OR. call_psolver_at_all_substeps ) THEN
204	CALL pres
205	ENDIF
206
207	!
208	!-- In case of a non-cyclic lateral wall, set the boundary conditions for
209	!-- the velocities at the outflow
210	IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN
211	CALL boundary_conds( 'outflow_uvw' )
212	ENDIF
213
214	!
215	!-- If required, compute virtuell potential temperature
216	IF ( moisture ) CALL compute_vpt
217
218	!
219	!-- If required, compute liquid water content
220	IF ( cloud_physics ) CALL calc_liquid_water_content
221
222	!
223	!-- Compute the diffusion quantities
224	IF ( .NOT. constant_diffusion ) THEN
225
226	!
227	!-- First the vertical fluxes in the Prandtl layer are being computed
228	IF ( prandtl_layer ) THEN
229	CALL cpu_log( log_point(19), 'prandtl_fluxes', 'start' )
230	CALL prandtl_fluxes
231	CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
232	ENDIF
233
234	!
235	!-- Compute the diffusion coefficients
236	CALL cpu_log( log_point(17), 'diffusivities', 'start' )
237	IF ( .NOT. moisture ) THEN
238	CALL diffusivities( pt )
239	ELSE
240	CALL diffusivities( vpt )
241	ENDIF
242	CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
243
244	ENDIF
245
246	ENDDO ! Intermediate step loop
247
248	!
249	!-- Increase simulation time and output times
250	current_timestep_number = current_timestep_number + 1
251	simulated_time = simulated_time + dt_3d
252	simulated_time_chr = time_to_string( simulated_time )
253	IF ( simulated_time >= skip_time_data_output_av ) THEN
254	time_do_av = time_do_av + dt_3d
255	ENDIF
256	IF ( simulated_time >= skip_time_do2d_xy ) THEN
257	time_do2d_xy = time_do2d_xy + dt_3d
258	ENDIF
259	IF ( simulated_time >= skip_time_do2d_xz ) THEN
260	time_do2d_xz = time_do2d_xz + dt_3d
261	ENDIF
262	IF ( simulated_time >= skip_time_do2d_yz ) THEN
263	time_do2d_yz = time_do2d_yz + dt_3d
264	ENDIF
265	IF ( simulated_time >= skip_time_do3d ) THEN
266	time_do3d = time_do3d + dt_3d
267	ENDIF
268	time_dvrp = time_dvrp + dt_3d
269	IF ( simulated_time >= skip_time_dosp ) THEN
270	time_dosp = time_dosp + dt_3d
271	ENDIF
272	time_dots = time_dots + dt_3d
273	IF ( .NOT. first_call_advec_particles ) THEN
274	time_dopts = time_dopts + dt_3d
275	ENDIF
276	IF ( simulated_time >= skip_time_dopr ) THEN
277	time_dopr = time_dopr + dt_3d
278	ENDIF
279	time_dopr_listing = time_dopr_listing + dt_3d
280	time_run_control = time_run_control + dt_3d
281
282	!
283	!-- Execute user-defined actions
284	CALL user_actions( 'after_integration' )
285
286	!
287	!-- If Galilei transformation is used, determine the distance that the
288	!-- model has moved so far
289	IF ( galilei_transformation ) THEN
290	advected_distance_x = advected_distance_x + u_gtrans * dt_3d
291	advected_distance_y = advected_distance_y + v_gtrans * dt_3d
292	ENDIF
293
294	!
295	!-- Check, if restart is necessary (because cpu-time is expiring or
296	!-- because it is forced by user) and set stop flag
297	CALL check_for_restart
298
299	!
300	!-- Carry out statistical analysis and output at the requested output times.
301	!-- The MOD function is used for calculating the output time counters (like
302	!-- time_dopr) in order to regard a possible decrease of the output time
303	!-- interval in case of restart runs
304
305	!
306	!-- Set a flag indicating that so far no statistics have been created
307	!-- for this time step
308	flow_statistics_called = .FALSE.
309
310	!
311	!-- If required, call flow_statistics for averaging in time
312	IF ( averaging_interval_pr /= 0.0 .AND. &
313	( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. &
314	simulated_time >= skip_time_dopr ) THEN
315	time_dopr_av = time_dopr_av + dt_3d
316	IF ( time_dopr_av >= dt_averaging_input_pr ) THEN
317	do_sum = .TRUE.
318	time_dopr_av = MOD( time_dopr_av, &
319	MAX( dt_averaging_input_pr, dt_3d ) )
320	ENDIF
321	ENDIF
322	IF ( do_sum ) CALL flow_statistics
323
324	!
325	!-- Sum-up 3d-arrays for later output of time-averaged data
326	IF ( averaging_interval /= 0.0 .AND. &
327	( dt_data_output_av - time_do_av ) <= averaging_interval .AND. &
328	simulated_time >= skip_time_data_output_av ) &
329	THEN
330	time_do_sla = time_do_sla + dt_3d
331	IF ( time_do_sla >= dt_averaging_input ) THEN
332	CALL sum_up_3d_data
333	average_count_3d = average_count_3d + 1
334	time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
335	ENDIF
336	ENDIF
337
338	!
339	!-- Calculate spectra for time averaging
340	IF ( averaging_interval_sp /= 0.0 .AND. &
341	( dt_dosp - time_dosp ) <= averaging_interval_sp .AND. &
342	simulated_time >= skip_time_dosp ) THEN
343	time_dosp_av = time_dosp_av + dt_3d
344	IF ( time_dosp_av >= dt_averaging_input_pr ) THEN
345	CALL calc_spectra
346	time_dosp_av = MOD( time_dosp_av, &
347	MAX( dt_averaging_input_pr, dt_3d ) )
348	ENDIF
349	ENDIF
350
351	!
352	!-- Computation and output of run control parameters.
353	!-- This is also done whenever the time step has changed or perturbations
354	!-- have been imposed
355	IF ( time_run_control >= dt_run_control .OR. &
356	( dt_changed .AND. timestep_scheme(1:5) /= 'runge' ) .OR. &
357	disturbance_created ) &
358	THEN
359	CALL run_control
360	IF ( time_run_control >= dt_run_control ) THEN
361	time_run_control = MOD( time_run_control, &
362	MAX( dt_run_control, dt_3d ) )
363	ENDIF
364	ENDIF
365
366	!
367	!-- Profile output (ASCII) on file
368	IF ( time_dopr_listing >= dt_dopr_listing ) THEN
369	CALL print_1d
370	time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
371	dt_3d ) )
372	ENDIF
373
374	!
375	!-- Graphic output for PROFIL
376	IF ( time_dopr >= dt_dopr ) THEN
377	IF ( dopr_n /= 0 ) CALL data_output_profiles
378	time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
379	time_dopr_av = 0.0 ! due to averaging (see above)
380	ENDIF
381
382	!
383	!-- Graphic output for time series
384	IF ( time_dots >= dt_dots ) THEN
385	CALL data_output_tseries
386	time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
387	ENDIF
388
389	!
390	!-- Output of spectra (formatted for use with PROFIL), in case of no
391	!-- time averaging, spectra has to be calculated before
392	IF ( time_dosp >= dt_dosp ) THEN
393	IF ( average_count_sp == 0 ) CALL calc_spectra
394	CALL data_output_spectra
395	time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
396	ENDIF
397
398	!
399	!-- 2d-data output (cross-sections)
400	IF ( time_do2d_xy >= dt_do2d_xy ) THEN
401	CALL data_output_2d( 'xy', 0 )
402	time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
403	ENDIF
404	IF ( time_do2d_xz >= dt_do2d_xz ) THEN
405	CALL data_output_2d( 'xz', 0 )
406	time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
407	ENDIF
408	IF ( time_do2d_yz >= dt_do2d_yz ) THEN
409	CALL data_output_2d( 'yz', 0 )
410	time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
411	ENDIF
412
413	!
414	!-- 3d-data output (volume data)
415	IF ( time_do3d >= dt_do3d ) THEN
416	CALL data_output_3d( 0 )
417	time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
418	ENDIF
419
420	!
421	!-- Output of time-averaged 2d/3d-data
422	IF ( time_do_av >= dt_data_output_av ) THEN
423	CALL average_3d_data
424	CALL data_output_2d( 'xy', 1 )
425	CALL data_output_2d( 'xz', 1 )
426	CALL data_output_2d( 'yz', 1 )
427	CALL data_output_3d( 1 )
428	time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
429	ENDIF
430
431	!
432	!-- Output of particle time series
433	IF ( time_dopts >= dt_dopts .OR. &
434	( simulated_time >= particle_advection_start .AND. &
435	first_call_advec_particles ) ) THEN
436	CALL data_output_ptseries
437	time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
438	ENDIF
439
440	!
441	!-- Output of dvrp-graphics (isosurface, particles, slicer)
442	#if defined( __dvrp_graphics )
443	CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
444	#endif
445	IF ( time_dvrp >= dt_dvrp ) THEN
446	CALL data_output_dvrp
447	time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
448	ENDIF
449	#if defined( __dvrp_graphics )
450	CALL DVRP_LOG_EVENT( 2, current_timestep_number )
451	#endif
452
453	!
454	!-- If required, set the heat flux for the next time step at a random value
455	IF ( constant_heatflux .AND. random_heatflux ) CALL disturb_heatflux
456
457	!
458	!-- Execute user-defined actions
459	CALL user_actions( 'after_timestep' )
460
461	CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
462
463	ENDDO ! time loop
464
465	#if defined( __dvrp_graphics )
466	CALL DVRP_LOG_EVENT( -2, current_timestep_number )
467	#endif
468
469	END SUBROUTINE time_integration

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |