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

Last change on this file since 1941 was 1933, checked in by hellstea, 8 years ago

last commit documented

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