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

Last change on this file since 2032 was 2032, checked in by knoop, 5 years ago

last commit documented

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