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

Last change on this file since 2020 was 2012, checked in by kanani, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 47.2 KB
Line 
1!> @file time_integration.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2016 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration.f90 2012 2016-09-19 17:31:38Z hellstea $
27!
28! 2011 2016-09-19 17:29:57Z kanani
29! Flag urban_surface is now defined in module control_parameters,
30! removed commented CALLs of global_min_max.
31!
32! 2007 2016-08-24 15:47:17Z kanani
33! Added CALLs for new urban surface model
34!
35! 2000 2016-08-20 18:09:15Z knoop
36! Forced header and separation lines into 80 columns
37!
38! 1976 2016-07-27 13:28:04Z maronga
39! Simplified calls to radiation model
40!
41! 1960 2016-07-12 16:34:24Z suehring
42! Separate humidity and passive scalar
43!
44! 1957 2016-07-07 10:43:48Z suehring
45! flight module added
46!
47! 1919 2016-05-27 14:51:23Z raasch
48! Initial version of purely vertical nesting introduced.
49!
50! 1918 2016-05-27 14:35:57Z raasch
51! determination of time step moved to the end of the time step loop,
52! the first time step is now always calculated before the time step loop (i.e.
53! also in case of restart runs)
54!
55! 1914 2016-05-26 14:44:07Z witha
56! Added call for wind turbine model
57!
58! 1878 2016-04-19 12:30:36Z hellstea
59! Synchronization for nested runs rewritten
60!
61! 1853 2016-04-11 09:00:35Z maronga
62! Adjusted for use with radiation_scheme = constant
63!
64! 1849 2016-04-08 11:33:18Z hoffmann
65! Adapted for modularization of microphysics
66!
67! 1833 2016-04-07 14:23:03Z raasch
68! spectrum renamed spectra_mod, spectra related variables moved to spectra_mod
69!
70! 1831 2016-04-07 13:15:51Z hoffmann
71! turbulence renamed collision_turbulence
72!
73! 1822 2016-04-07 07:49:42Z hoffmann
74! icloud_scheme replaced by microphysics_*
75!
76! 1808 2016-04-05 19:44:00Z raasch
77! output message in case unscheduled radiation calls removed
78!
79! 1797 2016-03-21 16:50:28Z raasch
80! introduction of different datatransfer modes
81!
82! 1791 2016-03-11 10:41:25Z raasch
83! call of pmci_update_new removed
84!
85! 1786 2016-03-08 05:49:27Z raasch
86! +module spectrum
87!
88! 1783 2016-03-06 18:36:17Z raasch
89! switch back of netcdf data format for mask output moved to the mask output
90! routine
91!
92! 1781 2016-03-03 15:12:23Z raasch
93! some pmc calls removed at the beginning (before timeloop),
94! pmc initialization moved to the main program
95!
96! 1764 2016-02-28 12:45:19Z raasch
97! PMC_ACTIVE flags removed,
98! bugfix: nest synchronization after first call of timestep
99!
100! 1762 2016-02-25 12:31:13Z hellstea
101! Introduction of nested domain feature
102!
103! 1736 2015-12-04 08:56:33Z raasch
104! no perturbations added to total domain if energy limit has been set zero
105!
106! 1691 2015-10-26 16:17:44Z maronga
107! Added option for spin-ups without land surface and radiation models. Moved calls
108! for radiation and lan surface schemes.
109!
110! 1682 2015-10-07 23:56:08Z knoop
111! Code annotations made doxygen readable
112!
113! 1671 2015-09-25 03:29:37Z raasch
114! bugfix: ghostpoint exchange for array diss in case that sgs velocities are used
115! for particles
116!
117! 1585 2015-04-30 07:05:52Z maronga
118! Moved call of radiation scheme. Added support for RRTM
119!
120! 1551 2015-03-03 14:18:16Z maronga
121! Added interface for different radiation schemes.
122!
123! 1496 2014-12-02 17:25:50Z maronga
124! Added calls for the land surface model and radiation scheme
125!
126! 1402 2014-05-09 14:25:13Z raasch
127! location messages modified
128!
129! 1384 2014-05-02 14:31:06Z raasch
130! location messages added
131!
132! 1380 2014-04-28 12:40:45Z heinze
133! CALL of nudge_ref added
134! bc_pt_t_val and bc_q_t_val are updated in case nudging is used
135!
136! 1365 2014-04-22 15:03:56Z boeske
137! Reset sums_ls_l to zero at each timestep
138! +sums_ls_l
139! Calculation of reference state (previously in subroutine calc_mean_profile)
140
141! 1342 2014-03-26 17:04:47Z kanani
142! REAL constants defined as wp-kind
143!
144! 1320 2014-03-20 08:40:49Z raasch
145! ONLY-attribute added to USE-statements,
146! kind-parameters added to all INTEGER and REAL declaration statements,
147! kinds are defined in new module kinds,
148! old module precision_kind is removed,
149! revision history before 2012 removed,
150! comment fields (!:) to be used for variable explanations added to
151! all variable declaration statements
152! 1318 2014-03-17 13:35:16Z raasch
153! module interfaces removed
154!
155! 1308 2014-03-13 14:58:42Z fricke
156! +netcdf_data_format_save
157! For masked data, parallel netcdf output is not tested so far, hence
158! netcdf_data_format is switched back to non-paralell output.
159!
160! 1276 2014-01-15 13:40:41Z heinze
161! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
162!
163! 1257 2013-11-08 15:18:40Z raasch
164! acc-update-host directive for timestep removed
165!
166! 1241 2013-10-30 11:36:58Z heinze
167! Generalize calc_mean_profile for wider use
168! Determine shf and qsws in dependence on data from LSF_DATA
169! Determine ug and vg in dependence on data from LSF_DATA
170! 1221 2013-09-10 08:59:13Z raasch
171! host update of arrays before timestep is called
172!
173! 1179 2013-06-14 05:57:58Z raasch
174! mean profiles for reference state are only calculated if required,
175! small bugfix for background communication
176!
177! 1171 2013-05-30 11:27:45Z raasch
178! split of prognostic_equations deactivated (comment lines), for the time being
179!
180! 1128 2013-04-12 06:19:32Z raasch
181! asynchronous transfer of ghost point data realized for acc-optimized version:
182! prognostic_equations are first called two times for those points required for
183! the left-right and north-south exchange, respectively, and then for the
184! remaining points,
185! those parts requiring global communication moved from prognostic_equations to
186! here
187!
188! 1115 2013-03-26 18:16:16Z hoffmann
189! calculation of qr and nr is restricted to precipitation
190!
191! 1113 2013-03-10 02:48:14Z raasch
192! GPU-porting of boundary conditions,
193! openACC directives updated
194! formal parameter removed from routine boundary_conds
195!
196! 1111 2013-03-08 23:54:10Z raasch
197! +internal timestep counter for cpu statistics added,
198! openACC directives updated
199!
200! 1092 2013-02-02 11:24:22Z raasch
201! unused variables removed
202!
203! 1065 2012-11-22 17:42:36Z hoffmann
204! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
205!
206! 1053 2012-11-13 17:11:03Z hoffmann
207! exchange of ghost points for nr, qr added
208!
209! 1036 2012-10-22 13:43:42Z raasch
210! code put under GPL (PALM 3.9)
211!
212! 1019 2012-09-28 06:46:45Z raasch
213! non-optimized version of prognostic_equations removed
214!
215! 1015 2012-09-27 09:23:24Z raasch
216! +call of prognostic_equations_acc
217!
218! 1001 2012-09-13 14:08:46Z raasch
219! all actions concerning leapfrog- and upstream-spline-scheme removed
220!
221! 849 2012-03-15 10:35:09Z raasch
222! advec_particles renamed lpm, first_call_advec_particles renamed first_call_lpm
223!
224! 825 2012-02-19 03:03:44Z raasch
225! wang_collision_kernel renamed wang_kernel
226!
227! Revision 1.1  1997/08/11 06:19:04  raasch
228! Initial revision
229!
230!
231! Description:
232! ------------
233!> Integration in time of the model equations, statistical analysis and graphic
234!> output
235!------------------------------------------------------------------------------!
236 SUBROUTINE time_integration
237 
238
239    USE advec_ws,                                                              &
240        ONLY:  ws_statistics
241
242    USE arrays_3d,                                                             &
243        ONLY:  diss, dzu, e, e_p, nr_p, prho, pt, pt_p, pt_init, q_init, q,    &
244               ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho, s, s_p, sa_p, &
245               tend, u, u_p, v, vpt, v_p, w, w_p
246
247    USE calc_mean_profile_mod,                                                 &
248        ONLY:  calc_mean_profile
249
250    USE control_parameters,                                                    &
251        ONLY:  advected_distance_x, advected_distance_y, average_count_3d,     &
252               averaging_interval, averaging_interval_pr,                      &
253               bc_lr_cyc, bc_ns_cyc, bc_pt_t_val,                              &
254               bc_q_t_val, call_psolver_at_all_substeps, cloud_droplets,       &
255               cloud_physics, constant_flux_layer, constant_heatflux,          &
256               create_disturbances, dopr_n, constant_diffusion, coupling_mode, &
257               coupling_start_time, current_timestep_number,                   &
258               disturbance_created, disturbance_energy_limit, dist_range,      &
259               do_sum, dt_3d, dt_averaging_input, dt_averaging_input_pr,       &
260               dt_coupling, dt_data_output_av, dt_disturb, dt_do2d_xy,         &
261               dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_domask,dt_dopts, dt_dopr,   &
262               dt_dopr_listing, dt_dots, dt_dvrp, dt_run_control,              &
263               end_time, first_call_lpm, galilei_transformation, humidity,     &
264               intermediate_timestep_count,                                    &
265               intermediate_timestep_count_max, large_scale_forcing,           &
266               loop_optimization, lsf_surf, lsf_vert, masks,                   &
267               microphysics_seifert, mid, nest_domain,                         &
268               neutral, nr_timesteps_this_run, nudging,                        &
269               ocean, on_device, passive_scalar,                               &
270               prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
271               run_coupled, simulated_time, simulated_time_chr,                &
272               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
273               skip_time_do3d, skip_time_domask, skip_time_dopr,               &
274               skip_time_data_output_av, sloping_surface,                      &
275               stop_dt, terminate_coupled, terminate_run, timestep_scheme,     &
276               time_coupling, time_do2d_xy, time_do2d_xz, time_do2d_yz,        &
277               time_do3d, time_domask, time_dopr, time_dopr_av,                &
278               time_dopr_listing, time_dopts, time_dosp, time_dosp_av,         &
279               time_dots, time_do_av, time_do_sla, time_disturb, time_dvrp,    &
280               time_run_control, time_since_reference_point,                   &
281               turbulent_inflow, urban_surface,                                &
282               use_initial_profile_as_reference,                               &
283               use_single_reference_value, u_gtrans, v_gtrans, virtual_flight, &
284               ws_scheme_mom, ws_scheme_sca
285
286    USE cpulog,                                                                &
287        ONLY:  cpu_log, log_point, log_point_s
288
289    USE flight_mod,                                                            &
290        ONLY:  flight_measurement
291
292
293    USE indices,                                                               &
294        ONLY:  i_left, i_right, j_north, j_south, nbgp, nx, nxl, nxlg, nxr,    &
295               nxrg, nyn, nyng, nys, nysg, nzb, nzt, nzb_u_inner, nzb_v_inner
296
297    USE interaction_droplets_ptq_mod,                                          &
298        ONLY:  interaction_droplets_ptq
299
300    USE interfaces
301
302    USE kinds
303
304    USE land_surface_model_mod,                                                &
305        ONLY:  land_surface, lsm_energy_balance, lsm_soil_model,               &
306               skip_time_do_lsm
307
308    USE ls_forcing_mod,                                                        &
309        ONLY:  ls_forcing_surf, ls_forcing_vert
310
311    USE microphysics_mod,                                                      &
312        ONLY: collision_turbulence
313
314    USE nudge_mod,                                                             &
315        ONLY:  calc_tnudge, nudge_ref
316
317    USE particle_attributes,                                                   &
318        ONLY:  particle_advection, particle_advection_start,                   &
319               use_sgs_for_particles, wang_kernel
320
321    USE pegrid
322
323    USE pmc_interface,                                                         &
324        ONLY:  nested_run, nesting_mode, pmci_datatrans,                       &
325               pmci_ensure_nest_mass_conservation, pmci_synchronize
326
327    USE production_e_mod,                                                      &
328        ONLY:  production_e_init
329
330    USE progress_bar,                                                          &
331        ONLY:  finish_progress_bar, output_progress_bar
332
333    USE prognostic_equations_mod,                                              &
334        ONLY:  prognostic_equations_acc, prognostic_equations_cache,           &
335               prognostic_equations_vector
336
337    USE radiation_model_mod,                                                   &
338        ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,&
339              skip_time_do_radiation, time_radiation
340
341    USE spectra_mod,                                                           &
342        ONLY: average_count_sp, averaging_interval_sp, calc_spectra, dt_dosp,  &
343              skip_time_dosp
344
345    USE statistics,                                                            &
346        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l, u_max,         &
347               u_max_ijk, v_max, v_max_ijk, w_max, w_max_ijk
348
349    USE surface_layer_fluxes_mod,                                              &
350        ONLY:  surface_layer_fluxes
351
352    USE urban_surface_mod,                                                     &
353        ONLY:  usm_material_heat_model, usm_material_model,     &
354               usm_radiation, usm_surface_energy_balance             
355
356    USE user_actions_mod,                                                      &
357        ONLY:  user_actions
358
359    USE wind_turbine_model_mod,                                                &
360        ONLY:  wind_turbine, wtm_forces
361
362    IMPLICIT NONE
363
364    CHARACTER (LEN=9) ::  time_to_string          !<
365
366    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
367                            !< steering of run control output interval
368
369!
370!-- At beginning determine the first time step
371    CALL timestep
372
373!
374!-- Synchronize the timestep in case of nested run.
375    IF ( nested_run )  THEN
376!
377!--    Synchronization by unifying the time step.
378!--    Global minimum of all time-steps is used for all.
379       CALL pmci_synchronize
380    ENDIF
381
382!
383!-- Determine and print out the run control quantities before the first time
384!-- step of this run. For the initial run, some statistics (e.g. divergence)
385!-- need to be determined first.
386    IF ( simulated_time == 0.0_wp )  CALL flow_statistics
387    CALL run_control
388
389!
390!-- Data exchange between coupled models in case that a call has been omitted
391!-- at the end of the previous run of a job chain.
392    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
393!
394!--    In case of model termination initiated by the local model the coupler
395!--    must not be called because this would again cause an MPI hang.
396       DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
397          CALL surface_coupler
398          time_coupling = time_coupling - dt_coupling
399       ENDDO
400       IF (time_coupling == 0.0_wp  .AND.                                      &
401           time_since_reference_point < dt_coupling )                          &
402       THEN
403          time_coupling = time_since_reference_point
404       ENDIF
405    ENDIF
406
407#if defined( __dvrp_graphics )
408!
409!-- Time measurement with dvrp software 
410    CALL DVRP_LOG_EVENT( 2, current_timestep_number )
411#endif
412
413    CALL location_message( 'start with time-stepping', .TRUE. )
414!
415!-- Start of the time loop
416    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. &
417                .NOT. terminate_run )
418
419       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
420
421!
422!--    Determine ug, vg and w_subs in dependence on data from external file
423!--    LSF_DATA
424       IF ( large_scale_forcing .AND. lsf_vert )  THEN
425           CALL ls_forcing_vert ( simulated_time )
426           sums_ls_l = 0.0_wp
427       ENDIF
428
429!
430!--    Set pt_init and q_init to the current profiles taken from
431!--    NUDGING_DATA
432       IF ( nudging )  THEN
433           CALL nudge_ref ( simulated_time )
434!
435!--        Store temperature gradient at the top boundary for possible Neumann
436!--        boundary condition
437           bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
438           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
439       ENDIF
440
441!
442!--    Execute the user-defined actions
443       CALL user_actions( 'before_timestep' )
444
445!
446!--    Calculate forces by wind turbines
447       IF ( wind_turbine )  THEN
448
449          CALL cpu_log( log_point(55), 'wind_turbine', 'start' )
450
451          CALL wtm_forces
452
453          CALL cpu_log( log_point(55), 'wind_turbine', 'stop' )
454
455       ENDIF       
456       
457!
458!--    Start of intermediate step loop
459       intermediate_timestep_count = 0
460       DO  WHILE ( intermediate_timestep_count < &
461                   intermediate_timestep_count_max )
462
463          intermediate_timestep_count = intermediate_timestep_count + 1
464
465!
466!--       Set the steering factors for the prognostic equations which depend
467!--       on the timestep scheme
468          CALL timestep_scheme_steering
469
470!
471!--       Calculate those variables needed in the tendency terms which need
472!--       global communication
473          IF ( .NOT. use_single_reference_value  .AND. &
474               .NOT. use_initial_profile_as_reference )  THEN
475!
476!--          Horizontally averaged profiles to be used as reference state in
477!--          buoyancy terms (WARNING: only the respective last call of
478!--          calc_mean_profile defines the reference state!)
479             IF ( .NOT. neutral )  THEN
480                CALL calc_mean_profile( pt, 4 )
481                ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
482             ENDIF
483             IF ( ocean )  THEN
484                CALL calc_mean_profile( rho, 64 )
485                ref_state(:)  = hom(:,1,64,0)
486             ENDIF
487             IF ( humidity )  THEN
488                CALL calc_mean_profile( vpt, 44 )
489                ref_state(:)  = hom(:,1,44,0)
490             ENDIF
491
492          ENDIF
493
494          IF ( .NOT. constant_diffusion )  CALL production_e_init
495          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
496               intermediate_timestep_count == 1 )  CALL ws_statistics
497!
498!--       In case of nudging calculate current nudging time scale and horizontal
499!--       means of u, v, pt and q
500          IF ( nudging )  THEN
501             CALL calc_tnudge( simulated_time )
502             CALL calc_mean_profile( u, 1 )
503             CALL calc_mean_profile( v, 2 )
504             CALL calc_mean_profile( pt, 4 )
505             CALL calc_mean_profile( q, 41 )
506          ENDIF
507
508!
509!--       Solve the prognostic equations. A fast cache optimized version with
510!--       only one single loop is used in case of Piascek-Williams advection
511!--       scheme. NEC vector machines use a different version, because
512!--       in the other versions a good vectorization is prohibited due to
513!--       inlining problems.
514          IF ( loop_optimization == 'cache' )  THEN
515             CALL prognostic_equations_cache
516          ELSEIF ( loop_optimization == 'vector' )  THEN
517             CALL prognostic_equations_vector
518          ELSEIF ( loop_optimization == 'acc' )  THEN
519             i_left  = nxl;         i_right = nxr
520             j_south = nys;         j_north = nyn
521             CALL prognostic_equations_acc
522
523!             i_left  = nxl;         i_right = nxl+nbgp-1
524!             j_south = nys;         j_north = nyn
525!             CALL prognostic_equations_acc
526!             i_left  = nxr-nbgp+1;  i_right = nxr
527!             j_south = nys;         j_north = nyn
528!             CALL prognostic_equations_acc
529
530!
531!--          Exchange of ghost points (lateral boundary conditions)
532             IF ( background_communication )  THEN
533
534                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
535               
536                send_receive = 'lr'
537                sendrecv_in_background = .TRUE.
538                req          = 0
539                req_count    = 0
540
541                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
542                   on_device = .TRUE.         ! to be removed after complete porting
543                ELSE                          ! of ghost point exchange
544                   !$acc update host( e_p, pt_p, u_p, v_p, w_p )
545                ENDIF
546
547                CALL exchange_horiz( u_p, nbgp )
548                CALL exchange_horiz( v_p, nbgp )
549                CALL exchange_horiz( w_p, nbgp )
550                CALL exchange_horiz( pt_p, nbgp )
551                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
552                IF ( ocean )  THEN
553                   CALL exchange_horiz( sa_p, nbgp )
554                   CALL exchange_horiz( rho, nbgp )
555                   CALL exchange_horiz( prho, nbgp )
556                ENDIF
557                IF ( humidity )  THEN
558                   CALL exchange_horiz( q_p, nbgp )
559                   IF ( cloud_physics .AND. microphysics_seifert )  THEN
560                      CALL exchange_horiz( qr_p, nbgp )
561                      CALL exchange_horiz( nr_p, nbgp )
562                   ENDIF
563                ENDIF
564                IF ( cloud_droplets )  THEN
565                   CALL exchange_horiz( ql, nbgp )
566                   CALL exchange_horiz( ql_c, nbgp )
567                   CALL exchange_horiz( ql_v, nbgp )
568                   CALL exchange_horiz( ql_vp, nbgp )
569                ENDIF
570                IF ( wang_kernel  .OR.  collision_turbulence  .OR.             &
571                     use_sgs_for_particles )  THEN
572                   CALL exchange_horiz( diss, nbgp )
573                ENDIF
574                IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
575
576                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
577                   on_device = .FALSE.        ! to be removed after complete porting
578                ELSE                          ! of ghost point exchange
579                   !$acc update device( e_p, pt_p, u_p, v_p, w_p )
580                ENDIF
581
582                sendrecv_in_background = .FALSE.
583
584                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'pause' )
585
586             ENDIF
587
588!             i_left  = nxl+nbgp;    i_right = nxr-nbgp
589!             j_south = nys;         j_north = nys+nbgp-1
590!             CALL prognostic_equations_acc
591!             i_left  = nxl+nbgp;    i_right = nxr-nbgp
592!             j_south = nyn-nbgp+1;  j_north = nyn
593!             CALL prognostic_equations_acc
594
595             IF ( background_communication )  THEN
596                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'start' )
597#if defined( __parallel )
598                CALL MPI_WAITALL( req_count, req, wait_stat, ierr )
599#endif
600                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'pause' )
601
602                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'continue' )
603
604                send_receive = 'ns'
605                sendrecv_in_background = .TRUE.
606                req          = 0
607                req_count    = 0
608
609                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
610                   on_device = .TRUE.         ! to be removed after complete porting
611                ELSE                          ! of ghost point exchange
612                   !$acc update host( e_p, pt_p, u_p, v_p, w_p )
613                ENDIF
614
615                CALL exchange_horiz( u_p, nbgp )
616                CALL exchange_horiz( v_p, nbgp )
617                CALL exchange_horiz( w_p, nbgp )
618                CALL exchange_horiz( pt_p, nbgp )
619                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
620                IF ( ocean )  THEN
621                   CALL exchange_horiz( sa_p, nbgp )
622                   CALL exchange_horiz( rho, nbgp )
623                  CALL exchange_horiz( prho, nbgp )
624                ENDIF
625                IF ( humidity )  THEN
626                   CALL exchange_horiz( q_p, nbgp )
627                   IF ( cloud_physics .AND. microphysics_seifert )  THEN
628                      CALL exchange_horiz( qr_p, nbgp )
629                      CALL exchange_horiz( nr_p, nbgp )
630                   ENDIF
631                ENDIF
632                IF ( cloud_droplets )  THEN
633                   CALL exchange_horiz( ql, nbgp )
634                   CALL exchange_horiz( ql_c, nbgp )
635                   CALL exchange_horiz( ql_v, nbgp )
636                   CALL exchange_horiz( ql_vp, nbgp )
637                ENDIF
638                IF ( wang_kernel  .OR.  collision_turbulence  .OR.             &
639                     use_sgs_for_particles )  THEN
640                   CALL exchange_horiz( diss, nbgp )
641                ENDIF
642                IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
643
644                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
645                   on_device = .FALSE.        ! to be removed after complete porting
646                ELSE                          ! of ghost point exchange
647                   !$acc update device( e_p, pt_p, u_p, v_p, w_p )
648                ENDIF
649
650                sendrecv_in_background = .FALSE.
651
652                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
653
654             ENDIF
655
656!             i_left  = nxl+nbgp;    i_right = nxr-nbgp
657!             j_south = nys+nbgp;    j_north = nyn-nbgp
658!             CALL prognostic_equations_acc
659
660             IF ( background_communication )  THEN
661                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'continue' )
662#if defined( __parallel )
663                CALL MPI_WAITALL( req_count, req, wait_stat, ierr )
664#endif
665                send_receive = 'al'
666                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'stop' )
667             ENDIF
668
669          ENDIF
670
671!
672!--       Particle transport/physics with the Lagrangian particle model
673!--       (only once during intermediate steps, because it uses an Euler-step)
674!--       ### particle model should be moved before prognostic_equations, in order
675!--       to regard droplet interactions directly
676          IF ( particle_advection  .AND.                         &
677               simulated_time >= particle_advection_start  .AND. &
678               intermediate_timestep_count == 1 )  THEN
679             CALL lpm
680             first_call_lpm = .FALSE.
681          ENDIF
682
683!
684!--       Interaction of droplets with temperature and specific humidity.
685!--       Droplet condensation and evaporation is calculated within
686!--       advec_particles.
687          IF ( cloud_droplets  .AND.  &
688               intermediate_timestep_count == intermediate_timestep_count_max )&
689          THEN
690             CALL interaction_droplets_ptq
691          ENDIF
692
693!
694!--       Exchange of ghost points (lateral boundary conditions)
695          IF ( .NOT. background_communication )  THEN
696
697             CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
698
699             IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
700                on_device = .TRUE.         ! to be removed after complete porting
701             ELSE                          ! of ghost point exchange
702                !$acc update host( e_p, pt_p, u_p, v_p, w_p )
703             ENDIF
704
705             CALL exchange_horiz( u_p, nbgp )
706             CALL exchange_horiz( v_p, nbgp )
707             CALL exchange_horiz( w_p, nbgp )
708             CALL exchange_horiz( pt_p, nbgp )
709             IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
710             IF ( ocean )  THEN
711                CALL exchange_horiz( sa_p, nbgp )
712                CALL exchange_horiz( rho, nbgp )
713                CALL exchange_horiz( prho, nbgp )
714             ENDIF
715             IF ( humidity )  THEN
716                CALL exchange_horiz( q_p, nbgp )
717                IF ( cloud_physics .AND. microphysics_seifert )  THEN
718                   CALL exchange_horiz( qr_p, nbgp )
719                   CALL exchange_horiz( nr_p, nbgp )
720                ENDIF
721             ENDIF
722             IF ( cloud_droplets )  THEN
723                CALL exchange_horiz( ql, nbgp )
724                CALL exchange_horiz( ql_c, nbgp )
725                CALL exchange_horiz( ql_v, nbgp )
726                CALL exchange_horiz( ql_vp, nbgp )
727             ENDIF
728             IF ( wang_kernel  .OR.  collision_turbulence  .OR.                &
729                  use_sgs_for_particles )  THEN
730                CALL exchange_horiz( diss, nbgp )
731             ENDIF
732             IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )     
733
734             IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
735                on_device = .FALSE.        ! to be removed after complete porting
736             ELSE                          ! of ghost point exchange
737                !$acc update device( e_p, pt_p, u_p, v_p, w_p )
738             ENDIF
739
740             CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
741
742          ENDIF
743
744!
745!--       Boundary conditions for the prognostic quantities (except of the
746!--       velocities at the outflow in case of a non-cyclic lateral wall)
747          CALL boundary_conds
748
749!
750!--       Swap the time levels in preparation for the next time step.
751          CALL swap_timelevel
752
753          IF ( nested_run )  THEN
754
755             CALL cpu_log( log_point(60), 'nesting', 'start' )
756!
757!--          Domain nesting. The data transfer subroutines pmci_parent_datatrans
758!--          and pmci_child_datatrans are called inside the wrapper
759!--          subroutine pmci_datatrans according to the control parameters
760!--          nesting_mode and nesting_datatransfer_mode.
761!--          TO_DO: why is nesting_mode given as a parameter here?
762             CALL pmci_datatrans( nesting_mode )
763
764             IF ( TRIM( nesting_mode ) == 'two-way' .OR.                               &
765                  nesting_mode == 'vertical' )  THEN
766!
767!--             Exchange_horiz is needed for all parent-domains after the
768!--             anterpolation
769                CALL exchange_horiz( u, nbgp )
770                CALL exchange_horiz( v, nbgp )
771                CALL exchange_horiz( w, nbgp )
772                IF ( .NOT. neutral            )  CALL exchange_horiz( pt, nbgp )
773                IF ( humidity                 )  CALL exchange_horiz( q, nbgp  )
774                IF ( passive_scalar           )  CALL exchange_horiz( s, nbgp  )
775                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp  )
776             ENDIF
777!
778!--          Correct the w top-BC in nest domains to ensure mass conservation.
779!--          This action must never be done for the root domain. Vertical
780!--          nesting implies mass conservation.
781             IF ( nest_domain )  THEN
782                CALL pmci_ensure_nest_mass_conservation
783             ENDIF
784
785             CALL cpu_log( log_point(60), 'nesting', 'stop' )
786
787          ENDIF
788
789!
790!--       Temperature offset must be imposed at cyclic boundaries in x-direction
791!--       when a sloping surface is used
792          IF ( sloping_surface )  THEN
793             IF ( nxl ==  0 )  pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - &
794                                                    pt_slope_offset
795             IF ( nxr == nx )  pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + &
796                                                    pt_slope_offset
797          ENDIF
798
799!
800!--       Impose a turbulent inflow using the recycling method
801          IF ( turbulent_inflow )  CALL  inflow_turbulence
802
803!
804!--       Impose a random perturbation on the horizontal velocity field
805          IF ( create_disturbances  .AND.                                      &
806               ( call_psolver_at_all_substeps  .AND.                           &
807               intermediate_timestep_count == intermediate_timestep_count_max )&
808          .OR. ( .NOT. call_psolver_at_all_substeps  .AND.                     &
809               intermediate_timestep_count == 1 ) )                            &
810          THEN
811             time_disturb = time_disturb + dt_3d
812             IF ( time_disturb >= dt_disturb )  THEN
813                !$acc update host( u, v )
814                IF ( numprocs == 1 )  on_device = .FALSE.  ! workaround, remove later
815                IF ( disturbance_energy_limit /= 0.0_wp  .AND.                 &
816                     hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
817                   CALL disturb_field( nzb_u_inner, tend, u )
818                   CALL disturb_field( nzb_v_inner, tend, v )
819                ELSEIF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
820!
821!--                Runs with a non-cyclic lateral wall need perturbations
822!--                near the inflow throughout the whole simulation
823                   dist_range = 1
824                   CALL disturb_field( nzb_u_inner, tend, u )
825                   CALL disturb_field( nzb_v_inner, tend, v )
826                   dist_range = 0
827                ENDIF
828                IF ( numprocs == 1 )  on_device = .TRUE.  ! workaround, remove later
829                !$acc update device( u, v )
830                time_disturb = time_disturb - dt_disturb
831             ENDIF
832          ENDIF
833
834!
835!--       Reduce the velocity divergence via the equation for perturbation
836!--       pressure.
837          IF ( intermediate_timestep_count == 1  .OR. &
838                call_psolver_at_all_substeps )  THEN
839             CALL pres
840          ENDIF
841
842!
843!--       If required, compute liquid water content
844          IF ( cloud_physics )  THEN
845             CALL calc_liquid_water_content
846             !$acc update device( ql )
847          ENDIF
848!
849!--       If required, compute virtual potential temperature
850          IF ( humidity )  THEN
851             CALL compute_vpt 
852             !$acc update device( vpt )
853          ENDIF 
854
855!
856!--       Compute the diffusion quantities
857          IF ( .NOT. constant_diffusion )  THEN
858
859!
860!--          Determine surface fluxes shf and qsws and surface values
861!--          pt_surface and q_surface in dependence on data from external
862!--          file LSF_DATA respectively
863             IF ( ( large_scale_forcing .AND. lsf_surf ) .AND. &
864                 intermediate_timestep_count == intermediate_timestep_count_max )&
865             THEN
866                CALL ls_forcing_surf ( simulated_time )
867             ENDIF
868
869!
870!--          First the vertical fluxes in the surface (constant flux) layer are computed
871             IF ( constant_flux_layer )  THEN
872                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
873                CALL surface_layer_fluxes
874                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
875             ENDIF
876
877!
878!--          If required, solve the energy balance for the surface and run soil
879!--          model
880             IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
881
882                CALL cpu_log( log_point(54), 'land_surface', 'start' )
883                CALL lsm_energy_balance
884                CALL lsm_soil_model
885                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
886             ENDIF
887
888!
889!--          If required, solve the energy balance for urban surfaces and run
890!--          the material heat model
891             IF (urban_surface) THEN
892                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
893                CALL usm_surface_energy_balance
894                IF ( usm_material_model )  THEN
895                   CALL usm_material_heat_model
896                ENDIF
897                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
898             ENDIF
899
900!
901!--          Compute the diffusion coefficients
902             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
903             IF ( .NOT. humidity ) THEN
904                IF ( ocean )  THEN
905                   CALL diffusivities( prho, prho_reference )
906                ELSE
907                   CALL diffusivities( pt, pt_reference )
908                ENDIF
909             ELSE
910                CALL diffusivities( vpt, pt_reference )
911             ENDIF
912             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
913
914          ENDIF
915
916!
917!--       If required, calculate radiative fluxes and heating rates
918          IF ( radiation .AND. intermediate_timestep_count                     &
919               == intermediate_timestep_count_max .AND. simulated_time >    &
920               skip_time_do_radiation )  THEN
921
922               time_radiation = time_radiation + dt_3d
923
924             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
925             THEN
926
927                CALL cpu_log( log_point(50), 'radiation', 'start' )
928
929                IF ( .NOT. force_radiation_call )  THEN
930                   time_radiation = time_radiation - dt_radiation
931                ENDIF
932
933                CALL radiation_control
934
935                CALL cpu_log( log_point(50), 'radiation', 'stop' )
936
937                IF (urban_surface)  THEN
938                   CALL cpu_log( log_point(75), 'usm_radiation', 'start' )
939                   CALL usm_radiation
940                   CALL cpu_log( log_point(75), 'usm_radiation', 'stop' )
941                ENDIF
942
943             ENDIF
944          ENDIF
945
946       ENDDO   ! Intermediate step loop
947
948!
949!--    Increase simulation time and output times
950       nr_timesteps_this_run      = nr_timesteps_this_run + 1
951       current_timestep_number    = current_timestep_number + 1
952       simulated_time             = simulated_time   + dt_3d
953       simulated_time_chr         = time_to_string( simulated_time )
954       time_since_reference_point = simulated_time - coupling_start_time
955
956
957
958       IF ( simulated_time >= skip_time_data_output_av )  THEN
959          time_do_av         = time_do_av       + dt_3d
960       ENDIF
961       IF ( simulated_time >= skip_time_do2d_xy )  THEN
962          time_do2d_xy       = time_do2d_xy     + dt_3d
963       ENDIF
964       IF ( simulated_time >= skip_time_do2d_xz )  THEN
965          time_do2d_xz       = time_do2d_xz     + dt_3d
966       ENDIF
967       IF ( simulated_time >= skip_time_do2d_yz )  THEN
968          time_do2d_yz       = time_do2d_yz     + dt_3d
969       ENDIF
970       IF ( simulated_time >= skip_time_do3d    )  THEN
971          time_do3d          = time_do3d        + dt_3d
972       ENDIF
973       DO  mid = 1, masks
974          IF ( simulated_time >= skip_time_domask(mid) )  THEN
975             time_domask(mid)= time_domask(mid) + dt_3d
976          ENDIF
977       ENDDO
978       time_dvrp          = time_dvrp        + dt_3d
979       IF ( simulated_time >= skip_time_dosp )  THEN
980          time_dosp       = time_dosp        + dt_3d
981       ENDIF
982       time_dots          = time_dots        + dt_3d
983       IF ( .NOT. first_call_lpm )  THEN
984          time_dopts      = time_dopts       + dt_3d
985       ENDIF
986       IF ( simulated_time >= skip_time_dopr )  THEN
987          time_dopr       = time_dopr        + dt_3d
988       ENDIF
989       time_dopr_listing          = time_dopr_listing        + dt_3d
990       time_run_control   = time_run_control + dt_3d
991
992!
993!--    Data exchange between coupled models
994       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
995          time_coupling = time_coupling + dt_3d
996
997!
998!--       In case of model termination initiated by the local model
999!--       (terminate_coupled > 0), the coupler must be skipped because it would
1000!--       cause an MPI intercomminucation hang.
1001!--       If necessary, the coupler will be called at the beginning of the
1002!--       next restart run.
1003          DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
1004             CALL surface_coupler
1005             time_coupling = time_coupling - dt_coupling
1006          ENDDO
1007       ENDIF
1008
1009!
1010!--    Execute user-defined actions
1011       CALL user_actions( 'after_integration' )
1012
1013!
1014!--    If Galilei transformation is used, determine the distance that the
1015!--    model has moved so far
1016       IF ( galilei_transformation )  THEN
1017          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
1018          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
1019       ENDIF
1020
1021!
1022!--    Check, if restart is necessary (because cpu-time is expiring or
1023!--    because it is forced by user) and set stop flag
1024!--    This call is skipped if the remote model has already initiated a restart.
1025       IF ( .NOT. terminate_run )  CALL check_for_restart
1026
1027!
1028!--    Carry out statistical analysis and output at the requested output times.
1029!--    The MOD function is used for calculating the output time counters (like
1030!--    time_dopr) in order to regard a possible decrease of the output time
1031!--    interval in case of restart runs
1032
1033!
1034!--    Set a flag indicating that so far no statistics have been created
1035!--    for this time step
1036       flow_statistics_called = .FALSE.
1037
1038!
1039!--    If required, call flow_statistics for averaging in time
1040       IF ( averaging_interval_pr /= 0.0_wp  .AND.  &
1041            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
1042            simulated_time >= skip_time_dopr )  THEN
1043          time_dopr_av = time_dopr_av + dt_3d
1044          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
1045             do_sum = .TRUE.
1046             time_dopr_av = MOD( time_dopr_av, &
1047                                    MAX( dt_averaging_input_pr, dt_3d ) )
1048          ENDIF
1049       ENDIF
1050       IF ( do_sum )  CALL flow_statistics
1051
1052!
1053!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
1054       IF ( averaging_interval /= 0.0_wp  .AND.                                &
1055            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
1056            simulated_time >= skip_time_data_output_av )                    &
1057       THEN
1058          time_do_sla = time_do_sla + dt_3d
1059          IF ( time_do_sla >= dt_averaging_input )  THEN
1060             CALL sum_up_3d_data
1061             average_count_3d = average_count_3d + 1
1062             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
1063          ENDIF
1064       ENDIF
1065
1066!
1067!--    Calculate spectra for time averaging
1068       IF ( averaging_interval_sp /= 0.0_wp  .AND.  &
1069            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
1070            simulated_time >= skip_time_dosp )  THEN
1071          time_dosp_av = time_dosp_av + dt_3d
1072          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
1073             CALL calc_spectra
1074             time_dosp_av = MOD( time_dosp_av, &
1075                                 MAX( dt_averaging_input_pr, dt_3d ) )
1076          ENDIF
1077       ENDIF
1078
1079!
1080!--    Call flight module and output data
1081       IF ( virtual_flight )  THEN
1082          CALL flight_measurement
1083          CALL data_output_flight
1084       ENDIF
1085
1086!
1087!--    Profile output (ASCII) on file
1088       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
1089          CALL print_1d
1090          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
1091                                                           dt_3d ) )
1092       ENDIF
1093
1094!
1095!--    Graphic output for PROFIL
1096       IF ( time_dopr >= dt_dopr )  THEN
1097          IF ( dopr_n /= 0 )  CALL data_output_profiles
1098          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
1099          time_dopr_av = 0.0_wp    ! due to averaging (see above)
1100       ENDIF
1101
1102!
1103!--    Graphic output for time series
1104       IF ( time_dots >= dt_dots )  THEN
1105          CALL data_output_tseries
1106          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
1107       ENDIF
1108
1109!
1110!--    Output of spectra (formatted for use with PROFIL), in case of no
1111!--    time averaging, spectra has to be calculated before
1112       IF ( time_dosp >= dt_dosp )  THEN
1113          IF ( average_count_sp == 0 )  CALL calc_spectra
1114          CALL data_output_spectra
1115          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
1116       ENDIF
1117
1118!
1119!--    2d-data output (cross-sections)
1120       IF ( time_do2d_xy >= dt_do2d_xy )  THEN
1121          CALL data_output_2d( 'xy', 0 )
1122          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
1123       ENDIF
1124       IF ( time_do2d_xz >= dt_do2d_xz )  THEN
1125          CALL data_output_2d( 'xz', 0 )
1126          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
1127       ENDIF
1128       IF ( time_do2d_yz >= dt_do2d_yz )  THEN
1129          CALL data_output_2d( 'yz', 0 )
1130          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
1131       ENDIF
1132
1133!
1134!--    3d-data output (volume data)
1135       IF ( time_do3d >= dt_do3d )  THEN
1136          CALL data_output_3d( 0 )
1137          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
1138       ENDIF
1139
1140!
1141!--    Masked data output
1142       DO  mid = 1, masks
1143          IF ( time_domask(mid) >= dt_domask(mid) )  THEN
1144             CALL data_output_mask( 0 )
1145             time_domask(mid) = MOD( time_domask(mid),  &
1146                                     MAX( dt_domask(mid), dt_3d ) )
1147          ENDIF
1148       ENDDO
1149
1150!
1151!--    Output of time-averaged 2d/3d/masked data
1152       IF ( time_do_av >= dt_data_output_av )  THEN
1153          CALL average_3d_data
1154          CALL data_output_2d( 'xy', 1 )
1155          CALL data_output_2d( 'xz', 1 )
1156          CALL data_output_2d( 'yz', 1 )
1157          CALL data_output_3d( 1 )
1158          DO  mid = 1, masks
1159             CALL data_output_mask( 1 )
1160          ENDDO
1161          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
1162       ENDIF
1163
1164!
1165!--    Output of particle time series
1166       IF ( particle_advection )  THEN
1167          IF ( time_dopts >= dt_dopts  .OR. &
1168               ( simulated_time >= particle_advection_start  .AND. &
1169                 first_call_lpm ) )  THEN
1170             CALL data_output_ptseries
1171             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
1172          ENDIF
1173       ENDIF
1174
1175!
1176!--    Output of dvrp-graphics (isosurface, particles, slicer)
1177#if defined( __dvrp_graphics )
1178       CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
1179#endif
1180       IF ( time_dvrp >= dt_dvrp )  THEN
1181          CALL data_output_dvrp
1182          time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
1183       ENDIF
1184#if defined( __dvrp_graphics )
1185       CALL DVRP_LOG_EVENT( 2, current_timestep_number )
1186#endif
1187
1188!
1189!--    If required, set the heat flux for the next time step at a random value
1190       IF ( constant_heatflux  .AND.  random_heatflux )  CALL disturb_heatflux
1191
1192!
1193!--    Execute user-defined actions
1194       CALL user_actions( 'after_timestep' )
1195
1196!
1197!--    Determine size of next time step. Save timestep dt_3d because it is
1198!--    newly calculated in routine timestep, but required further below for
1199!--    steering the run control output interval
1200       dt_3d_old = dt_3d
1201       CALL timestep
1202
1203!
1204!--    Synchronize the timestep in case of nested run.
1205       IF ( nested_run )  THEN
1206!
1207!--       Synchronize by unifying the time step.
1208!--       Global minimum of all time-steps is used for all.
1209          CALL pmci_synchronize
1210       ENDIF
1211
1212!
1213!--    Computation and output of run control parameters.
1214!--    This is also done whenever perturbations have been imposed
1215       IF ( time_run_control >= dt_run_control  .OR.                     &
1216            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
1217       THEN
1218          CALL run_control
1219          IF ( time_run_control >= dt_run_control )  THEN
1220             time_run_control = MOD( time_run_control, &
1221                                     MAX( dt_run_control, dt_3d_old ) )
1222          ENDIF
1223       ENDIF
1224
1225!
1226!--    Output elapsed simulated time in form of a progress bar on stdout
1227       IF ( myid == 0 )  CALL output_progress_bar
1228
1229       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
1230
1231
1232    ENDDO   ! time loop
1233
1234    IF ( myid == 0 )  CALL finish_progress_bar
1235
1236#if defined( __dvrp_graphics )
1237    CALL DVRP_LOG_EVENT( -2, current_timestep_number )
1238#endif
1239
1240    CALL location_message( 'finished time-stepping', .TRUE. )
1241
1242 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.