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

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

changes related to steering and formating of urban surface model

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