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

Last change on this file since 2053 was 2051, checked in by gronemeier, 8 years ago

last commit documented

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