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

Last change on this file since 707 was 707, checked in by raasch, 10 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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