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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

File size: 20.1 KB
Line 
1 SUBROUTINE time_integration
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: time_integration.f90,v $
11! Revision 1.8  2006/08/22 14:16:05  raasch
12! Disturbances are imposed only for the last Runge-Kutta-substep
13!
14! Revision 1.7  2006/08/04 15:04:39  raasch
15! Output of particle timeseries
16!
17! Revision 1.6  2006/02/23 12:58:38  raasch
18! Data output (1d, 2d, 3d, sp) can be skipped for a given time interval from
19! simulation start, additional argument nzb_.._inner added to call of
20! disturb_field, nt_anz renamed current_timestep_number, pl.. renamed do..,
21! ..anz renamed ..n, dt_average renamed dt_averaging_interval_pr,
22! time_average_.. renamed time_do.._av, routines plot_.. renamed data_output_..,
23! additional call of user_actions (after_integration)
24!
25! Revision 1.5  2005/06/26 20:10:20  raasch
26! +call of interaction_droplets_ptq (+module of the same name)
27!
28! Revision 1.4  2005/04/23 09:39:01  raasch
29! Revised calculation of output time counters regarding a possible decrease of
30! the output time interval in case of restart runs
31!
32! Revision 1.3  2005/03/26 21:10:26  raasch
33! Pres is called by default only at the last Runge-Kutta-substep.
34! Determination of the timestep moved from end to the beginning of the
35! time loop in order to assure consistent output in run_control (first line
36! in restart runs must be equal to the last line in the corresponding
37! previous run). Exchange of ghost points for km and kh removed.
38! Changes regarding non-cyclic boundary conditions:
39! arguments are added to call of routines disturb_field, exchange_horiz,
40! and boundary_conds, which is now called for a second time after pres in
41! order to set the velocity boundary conditions at the outflow. Velocity
42! perturbations are imposed near the inflow throughout the whole run.
43!
44! Revision 1.2  2004/04/30 13:03:40  raasch
45! decalpha-specific warning removed, routine name changed to time_integration,
46! particle advection is carried out only once during the intermediate steps,
47! impulse_advec renamed momentum_advec
48!
49! Revision 1.1  2004/04/30 12:58:21  raasch
50! Initial revision
51!
52! Revision 1.44  2004/01/28 15:17:26  raasch
53! Additional intermediate timesteps for the Runge-Kutta schemes, setting
54! of the steering factors for prognostic equations within the new
55! routine timestep_scheme_steering at the beginning of each step (former
56! part of routine timestep called at the end of each step), routine
57! run_control is not called for the Runge-Kutta schemes automatically in case
58! of timestep changes.
59!
60! Revision 1.43  2003/03/16 09:40:12  raasch
61! Two underscores (_) are placed in front of all define-strings
62!
63! Revision 1.42  2003/03/12 16:33:23  raasch
64! New routine prognostic_equations_vec is used on NEC machines
65!
66! Revision 1.41  2003/03/04 11:31:39  raasch
67! Module DVRP added in case that dvrp graphics package is switched on
68!
69! Revision 1.40  2002/12/19 15:45:13  raasch
70! Routine check_cpu_time renamed to check_for_restart
71!
72! Revision 1.39  2002/06/11 13:11:44  raasch
73! +USE user_actions_mod, prognostic_equations_mod,
74! temperature offset at cyclic boundaries in case of sloping surface moved
75! from prognostic_equations to this subroutine. Particle advection moved from
76! prognostic_equations to this routine. A fast version of prognostic_equations
77! is called in case that the Piascek-Williams advection scheme is used.
78!
79! Revision 1.38  2002/05/02 18:52:19  raasch
80! Application of time filter and exchange of ghost points moved from
81! prognostic_equations to this routine.
82! Error in calculating output time intervals for various quantities removed.
83!
84! Revision 1.37  2001/08/21  09:51:26  09:51:26  raasch (Siegfried Raasch)
85! Time measurements for user_actions moved to that subroutine
86!
87! Revision 1.36  2001/03/30 07:34:24  raasch
88! Timelevel t+dt eliminated, Asselin filter now included in
89! prognostic_equations, arguments removed from pres,
90! Translation of remaining German identifiers (variables, subroutines, etc.)
91!
92! Revision 1.35  2001/01/25 07:06:43  raasch
93! Time measurements with dvrp software included
94!
95! Revision 1.34  2001/01/05 15:12:50  raasch
96! Calculation and plot of spectra
97!
98! Revision 1.33  2000/04/27 06:40:49  raasch
99! Old revision remarks deleted, output of isosurfaces and particles now
100! combined into new routine data_output_dvrp (includes former routines
101! plot_particles and plot_isosurface)
102!
103! Revision 1.32  2000/04/13 13:39:30  schroeter
104! computation of prognostic equations moved to the new routine
105! prognostic_equation, new routines calc_liquid_water_content and
106! compute_vpt
107!
108! Revision 1.31  2000/01/25  16:25:24  16:25:24  letzel (Marcus Letzel)
109! All comments translated into English
110!
111! Revision 1.1  1997/08/11 06:19:04  raasch
112! Initial revision
113!
114!
115! Description:
116! ------------
117! Integration in time of the model equations, statistical analysis and graphic
118! output
119!------------------------------------------------------------------------------!
120
121    USE arrays_3d
122    USE averaging
123    USE control_parameters
124    USE cpulog
125#if defined( __dvrp_graphics )
126    USE DVRP
127#endif
128    USE grid_variables
129    USE indices
130    USE interaction_droplets_ptq_mod
131    USE interfaces
132    USE particle_attributes
133    USE pegrid
134    USE prognostic_equations_mod
135    USE statistics
136    USE user_actions_mod
137
138    IMPLICIT NONE
139
140    CHARACTER (LEN=9) ::  time_to_string
141    INTEGER ::  i, j, k
142
143!
144!-- At the beginning of a simulation determine the time step as well as
145!-- determine and print out the run control parameters
146    IF ( simulated_time == 0.0 )  CALL timestep
147    CALL run_control
148
149#if defined( __dvrp_graphics )
150!
151!-- Time measurement with dvrp software 
152    CALL DVRP_LOG_EVENT( 2, current_timestep_number )
153#endif
154
155!
156!-- Start of the time loop
157    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. &
158                .NOT. terminate_run )
159
160       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
161
162!
163!--    Determine size of next time step
164       IF ( simulated_time /= 0.0 )  CALL timestep
165
166!
167!--    Execute the user-defined actions
168       CALL user_actions( 'before_timestep' )
169
170!
171!--    Start of intermediate step loop
172       intermediate_timestep_count = 0
173       DO  WHILE ( intermediate_timestep_count < &
174                   intermediate_timestep_count_max )
175
176          intermediate_timestep_count = intermediate_timestep_count + 1
177
178!
179!--       Set the steering factors for the prognostic equations which depend
180!--       on the timestep scheme
181          CALL timestep_scheme_steering
182
183!
184!--       Solve the prognostic equations. A fast cache optimized version with
185!--       only one single loop is used in case of Piascek-Williams advection
186!--       scheme. NEC vector machines use a different version, because
187!--       in the other versions a good vectorization is prohibited due to
188!--       inlining problems.
189          IF ( host(1:3) == 'nec' )  THEN
190             CALL prognostic_equations_vec
191          ELSE
192             IF ( momentum_advec == 'ups-scheme'  .OR.  &
193                  scalar_advec == 'ups-scheme'   .OR.  &
194                  scalar_advec == 'bc-scheme' )        &
195             THEN
196                CALL prognostic_equations
197             ELSE
198                CALL prognostic_equations_fast
199             ENDIF
200          ENDIF
201
202!
203!--       Particle advection (only once during intermediate steps, because
204!--       it uses an Euler-step)
205          IF ( simulated_time >= particle_advection_start  .AND. &
206               intermediate_timestep_count == 1 )  THEN
207             CALL advec_particles
208             first_call_advec_particles = .FALSE.
209          ENDIF
210
211!
212!--       Interaction of droplets with temperature and specific humidity.
213!--       Droplet condensation and evaporation is calculated within
214!--       advec_particles.
215          IF ( cloud_droplets  .AND.  &
216               intermediate_timestep_count == intermediate_timestep_count_max )&
217          THEN
218             CALL interaction_droplets_ptq
219          ENDIF
220
221!
222!--       Exchange of ghost points (lateral boundary conditions)
223          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
224          CALL exchange_horiz( u_p, uxrp,    0 )
225          CALL exchange_horiz( v_p,    0, vynp )
226          CALL exchange_horiz( w_p,    0,    0 )
227          CALL exchange_horiz( pt_p,   0,    0 )
228          IF ( .NOT. constant_diffusion       )  CALL exchange_horiz( e_p, 0, 0)
229          IF ( moisture  .OR.  passive_scalar )  CALL exchange_horiz( q_p, 0, 0)
230          IF ( cloud_droplets )  THEN
231             CALL exchange_horiz( ql,    0, 0 )
232             CALL exchange_horiz( ql_c,  0, 0 )
233             CALL exchange_horiz( ql_v,  0, 0 )
234             CALL exchange_horiz( ql_vp, 0, 0 )
235          ENDIF
236
237          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
238
239!
240!--       Apply time filter in case of leap-frog timestep
241          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
242             CALL asselin_filter
243          ENDIF
244
245!
246!--       Swap the time levels in preparation for the next time step.
247          CALL swap_timelevel
248
249!
250!--       Boundary conditions for the prognostic quantities (except of the
251!--       velocities at the outflow in case of a non-cyclic lateral wall)
252          CALL boundary_conds( 'main' )
253
254!
255!--       Temperature offset must be imposed at cyclic boundaries in x-direction
256!--       when a sloping surface is used
257          IF ( sloping_surface )  THEN
258             IF ( nxl ==  0 )  pt(:,:,nxl-1) = pt(:,:,nxl-1) - pt_slope_offset
259             IF ( nxr == nx )  pt(:,:,nxr+1) = pt(:,:,nxr+1) + pt_slope_offset
260          ENDIF
261
262!
263!--       Impose a random perturbation on the horizontal velocity field
264          IF ( create_disturbances  .AND.  &
265               intermediate_timestep_count == intermediate_timestep_count_max )&
266          THEN
267             time_disturb = time_disturb + dt_3d
268             IF ( time_disturb >= dt_disturb )  THEN
269                IF ( hom(nzb+5,1,var_hom,0) < disturbance_energy_limit )  THEN
270                   CALL disturb_field( nzb_u_inner, tend, u, uxrp, 0    )
271                   CALL disturb_field( nzb_v_inner, tend, v ,   0, vynp )
272                ELSEIF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
273!
274!--                Runs with a non-cyclic lateral wall need perturbations
275!--                near the inflow throughout the whole simulation
276                   dist_range = 1
277                   CALL disturb_field( nzb_u_inner, tend, u, uxrp, 0    )
278                   CALL disturb_field( nzb_v_inner, tend, v ,   0, vynp )
279                   dist_range = 0
280                ENDIF
281                time_disturb = time_disturb - dt_disturb
282             ENDIF
283          ENDIF
284
285!
286!--       Reduce the velocity divergence via the equation for perturbation
287!--       pressure.
288          IF ( intermediate_timestep_count == intermediate_timestep_count_max &
289               .OR.  call_psolver_at_all_substeps )  THEN
290             CALL pres
291          ENDIF
292
293!
294!--       In case of a non-cyclic lateral wall, set the boundary conditions for
295!--       the velocities at the outflow
296          IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
297             CALL boundary_conds( 'outflow_uvw' )
298          ENDIF
299
300!
301!--       If required, compute virtuell potential temperature
302          IF ( moisture ) CALL compute_vpt
303
304!
305!--       If required, compute liquid water content
306          IF ( cloud_physics ) CALL calc_liquid_water_content
307
308!
309!--       Compute the diffusion quantities
310          IF ( .NOT. constant_diffusion )  THEN
311
312!
313!--          First the vertical fluxes in the Prandtl layer are being computed
314             IF ( prandtl_layer )  THEN
315                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'start' )
316                CALL prandtl_fluxes
317                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
318             ENDIF
319
320!
321!--          Compute the diffusion coefficients
322             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
323             IF ( .NOT. moisture ) THEN
324                CALL diffusivities( pt )
325             ELSE
326                CALL diffusivities( vpt )
327             ENDIF
328             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
329
330          ENDIF
331
332       ENDDO   ! Intermediate step loop
333
334!
335!--    Execute user-defined actions
336       CALL user_actions( 'after_integration' )
337
338!
339!--    Increase simulation time and output times
340       current_timestep_number = current_timestep_number + 1
341       simulated_time     = simulated_time   + dt_3d
342       simulated_time_chr = time_to_string( simulated_time )
343       IF ( simulated_time >= skip_time_data_output_av )  THEN
344          time_do_av         = time_do_av       + dt_3d
345       ENDIF
346       IF ( simulated_time >= skip_time_do2d_xy )  THEN
347          time_do2d_xy       = time_do2d_xy     + dt_3d
348       ENDIF
349       IF ( simulated_time >= skip_time_do2d_xz )  THEN
350          time_do2d_xz       = time_do2d_xz     + dt_3d
351       ENDIF
352       IF ( simulated_time >= skip_time_do2d_yz )  THEN
353          time_do2d_yz       = time_do2d_yz     + dt_3d
354       ENDIF
355       IF ( simulated_time >= skip_time_do3d    )  THEN
356          time_do3d          = time_do3d        + dt_3d
357       ENDIF
358       time_dvrp          = time_dvrp        + dt_3d
359       IF ( simulated_time >= skip_time_dosp )  THEN
360          time_dosp       = time_dosp        + dt_3d
361       ENDIF
362       time_dots          = time_dots        + dt_3d
363       IF ( .NOT. first_call_advec_particles )  THEN
364          time_dopts      = time_dopts       + dt_3d
365       ENDIF
366       IF ( simulated_time >= skip_time_dopr )  THEN
367          time_dopr       = time_dopr        + dt_3d
368       ENDIF
369       time_dopr_listing          = time_dopr_listing        + dt_3d
370       time_run_control   = time_run_control + dt_3d
371
372!
373!--    If Galilei transformation is used, determine the distance that the
374!--    model has moved so far
375       IF ( galilei_transformation )  THEN
376          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
377          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
378       ENDIF
379
380!
381!--    Check, if restart is necessary (because cpu-time is expiring or
382!--    because it is forced by user) and set stop flag
383       CALL check_for_restart
384
385!
386!--    Carry out statistical analysis and output at the requested output times.
387!--    The MOD function is used for calculating the output time counters (like
388!--    time_dopr) in order to regard a possible decrease of the output time
389!--    interval in case of restart runs
390
391!
392!--    Set a flag indicating that so far no statistics have been created
393!--    for this time step
394       flow_statistics_called = .FALSE.
395
396!
397!--    If required, call flow_statistics for averaging in time
398       IF ( averaging_interval_pr /= 0.0  .AND.  &
399            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
400            simulated_time >= skip_time_dopr )  THEN
401          time_dopr_av = time_dopr_av + dt_3d
402          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
403             do_sum = .TRUE.
404             time_dopr_av = MOD( time_dopr_av, &
405                                    MAX( dt_averaging_input_pr, dt_3d ) )
406          ENDIF
407       ENDIF
408       IF ( do_sum )  CALL flow_statistics
409
410!
411!--    Sum-up 3d-arrays for later output of time-averaged data
412       IF ( averaging_interval /= 0.0  .AND.                                &
413            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
414            simulated_time >= skip_time_data_output_av )                    &
415       THEN
416          time_do_sla = time_do_sla + dt_3d
417          IF ( time_do_sla >= dt_averaging_input )  THEN
418             CALL sum_up_3d_data
419             average_count_3d = average_count_3d + 1
420             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
421          ENDIF
422       ENDIF
423
424!
425!--    Calculate spectra for time averaging
426       IF ( averaging_interval_sp /= 0.0  .AND.  &
427            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
428            simulated_time >= skip_time_dosp )  THEN
429          time_dosp_av = time_dosp_av + dt_3d
430          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
431             CALL calc_spectra
432             time_dosp_av = MOD( time_dosp_av, &
433                                 MAX( dt_averaging_input_pr, dt_3d ) )
434          ENDIF
435       ENDIF
436
437!
438!--    Computation and output of run control parameters.
439!--    This is also done whenever the time step has changed or perturbations
440!--    have been imposed
441       IF ( time_run_control >= dt_run_control  .OR.                     &
442            ( dt_changed  .AND.  timestep_scheme(1:5) /= 'runge' )  .OR. &
443            disturbance_created )                                        &
444       THEN
445          CALL run_control
446          IF ( time_run_control >= dt_run_control )  THEN
447             time_run_control = MOD( time_run_control, &
448                                     MAX( dt_run_control, dt_3d ) )
449          ENDIF
450       ENDIF
451
452!
453!--    Profile output (ASCII) on file
454       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
455          CALL print_1d
456          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
457                                                           dt_3d ) )
458       ENDIF
459
460!
461!--    Graphic output for PROFIL
462       IF ( time_dopr >= dt_dopr )  THEN
463          IF ( dopr_n /= 0 )  CALL data_output_profiles
464          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
465          time_dopr_av = 0.0    ! due to averaging (see above)
466       ENDIF
467
468!
469!--    Graphic output for time series
470       IF ( time_dots >= dt_dots )  THEN
471          IF ( dots_n /= 0 )  CALL data_output_tseries
472          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
473       ENDIF
474
475!
476!--    Output of spectra (formatted for use with PROFIL), in case of no
477!--    time averaging, spectra has to be calculated before
478       IF ( time_dosp >= dt_dosp )  THEN
479          IF ( average_count_sp == 0 )  CALL calc_spectra
480          CALL data_output_spectra
481          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
482       ENDIF
483
484!
485!--    2d-data output (cross-sections)
486       IF ( time_do2d_xy >= dt_do2d_xy )  THEN
487          CALL data_output_2d( 'xy', 0 )
488          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
489       ENDIF
490       IF ( time_do2d_xz >= dt_do2d_xz )  THEN
491          CALL data_output_2d( 'xz', 0 )
492          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
493       ENDIF
494       IF ( time_do2d_yz >= dt_do2d_yz )  THEN
495          CALL data_output_2d( 'yz', 0 )
496          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
497       ENDIF
498
499!
500!--    3d-data output (volume data)
501       IF ( time_do3d >= dt_do3d )  THEN
502          CALL data_output_3d( 0 )
503          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
504       ENDIF
505
506!
507!--    Output of time-averaged 2d/3d-data
508       IF ( time_do_av >= dt_data_output_av )  THEN
509          CALL average_3d_data
510          CALL data_output_2d( 'xy', 1 )
511          CALL data_output_2d( 'xz', 1 )
512          CALL data_output_2d( 'yz', 1 )
513          CALL data_output_3d( 1 )
514          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
515       ENDIF
516
517!
518!--    Output of particle time series
519       IF ( time_dopts >= dt_dopts  .OR. &
520            ( simulated_time >= particle_advection_start  .AND. &
521              first_call_advec_particles ) )  THEN
522          CALL data_output_ptseries
523          time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
524       ENDIF
525
526!
527!--    Output of dvrp-graphics (isosurface, particles, slicer)
528#if defined( __dvrp_graphics )
529       CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
530#endif
531       IF ( time_dvrp >= dt_dvrp )  THEN
532          CALL data_output_dvrp
533          time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
534       ENDIF
535#if defined( __dvrp_graphics )
536       CALL DVRP_LOG_EVENT( 2, current_timestep_number )
537#endif
538
539!
540!--    If required, set the heat flux for the next time step at a random value
541       IF ( constant_heatflux  .AND.  random_heatflux )  CALL disturb_heatflux
542
543!
544!--    Execute user-defined actions
545       CALL user_actions( 'after_timestep' )
546
547       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
548
549    ENDDO   ! time loop
550
551#if defined( __dvrp_graphics )
552    CALL DVRP_LOG_EVENT( -2, current_timestep_number )
553#endif
554
555 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.