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

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

Id keyword set as property for all *.f90 files

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