Home

Context Navigation

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

Last change on this file since 70 was 63, checked in by raasch, 18 years ago
preliminary changes concerning update of BC-scheme
Property svn:keywords set to `Id`
File size: 16.3 KB

Rev	Line
[1]	1	SUBROUTINE time_integration
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[46]	6	! Move call of user_actions( 'after_integration' ) below increment of times
[63]	7	! and counters,
	8	! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
	9	! .._vector, these calls are now controlled by switch loop_optimization
[1]	10	!
	11	! Former revisions:
	12	! -----------------
[3]	13	! $Id: time_integration.f90 63 2007-03-13 03:52:49Z raasch $
	14	! RCS Log replace by Id keyword, revision history cleaned up
	15	!
[1]	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.
[63]	102	IF ( loop_optimization == 'vector' ) THEN
	103	CALL prognostic_equations_vector
[1]	104	ELSE
	105	IF ( momentum_advec == 'ups-scheme' .OR. &
	106	scalar_advec == 'ups-scheme' .OR. &
	107	scalar_advec == 'bc-scheme' ) &
	108	THEN
[63]	109	CALL prognostic_equations_noopt
[1]	110	ELSE
[63]	111	CALL prognostic_equations_cache
[1]	112	ENDIF
	113	ENDIF
	114
	115	!
	116	!-- Particle advection (only once during intermediate steps, because
	117	!-- it uses an Euler-step)
[63]	118	IF ( particle_advection .AND. &
	119	simulated_time >= particle_advection_start .AND. &
[1]	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	!-- Swap the time levels in preparation for the next time step.
	161	CALL swap_timelevel
	162
	163	!
	164	!-- Boundary conditions for the prognostic quantities (except of the
	165	!-- velocities at the outflow in case of a non-cyclic lateral wall)
	166	CALL boundary_conds( 'main' )
	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	!
[46]	283	!-- Execute user-defined actions
	284	CALL user_actions( 'after_integration' )
	285
	286	!
[1]	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
[48]	385	CALL data_output_tseries
[1]	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 |