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

Last change on this file since 108 was 108, checked in by letzel, 14 years ago
  • Improved coupler: evaporation - salinity-flux coupling for humidity = .T.,

avoid MPI hangs when coupled runs terminate, add DOC/app/chapter_3.8;

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