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

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