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

Last change on this file since 2050 was 2050, checked in by gronemeier, 5 years ago

Implement turbulent outflow condition

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