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

Last change on this file since 1958 was 1958, checked in by suehring, 8 years ago

last commit documented

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