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

Last change on this file since 392 was 392, checked in by raasch, 12 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

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